Observation of forward neutron multiplicity dependence of dimuon acoplanarity in ultraperipheral Pb-Pb collisions at $\sqrt{s_\mathrm{NN}} =$ 5.02 TeV

The first measurement of the dependence of $\gamma\gamma$ $\to$ $\mu^{+}\mu^{-}$ production on the multiplicity of neutrons emitted very close to the beam direction in ultraperipheral heavy ion collisions is reported. Data for lead-lead interactions at $\sqrt{s_\mathrm{NN}} =$ 5.02 TeV, with an integrated luminosity of approximately 1.5 nb$^{-1}$, were collected using the CMS detector at the LHC. The azimuthal correlations between the two muons in the invariant mass region 8 $\lt$ $m_{\mu\mu}$ $\lt$ 60 GeV are extracted for events including 0, 1, or at least 2 neutrons detected in the forward pseudorapidity range $|\eta|$ $\gt$ 8.3. The back-to-back correlation structure from leading-order photon-photon scattering is found to be significantly broader for events with a larger number of emitted neutrons from each nucleus, corresponding to interactions with a smaller impact parameter. This observation provides a data-driven demonstration that the average transverse momentum of photons emitted from relativistic heavy ions has an impact parameter dependence. These results provide new constraints on models of photon-induced interactions in ultraperipheral collisions. They also provide a baseline to search for possible final-state effects on lepton pairs caused by traversing a quark-gluon plasma produced in hadronic heavy ion collisions.

The momentum of emitted quasireal photons is predominantly along the beam direction and the transverse momentum (p T ) is small, typically less than 30 MeV [5,6].Therefore, the lepton pairs produced from leading-order photon-photon scattering (γγ → + − ) possess small pair p T and are nearly back-to-back in the azimuthal angle (φ).Recently, photon-photon [26,27] and photon-nucleus [28, 29] processes have been observed at very low p T in hadronic (b < 2R A ) heavy ion collisions.Interestingly, a broadening of lepton pair azimuthal angle correlations (or, equivalently, an increase of lepton pair p T ) is observed in hadronic collisions compared to that from UPCs [26,27].In hadronic events, a deconfined state of partonic matter, known as the quark-gluon plasma (QGP), can be formed.Therefore, final-state EM modifications of lepton pairs inside a QGP medium have been proposed as possible interpretations of the broadening effect [26,27,30].The initial p T of the lepton pairs depends on the overlap integral of the photon fluxes produced by the two nuclei, and as a result, the average pair p T ( p T ) could depend on the b between the two colliding ions.Although models of the flux of photons integrated over a given b range have large uncertainties [7,31,32], a QED calculation [32] predicts larger p T for smaller b values.Such a larger p T in the initial state would broaden the pair angular correlation, which could explain the effects observed in more central hadronic collisions.
To disentangle possible contributions from initial-and final-state effects to the modifications observed in hadronic heavy ion collisions, an experimental handle on the b dependence of lepton pair production in UPCs is essential.The photon-photon interactions can occur in conjunction with the excitation of one or both of the ions via photon absorption into giant dipole resonances or higher excited states [5-8, 33, 34].The giant dipole resonances typically decay by emitting a single neutron, while higher excited states may emit two or more neutrons.These forward neutrons have very low relative momentum with respect to their parent ions, and therefore approximately retain the beam rapidity.The contribution of higher excitations becomes larger as b gets smaller [5][6][7][8].Therefore, the number of emitted neutrons detected in the forward region can be used to classify UPC events into different b ranges.This Letter reports the first measurement of the forward neutron multiplicity dependence of γγ → µ + µ − production in the muon pair invariant mass region 8 < m µµ < 60 GeV in lead-lead (Pb-Pb) UPCs at a nucleon-nucleon center-of-mass energy √ s NN = 5.02 TeV, using data collected with the CMS detector during the 2018 LHC run.The Pb-Pb sample that includes information about forward neutrons corresponds to an integrated luminosity of approximately 1.5 nb −1 .Azimuthal correlations of muon pairs, quantified by the acoplanarity, α = 1 − |φ + − φ − |/π, are presented for several different classes of neutron multiplicity detected in the forward pseudorapidity range |η| > 8.3.Here, φ ± represent the azimuthal angle of each muon in the lab frame.A larger average α for lepton pairs from leading-order γγ scatterings corresponds to fewer back-to-back azimuthal correlations, and thus larger initial p T of the interacting photons.The muon azimuthal angle is used instead of p T because of its superior experimental resolution.The average invariant mass of muon pairs in various neutron multiplicity classes is also presented as a probe of the initial photon energy and its b dependence.Tabulated results are provided in the HEPData record for this analysis [35].
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, there are four subdetectors, including a silicon pixel and strip tracker detector, a lead-tungstate crystal electromagnetic calorimeter, and a brass and scintillator hadron calorimeter, each composed of a barrel and two end cap sections.Muons are detected in the range |η| < 2.4 in gas-ionization detectors embedded in the steel flux-return yoke outside the solenoid, with detection planes made using three technologies: drift tubes, cathode strip chambers, and resistive-plate chambers.Matching muons to tracks measured in the silicon tracker leads to a relative p T resolution around 1% [36] and an azimuthal angle resolution better than 7 × 10 −4 rad for a typical muon in this analysis.The CMS experiment has extensive forward calorimetry, including two steel and quartz-fiber Cherenkov hadron forward (HF) calorimeters that cover the range of 2.9 < |η| < 5.2, which are used to reject hadronic Pb-Pb collision events.Two zero degree calorimeters (ZDC) [37], made of quartz fibers and plates embedded in tungsten absorbers, are used to detect neutrons from nuclear dissociation events in the range |η| > 8.3.A 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. [38].
Events used in this study were selected online using a hardware-based trigger system that requires at least one muon candidate coincident with a Pb-Pb bunch crossing [39].On the trigger level, there is no explicit selection on the minimum muon p T and events with an energy deposit above the noise threshold in both HF calorimeters are vetoed.For the off-line analysis, events have to pass a set of selection criteria designed to reject beam-related background processes (beam-gas collisions and beam scraping events) and hadronic collisions.Events are required to have a primary interaction vertex, formed by two or more tracks, within 20 cm from the CMS detector center along the beam axis.The cluster shapes in the pixel detector must be compatible with those expected from particles produced by a Pb-Pb collision [40].To suppress hadronic Pb-Pb collisions, the largest energy deposits in the HF calorimeters are required to be below 7.3 and 7.6 GeV in the positive and negative rapidity sides, respectively, where these noise thresholds are determined from empty bunch crossing events.In addition, events must contain exactly two muon candidates and no additional tracks in the range |η| < 2.4.Selected events are then classified by neutron multiplicity, which is determined by the energies deposited in the ZDCs.For single neutrons, the relative energy resolution of the ZDCs is ∼22%-26%, while the detection efficiency is close to 100% in simulated events [37].Based on neutron peaks observed in the total ZDC energy distribution (see Appendix A), events are divided into three neutron multiplicity classes (0n, 1n, and Xn with X ≥ 2) on each side.The corresponding purities of selected neutron multiplicity classes are estimated by a multi-Gaussian function fit to the energy distribution.The purities are nearly 100% for the 0n and Xn classes, but only ∼93%-95% for the 1n class because of detector resolution effects.From the combinations of the number of neutrons in each ZDC separately, a total of six neutron multiplicity classes, labeled as 0n0n, 0n1n, 0nXn, 1n1n, 1nXn, and XnXn, are used in this study.The 0n0n class corresponds to no Coulomb breakup of either nucleus and the 1nXn class corresponds to one neutron emitted from one nucleus and at least two neutrons emitted from the other nucleus.
Muons are selected in the kinematic range of p µ T > 3.5 GeV and |η µ | < 2.4.They are reconstructed using the combined information of the tracker and muon detectors (so-called "soft muons" defined in Ref. [36]).The opposite-sign distribution (signal and background) is reconstructed by combining µ + and µ − candidates, while the combinatorial background is estimated using events containing same-sign muons.One of the muon candidates in the oppositeor same-sign pair is required to match a trigger muon.The studied dimuon kinematic range is 8 < m µµ < 60 GeV and |y µµ | < 2.4 to ensure high efficiency and also to suppress the contribu-tion from photoproduced resonances (charmonia and Z bosons).
The detector reconstruction efficiency is estimated using a dedicated γγ → µ + µ − Monte Carlo simulation sample produced by the STARLIGHT (v3.0) event generator [43] without restriction on the Coulomb breakup of either nucleus.Only + − pairs from the leading-order γγ scattering are generated, and the calculation is performed by integrating over the entire b space for UPC events.No differential b dependence of the initial photon p T is considered in STARLIGHT.The CMS detector response is simulated further using GEANT4 with these STARLIGHT generated events [44].The muon trigger (ε µ trig ) and reconstruction (ε µ reco ) efficiencies are estimated as functions of p µ T , η µ , and φ µ .To correct for detector inefficiencies, each muon pair event is scaled by (ε trig ε reco ) −1 , where The cross section of single electromagnetic dissociation (EMD) [45,46] of Pb nuclei in Pb-Pb collisions was measured to be 187.4± 0.2 (stat) +13.2 −11.2 (syst) b at 76 TeV [47].It is expected to be even larger at √ s NN = 5.02 TeV given the stronger EM fields.Because of the large single-EMD cross section, a single measured γγ → µ + µ − event may contain concurrent EMD Pb-Pb events in the same bunch crossing.These concurrent events can emit neutrons and migrate the neutron multiplicity of a single γγ → µ + µ − interaction to higher values.This EMD pileup effect is quantified by measuring the ZDC energy distributions from "zero-bias" triggered events that require only the presence of both beams in the same bunch crossing.No valid collision vertex or track is allowed to be present in the event.The same HF veto thresholds as for the γγ → µ + µ − events are applied.The neutron multiplicity classes in these selected zero-bias events are used to estimate the probability of a γγ → µ + µ − event being assigned an incorrect neutron multiplicity because of pileup effects.By inverting a matrix of these migration probabilities, the true observable distributions are extracted from the measured data.In this study, about 11% of measured γγ → µ + µ − events have neutron multiplicity migration caused by EMD pileup.
Figure 1 shows the corrected α distributions of µ + µ − pairs in Pb-Pb collisions within the kinematic range (p µ T > 3.5 GeV, |η µ | < 2.4, and |y µµ | < 2.4) for different neutron multiplicity classes.The α distributions are normalized to unit integral over their measured range [(1/N s )dN s /dα, where N s represents the signal yield].Each α spectrum is characterized by a narrow core close to zero and a long tail.The core component mostly originates from the leading-order γγ scattering, while in the tail component, higher-order γγ processes dominate.These higher-order processes include, e.g., extra photon radiation from the produced lepton(s), multiple-photon interactions, or scattering of (one or both) photons emitted from one of the protons inside the nucleus [9,30].The tail contribution in the XnXn class is larger than that in the 0n0n class.This is consistent with the expectation of larger contributions of higher-order γγ processes in UPC events that have smaller b and produce more neutrons in the forward region.
To investigate a possible b dependence of the initial photon p T , the core contribution to the α distribution is decoupled from the tail contribution using a two-component empirical fit function (where c i and t i are the fit parameters), as shown in Fig. 1, core : c 1 e −α/c 2 +c 3 α 0.25 , tail :  1).The vertical lines on data points depict the statistical uncertainties, while the systematic uncertainties and horizontal bin widths are shown as gray boxes.
except for the case of 1n1n, where a simple exponential function is used for the tail component, given the limited number of events.The core component is largely modeled by an exponential function with a correction term (c 3 ) to account for the small depletion in the very small α (e.g., < 5 × 10 −4 ) region, which tends to become more evident as the neutron multiplicity increases.This core functional form is validated by the STARLIGHT event generator and leading-order QED calculations, resulting in a <0.3% discrepancy on the average acoplanarity from the fit and theoretical predictions.A binned χ 2 goodness-of-fit minimization is performed using the integral of the function across each bin to account for the finite binning effect of the histogram.The average acoplanarity of µ + µ − pairs from the core component ( α core ) is then calculated using the fit function.
The measured α distribution and α core of µ + µ − pairs have several sources of systematic uncertainty arising from the contamination of hadronic collisions, EMD pileup correction, neutron multiplicity classification, and fit procedure.The uncertainty of the hadronic contamination is estimated by removing the requirement that selected events only contain two muons and is found to be <1.1%.To estimate the systematic uncertainty associated with the HF noise threshold, the threshold to define the hadronic contamination is tightened to 5 GeV for both UPCs and zero-bias triggered events.The difference from the nominal result is quoted as the systematic uncertainty and contributes <2.7%.The uncertainty arising from impure 1n class selection (<0.7%) is estimated by subtracting the contributions of 2n events selected with tight energy requirements, according to the 2n contamination probability.The systematic uncertainty associated with contamination of photoproduced Υ mesons (∼0.6%) is estimated by comparing α distributions from STARLIGHT between pure γγ → µ + µ − and γγ → µ + µ − mixed with photoproduced coherent Υ(1S), with the relative yield ratio of Υ(1S) over γγ → µ + µ − estimated by fitting the invariant mass distribution.The systematic uncertainty in α core associated with the binned χ 2 fit procedure is estimated by varying the bin width of α distributions, and is found to be less than 4%.The total systematic uncertainties are derived from a quadratic sum of all systematic sources and are found to be at most 5.1% in α core .To measure m µµ , a second-order polynomial function is fit to the mass spectrum (see Appendix A), to interpolate the contribution of γγ scattering to dimuon pair production over the Υ mass region.The systematic uncertainty related to this procedure is estimated by comparing the nominal result to the one obtained by a third-order polynomial function fit.Together with the aforementioned systematic sources, the total systematic uncertainty in m µµ is below 1.8%, across all neutron multiplicity classes.The vertical lines on data points depict the statistical uncertainties, while the systematic uncertainties of the data are shown as shaded areas.The dot-dashed line shows the STARLIGHT prediction, and the dashed line corresponds to the leading-order QED calculation of Ref. [48].
The neutron multiplicity dependence of α core for µ + µ − pairs in ultraperipheral Pb-Pb collisions at √ s NN = 5.02 TeV is shown in Fig. 2 (upper), in the mass region 8 < m µµ < 60 GeV.A strong neutron multiplicity dependence of α core is clearly observed, while the α core predicted by STARLIGHT is almost constant at a value of about 1.35 × 10 −3 , shown as the dotdashed line in Fig. 2 (upper).The α core for inclusive UPCs is measured to be [1227 ± 7 (stat) ± 8 (syst)] × 10 −6 , about 10% lower than the STARLIGHT prediction.In general, the α core in data becomes larger as the emitted neutron multiplicity increases.A fit to the dependence of α core on the neutron multiplicity with a constant value is rejected with a p value corresponding to 5.7 standard deviations.This observation demonstrates that initial photons producing µ + µ − pairs have a significant b dependence of their p T , which impacts the p T and acoplanarity of muon pairs in the final state.This initial-state contribution must be properly taken into account when exploring possible final-state EM effects arising from a hot QGP medium formed in hadronic heavy ion collisions [26,27].A recent leading-order QED calculation [48], incorporating a b dependence of the initial photon p T [32], has provided results for all the reported neutron mul-tiplicity classes.The average b values estimated in Ref. [48] range from about 112 to 22 fm for the 0n0n to XnXn neutron multiplicity classes, respectively.The model calculation can qualitatively describe the increasing trend of α core data with the neutron multiplicity, shown as the dashed line in Fig. 2 (upper).However, the data are systematically higher than the model calculation (plotted without uncertainties) by about 5%, which may be related to the presence in data of soft photon radiation from the muons [30].
Figure 3 shows a direct comparison between data and model calculations for the α distributions with α < 0.008 in all the reported neutron multiplicity classes.The α distribution is clearly observed to broaden in the high neutron multiplicity class.The QED calculation, incorporating a b dependence of the initial photon p T , can describe the α distributions reasonably well while STARLIGHT fails to describe the data.A direct comparison between data and model calculations for the α distributions with α < 0.008.The dot-dashed line shows the STARLIGHT prediction and the dashed line corresponds to the leading-order QED calculation of Ref. [48].The dot-dot-dashed and dotted lines indicate the core and tail contributions, respectively.The vertical lines and shaded areas represent the statistical and systematic uncertainties, respectively.
A rapidity dependence of the α distribution is also investigated for 0n1n, 0nXn and 1nXn classes (see Appendix A) for dimuon rapidity in the hemisphere containing larger (smaller) neutron multiplicity.In the 0nXn class, the tail contribution in the rapidity hemisphere with Xn is significantly larger than that in the rapidity hemisphere with zero neutrons, suggesting contributions from different higher-order processes that correlate with the dimuon pair production.However, no rapidity dependence is observed for the α core values extracted from the fits using Eq. ( 1).This is consistent with the expectation that the α core is dominated by leading-order γγ → µ + µ − scatterings and illustrates that the employed core functional form is robust.
In Fig. 2 (lower), the average invariant mass m µµ of all muon pairs passing the selection criteria, is shown as a function of the neutron multiplicity.A clear neutron multiplicity dependence of m µµ is observed, with the m µµ value measured in XnXn events being larger than that in 0n0n events with a significance exceeding 5 standard deviations.This trend of m µµ can be qualitatively described by both model calculations.As the muon pair invariant mass is largely determined by the initial photon energy, this observation suggests that the energy of the photons involved in UPCs is, on average, larger in collisions with smaller b, a conclusion similar to that previously drawn for the initial photon p T .
In summary, the first measurements of γγ → µ + µ − production as a function of forward neutron multiplicity in ultraperipheral lead-lead collisions at a nucleon-nucleon center-of-mass energy of 5.02 TeV are reported.A significant broadening of back-to-back azimuthal correlations is seen, with respect to the leading-order γγ → µ + µ − process, for increasing multiplicities of emitted forward neutrons.This observed trend is qualitatively reproduced by a leading-order quantum electrodynamics calculation, demonstrating the importance of an impact-parameterdependent photon p T .A similar trend of increasing average invariant mass of muon pairs with neutron multiplicity is also observed.These measurements provide the first experimental demonstration that the initial energy and transverse momentum of photons exchanged in ultraperipheral heavy ion collisions depend on the impact parameter of the interaction.These results call for theoretical efforts to improve the precision in modeling photon-induced interactions.Future searches for electromagnetic interactions of leptons inside the quark-gluon plasma created in heavy ion collisions should incorporate a baseline where the initial broadening effects presented in this Letter are properly taken into account.
[17] R.  .1 (left) shows the correlation between energy distributions of the ZDC detectors, located on the positive (Plus) and negative (Minus) directions with respect to the CMS interaction point, for events selected in the analysis.For the measured neutron multiplicity class with asymmetric neutron numbers, the dimuon rapidity is divided into two hemispheres using the plane defined by y = 0.In each rapidity hemisphere, the α distribution from γγ → µ + µ − is normalized by the total yields in this neutron multiplicity class ((1/N s )dN to values found by multiplying those for Υ(1S) by the ratio of the published masses of the states [42].The contribution of γγ scattering to dimuon pair production in the Υ mass region is extracted using a second order polynomial function.

A ZDC energy distributions and rapidity dependence of acoplanarity distributions and subtraction of Υ mesons
The reconstruction and trigger efficiencies rapidly reach a plateau as functions of p T with values of ∼95%-99% above p µ T ≈ 6 GeV for |η µ | < 1.2 and above p µ T ≈ 4 GeV for 1.2 < |η µ | < 2.4.Systematic uncertainties associated with the efficiency corrections are negligible since they largely cancel out in the final observables, which are normalized by the total yield.

Figure 1 :
Figure 1: Neutron multiplicity dependence of acoplanarity distributions from γγ → µ + µ − for p µ T > 3.5 GeV, |η µ | < 2.4, |y µµ | < 2.4, and 8 < m µµ < 60 GeV in ultraperipheral Pb-Pb collisions at √ s NN = 5.02 TeV.The α distributions are normalized to unit integral over their measured range.The dot-dot-dashed and dotted lines indicate the core and tail contributions, respectively, found using a fit to Eq. (1).The vertical lines on data points depict the statistical uncertainties, while the systematic uncertainties and horizontal bin widths are shown as gray boxes.

Figure 2 :
Figure 2: Neutron multiplicity dependence of (upper) α core and (lower) m µµ of µ + µ − pairs in ultraperipheral Pb-Pb collisions at √ s NN = 5.02 TeV.The vertical lines on data points depict the statistical uncertainties, while the systematic uncertainties of the data are shown as shaded areas.The dot-dashed line shows the STARLIGHT prediction, and the dashed line corresponds to the leading-order QED calculation of Ref. [48].
Figure3: A direct comparison between data and model calculations for the α distributions with α < 0.008.The dot-dashed line shows the STARLIGHT prediction and the dashed line corresponds to the leading-order QED calculation of Ref.[48].The dot-dot-dashed and dotted lines indicate the core and tail contributions, respectively.The vertical lines and shaded areas represent the statistical and systematic uncertainties, respectively.

Figure A. 1 (Figure A. 1 :
Figure A.1 (left)  shows the correlation between energy distributions of the ZDC detectors, located on the positive (Plus) and negative (Minus) directions with respect to the CMS interaction point, for events selected in the analysis.Figure A.1 (right) shows the measured Minus ZDC energy distribution together with a multi-Gaussian function fit.

Figure A. 2 :
Figure A.2: Acoplanarity distributions of γγ → µ + µ − events for three different neutron multiplicity classes with asymmetric neutron numbers.The solid red (open blue) symbols correspond to events where the dimuon rapidity is in the hemisphere containing larger (smaller) neutron multiplicity.The vertical lines on data points depict the statistical uncertainties while the systematic uncertainties are shown as shaded areas.The yields of muon pairs from γγ scattering in the Υ mass region (9 < m µµ < 11 GeV) are extracted by a binned χ 2 fit to the invariant mass spectrum, as shown in Fig. A.3.Each Υ state is modeled by a Gaussian function.All the parameters of the Υ(1S) fit are left free.For the Υ(2S) and Υ(3S) states, the yields are allowed to vary while the mean and width are fixed

Figure A. 3 :
Figure A.3: The efficiency corrected invariant mass distribution of muon pairs in inclusive ultraperipheral Pb-Pb collisions, for the kinematic range p µ T > 3.5 GeV, |η µ | < 2.4, and |y µµ | < 2.4.The result of the fit to the data is shown as solid blue line.The yields of muon pairs from γγ scattering in the Υ mass region are shown as dashed red line.The separate yields for each Υ state are shown as dotted violet lines.