Exclusive dimuon production in ultraperipheral Pb+Pb collisions at $\sqrt{s_{\mathrm{NN}}} = 5.02$ TeV with ATLAS

Exclusive dimuon production in ultraperipheral collisions (UPC), resulting from photon-photon interactions in the strong electromagnetic fields of colliding high-energy lead nuclei, $\mathrm{PbPb}(\gamma\gamma) \rightarrow \mu^+\mu^- (\mathrm{Pb}^{(\star)}\mathrm{Pb}^{(\star)} )$, is studied using $\mathcal{L}_{\mathrm{int}} = 0.48$ nb$^{-1}$ of $\sqrt{s_\mathrm{NN}}=5.02$ TeV lead-lead collision data at the LHC with the ATLAS detector. Dimuon pairs are measured in the fiducial region $p_{\mathrm{T}\mu}>4$ GeV, $|\eta_{\mu}|<2.4$, invariant mass $m_{\mu\mu}>10$ GeV, and $p_{\mathrm{T,\mu\mu}}<2$ GeV. The primary background from single-dissociative processes is extracted from the data using a template fitting technique. Differential cross sections are presented as a function of $m_{\mu\mu}$, absolute pair rapidity ($|y_{\mu\mu}|$), scattering angle in the dimuon rest frame ($|\cos \vartheta^{\star}_{\mu\mu}|$) and the colliding photon energies. The total cross section of the UPC $\gamma \gamma \rightarrow \mu^{+}\mu^{-}$ process in the fiducial volume is measured to be $\sigma_{\mathrm{fid}}^{\mu\mu} = 34.1 \! \pm \! 0.3 \mathrm{(stat.)} \! \pm \! 0.7 \mathrm{(syst.)}$ $\mu\mathrm{b}$. Generally good agreement is found with calculations from STARlight, which incorporate the leading-order Breit-Wheeler process with no final-state effects, albeit differences between the measurements and theoretical expectations are observed. In particular, the measured cross sections at larger $|y_{\mu\mu}|$ are found to be about 10-20% larger in data than in the calculations, suggesting the presence of larger fluxes of photons in the initial state. Modification of the dimuon cross sections in the presence of forward and/or backward neutron production is also studied and is found to be associated with a harder incoming photon spectrum, consistent with expectations.


Introduction
example of which is shown in Figure 1(b), where the muons are accompanied by additional resolved soft photons in the final state. Dissociative processes, where one photon is emitted by charged constituents of a nucleon, as shown in Figure 1(c), are also neglected by most models, in part due to the fact that these processes are not coherently enhanced.
The study of exclusive dimuon cross sections, conditional on observations of forward neutron production in the direction of one or both incoming nuclei, provides an additional experimental handle on the impact parameter range sampled in the observed events [12,[18][19][20]. In any particular collision, soft photons emitted by one lead nucleus (Pb) can excite the other (Pb ★ ), typically through the giant dipole resonance [21], and induce the emission of one or more neutrons, each of which carry, on average, the full per-nucleon beam energy. Since the probability of these excitations, as well as the overall hardness of the photon spectrum, is correlated with the nucleus-nucleus impact parameter [12], events with neutron excitation are typically correlated with harder photon collisions. In STARlight, dilepton cross sections associated with forward neutron production are calculated by convolving differential cross sections for low-energy photonuclear neutron production with the expected photon fluxes, thus in principle providing an essentially parameter-free prediction. Of course, the contribution from nucleonic dissociative processes must be subtracted before comparisons with data. detail in Section 9, and these variables are perhaps the most transparent way to compare the experimental data with yields determined by the nuclear photon fluxes assumed by the models.
The ALICE and CMS experiments have performed several measurements of exclusive vector-meson production in Pb+Pb collisions, for which the Breit-Wheeler process is one of the backgrounds. Reference [22] presents the measured contribution from the dilepton continuum out to invariant mass = 10 GeV. They found no observable contribution from the Υ, which is expected to have a cross section at least an order of magnitude below the continuum at LHC energies. ALICE has also measured / in UPC in +Pb collisions [23]. ATLAS [24,25] and CMS [26] have both measured exclusive high-mass dileptons in proton-proton ( ) collisions, and LHCb has measured exclusive Υ production [27]. At RHIC, where the maximum photon energy is about 3 GeV, STAR has measured low-mass continuum electron pairs [28,29], 0 photoproduction [30]. PHENIX has measured the cross section for / and the continuum below the / peak [31]. CDF has observed exclusive electron pair production out to = 30 GeV [32] and exclusive muons [33] in 3 < < 4 GeV. CMS has measured modifications of distributions, but not cross sections, as a function of the forward neutron topology [34]. In no case has a heavy-ion cross section measurement of the dilepton pair continuum been performed above 10 GeV, which impedes precise predictions for higher-mass processes, particularly light-by-light scattering [35][36][37][38].
This paper presents measurements of cross sections for exclusive dimuon production in UPC with lead nuclei, PbPb( ) → + − (Pb (★) Pb (★) ), using L int = 0.48 nb −1 of data taken in 2015. The measurements are performed for a fiducial region matched to the optimal performance for measuring muons in ATLAS. This requires single muons with transverse momentum T, > 4 GeV and pseudorapidity | | < 2.4, a pair invariant mass range 10 < < 200 GeV, and a pair T, < 2 GeV. The cross sections are measured by selecting events using a single-muon trigger, in association with an otherwise low-multiplicity event.
The events are required to consist only of two oppositely charged muons, with no other measured activity in the ATLAS inner detector. Events are then corrected for muon trigger and reconstruction efficiencies, migration relative to the fiducial region, and for dissociative backgrounds. The cross sections are presented as a function of the final-state dimuon pair mass ( ), pair rapidity ( ), center-of-mass scattering angle (| cos ★ |), and the dimuon acoplanarity. The two muon four-momenta can be combined to present cross sections as a function of the maximum and minimum energies of the incoming photons. Detailed comparisons of these data with predictions from STARlight provide the first detailed comparison at LHC energies of the various aspects of STARlight -from the nuclear photon fluxes to the implementation of the exclusivity condition -with experimental data over a very wide range in initial photon energies (up to 100 GeV). To vary the impact parameter range sampled by the dimuon events, events are classified according to their forward neutron topology, i.e. whether or not the neutrons are emitted in the forward or backward direction, or both, and corrected for the presence of EM pileup processes. Finally, cross sections for events with no forward neutron activity, which removes dissociative events, are presented as a function of dimuon acoplanarity. These data can be used to study the impact of QED showering on the final-state muons in a environment with very low backgrounds.

ATLAS detector
The ATLAS detector [39] at the LHC covers nearly the entire solid angle around the collision point. It consists of an inner tracking detector surrounded by a thin superconducting solenoid, electromagnetic and hadronic calorimeters, and a muon spectrometer incorporating three large superconducting toroidal magnets with eight coils each.
The inner detector (ID) is immersed in a 2 T axial magnetic field and provides charged-particle tracking in the range | | < 2.5. The ID is composed of three major subsystems. The high-granularity silicon pixel detector covers the vertex region and typically provides four measurement points per track. It is followed by the silicon microstrip tracker, which usually provides four two-dimensional measurement points per track. These silicon detectors are complemented by the transition radiation tracker, which covers the pseudorapidity 1 interval | | < 2 and thus enables radially extended track reconstruction in that range.
The calorimeter system covers the pseudorapidity range | | < 4.9. Within the region | | < 3.2, electromagnetic calorimetry is provided by barrel and endcap high-granularity lead/liquid-argon (LAr) electromagnetic calorimeters, with an additional thin LAr presampler covering | | < 1.8 to correct for energy loss in material upstream of the calorimeters. Hadronic calorimetry is provided by the steel/scintillator-tile calorimeter, segmented into three barrel structures within | | < 1.7, and two copper/LAr hadronic endcap calorimeters. The solid angle coverage is completed with forward copper/LAr and tungsten/LAr calorimeter modules optimized for electromagnetic and hadronic measurements respectively.
The muon spectrometer (MS) comprises separate trigger and high-precision tracking chambers measuring the deflection of muons in a magnetic field generated by the superconducting air-core toroids. The field integral of the toroids ranges between 2.0 and 6.0 T m across most of the detector. A set of precision chambers covers the region | | < 2.7 with three layers of monitored drift tubes, with the innermost chambers replaced by cathode strip chambers in the forward region (| | > 2), where the background is highest. The muon trigger system covers the range | | < 2.4 with resistive plate chambers in the barrel, and thin gap chambers in the endcap regions.
Muon reconstruction is first performed independently in the ID and MS. The information from the individual subdetectors is then combined to form the muon tracks that are used in the analysis [40].
A multi-level trigger system is used to select interesting events [41]. The first-level (L1) trigger is implemented in hardware and uses a subset of detector information to reduce the event rate to a design value of at most 75 kHz. This is followed by a software-based high-level trigger (HLT) which reduces the event rate to several kHz.
To reject hadronic interactions, this measurement makes use of the minimum-bias trigger scintillators (MBTS), which cover 2.08 < | | < 3.75 with two rings of scintillator counters positioned at = ±3.28 m from the ATLAS IP. The inner MBTS ring covers 2.78 < | | < 3.75 with eight slats, each subtending approximately 45 • in azimuth.
During heavy-ion running, two zero-degree calorimeters (ZDC) are installed at = ±140 m, upstream and downstream of the IP in both directions, each covering approximately | | > 8.3 for neutral particles, which are primarily neutrons and photons. Each detector side (labeled ZDC+ for > 8.3 and ZDC− for < −8.3) consists of four modules, each with approximately one interaction length of tungsten absorber surrounding layers of vertical quartz rods. Hadronic showers induce Cherenkov photons in the quartz which are transported to a photomultiplier tube about 35 cm above the top edge of the absorber. The ZDC calibration is performed in each set of four modules using photonuclear processes that deposit one or more neutrons on one side, used for triggering, and a single neutron, carrying the full per-nucleon beam energy, on the other. Time-dependent module weights are determined in short time intervals to minimize the 1 ATLAS uses a right-handed coordinate system with its origin at the nominal interaction point (IP) in the center of the detector and the -axis along the beam pipe. The -axis points from the IP to the centre of the LHC ring, and the -axis points upwards. Cylindrical coordinates ( , ) are used in the transverse plane, being the azimuthal angle around the -axis. The pseudorapidity is defined in terms of the polar angle as = − ln tan( /2). Angular distance is measured in units of variance around the nominal per-nucleon beam energy. Energy resolutions achieved are typically around Δ / ≈ 16%.

Trigger and luminosity
The data analyzed in this paper are from the 2015 Pb+Pb run at the LHC. The primary event trigger used in this work utilized nearly the full ATLAS detector system: the MS trigger, the full calorimeter system, the ID, and the MBTS.

Event trigger
To collect a large sample of exclusive dimuon events, the L1 trigger signature was a coincidence of a muon track in the MS, with no specific transverse momentum selection, for events with a maximum transverse energy of 50 GeV registered in the entire calorimeter system. The MBTS detectors were used in the HLT to reject events with more than one hit in the inner ring in either the forward or backward direction. Finally, the ID was used in the HLT to select events with a well-reconstructed track with T > 200 MeV. The efficiency of the HLT track trigger is estimated to be approximately 100%. The MBTS veto in the trigger could potentially lead to inefficiencies if the noise rate is sufficiently high. Using events triggered on empty bunch crossings, it is found that the fraction of events failing the MBTS selection is negligible. Thus, the efficiency of the trigger is driven primarily by the efficiency of the L1 MS trigger requirement for single muons. This was determined using minimum-bias data, and cross-checked using exclusive dimuon events, as described below in Section 7.2.
All selected events are also required to pass standard quality preselections. Events are required to be taken during periods where ATLAS is fully operational. Furthermore, events with identified problems associated with the ID, calorimeter, or data acquisition are removed. This only occurs in 0.15% of the total number of events, so no correction for this procedure is applied.

Luminosity and electromagnetic (EM) pileup
The primary trigger provided rejection against hadronic events using the MBTS requirement. Thus, it sampled the full luminosity of the 2015 heavy-ion run. The total integrated luminosity used for this analysis is found to be L int = 0. 48 Similarly, the mutual interaction probability per event based on rescaling using the RELDIS ratio is M = (0.174 ± 0.011) · 10 −2 . These quantities are used in Section 7.6 to correct the observed fractions of events with a particular forward neutron topology.

Theoretical background
Analytic calculations of dilepton rates, which are used as the basis for Monte Carlo (MC) generators, typically involve an integral over the factorized product of a two-photon flux (which acts as an effective luminosity) and the → + − cross section from perturbative QED. An event-by-event sampling of the scattering angle, as well as the transverse momenta of the scattering photons, is performed to generate a sample of events satisfying the fiducial selection.
Primarily following the discussion in Ref. [12], the two-photon flux is the joint distribution of the incoming photon energies 1 and 2 where the functions ( , ) are single-nucleus photon differential fluxes 3 / 2 , which depend on the magnitudes of ì 1 or ì 2 , the position vectors relative to each nucleus, with the origin at the nuclear barycenter. The function ( ) encodes the interaction probability as a function of the nucleusnucleus impact parameter , which is the magnitude of the impact parameter vector ì = ì 2 − ì 1 . This function is a product of different contributions, the two most important being the probability of having no hadronic interaction and the probability of a particular configuration of forward neutron production. For ultrarelativistic nuclei, the photon flux for each nucleus ( , ) can be approximated as To present the flux in a form that best matches the dilepton production cross section, the photon energies need to be presented in terms of their invariant mass ( ) and pair rapidity ( ), defined as: 1,2 = ( /2) exp(± ) and = (1/2) ln ( 1 / 2 ).
The observed rate of dimuon events is then the product of the two-photon flux and the dilepton cross section where ( → + − ) is the dimuon production cross section calculated from perturbative QED. Most UPC calculations use the lowest-order version of this: the Breit-Wheeler process [4]. The rapidities of the final-state muons also depend on the distribution of the center-of-mass scattering angle ( ★ ), calculated using lowest-order QED [45].
While the formalism shown here is shared between most theoretical groups, there are differences both in how the integrals are restricted to exclude hadronic events, and where dimuon pairs are allowed to form. This affects the integration limits over ì 1 and ì 2 , as well as the impact-parameter-dependent function ( ): • The photon flux from each nucleus always scales with 2 but its dependence on the radial coordinate ( ) has to be regulated to account for the extent of the nuclear charge distribution. In some calculations [16], this is done by convolving the point charge with a form factor derived from the measured nuclear densities. In others, the nuclear radius is imposed as a minimum radial distance for each nucleus, such that a dilepton pair cannot be produced inside it [12].
• At small impact parameters, hadronic processes occur simultaneously with a process and contaminate the exclusive process with additional hadrons. To avoid needing to disentangle these two processes, most calculations exclude the impact parameter range where the nuclei overlap. This is sometimes approximated as a step function imposing a condition on , the relative impact parameter between the two nuclei, > 1 + 2 , where 1 and 2 are the radii of the two nuclei. For realistic nuclear densities, the condition is imposed as the impact-parameter-dependent probability of not interacting, so ( ) ∝ 1 − H ( ), where H ( ) is the probability of a hadronic interaction at a given : where the nuclear thickness function (|ì|) is the integral over (the longitudinal coordinate along the beam axis) of the nuclear density ( ì , ). This function is essentially unity for impact parameters less than the sum of the radii of the colliding objects, e.g. ≤ 2 A for symmetric nuclei, and drops to zero rapidly above it.
• Until recently (e.g. Ref. [46,47]), most calculations did not account for the final-state radiation (FSR) from the outgoing charged leptons, and none include contributions from dissociative processes. This process, shown diagrammatically in Figure 1(c), also produces dileptons, but one of the incoming photons is radiated from a charged constituent of a nucleon, and this radiative process is associated with a momentum scale large enough to break up the nucleus. Both of these effects tend to populate the dimuon acoplanarity distribution in the region > 0.01. However, dissociative processes also induce nuclear breakup, allowing a systematic disentangling of these two processes in the experimental data by the comparison of different ZDC selections.
• The presence or absence of forward neutrons also affects the impact parameter dependence of the two-photon flux. STARlight implements this using fn ( ), an additional contribution to ( ). The function fn ( ) is the probability of a dissociative photonuclear interaction, where the precise form depends on the desired forward neutron topology. This can be nondissociative ("0 0 ", with no neutrons in either direction); single dissociative ("X 0 ", with neutrons in only one direction and not the other); or double dissociative ("X X ", with neutrons emitted in both directions). The primary quantity which determines fn ( ) is 1 X ( ), the rate of photodissociation of one nucleus ( → * ) as a function of the distance between the two nuclei: where → * ( ) are measured photonuclear cross sections [48], parameterized by , the incoming photon energy measured in the target rest frame. Using this, the probability of neutrons being emitted toward each ZDC arm is X ( ) = 1 − exp (− 1 X ( )). In STARlight the probability for exclusive single-arm X 0 events is 2 X (1 − X ); two-arm X X events occur with probablity 2 X ; and 0 0 events without neutrons have probability (1 − X ) 2 . For the beam energy and kinematic selections studied in this paper, typical probabilities (for example, from STARlight) for zero, one or two ZDC arms are 60%, 35% and 5% respectively.
• The transverse dimensions of the nuclear charge distributions are reflected in the transverse momenta of the incoming photons. There are different ways to implement this, but the most common approach is to provide a transverse momentum kick to each photon with a distribution determined primarily by the nuclear form factor. These are then summed vectorially to give the intrinsic photon-pair T that is imparted to the lepton pair. While STARlight does not implement an impact-parameter dependence of the photon-pair T , recent calculations suggest that this may be an observable quantity [49].
These differences are relevant to all cross-section calculations. However, it should be noted that while the most widely used framework for this, the STARlight model, has the most sophisticated handling of the nuclear photon fluxes as well as the nuclear overlap, it has no higher-order QED or nucleon dissociative processes.

Monte Carlo samples
Several different models are used to simulate signal and background processes, as well as to add FSR. Fully simulated samples are used to calculate experimental corrections using single muons, and to test their accuracy on simulated dimuon events. Fast simulations are used to study the impact of QED showering on the measured dimuon final state. Finally generator-level samples of dissociative pair production are used to estimate the background contributions.

STARlight
STARlight [50] implements calculations of both the nuclear photon flux (including the nuclear form factor) and the lowest-order QED cross sections for exclusive processes [51], such as the process → + − studied in this paper. Reference [12] provides the most complete presentation of the STARlight formalism, including the modifications of the dimuon spectrum when selecting on the number of ZDC arms (zero, one or two) with forward neutrons. Each dimuon event is fully determined by , , | cos ★ | and T, .

Simulated samples: single muons and STARlight
A large sample of single muons, distributed uniformly in the ranges | | < 3, − < < and 2 < T, < 50 GeV, is utilized to evaluate the muon reconstruction efficiencies R ( , T, ) as a function of (the product of the muon charge and pseudorapidity) and T, . To perform detailed comparisons of data with simulated exclusive dimuon events, a sample of STARlight 1.1 events is passed through a detector simulation [52] based on G 4 [53] and is reconstructed with the standard ATLAS reconstruction software.

QED showering with STARlight+P 8
STARlight does not implement final-state radiation (FSR), which is observed in the acoplanarity distributions in data as a notable tail that is absent in the STARlight simulations. STARlight only provides initial-state transverse momentum from the nuclear form factors. Higher-order contributions have been calculated in Ref. [46] in a Sudakov formalism and are found to qualitatively compare well with the preliminary version of this analysis. To include this effect, STARlight dimuon events were provided as input into P 8 (version 8.226 [54,55]) for QED showering, with the hard scale set to one muon's T (the standard choice in P 8 for a -channel process). By construction, the showering process conserves the overall cross section, but the FSR can lead to a fraction of the events migrating in and out (but primarily out) of the selected fiducial region used here. The sensitivity of the showering to the assumptions used in P 8 was tested by varying the available phase space for the radiation, relative to the selected hard scale, and no change in the final distributions was observed for the kinematic regions explored in this paper.

LP
The LP 4.0 generator [56] is used to calculate the dissociative process ★ → → + − , a well-known background in proton-proton collisions [57]. LP models single-dissociative proton-proton collisions using an EPA spectrum for the nondissociated nucleon, and the Suri-Yennie form factor [58] for the dissociated nucleon. This leads to the production of dimuons, many of which are emitted into the central region, as well as forward hadron production, although the hadrons are typically not emitted into the region covered by the ATLAS ID. Since the only observable from LP used in this analysis is the acoplanarity, and this is determined primarily by the much larger initial-state T associated with the dissociated proton, no attempt is made to reweight the initial photon spectra. However, comparisons of LP with STARlight+P 8 are made only after identical selections in and .
6 Event characterization

Event selection
The signal process PbPb( ) → + − (Pb (★) Pb (★) ) that can be treated using the EPA is only straightforward to calculate when the two nuclei do not interact hadronically. In this case, there is no potential background from hadronic processes and the full nuclear charge contributes to the two-photon flux. To suppress hadronic processes, the ATLAS ID is used to require no activity beyond that associated with the two oppositely charged muons. Of the triggered events passing all quality selections, events are selected with two "tight" muons [40]. Only a single event has a third tight muon, so this event is discarded as potential background. One of the two muons is required to match the L1 trigger muon that caused the event to be recorded. Next the single muon and muon pair kinematics are required to satisfy the fiducial selection criteria: • The two muons must have opposite charges.
• Each muon must have T, > 4 GeV, | | < 2.4, and a transverse impact parameter relative to the beamspot of less than 1.5 mm.  • The muon pair must have invariant mass > 10 GeV, and a transverse momentum of the dimuon system T, < 2 GeV.
The fiducial criteria remove about two-thirds of the triggered events containing two good muons. Of the remaining events, none are observed to have more than one reconstructed vertex, so no vertex requirement (for either a single vertex or multiple vertices) is imposed. Finally, events are rejected if there are any other well-reconstructed tracks in the acceptance with T > 100 MeV; this removes only 2% of the events. As a cross-check, a more restrictive condition was also considered -to reject events with any additional track of any track quality -but this was found to only reject 3% of the selected events. For the selected events, the distributions of track transverse impact parameter are consistent with simulated STARlight events, implying that there is no residual contribution from heavy-flavor production.
In fully simulated STARlight events, the requirement of having a reconstructed vertex leads to a selection inefficiency of approximately 5% while providing no additional background rejection. This was demonstrated by noting that after the vertex requirement, the fraction of events with an additional good track (2%) is found to be the same as for the event selection without it. Thus, the vertex requirement is considered only as a cross-check and is not part of the primary event selection.
After all fiducial selections, 12132 candidate events passing fiducial selections are utilized in the analysis.

ZDC event class definition
Each selected exclusive dimuon event can be classified according to its forward neutron topology, utilizing the two ZDCs. The ZDC energies are reconstructed by extracting the signal amplitude for each of the four modules in each direction and then applying a time-dependent weight to each module. These weights are determined so as to minimize the width of the energy distribution in the single-neutron region. The energy resolution is sufficient to clearly separate the one-neutron peak from the distribution at lower energies, which is mainly a combination of electronic noise and forward photons. The presence of ZDC activity on either side is defined as an energy greater than 40% of the single-neutron peak position. This is indicated in the individual energy distributions in the left panel of Figure 2. The correlation of energies in the right panel of Figure 2 illustrates the three primary topologies available for these events: 1) the most probable configuration is no activity in either ZDC ("0 0 "), 2) the next mostly likely configuration is observing one or more forward neutrons in one ZDC, and none in the other ("X 0 "), and 3) finally, the rarest configuration is observing one or more forward neutrons in both ZDC arms ("X X ").
Due to EM pileup in the LHC (as discussed above in Section 3.2) the neutrons detected in one or both arms of the ZDC are not necessarily associated with the observed dimuon pair. However, the rate of additional neutrons can be predicted when the single and mutual dissociation probabilities are known. Corrections associated with this are calculated below in Section 7.6.

Data analysis 7.1 Cross-section definitions
After the exclusive dimuon event selection, cross sections are derived by scaling the dimuon yields by the Pb+Pb collision luminosity (L int ). This is done after correcting the observed number of events ( ) for: 1) muon pair trigger and reconstruction efficiencies ( T and R , respectively), 2) bin migration at the edge of the T, > 4 GeV fiducial region ( mig ) and 3) backgrounds from dissociative processes ( dis ). This is represented by the following formula: where is one of the dimuon kinematic variables, described in Section 1. These include , , | cos ★ | and the initial energies of the photons, inferred by combining and .
The efficiencies and background corrections are applied as event-level weights, while mig is applied to the final spectra. The trigger efficiency, first discussed in Section 3.1 is derived using Pb+Pb data, using both minimum-bias hadronic events with reconstructed muons and UPC dimuon events. The muon reconstruction efficiency is based on simulation, but corrected by an -dependent scale factor (SF) based on data.

Trigger efficiencies
The single-muon L1 trigger efficiencies used for final corrections are derived using the minimum-bias data as a function of the product of the muon charge and pseudorapidity, , and the muon transverse momentum T, . The trigger efficiency does not depend on the overall event activity, which was checked by comparing different centrality intervals in the data. The efficiency values were also cross-checked by an independent tag-and-probe method using the selected exclusive dimuon events. In these events, only one of the muons is required to be associated with a L1 muon object. Thus, one can use each muon associated with a trigger as a tag, and then measure the probability of the other muon being emitted into an angular region, opposite in azimuth, that also has a L1 trigger muon. This gives an independent estimate of the trigger efficiency, which is found to compare well with the minimum-bias data. For application to the data, the single-trigger efficiencies are smoothed using a Fermi function.
The total event-level trigger efficiency is the probability that at least one of the muons led to the recording of the event: where + T and − T are the single-trigger efficiencies for positive and negative muons, respectively. The typical trigger efficiency is 93% at smaller (less than 20 GeV) and | | (less than 1), and increases to 97% at larger values of (greater than 40 GeV) and | | (greater than 1.5).
To assess the overall level of agreement between data and simulation after trigger efficiency corrections, the normalized pseudorapidity and azimuthal distributions of single muons in data and simulation are compared separately for both positive and negative charges, for all events in the fiducial acceptance. While the data are corrected for the trigger efficiency, neither data nor simulation are corrected for the reconstruction efficiency. This is to provide a low-level comparison of data and MC simulation, as it should also be noted that the data are not yet corrected for dissociative backgrounds. Figure 3 shows that the distributions in and are in good agreement for both charges, with the overall shape and some smaller-scale features of the data also seen in the simulations.

Reconstruction efficiencies and scale factors
The muon reconstruction efficiency is defined as the probability for a generated muon to be fully reconstructed and characterized as a "tight" muon [40]. It is determined as a function of the muon pseudorapidity times its charge ( ) using single-muon MC samples as described in Section 5.2, with a per-muon scale factor determined using UPC events in a tag-and-probe approach. In events with fewer than three tracks, a "tag" muon is selected using the single-muon trigger, and is required to satisfy tight selection criteria. Two types of "probes" can then be selected, track probes (reconstructed ID tracks with a transverse impact parameter | 0 | < 1.5 mm and passing selections identical to those for muons) and extrapolated MS track probes. Both are required to have a reconstructed charge opposite to the tag, to have an acoplanarity of less than 0.02, and a pair T < 2 GeV. For each ID track probe, one searches for a muon passing tight selections, which defines ( |ID). For each extrapolated MS track probe, one searches for the closest track passing the same spectrometer and impact parameter selections, and (ID|MS) is the matching probability. In each case, the two objects are considered matched if they are separated by a relative -distance of Δ < 0.1. The full reconstruction efficiency is the product ( |ID) (ID|MS), and this can be compared directly between data and STARlight simulations. The ratio of the full reconstruction efficiency in data to that in simulation is defined as the SF, which is a multiplicative correction to the simulated efficiency to account for differences between data and MC simulation that depend on emission angle. The SF is measured as a function of the muon pseudorapidity, and for each range in is found to be constant with T, . Thus, only one SF is extracted for each range. It is found to be within −1% of unity for | | > 1.5 and to decrease systematically toward | | = 0, with a maximum deviation of around −4% near = 0.
Since the muons are well separated in emission angle, the pair efficiency is the product of the two individual efficiencies: where each individual muon efficiency is multiplied by the appropriate SF. Due to the turn-on of the efficiency at T, = 4 GeV, the impact of correcting for the reconstruction efficiency is about 40-50% for < 20 GeV and | | < 0.8, decreasing to 15% at larger values.
The combined weight, including both the trigger and reconstruction efficiency corrections, is applied as a per-event weight = 1/( T R ).

Outflow from fiducial region
Migration between bins in the dimuon variables was checked and found to be negligible, due to the high precision of the ATLAS MS. Despite the excellent reconstruction performance of the MS, the fact that the fiducial acceptance region requires a two-dimensional cut in the joint transverse momentum distribution of both + and − leads to some net loss owing to momentum resolution. For events within the fiducial volume, but just at the boundary with both muon T, > 4 GeV, if either muon's transverse momentum fluctuates even slightly below 4 GeV, due to finite momentum resolution, then the event is rejected, while both muon transverse momenta would have to fluctuate upward to remain in the fiducial volume. This effect is exacerbated by the strong T correlations induced by the UPC initial state, with nearly vanishing total transverse momentum. This can be seen clearly by reconstructing simulated STARlight events, and applying the reconstruction efficiency weights event by event. While most of the spectrum in pair rapidity and pair mass is reconstructed correctly, within 1%, a systematic depletion is observed when one or both of the muon T, values approach 4 GeV. This depletion, quantified as a function of the different kinematic variables, is derived using the STARlight simulation, and is found to lead to corrections of 2-3%, but only for | | < 1.6 and < 20 GeV. The edges of the acceptance in | cos ★ | vs. are affected somewhat more, with typical corrections of about 10%.
Since these corrections are calculated as a function of the dimuon kinematic variables, they are applied to the multidimensional yields (e.g. vs. or | cos ★ | vs. ) before the final results are calculated.

Dissociative background corrections
After all corrections are applied, there remains an irreducible background from the presence of dissociative events. This process is expected to contribute at large values of acoplanarity ( > 0.01). However, one cannot simply select only events with small acoplanarity, since the large-acoplanarity region also has a non-negligible contribution from QED radiative effects. For each dimuon kinematic selection, a binned maximum-likelihood fitting procedure applied to the dimuon acoplanarity distribution is used to estimate the fraction of dissociative events. Due to the limited size of the data sample, it was split into coarse intervals in (three intervals with boundaries at 10, 20, 40 and 80 GeV) and absolute pair rapidity (three intervals with boundaries at 0, 0.8, 1.6 and 2.4).
The fit model utilizes two templates, where the signal distribution is calculated using STARlight + P 8, and the background is from LP . The PDF of the acoplanarity distribution, in a particular selection of and , ( , , ), is a sum of the signal and LP contributions, parameterized by the fraction dis : For the nominal determination of the correction factors, no weights are applied to the events during the fitting, since the corrections do not vary strongly within each individual dimuon kinematic interval and they are independent of . The fits are applied to distributions both for a specific forward neutron topology (0 0 , X 0 , and X X ) and for ZDC-inclusive data. Figure 4 shows the outputs of the fits for the selection with the largest number of events: 10 < < 20 GeV and | | < 0.8. Two notable features are observed. First, the 0 0 selection is very well described by the STARlight+P 8 distribution, with no significant contribution from LP 4.0 needed. This suggests that events with no ZDC activity (0 0 ) have no dissociation, and the tails in the acoplanarity distributions derive only from higher-order QED contributions. By contrast, the X 0 and X X selections have a much less steep falloff at larger acoplanarity values. Inclusion of the LP contribution gives a good description in the region beyond > 0.02. The overall contributions to the X 0 and X X integrals from the dissociative contribution are 7% and 12%, respectively. Finally, the distribution inclusive in forward neutron topology has an overall contribution of about 3% for this kinematic selection. Intervals with much higher pair mass or pair rapidity have fewer events that can be utilized in the fit, giving larger statistical uncertainties to the extracted dis . These uncertainties are propagated into the final results as systematic uncertainties.
The dissociative contributions come primarily from the process → + − where one photon comes from a lead ion (and thus has T ≈ 25 MeV), while the other is emitted by the charged constituents of an individual nucleon with a GeV-scale transverse momentum. Of the two, the nucleon has the smaller radius, so its photon spectrum is much harder, and thus a strong asymmetry is expected in the dimuon rapidity distribution, with more pairs emitted in the nucleon-going direction. One can test this expectation after assigning the pair rapidity a sign based on the direction of the ZDC signal in X 0 events, as one expects to observe nuclear dissociation in the direction of the dissociated nucleon. The fractions of dissociative events are shown as a function of signed in Figure 5. It is observed that the extracted dissociative contribution increases with pair rapidity in the direction of the dissociated nucleus, consistent with being induced by a dissociated nucleon.

Forward neutron topology fractions
In Section 4, it is shown that the ZDC selection modifies the impact parameter profile ( ) found in the integral for the two-photon luminosity. Due to the increasing photon flux at smaller radial distances from the nuclear center ( ), events with a ZDC signal (either X 0 or X X ) typically have smaller impact parameters than for 0 0 . The smaller impact parameters lead to interactions involving higher photon energies from one or both nuclei. This tends to harden the observed diphoton (and thus dimuon) invariant mass spectrum for a fixed interval in pair rapidity. Conversely, for a fixed interval in dimuon mass, the rapidity distribution gets narrower. To measure this in each dimuon kinematic selection, the fractions of events with X 0 and X X ( X 0 and X X , respectively) are estimated with a single fit over all three ZDC topologies that includes the contribution from dissociative background, as discussed previously. Once the measured ZDC fractions are determined, they must be corrected for EM pileup processes and, in principle, for any ZDC inefficiency. However, while the ZDC is only 4.6 interaction lengths deep, which implies an approximately 2% inefficiency for single neutrons, more than half of the EMD events on a single side produce more than one neutron. This leads to an estimate of an overall efficiency of over 99%. Accounting for this sub-percent inefficiency is addressed in the systematic uncertainties. The combined fit model is formulated as a set of three equations for each kinematic region, with a total number of events in each ZDC category ( , where is an index corresponding to 0 0 , X 0 , or X X ) given by with parameters EPA , the overall yield of dimuons from the signal distributions, Z the measured fraction of events for each ZDC class (i.e. = 0 0 , X 0 , and X X ), and dis the background contribution for each ZDC selection. The fractions are subject to the constraint 0 0 + X 0 + X X = 1 since all events have to fall into one of the categories. The simultaneous fit incorporates the statistical correlations between the samples. The choice of using the same template PDFs EPA and dis for all three ZDC selections is motivated by direct comparisons of the acoplanarity distributions, and validated by the good quality of the fits.
EM pileup processes (as discussed in Section 3.2) generate neutrons that are detected in one or both arms of the ZDC but are not associated with the scattering process that generated the dimuon pair. This leads to an outflow from 0 0 and X 0 events to both the X 0 and X X events. The two fundamental parameters are S and M , the probabilities for exclusive single (for each arm) and mutual (two arm) dissociation per bunch crossing. The probability of pileup preserving a 0 0 configuration is (1 − S ) (1 − S ) (1 − M ), i.e. the joint probability of having no single pileup in each arm, as well as no mutual process occurring. Conversely, the probability of pileup turning 0 0 into X X is M + 2 S , i.e. the sum of probabilities of mutual breakup or two independent single breakup processes, one in each arm. The probability turning 0 0 into X 0 is just the complement of the other two quantities. Similarly, the probability of preserving X 0 is (1 − S ) (1 − M ), since one arm has already seen activity, and its complement is the probability for converting X 0 to X X . Finally, X X is unaffected by the presence of any pileup. Putting these together, one gets the migration matrix: Here the fractions are the measured ones extracted from the combined fits, and the fractions are the corrected values reported in the results. The terms which go as 2 S reflect the possibility that a pileup event can be composed of two separate pileup interactions in addition to the primary physics event. While not identical, as assumed by STARlight, 2 S is of similar order to the M estimated using measured cross sections, as discussed in Section 3.2. This matrix is inverted to transform the measured ZDC fractions into the corrected ones. The statistical uncertainty of the corrected values includes the full fit covariance information.

Unfolding of the acoplanarity distribution
In order to facilitate comparisons of QED calculations with the experimental acoplanarity ( ) distribution, a bin-by-bin unfolding procedure is applied to the measured values. To maximize the statistical precision, the entire fiducial region with > 10 GeV is included. The response used in the unfolding is derived from a large sample of events produced with STARlight + P 8, with each muon's azimuthal angle smeared in accord with parameterizations of the full G 4 simulations. The unfolding corrections range from an upward correction of about 3% for < 0.0001 to a downward correction of a bit less than 4% around = 0.005. Around = 0.01, the corrections converge rapidly to zero, as the underlying spectrum becomes flatter. The net effect is to slightly narrow the measured peak near = 0, which is smeared slightly by detector effects. The corrections were also compared with those for a sample in which the simulated spectrum is reweighted to agree with data, a procedure which is limited by the size of the data sample, and while they are generally similar to each other, they disagree in some intervals by up to 2%, which sets the overall scale for the systematic uncertainty. The corrections were checked against a fully simulated sample of STARlight events and the corrections agree to better than 0.4% for < 0.01. To eliminate the contribution from dissociative events, only the 0 0 selection is used, and the results are normalized to the pileup-corrected 0 0 cross section. This is done by scaling the measured total cross section (reported in the next section) by the 0 0 fraction 0 0 = (72.1 ± 0.6)% after correcting the measured fraction of 64% for the time-averaged EM pileup.

Systematic uncertainties 8.1 Trigger efficiency
The trigger selection is based on four contributions: the L1 muon trigger, the rejection of events with total transverse energy greater than 50 GeV, the MBTS veto in the HLT, and the track selection in the HLT. The muon trigger efficiency has been measured directly in peripheral minimum-bias data and is parameterized by a Fermi function. The assigned systematic error is the statistical uncertainty of the fits, propagated using 200 replicas, and is about 0.2%. The tag-and-probe results are consistent with the minimum-bias results, so either could be used in principle, but the minimum-bias data have better statistical precision and are thus preferred. The transverse energy veto has a negligible inefficiency, as there are no selected dimuon events with a L1 transverse energy above 15 GeV, far from the veto threshold. The HLT track selection is determined to have 100% efficiency, with an uncertainty of 0.4%, using a support trigger based on the MS accompanied by a total transverse energy requirement of T < 50 GeV. Finally, the false-positive rate for the MBTS veto was checked using empty events, where neither LHC ring had filled bunches colliding in a particular bunch crossing. In these events, the rate of false positives is at the 10 −5 level, indicating that this contribution is negligible.

Reconstruction scale factors
The uncertainty in the reconstruction SFs has two contributions. The statistical uncertainties of the SF are propagated by constructing 200 replicas, with the SF in each interval varied according to its statistical uncertainty, and the variance in the yields determines the uncertainty. The sensitivity to the details of the procedure was checked by changing the selection criteria for the events used in the SF calculation, as well as the matching between the ID and MS. Given that the dimuon selection is quite pure, these changes have only a small impact, and a constant 1% systematic uncertainty is assigned to cover this. This uncertainty also includes a small deviation of the reconstruction efficiency determined in the simulation using tag-and-probe from the one extracted directly from the simulated muons for 2.0 < | | < 2.4.

Muon momentum scale and resolution
The muon momentum scale and momentum resolution are studied using a set of weights that vary the transverse momentum of simulated muons to determine a set of modified corrections that encapsulate the known uncertainties in muon performance in collisions [40]. The extracted corrections have a negligible impact on this analysis, and so are neglected. The impact of possible bias on the momentum measurement due to relative detector misalignment between data and MC simulation is studied by comparing the transverse momentum of each muon as measured in the ID with that measured in the MS. The differences observed in data and MC simulation are found to agree to better than 1% and no specific uncertainty is assigned to this contribution.

Fiducial acceptance definition
The correction factors determined using the fully reconstructed STARlight events account for the net outflow of one or both reconstructed muons from the T, > 4 GeV fiducial acceptance. This correction is up to 3% for < 20 GeV and decreases with increasing pair mass and pair rapidity. This is sensitive to the photon spectral shape. Thus, the uncertainty in this correction was estimated by reweighting the MC spectrum to resemble the data, as a function of 1 and 2 , the energies of the incoming photons. The difference between these two sets of corrections is at the 0.1% level, and no systematic uncertainty is assigned.

Dissociative background
The background from dissociative processes is estimated directly from data. The uncertainties in this derive both from the fit statistical uncertainties, and from the sensitivity to variations of the fitting procedure. As a cross-check, the fits are also performed using an unbinned fitting procedure, instead of the binned procedure used by default, and no significant differences are found. In addition, several other variations are performed: 1. Signal templates are allowed to be rescaled, in the fitting procedure, to better agree with the distributions near = 0.
2. The measured distribution for the 0 0 selection is used as the signal distribution for X 0 and X X selections. This is justified by the observed dissociative fraction being consistent with zero, e.g. as shown in Figure 4.
3. Applying reconstruction and trigger efficiency corrections before the fit.
The impact of these changes is quite small, since the dissociative tail has a distinct shape, especially when requiring forward neutrons in either ZDC arm. The differences turn out to be about 1% (absolute), constant in and . A 1% uncertainty is applied to the yields to cover this.
The statistical contributions to the background correction are propagated by keeping track of the sum of background correction weights, i.e. (Δ(1 − dis )) 2 = (Δ dis ) 2 , for each kinematic selection. The final uncertainty associated with this is then √︁ (Δ dis ) 2 .

Luminosity
The uncertainty in the integrated luminosity for the 2015 Pb+Pb run at √ NN = 5.02 TeV is 1.5%. It is derived from the calibration of the luminosity scale using -beam-separation scans, following a methodology similar to that detailed in Ref. [42], and using the LUCID-2 detector for the baseline luminosity measurements [59].

ZDC neutron topology fractions
The analysis of forward neutron topology fractions involves ratios of selected distributions, classifying them by the presence or absence of a signal in one or both ZDCs. As neither the event trigger nor the muon reconstruction is affected by the ZDC event class (0 0 , X 0 or X X ), all of these uncertainties cancel out. The systematic uncertainties related to the fitting procedure itself are the same as in the previous subsection, and are at most 1% on an absolute scale. The impact of EM pileup on the extracted ZDC fractions depends on our knowledge of the measured EMD cross sections, which have not been measured for √ NN = 5.02 TeV, but only for √ NN = 2.76 TeV. Thus, the cross sections used to extrapolate to the

Unfolded acoplanarity distributions
The primary systematic uncertainties applicable to the unfolded differential cross sections are those that apply to all cross sections, but now integrated over the full fiducial region > 10 GeV. Besides the luminosity uncertainty (1.5%), the 0 0 cross section has an uncertainty of 0.8%, dominated by the extrapolation of the EM pileup expectations to 5.02 TeV. Finally, there is an overall 2% uncertainty, constant in , which contains the variations testing the sensitivity to the input spectrum (nominal compared with reweighted) as well as the resolution model (full simulation compared with fast simulation). Combining these uncertainties in quadrature, the total systematic uncertainty is estimated to be 2.6%, which is applied uniformly to the unfolded acoplanarity distribution.

Results
This section presents comparisons of data with STARlight 2.0, which implements all of the primary physics mechanisms, except dissociative processes and QED FSR. All of the systematic uncertainties presented in Section 8 are combined in quadrature, and are typically about 2.2%. The dissociative processes are explicitly corrected for by the fitting procedure described in Section 7.5, but the FSR is necessarily included in the signal yields. It should be noted that only 1.5% of the data are removed by the T, < 2 GeV selection, which was chosen to limit the dissociative contributions. The same T selection only removes a small fraction of STARlight events even including FSR effects. However, the FSR separately has the effect of reducing the STARlight cross sections by about 4%, due to the reduction in T of one or both muons from the extra radiation. This effect is not included in the nominal STARlight comparisons, except in comparisons of data with the full distribution.
After integrating over the full fiducial phase space ( T, > 4 GeV, | | < 2.4, > 10 GeV, T, < 2 GeV), the measured total cross section is fid = 34.1±0.3(stat.) ±0.7(syst.) b. This should be compared with 32.1 b for STARlightand 30.8 b for STARlight + P 8. Figure 6 shows the absolute dimuon rapidity distributions for three mass intervals (10 < < 20 GeV, 20 < < 40 GeV, 40 < < 80 GeV). It is observed that while the measured cross sections are consistent with predictions near | | = 0, there is an excess in the data which increases monotonically with increasing | |. There are too few events in the higher bins to allow any strong statements to be made about whether the increase remains the same as in the lowest mass bin.
The mass distributions are presented in three pair rapidity bins (| | < 0.8, 0.8 < | | < 1.6, and 1.6 < | | < 2.4) in Figure 7. It is observed that the overall shape of the spectra is well described out to the highest masses in the available event sample (which limits the mass range for higher values of | |). The invariant mass interval extending from 100 to 200 GeV contains only four events, including the highest-mass candidate with = 173 GeV. While the pair mass and rapidity distributions are clearly important for characterizing the interaction between the two photons, the scattering angle in the initial center-of-mass system is particularly important for assessing the compatibility of the data with QED. Without the fiducial selection, these distributions would have a steep rise at large | cos ★ |. However, the restriction on muon pseudorapidity (| | < 2.4) strongly cuts into the acceptance for | cos ★ | near unity, i.e. at small scattering angles. Since the differential cross section calculated in leading-order QED is only a function of [45], the disagreement already observed as an enhancement of the cross section at forward | | makes comparisons of data and STARlight nontrivial. If one compares the differential cross sections / | cos ★ | for the same selections as in Figure 6, one arrives at the distributions shown in Figure 8. STARlight calculations show that events with smaller values of | cos ★ | also span a large range in | |, while these large values of | cos ★ | are only accessible in a more limited range, close to | | = 0. The cross sections are recalculated for | | < 0.8, the region where data and STARlight agree well. With this selection, there is a better description of the | cos ★ | distributions in the lowest mass range, as is observed in Figure 9. The systematic enhancement at smaller | cos ★ | can thus be attributed to the residual disagreement in the pair rapidity distributions.
Since the observed final state consists of just the final-state muons, and two-photon invariant mass and rapidity are conserved in the final-state dimuon system, the muon kinematics can be used to estimate the initial photon energies, 1 and 2 (as discussed in Section 4): 1,2 = (1/2) exp(± ). No corrections are applied for the presence of a soft final-state photon, but the impact on the extracted mass and rapidity is generally found to be quite small. Since the two photons are emitted independently, each event can be characterized by the maximum and minimum photon energies max and min , where max is the larger of the two photon energies. The corrected differential cross sections in max and min are presented in Figure 10, for both data and STARlight in the top panel, and their ratio is provided in the bottom panel.
It is observed that the ratio is near unity for photon energies (both minimum and maximum) around 10-20 GeV, but is significantly higher at both lower and higher energies. The disagreement for 75 < max < 100 GeV is approximately 40%, while an enhancement of about 15% is observed for min < 2 GeV. In the region between 5 and 20 GeV, the distributions for max and min are also observed to overlap, within uncertainties, suggesting that the distributions factorize (as represented in Eq. (1), i.e. it does not matter whether a 10 GeV photon has the larger or smaller energy). They also have the smallest difference relative to STARlight, with an excess of only around 5%. Combined, these suggest that a systematic modification of the initial energy spectrum is necessary to explain the differences between the data and STARlight calculations.
One possible avenue to explore is whether the integration limits in Eq. (1) could be relaxed without violating any important physics constraints. In particular, the authors of Ref. [12] note in their derivation that this restriction against production inside the two nuclei is to avoid the produced particles scattering off the nuclei. Dileptons have a very small probability of interacting with the nuclei, and this is even less likely within the nuclear skin. Thus, this condition could be weakened somewhat, so long as 1 − H (as defined in Section 4 ) is not too small. If the integration limits are decreased slightly, by one nuclear skin depth, it is observed that, relative to STARlight, the distribution of min is enhanced below 2 GeV, and the distribution of max is enhanced above 20 GeV.
All of the results shown so far are inclusive of the different ZDC topologies, which indicate either no activity (0 0 ), activity in either the forward or backward side (X 0 ), or activity on both sides (X X ). Events with smaller impact parameters, where the nuclei are closer together, are more likely to be accompanied by neutron dissociation in one or both arms and to have photons with higher energies. The procedure to extract the ZDC event class fractions is discussed in detail in Section 7.6, and is based on a simultaneous fit to the acoplanarity distributions for all three ZDC selections, assuming that all events arise from partitioning the original selection of signal events, along with backgrounds from dissociative processes that can be different for each forward neutron topology (X 0 and X X ). The EM-pileup-corrected results are shown in Figure 11, which displays X 0 and X X , the fraction of events with X 0 and X X , as functions of and | |. It should be noted that the two sets of results are different representations of the same data, but it is useful to see them plotted separately vs. and | |. The uncorrected data are also shown, to indicate the size of the EM pileup corrections. The data compared with STARlight suggest overall agreement, but STARlight generally tends to predict too large a fraction of events with forward neutrons, which is consistent with previously reported ALICE data at 2.76 TeV [60]. The unfolded differential cross sections / for the 0 0 topology are shown in Figure 12. The measured cross sections are compared with predictions of both generator-level STARlight and generatorlevel STARlight + P 8, both for an inclusive ZDC selection, but scaled by the same 0 0 fraction as observed in data. The STARlight distributions have no dependence on impact parameter, and thus no dependence on the ZDC selection. It is observed that the shape of the spectrum at large > 0.01, which reflects only the QED showering after applying the 0 0 condition, agrees well with STARlight + P 8. However, a difference in shape is observed for smaller values < 0.01. This could be explained by a small change in the T spectrum assumed by STARlight, which controls the width of the distribution.    Figure 12: (a) Fully corrected differential cross sections / for 0 0 -selected data. Data are compared with absolute cross sections from STARlight with, and without, P 8 QED showering. Statistical uncertainties are shown as error bars. (b) Ratios of STARlight + P 8 cross sections (black circles) and STARlight cross sections (magenta circles) to the data. The STARlight ratios do not extend beyond = 0.01 due to the absence of higher-order QED effects. The blue band around unity indicates the overall systematic uncertainty, while the gray bands around the data points reflect the uncertainties associated with the bin-by-bin unfolding.

Conclusion
The ATLAS experiment at the LHC has performed a measurement of cross sections for exclusive dimuon production in the process PbPb( ) → + − (Pb (★) Pb (★) ) using 0.48 nb −1 of Pb+Pb collision data taken at √ NN = 5.02 TeV. This reaction directly probes the quasi-real photon fluxes surrounding the highly boosted nuclei. The cross section of dimuons is corrected for detector effects, as well as for backgrounds from dissociative processes. The acoplanarity distributions are used to subtract dissociative backgrounds using fits incorporating contributions from signal and background. The signal is modeled by STARlight 2.0 with QED showering using P 8, and backgrounds are modeled using the single dissociative process from LP 4.0. After all corrections, differential cross sections in a fiducial acceptance ( T, > 4 GeV, | | < 2.4, > 10 GeV, T, < 2 GeV) are presented as a function of , | |, | cos ★ |, max and min , and compared with STARlight 2.0 calculations. Generally, good agreement is found but some systematic differences are seen, which may be explained by deficiencies in the modeling of the incoming photon flux. In particular, allowing dilepton pairs to be produced deeper within the nuclear skin may be sufficient to explain the observed differences, something which could be addressed systematically within the currently available models. Progress in modeling this process, using the data presented here, will be important in reducing uncertainties in the photon fluxes. These reduced uncertainties will be needed for precision studies of QED and QCD in nuclear collisions, as well as to probe physics beyond the Standard Model, both at the LHC, especially with the increased luminosity expected, and at future machines.

Acknowledgments
We thank CERN for the very successful operation of the LHC, as well as the support staff from our institutions without whom ATLAS could not be operated efficiently.