Measurement of muon pairs produced via $\gamma\gamma$ scattering in non-ultraperipheral Pb+Pb collisions at $\sqrt{s_{_\mathrm{NN}}} = 5.02$ TeV with the ATLAS detector

Results of a measurement of dimuon photoproduction in non-ultraperipheral Pb+Pb collisions at $\sqrt{s_{_\mathrm{NN}}} = 5.02$ TeV are presented. The measurement uses ATLAS data from the 2015 and 2018 Pb+Pb data-taking periods at the LHC with an integrated luminosity of $1.94~\mathrm{nb}^{-1}$. The $\gamma\gamma \rightarrow \mu^+\mu^-$ pairs are identified via selections on pair momentum asymmetry and acoplanarity. Differential cross-sections for dimuon production are measured in different centrality, average muon momentum and pair rapidity intervals as functions of acoplanarity and $k_\perp$, the transverse momentum kick of one muon relative to the other. Measurements are also made as a function of the rapidity separation of the muons and the angle of the muon pair relative to the second-order event plane to test whether magnetic fields generated in the quark-gluon plasma affect the measured muons. A prior observation of a centrality-dependent broadening of the acoplanarity distribution is confirmed. Furthermore, the improved precision of the measurement reveals a depletion in the number of pairs having small acoplanarity or $k_{\perp}$ values in more central collisions. The acoplanarity distributions in a given centrality interval are observed to vary with the mean-$p_{\mathrm{T}}$ of the muons in the pair, but the $k_{\perp}$ distributions do not. Comparisons with recent theoretical predictions are made. The predicted trends associated with effects of magnetic fields on the dimuons are not observed.


Introduction
In relativistic heavy ion collisions, the Lorentz contracted electromagnetic fields of the ions act as a source of high energy quasi-real photons [1][2][3]. In collisions involving Pb ions at the LHC, the photons emitted coherently by the incoming nuclei can have energies of up to 80 GeV in the center of mass frame [4]. The processes involving the interactions of such coherently produced photons are typically studied in ultra-peripheral collisions (UPC); events where the contributions from strong interactions are absent as they typically occur at impact parameters larger than twice the nuclear radius. Such UPC events can lead to exclusive final states with dileptons [4][5][6]. The induced processes are also present in heavy-ion collisions with hadronic interactions, but are more difficult to observe due to the large amounts of background particles produced in such collisions.
Recent measurements of dilepton production via scattering in non-ultraperipheral (non-UPC) heavy-ion collisions [7,8] have stimulated significant interest due to the possibility that such pairs can be used as electromagnetic probes of the quark-gluon plasma created in such collisions [7,9,10]. In particular, the alignment of the two leptons in the transverse plane and/or the lepton-pair transverse momenta may be modified by the scattering of the leptons off constituents of the plasma or through interaction with long-range magnetic fields generated by the charge flow in the plasma. 1 Previous measurements [7] of the centrality dependence of dimuon acoplanarity in non-UPC Pb+Pb collisions at the Large Hadron Collider (LHC) have shown a systematic broadening of the acoplanarity distribution that is consistent with expectations of electromagnetic scattering of one or both of the muons in the quark-gluon plasma [9]. Similarly, a measurement of electron-pair transverse momentum distributions in non-UPC Au+Au collisions at the Relativistic Heavy Ion Collider showed a broadening compared to theoretical expectations for vacuum production [8]. That broadening was found to be consistent with the expected deflection of the electrons in magnetic fields produced during the collision [10]. However, it has also been argued that the effects observed in non-UPC heavy-ion collisions may result from physics associated with the initial state, particularly a variation of the transverse momenta of the initial photons [13] with the impact parameter of the nuclear collision. That idea led to a prediction [14] that the acoplanarity of dilepton pairs produced in ultraperipheral collisions may depend on the breakup of the incoming nuclei, due to the dependence of the breakup probability on the impact parameter. 2 Such an effect was observed recently [6].
A recent alternative analysis [11] proposes a quantum-mechanical extension of the equivalent photon approximation using the Wigner functions to describe, simultaneously, the initial-state photon transverse coordinate and transverse momentum distributions. Calculations using these photon Wigner functions [11,15] show effects similar to those seen in the data due to the modulation of the photon Wigner distributions as a function of transverse momentum and/or impact parameter with respect to the parent nucleus.
Specific tests, not possible with previous data, have been suggested for evaluating the possible role of magnetic effects [9,11]. In particular, it was proposed to study the dependence of the dilepton acoplanarity and/or transverse momentum imbalance on the rapidity separation of the leptons and on the dilepton momentum direction relative to the impact parameter vector. Magnetic fields are predicted to affect the leptons in a rapidly increasing way with increasing rapidity separation between the leptons and to be largest for leptons emitted perpendicular to the magnetic field, since their direction in the transverse plane is correlated with the direction of the impact parameter vector [16,17]. In noncentral heavy-ion collisions, the "elliptic" second-order event plane angle is well-correlated with the direction of the impact parameter [18]. Thus, a measurement of dileptons as a function of the angle between the dilepton axis and the second-order event plane may provide independent sensitivity to magnetic fields [19] that is crucial to the understanding of possible chiral magnetic effects [20,21] in nucleus-nucleus collisions.
In order to experimentally distinguish between the various explanations of the observed effects in dilepton production from scattering in non-UPC heavy-ion collisions, a larger data set, as well as a broader set of measurements, is clearly needed. To this end, this paper repeats and extends the measurements presented in Ref.
[7], taking advantage of the increased integrated luminosity of the 2018 Pb+Pb run at the LHC and improved trigger selections. It presents results from an ATLAS measurement of non-UPC production of → + − pairs in √ NN = 5.02 TeV Pb+Pb collisions at the LHC using a combination of 2015 and 2018 Pb+Pb data sets with a total integrated luminosity of 1.93 nb −1 . This integrated luminosity is about four times larger than that available for the previous measurement. The results presented here complement a recent ATLAS measurement [4] of exclusive UPC dimuon production performed over a more restrictive mass range but a factor of eight times larger range in acoplanarity. Small contributions from QED final-state radiation and dissociative processes studied in that measurement are both strongly suppressed in this analysis by kinematic selections (requirements on the , ⊥ and variables discussed below) and obscured by backgrounds generated in non-UPC Pb+Pb collisions. The centrality of those collisions is characterized by the total transverse energy within the acceptance of the ATLAS forward calorimeters. The → + − measurement is performed over a kinematic fiducial region: | | < 2.4, T > 3.7 GeV, and + − < 45 GeV, where and T are single-muon rapidities and transverse momenta, respectively, and + − is the dimuon invariant mass. Candidate muon pairs are identified using selections on pair acoplanarity, , and asymmetry, , defined as 3 where 1,2 and T 1,2 represent the azimuthal angles and the transverse momenta of each of the two muons, respectively.
The dominant background in this measurement, namely pairs of muons resulting from heavy-flavor (HF) decays, is suppressed through requirements on the pointing of the muons to the primary vertex. The remaining background is estimated using a template-fitting procedure based on the combined distance of closest approach of the two muons to the collision vertex. Potential non-negligible backgrounds from Drell-Yan (DY) processes are estimated using Monte Carlo (MC) simulations.
Distributions of , and the associated transverse momentum scale ⊥ , are measured as a function of Pb+Pb collision centrality and the average of the transverse momenta of the two muons,¯T ≡ ( T 1 + T 2 ) /2. It is shown later that the distributions vary significantly with¯T, while the ⊥ distributions do not. Thus, the ⊥ distributions are better suited for assessing the centrality-dependent modifications of the dimuon alignment. However, some of the theoretical calculations are only available for acoplanarity, so results are presented using both variables. Moments of the ⊥ distributions are used to quantify the centrality-dependent modifications of the dimuon alignment.
It was shown in Ref.
[7] that the dimuon asymmetry distributions are not sensitive to the small transverse momentum scales associated with the observed modifications of the dimuon alignment. However, the → + − asymmetry distributions are much narrower than those from background QCD processes. To further strengthen an already robust demonstration that the observed signal represents → + − pairs, the asymmetry distributions are measured in different centrality bins and compared with the results obtained from the STAR [9] event generator.
To test predictions [9,11] that magnetic broadening of the dimuon or ⊥ distributions should depend on the rapidity separation of the two muons, the ⊥ distributions are measured as a function of |Δ | ≡ | 1 − 2 |, where 1 and 2 represent the rapidities of the two muons. To test whether magnetic broadening effects have a directional dependence in the transverse plane, the ⊥ distributions are also measured as a function of 2Δ ≡ 2 − Ψ 2 , where , defined in Section 6.5, represents the orientation of the dimuons in the transverse plane, and Ψ 2 is the second-order event plane angle.
To allow tests of theoretical calculations, cross-sections for the production of muon pairs are measured in different centrality and¯T intervals. The total cross-section for → + − production, including UPC contributions, within the fiducial constraints of the measurement is also obtained. Separate "normalized pair yields" representing the fraction of the total Pb+Pb → + − yield measured in 5 a given centrality and kinematic interval are also presented. In these relative yields, some systematic uncertainties in the measurement cancel out.
The remainder of this paper is structured as follows: Section 2 describes the ATLAS detector; Section 3 describes the data and Monte Carlo samples used in the measurement and the applied event and dimuon selections; Section 4 describes the corrections for trigger and reconstruction inefficiency and the estimation and subtraction of backgrounds; Section 5 describes the systematic uncertainties in the measurement; Section 6 presents the results; and Section 7 presents a summary of the results and conclusions.

ATLAS detector
The ATLAS detector [22] is composed of an inner tracking detector (ID) inside a superconducting solenoid magnet, electromagnetic and hadronic calorimeters, and a muon spectrometer (MS) with superconducting toroid magnets, and has a high-speed trigger and data-acquisition system.
The ID, consisting of a silicon pixel detector, a silicon microstrip tracker, and a transition radiation tracker, is immersed in a 2 T axial magnetic field [23]. The ID provides charged-particle tracking in the pseudorapidity range | | < 2.5. The high-granularity silicon pixel detector covers the interaction region and typically provides four measurements per track. The pixel detector is followed by the silicon microstrip tracker, which typically provides measurements of four two-dimensional space points per track. These silicon detectors are complemented by the transition radiation tracker, which enables radially extended track reconstruction up to | | = 2.0, providing around 30 hits per track.
The MS comprises separate trigger and high-precision tracking chambers that measure the deflection of muons in the magnetic field of the superconducting air-core toroids. The precision chamber system covers the region | | < 2.7 with three layers of monitored drift tubes complemented by cathode strip chambers in the forward region. The muon trigger system covers the range | | < 2.4 with resistive plate chambers in the barrel and thin gap chambers in the endcap regions.
Two zero-degree calorimeters (ZDCs), which measure neutrons emitted from the incident nuclei at large absolute rapidities, are used for triggering and for offline event selection. The ZDCs are located symmetrically at a distance of ±140 m from the nominal interaction point and cover | | > 8. 3. Each of the ZDCs consists of four modules, each containing slightly more than one interaction length of tungsten absorber.
The ATLAS trigger system [24] consists of a first-level (L1) trigger implemented using a combination of dedicated electronics and programmable logic, and a software-based high-level trigger (HLT).
Muon triggers are formed using a combination of L1 triggers that find candidate muons from the MS trigger chambers and HLT triggers that combine ID and MS tracks. The L1 muon trigger has lower acceptance (| | < 2.4) than the full MS (| | < 2.7).
An extensive software suite [25] is used in the reconstruction and analysis of real and simulated data, in detector operations, and in the trigger and data acquisition systems of the experiment.
3 Data sets, reconstruction, and event selection

Data Sets
The Pb+Pb data used in this measurement were recorded in 2015 and 2018 with integrated luminosities of 0.49 nb −1 and 1.44 nb −1 , respectively. Two dimuon triggers were used for this analysis. The first trigger (L1Single) required a single muon with T > 4 GeV at L1, and two muons with T > 4 GeV at the HLT. The second trigger (L1Pair) required two muons with T > 4 GeV at L1 and also at the HLT. By requiring only one muon at L1, the first trigger has greater L1 trigger efficiency, but it was prescaled for a small portion of the data-taking and thus did not sample the full luminosity. A total of 4.7 million and 12.2 million events were obtained from these triggers in 2015 and 2018, respectively.
Separate samples of minimum-bias Pb+Pb events and Pb+Pb events selected by a combination of single-muon triggers are used to evaluate the efficiency of the muon triggers used in this analysis.
The minimum-bias sample is built using a combination of three mutually exclusive triggers. The first (second) trigger required the total transverse energy in the calorimeters at L1, L1 T , to be greater than 600 GeV (between 50 GeV and 600 GeV) without any additional requirements at the HLT. The third trigger required L1 T to be less than 50 GeV with the additional requirement of an energy deposition above the single-neutron threshold in either one of the ZDCs. At the HLT, this trigger additionally required a reconstructed track having T > 0.2 GeV. The single-muon triggered sample is built using a combination of three triggers that require a muon with T > 4, 6 and 8 GeV at the HLT.
The performance of the ATLAS detector in reconstructing muon pairs is evaluated using an MC sample obtained by overlaying → + − events produced with the STAR event generator onto minimum-bias Pb+Pb events simulated using the H v1.383 [26] event generator. The detector response in the MC samples was simulated using G 4 [27], and the resulting events are reconstructed using the same algorithms that are applied to the data [28]. A total of 4 million such events are analyzed using the same methods as applied in the data analysis. The STAR MC sample is also used for comparison with the measured observables.
Potential backgrounds from DY processes are estimated using the P B v2 [29][30][31] generator interfaced to P 8 configured using parameter values set to the AZNLO tune [32] and CTEQ6L1 [33] parton distribution functions (PDFs). Separate samples of , , and events were generated and combined with appropriate isospin weights. The P B generator was configured to provide per-event weights for five different nuclear PDF sets: nCTEQ15 [34], EPPS16 NLO [35], nNNPDF1.0 NNLO [36], nNNPDF2.0 NLO [37], and TUJU19 NNLO [38]. Thus, separate evaluations of the DY background in this measurement are obtained for all five PDF sets. Variations of the renormalization and factorization scales are performed using the nCTEQ15 set and are compared with a similar set of variations performed in P +P 8 using the nucleon CT14 set [39].

Event and muon-pair selections
Events used in the analysis are required to have been recorded during stable running conditions of the LHC, to have no detector hardware or readout error, and to have a reconstructed collision vertex. Charged-particle tracks and collision vertices are reconstructed using standard methods [40] tuned for the conditions of Pb+Pb collisions and assuming a single collision vertex per event. In addition to track kinematic parameters, the ID reconstruction also provides information about the minimum distances 0 and 0 between the projected track and the reconstructed vertex in the transverse and longitudinal planes, respectively. Muons are reconstructed by combining ID tracks with tracks reconstructed in the muon spectrometer. The muons are required to pass the "medium" muon selection requirements described in Ref. [41].
Opposite-sign muon pairs passing the following preselections are used for the analysis: each muon has T > 3.7 GeV 4 and | | < 2.4; the pair has a dimuon invariant mass less than 45 GeV; and both muons must be matched in angular space to HLT-reconstructed muons. These kinematic selections are largely determined by the acceptance of the MS; the mass restriction is applied to avoid contamination from boson decays.
To reduce the background from semileptonic decays of heavy-flavor hadrons, requirements are imposed on the pointing of the muons to the vertex using a combination of the single-muon 0 and 0 sin values (where sin is the polar angle of the muon track): Distributions of 0pair and ( 0 sin ) pair for pairs passing the above preselections are shown in Figure 1 in the top and bottom panels, respectively. Also shown for comparison are 0pair and ( 0 sin ) pair distributions for pairs passing a kinematic fiducial selection, described below, that suppresses the HF decay contribution. The fiducial selection strongly suppresses the yield of pairs with large 0pair  Figure 1: Distributions of 0pair (top) and ( 0 sin ) pair (bottom) for muon pairs passing the preselections (black), and for pairs additionally passing the Fid-( < 0.06, < 0.012) selection (red). The error bars in both panels correspond to the statistical uncertainties, and are typically too small to be seen. The pairs passing the Fid-selection have much smaller 0pair and ( 0 sin ) pair values than pairs passing the preselections, due to the large HF background in preselected pairs. and ( 0 sin ) pair values that predominantly result from HF-decay background pairs. The following selections are imposed on muon pairs used in the measurement: These requirements reduce the yield of HF-decay pairs by a factor of ≈2 while introducing an inefficiency for → + − pairs of 2%.
Following the methods of Ref.
[7], candidate → + − pairs are obtained from those passing the preselection and the 0pair and ( 0 sin ) pair requirements by imposing stringent requirements on the pair asymmetry and either the acoplanarity or the ⊥ value. For this paper, two different fiducial selections are defined: < 0.06 ∧ < 0.012 or < 0.06 ∧ ⊥ < 150 MeV, labeled Fid-and Fid-⊥ , respectively. Both fiducial selections include the muon pseudorapidity and T requirements and the pair mass constraints included in the preselections. The separate fiducial selections are motivated by the HF and DY subtraction that is discussed later in Section 4.3. In particular, the backgrounds are observed to be uniform as a function of and ⊥ as long as no requirement is imposed on the other variable. However, because of the direct relationship between and ⊥ made explicit in Eq. (1), a selection on introduces a¯T-dependent constraint on ⊥ and vice versa. A single fiducial selection would thus distort the shapes of the HF backgrounds and make those shapes sensitive to an applied¯T selection. For that reason, separate fiducial regions are used. Specifically, the Fid-selection is used for measurements of acoplanarity distributions, and the Fid-⊥ selection is applied in measurements of ⊥ distributions.
Otherwise, for consistency with the previous measurement [7], the Fid-selection is applied for many of the plots in Section 4 that document technical details of the analysis, while measurements of the production cross-sections and related quantities are presented using only the Fid-⊥ selection.

Centrality
In ATLAS heavy-ion measurements, the Pb+Pb collision centrality is characterized by the total transverse energy, Σ FCal T , measured in the ATLAS forward calorimeters. The relationship between Σ FCal T and the geometry of the Pb+Pb collisions is evaluated using a Glauber model analysis [42,43] following standard methods (see details presented in Ref. [44]). That analysis also provides values for the nuclear overlap parameter, AA , which describes the effective nucleon-nucleon luminosity of a Pb+Pb collision. While the → + − process may have a very different dependence on Pb+Pb collision geometry than soft or hard QCD scattering processes, the use of standard Pb+Pb centrality intervals nonetheless remains useful, for example, in estimating backgrounds from QCD processes. This analysis uses the following set of centrality intervals, defined in terms of percentiles of the minimum-bias Pb+Pb Σ FCal T distribution: ten-percent intervals spanning the range 10-90%, plus the two intervals 0-5% and 5-10%. A set of larger, combined centrality intervals are used in specific instances to reduce statistical uncertainties in some of the measured quantities or distributions. To cover the 10% most peripheral collisions, which have Σ FCal T < 24 GeV and, for which, centrality calibrations are not available, 5 a set of Σ FCal T intervals, labeled ET0-ET3, are used. These have lower Σ FCal T boundaries at integer multiples of 5 GeV except for ET0, which includes all events having Σ FCal T < 5 GeV. To suppress the contribution of ultraperipheral collisions for the measurement of non-UPC → + − production, events included in the above centrality intervals are required to have at least one neutron in each ZDC. A separate set of exclusive UPC events which have zero neutrons in one or both of the ZDCs and no additional charged particles beyond the two muons, is used to evaluate the effects of the 0pair and ( 0 sin ) pair selections and for comparison with the non-UPC results.
The inclusive cross-section for production of dimuons in √ NN = 5.02 TeV Pb+Pb collisions is measured for tests of theoretical calculations and for use in this paper. For this purpose, no event selections are imposed beyond those described in Section 3.2. The inclusive measurement includes pairs that pass neither the exclusive UPC requirement nor the ZDC coincidence requirement applied to the non-UPC intervals. Such pairs are referred to as "unassigned". Table 1 summarizes the Σ FCal T intervals used in this measurement, indicates the corresponding centrality range, where available, and for those with a centrality calibration, provides the corresponding values. Figure 2 shows, in the left panel, the Σ FCal T distribution for events with muon pairs passing different stages of pair selection: preselections only; also passing the 0pair < 0.1 mm and ( 0 sin ) pair < 0.2 mm requirements; and passing these plus the Fid-selection. Also shown is the Σ FCal  the yield increases with increasing Σ FCal T relative to that for minimum-bias Pb+Pb events and the Fid-selection. This behavior results from the geometric enhancement of QCD hard-scattering processes -particularly the production of heavy flavor -the rates of which are proportional to AA . The application of the Fid-selection yields a Σ FCal T distribution that is almost flat at large transverse energies. However, even with the fiducial selection, there remains a non-negligible and centrality-dependent HF-decay background that must be subtracted. Thus, the apparent flatness of the Σ FCal T distribution is at least partially accidental. However, as is seen below, even in the most central collisions, the HF backgrounds are, at most, comparable to the signal pair yields, so → + − production is observed over the full range of Pb+Pb collision centralities, including the most central collisions.

Analysis
A total of 69490 muon pairs pass the combination of preselections, 0pair and ( 0 sin ) pair requirements, and Fid-selection. For comparison, 67789 pairs pass the full selections and the Fid-⊥ selection. The smaller Fid-⊥ yield results from the ⊥ selection being more restrictive than that for the acoplanarity. Both fiducial selections reject about 10% of the pairs in the UPC sample that pass the preselections and pair vertex requirements. The excluded events primarily result from hard QED radiation or dissociative photon-induced processes (see Ref. [45] and references therein) but may also include background, non-→ + − events. For non-UPC events the fiducial selections also Table 2: Numbers of muon pairs in each centrality interval for different¯T ranges and the two fiducial selections. Also listed are the total → + − yields. The interval labeled as Unassigned includes pairs that contribute to the inclusive yield but do not appear in any of the listed centrality intervals, as they do not satisfy the ZDC coincidence requirement imposed on the non-UPC centrality intervals, and also do not satisfy the exclusive UPC requirement. suppress QCD backgrounds which generate muon pairs that typically have much larger acoplanarity and ⊥ values than → + − pairs. Table 2 lists the number of pairs passing the Fid-and Fid-⊥ selections in each centrality interval for three intervals of¯T.
Dimuon mass and¯T distributions for all dimuons passing the Fid-selection are shown in Figure 3. The distributions for measured pairs are suppressed at low + − and¯T values owing to inefficiencies in the muon trigger and offline reconstruction. Corrections for the resulting losses are discussed in the following section. An additional suppression at low + − results from the acoplanarity and asymmetry requirements imposed as part of the fiducial selection. A decrease observed in the¯T distribution at the highest values results from the + − < 45 GeV requirement.

Trigger and reconstruction efficiency
In the analysis described below, muon pairs are corrected for trigger and reconstruction efficiency by application of a per-pair weight to data events: where trig is the pair trigger efficiency, rec 1 and rec 2 are the single-muon reconstruction efficiencies, and vtx is the vertex-pointing efficiency.
Since all dimuons used in this analysis are back-to-back in azimuth and are thus well-separated in the detector, the efficiency of the dimuon triggers is evaluated using independent, single-muon trigger efficiencies. The single-muon trigger efficiencies are evaluated by testing whether offlinereconstructed muons passing the preselections described in Section 3 are matched to a muon found by the trigger. These efficiencies are evaluated as a function of muon T and of the product of the muon charge and pseudorapidity, . The single-muon trigger efficiencies are combined -separately for the L1Single and L1Pair triggers -to produce per-pair trigger efficiencies. The trigger efficiencies are found to change only by a few percent over the full centrality range of the measurement, whereas the detector occupancy varies by orders of magnitude.
The efficiency for reconstructing dimuons in the offline analysis is evaluated using the STAR MC simulation sample. As with the trigger efficiencies, the pair reconstruction efficiency is taken as the product of the single-muon reconstruction efficiencies. The single-muon reconstruction efficiencies are evaluated as a function of muon T and . They are corrected for small data-MC differences observed in previous measurements [41]. The reconstruction efficiencies have negligible centrality dependence at mid-rapidity (| | < 1), while at forward rapidity (1 < | | < 2.4) they decrease by ≈10% between peripheral and central collisions.
The pair pointing requirements introduce an inefficiency of a few percent for → + − pairs. The corrections for this inefficiency are obtained from simulation and validated using the UPC data sample, where the 0pair and ( 0 sin ) pair selections remove 1.4% and 0.5% of the pairs, respectively. The 0pair and ( 0 sin ) pair distributions in the STAR MC sample generally agree well with those in data, although the simulation slightly underestimates the 0pair resolution for UPC collisions. To rectify this disagreement, additional T -dependent Gaussian smearing is applied to the MC single-muon 0 values in all events. With this adjustment, the MC-evaluated efficiency for → + − pairs to pass the 0pair and ( 0 sin ) pair vertex selections, vtx , agrees well with the efficiency in UPC events. For non-UPC collisions, the vertex selection efficiency decreases by a few percent between the most peripheral and the most central collisions.
The effect of correcting the dimuon kinematic distributions for the trigger and muon reconstruction efficiencies is illustrated in Figure 3, where the efficiency-corrected + − and¯T distributions for the Fid-selection are shown with the black points. The corrections for efficiency increase the yield relative to the measured distributions by a factor of ≈2 for + − = 8 GeV and¯T = 4 GeV. At higher + − and¯T, the corrections increase the dimuon yield by ≈30% for the most peripheral collisions (ET0-ET3 or UPC) and ≈50% in the most central (0-10%) collisions.

Heavy-flavor decay background estimation
As mentioned above, the dominant background in this measurement results from muons produced in semileptonic decays of heavy quarks. At lower¯T values, the background is dominated by combinatoric pairs produced in the decays of uncorrelated heavy quarks. However, with increasing¯T, the background becomes dominated by correlated pairs. The HF-decay muons can be distinguished from prompt muons by their displaced secondary vertices that are experimentally observed via the 0 and 0 sin parameters. The HF background surviving the 0pair and ( 0 sin ) pair selections is estimated using a template-fitting procedure applied to the 0pair distributions. These distributions are assumed to result from the combination of pure → + − signal and a background contribution. The shape of the signal 0pair distribution, called the signal template below, is obtained from the STAR +H MC sample. The shape of the background 0pair distribution, called the background template below, is obtained directly from data using muon pairs that pass the preselections, but have > 0.06 and > 0.012. These last selections effectively eliminate the → + − contribution, leaving only background pairs. While the background subtraction is primarily designed to remove HF-decay contributions, the background template contains all muon pairs without a strong back-to-back angular correlation.
The template fits are performed over an extended 0pair range, 0pair < 0. 3 Figure 4: Results of template fits to measured 0pair distributions for pairs passing the muon preselections, the ( 0 sin ) pair requirement and the Fid-selection ( < 0.06 and < 0.012). Each panel represents a different centrality interval. The error bars shown on the data and the templates represent statistical uncertainties only. For many of the points, the error bars are smaller than the size of the marker.

70-80%
of the signal and combined templates over 0pair < 0.1 mm: . (2) The above-described template fits are performed separately for the Fid-and Fid-⊥ selections. The resulting signal fractions are denoted by sig and ⊥ sig , respectively. To allow measurements as a function of¯T, separate template fits are performed in all centrality intervals for pairs in restricted T intervals: 4 <¯T < 5 GeV, 5 <¯T < 6 GeV, 6 <¯T < 8 GeV, and¯T > 8 GeV, as well as for an integrated interval,¯T > 4 GeV. Other template fits are performed in intervals of | |, |Δ |, and |2Δ |, to allow differential measurements as a function of these quantities. Figure 4 shows examples of the template fits applied to the 0pair distributions in several centrality intervals for pairs satisfying the preselections, the ( 0 sin ) pair requirement, and the Fid-selection. In all of the analyzed centrality intervals, the template fits reproduce the data well; the minimum 2 -equivalent values are consistent with expectations given the number of degrees of freedom. The  (2) as a function of centrality for different selections on¯T. Generally, the signal fractions become larger with increasing¯T and for more peripheral collisions. In particular, for¯T > 4 GeV, the signal fraction in the 0-5% centrality interval is ≈50% while, for the UPC selection, it is consistent with one. This decrease of the signal fractions between peripheral and central collisions results from the geometric enhancement of heavy-quark production relative to → + − processes in more central collisions. However, with increasing¯T, the centrality dependence becomes weaker such that for the highest¯T interval,¯T > 8 GeV, the signal fraction varies by only ≈12% as a function of centrality. This behavior is understood to result from the quenching of heavy quarks (see Refs. [48,49] and [50] and references therein) that suppresses the correlated HF-decay background in more central collisions.
To evaluate the sensitivity of the above results to the shape of the background template, three alternative selections for the background template were applied: > 0.1 ∧ > 0.012, 0.06 < < 0.3 ∧ > 0.012, and > 0.06 ∧ > 0.2. The first two significantly change the requirements on the muon momentum balance while the third places a tighter requirement on the angular alignment. The results obtained with these alternative background templates are consistent with the nominal results.
The numbers of signal and background pairs, fid sig and fid bkg for a given fiducial selection and in a given centrality and¯T interval are given by fid sig (cent,¯T) = fid sig (cent,¯T) fid (cent,¯T), where fid sig (cent,¯T) and fid (cent,¯T) are, respectively, the signal fraction and efficiency-corrected number of pairs satisfying a given fiducial selection in the specified centrality and¯T intervals.
Once the signal fractions fid sig (cent,¯T) are determined, the measured acoplanarity, ⊥ , and distributions can be corrected by subtracting the contribution of background pairs. The shapes of the background and ⊥ distributions are evaluated by selecting muon pairs passing the pair preselections and having > 0.06. This asymmetry requirement effectively removes contributions from → + − pairs leaving, ostensibly, the HF-decay background. To determine the shape of the background asymmetry distribution, pairs passing preselections and having > 0.012 are used. Signal fraction The and ⊥ distributions are obtained by applying > 0.06 to preselected pairs. The distribution is obtained by applying > 0.012. The error bars correspond to statistical uncertainties. The lines indicate fits to constant functions for and ⊥ and to a linear function for .
The resulting differential distributions in , ⊥ , and , normalized to unit integral, are shown in Figure 6. The and ⊥ distributions are found to be uniform within their statistical uncertainties, so they are taken to be constants, HF , with values given by such that the integral of each constant over the fiducial range of the variable yields the number of background counts obtained from the template fitting. The distribution has nonzero slope, so it is fitted with a linear function in each centrality interval.  Figures 7 and 8 show, for several centrality intervals, differential dimuon yields versus and ⊥ , respectively, for the production of dimuons satisfying the corresponding fiducial selection. The distributions are plotted over ranges that extend beyond the fiducial limits in order to emphasize the behavior of the HF-decay background. Also shown on the plots are the background levels represented by the constants, HF and ⊥ HF , together with statistical uncertainties of the constants resulting from the template fit. The backgrounds provide a good description of the behavior of the and ⊥ distributions at large and ⊥ , although the data show a slight, centrality-dependent 20 excess over those HF background estimates. When viewed over ranges of and ⊥ that extend well beyond the fiducial bounds, that excess is seen to be uniform versus and ⊥ , so it is unlikely to have significant contributions from actual → + − pairs.

Differential distributions
To evaluate the contribution from the observed excess, the measured and ⊥ distributions are fitted with constant functions over the ranges 0.012 < < 0.02 and 150 MeV < ⊥ < 500 MeV, respectively. The results of those constant fits -referred to as "asymptotic fits" in the remainder of this paper -are indicated with the blue, long-dashed lines in Figures 7 and 8.

Drell-Yan background estimation
As noted above, Drell-Yan production of prompt dimuons represents a potential background in this measurement. The relevant DY pairs are produced in the scattering of quarks and antiquarks with momentum fractions typically less than 10 −2 . Thus, the DY pair yield is sensitive to the degree of shadowing in nuclear PDFs.
The DY pairs obtained from the P +P 8 MC sample are required to pass the same kinematic preselections that are applied to the data. The resulting P +P 8 differential cross-sections as a function of and ⊥ are shown in Figure 9. They are plotted over both the full dynamic range populated by the MC sample and, in the insets, for pairs satisfying < 0.06 over intervals of and ⊥ that are twice the corresponding fiducial ranges. Results are shown for three of the five nuclear PDFs sets, but the results from all sets differ primarily in the overall normalization. For the sake of clarity, the effects of PDF uncertainty and renormalization and factorization scale variations are not shown in this figure. The and ⊥ distributions vary substantially as a function of and ⊥ but are constant, within statistical uncertainties, over twice the fiducial ranges of the corresponding variables. Only a small fraction of the DY pairs passing the preselections used in this measurement subsequently satisfy the severely restrictive fiducial selections. For example, using the nCTEQ15 PDF set, 1.8% and 1.2% of the preselected DY pairs satisfy the Fid-and Fid-⊥ requirements, respectively. Table 3 lists the effective nucleon-nucleon (NN) cross-sections, Fid DY,NN , for production of DY dimuons, within the Fid-and Fid-⊥ fiducial regions, obtained from P +P 8 for different nuclear PDF sets and for the CT14 nucleon PDF set. The table also shows uncertainties obtained by propagating PDF systematic variations through P +P 8 [51].
The P +P 8 DY cross-sections within the fiducial regions used in this measurement vary by ≈30% between the different nuclear PDF sets, with the nCTEQ15 PDFs yielding the smallest cross-sections and nNNPDF2.0 and TUJU19 producing the largest cross-sections. However, the cross-sections for even those two nuclear PDF sets are smaller than that obtained from the CT14 set by about 30% due to nuclear shadowing. Because the most recent of the implemented nuclear PDF sets, TUJU19 and nNNPDF2.0, yield consistent results and have smaller uncertainties, the P +P 8 simulations produced with the nNNPDF2.  Figure 9: Differential cross-sections versus (left) and ⊥ (right) for Drell-Yan production of dimuons in √ NN = 5.02 TeV nucleon-nucleon collisions from P +P 8 for three nuclear PDF sets. The error bars indicate statistical uncertainties. The insets show the differential cross-sections for pairs with < 0.06. They are plotted over ranges of or ⊥ that are twice as wide as the corresponding fiducial regions. The statistical uncertainties are correlated between the different nPDF sets because all results derive from the same set of P +P 8 events.
the DY background in this measurement.
The number of DY pairs produced in a given centrality interval within a fiducial region can be estimated using the effective NN cross-sections obtained above: where Pb+Pb had is the total Pb+Pb hadronic (i.e. excluding UPC) cross-section and Δcent represents the width of a centrality interval expressed as an absolute fraction. The terms in the first bracket in the equation express the total number of Pb+Pb collisions sampled for the luminosity used in Table 3: Effective nucleon-nucleon cross-sections obtained from P +P 8 for the production of Drell-Yan muon pairs in √ NN = 5.02 TeV collisions using different nuclear PDF sets. The systematic uncertainties of the fiducial cross-sections obtained by propagating PDF uncertainties through P +P 8 are also shown. A separate ±15% uncertainty in the cross-sections, due to factorization and renormalization scale uncertainties, is not included in the shown uncertainties. Because the AA values increase in more central collisions -they vary by more than two orders of magnitude between the most peripheral and most central intervals used in this measurement -the DY background will be largest in the most central collisions.

PDF set
Since the P +P 8 DY and ⊥ distributions are uniform within uncertainties over the fiducial ranges of the measurement, for the differential measurements of these quantities, the DY backgrounds are added to the HF backgrounds to yield constants that are subtracted from the measured distributions: For the measurement of the asymmetry (| |) distributions, the P +P 8 dimuon asymmetry (| |) differential cross-sections, scaled according to Eq. (4) to produce differential background yields, are subtracted from the data. Figure 10 shows, in the top panels, the centrality dependence of the estimated DY pair yields for the nCTEQ15 and nNNPDF2.0 PDF sets compared with the excess pair yields at large and ⊥ described in Section 4.3. The production of DY pairs can account for a substantial fraction of the observed excess, but the P +P 8 simulations systematically underestimate the excess yield even with the PDF set(s) that predict the largest DY rates. The differences between the observed background excess and the DY rates obtained using the nNNPDF2.0 PDF set, expressed as a fraction of the fiducial → + − yields, are shown in the bottom panels of Figure 10 for the two fiducial selections. No significant centrality dependence of the excess is evident within the relatively large uncertainties. Notably, the excess persists in peripheral events where HF and DY backgrounds are small or negligible. Averaged over centrality, the excess corresponds to roughly 4% of the fiducial yields.

Excess observed at large and ⊥
Dimuons produced in the decays of Υ states 6 represent a plausible potential background in this measurement. However, the imposed fiducial selections strongly suppress such a background by restricting the acceptance to Υ mesons having essentially zero T . The mass distribution shown in Figure 3 has an enhancement near the Υ mass, but that enhancement is purely a result of the fiducial selection. When examined with much finer binning, no evidence of an excess in the region around the Υ mass is seen for pairs passing the Fid-selection, while a significant enhancement is observed in events passing the preselection and vertex-pointing requirements but failing the acoplanarity and 6 There is zero acceptance for / decays due to the mass requirement implicit in the fiducial selections.  Table 3. The data-points for the DY calculations are shifted slightly along the -axis, for clarity. Bottom row: The differences between the observed background excess and the DY rates obtained using nNNPDF2.0 PDF set, expressed as a fraction of the fiducial → + − yields. The left and right panels correspond to the Fid-and Fid-⊥ selections, respectively. asymmetry selections.
To test for a possible contribution of events associated with dissociative photon emission by one of the nuclei, the rapidity distribution of dimuons in the ⊥ interval 150 < ⊥ < 300 MeV, where the coherent → + − pairs make little contribution, was compared with the corresponding distribution for dissociative processes obtained from SuperChic4 [52]. As observed in the recent ATLAS UPC dimuon measurement [4], dissociatively produced pairs lie preferentially at large rapidities because the dissociative photons tend to be significantly higher in energy than photons produced coherently by a nucleus. However, the shape of the measured rapidity distribution for dimuons within the above ⊥ interval was consistent, within uncertainties, with that for pairs having ⊥ < 150 MeV. Thus, the excess background does not appear to result from dissociative processes that involve breakup of a nucleon.
The absence of a clear dissociative component in non-UPC collisions does not contradict the observation in Ref. [4] that dissociative events comprise 15% of the so-called XnXn sample which requires at least one neutron in each ZDC. In UPC collisions, nuclei emit neutrons primarily as a result of Coulomb excitation processes [53] that have a low probability to break up one, or especially both, nuclei. For example, in the mass region used for the measurement in Ref. [4] the per-nucleus single-breakup probability was estimated to be 23% while the double-breakup probability is 5%. Since dissociative processes substantially increase the likelihood that the photon-emitting nucleus emits neutrons, dissociative events are significantly enhanced in UPC XnXn events. However, in non-UPC Pb+Pb collisions, neutrons are produced with unit probability for all except the most peripheral collisions. Thus, the XnXn requirement does not enhance the rate for dissociative processes, with the result that the dissociative contribution is, at most, a few percent of the pairs over the kinematic range considered in Ref. [4]. Of those, about 40% fall within the Fid-range of this measurement. Such a small contribution would not be visible in the presence of the much larger HF and DY backgrounds, except possibly in the very most peripheral collisions.
Since the origin of the excess observed in the data is not known, it is not subtracted as part of this measurement. However, to account for the possibility that it represents an unidentified background, it is treated as a systematic uncertainty of the background subtraction procedure.

Background-subtracted , ⊥ , and distributions
The HF+DY background-subtracted acoplanarity distributions for → + − pairs having¯T > 4 GeV are shown in Figure 11 for several centrality intervals. The error bars indicate statistical uncertainties only. The figures also show the distributions from the STAR MC sample for both the generated and reconstructed muon pairs. The generated and reconstructed distributions differ only slightly due to the excellent angular resolution of the inner detector. The data and reconstructed MC distributions agree well in the UPC and most peripheral non-UPC centrality intervals, but in more central collisions the data systematically deviate from the MC predictions, with the data having wider distributions and a suppression at the smallest values.
Results similar to those shown in Figure 11, but for the HF+DY-subtracted ⊥ distributions, are presented in Figure 12. Like the distributions, the ⊥ distributions broaden and show a depletion near ⊥ = 0 with decreasing centrality percentile. However, the depletion becomes more significant than that observed in the distributions and becomes apparent at larger centrality intervals. Similar to the results, the MC generated and reconstructed distributions barely differ while the reconstructed MC and data distributions have very different behavior except in the most peripheral centrality intervals. Figure 13 shows HF+DY background-subtracted distributions compared with the generated and reconstructed STAR MC distributions. The data agree well with the reconstructed STAR distributions in all centrality intervals. This provides further evidence that the muon pairs in this measurement result from → + − processes. The resolution of the measurement, demonstrated by the difference between the MC generated and reconstructed distributions, is poor relative to that of the and ⊥ measurements and is not sufficient to probe effects at transverse momentum scales of 100 MeV.    Figure 12: HF+DY background-subtracted ⊥ distributions for pairs satisfying the Fid-⊥ selection and havinḡ T > 4 GeV in different centrality intervals from the most central 0-5% (top left) to the UPC interval (bottom right). For a few panels the distributions are scaled to allow a common -axis range for the plots. The scale factors are stated on the panels. Also shown for comparison are the generated and reconstructed distributions obtained from the STAR simulation samples. The STAR generated and reconstructed distributions are scaled to match the the corresponding data distributions over the ⊥ < 150 MeV interval.    Figure 13: HF+DY background-subtracted asymmetry distributions for pairs having¯T > 4 GeV in different centrality intervals from the most central 0-5% (top left) to the UPC interval (bottom right). In a few of the panels, the distributions are scaled to allow a common -axis range for the plots. The scale factors are stated on the panels. Also shown for comparison are the generated and reconstructed distributions obtained from the STAR simulation samples. The STAR reconstructed distributions are scaled to match the the corresponding data distributions over the 0 < < 0.1 interval. The STAR generated distributions, which are much narrower than the measured distributions or the data, are scaled vertically to allow direct comparison with the other distributions.   Figure 14: The (top) and ⊥ (bottom) distributions of dimuons before and after unfolding. The left and right panels correspond to the 0-10% and UPC intervals, respectively.

Unfolding
Iterative Bayesian unfolding [54] implemented using the ROOUnfold [55] framework is used to unfold the and ⊥ distributions. Although the resolution in each of these variables is comparable to or smaller than the histogram bin widths used in Figures 11 and 12, the observed suppression at small and ⊥ values may be influenced by the finite resolution in these quantities. Also, the distributions in the most peripheral centrality intervals are sufficiently steep that the finite or ⊥ resolution noticeably affects their shape. Response matrices in both and ⊥ are produced from the STAR MC sample. Because the STAR sample contains few pairs at large or ⊥ , those regions of the response matrices are augmented using parameterizations of the pair and ⊥ response. Different response matrices are produced for all combinations of centrality and¯T intervals used in the measurement.
Although the response matrix for the unfolding is nearly diagonal, the steep derivative in the and ⊥ distributions at small values means that the unfolding takes a few iterations to converge.
Thus, a conservative choice to use three iterations was made. The statistical uncertainties of the unfolded results, obtained using pseudoexperiments on the data and on the response matrix, are essentially unchanged from the distributions before unfolding. The reweighting of the MC ⊥ and distributions to match the data has negligible impact on the unfolded distributions.
The and ⊥ distributions are unfolded prior to the HF+DY background subtraction. Since the background is flat in both and ⊥ , it is unaffected by the unfolding, so the HF+DY subtraction is the same for the unfolded distributions as it is in the measured distributions shown above. No attempt is made to unfold the asymmetry distributions, as the resolution precludes sensitivity to effects at the momentum scales observed in the and ⊥ distributions. However, an efficiency correction, obtained from the STAR MC sample, is applied to account for a loss of pairs resulting from the < 0.06 selection due to migration in . The use of the STAR asymmetry distribution is acceptable for this purpose because of the insensitivity of the asymmetry to the observed differences between the data and MC samples.

Observables
Cross-sections for → + − production are obtained in intervals of centrality and¯T using: where L represents the integrated luminosity, and sig represents the number of backgroundsubtracted pairs (see Eq. (3) and Eq. (4)). The cross-sections in Eq. (5) can be evaluated using the Fid-or Fid-⊥ fiducial selections. However, the two fiducial selections give results that agree within statistical and systematic uncertainties. For this reason, the cross-sections are presented in this paper only for the Fid-⊥ selection.
To reduce systematic uncertainties resulting from overall normalization factors, e.g. due to efficiency corrections and luminosity, "normalized yields" ( ) are calculated according to where cent (¯T) is defined in Eq. (5) and tot (¯T) is the total cross-section, obtained by summing cent (¯T) over all centrality intervals, including the UPC interval.
Differential cross-sections for → + − production are obtained from the unfolded and ⊥ distributions. In the ⊥ case, the differential cross-section is calculated as where Δ represents the number of unfolded counts in a given ⊥ interval of width Δ ⊥ , and other factors are defined in Section 4. A similar formula applies to the calculation of / . These cross-sections are obtained for each combination of centrality and¯T intervals included in the analysis using the corresponding fiducial selection.

Systematic uncertainties
This section discusses estimates of the systematic uncertainties in the cross-sections, normalized yields, and ⊥ distributions. The uncertainties in the distributions mirror those in the ⊥ distributions and are not discussed separately. The following sources of systematic uncertainty are considered: • Luminosity: The systematic uncertainty of the combined 2015 and 2018 Pb+Pb integrated luminosity is 1.5%, obtained using the LUCID-2 detector [56] for the primary luminosity measurement [57].
• Muon working point: The analysis is performed using the "medium" muon selection requirements [41]. The measurements are repeated using the "tight" working point [41], which results in an increase in the muon purity, but reduces the muon-pair yields by about 20%. After applying an efficiency correction appropriate for the alternative working point, the cross-sections and normalized yields differ by 7% and 2% respectively. This variation is included as a systematic uncertainty that covers possible contributions from misidentified muons. For the differential cross-sections, the change in the muon working point mainly introduces point-by-point statistical fluctuations with no systematic trend beyond the 7% change in the overall normalization. Thus, no additional systematic uncertainty is applied to the and ⊥ differential cross-sections or yields.
• Trigger efficiency: The trigger efficiencies are evaluated directly from a data sample large enough for the statistical uncertainties in the efficiencies to be negligible. However, a systematic uncertainty is applied to cover the effects of possible centrality-dependent differences between the composition of muons used for the efficiency measurement and that of muons used in the measurement. It is evaluated from the following difference: where the averages are taken over centrality intervals. The systematic uncertainty essentially compares the efficiencies obtained using "medium" and "tight" muon selections in a given centrality interval while removing average changes due to the different muon populations. This uncertainty is within 3% across most of the centrality range, but increases to 6% for the 0-10% most central collisions.
• Reconstruction efficiency: The effects of the muon reconstruction efficiency are accounted for by applying weights -the inverse product of the single-muon efficiencies -to each muon pair used in the analysis. The efficiencies are functions of T , , and collision centrality. Uncertainties in the efficiencies produce uncertainties in the measured distributions and in the cross-sections and yields. These uncertainties are evaluated by varying the reconstruction efficiencies within their own uncertainties [41] and by evaluating the resulting changes in the measurements. This procedure yields a ≈2% systematic uncertainty in the cross-sections and a 0.5% uncertainty in the normalized yields. For the differential measurements, the variation of the reconstruction efficiency introduces no systematic effects beyond a change in the normalization which is already accounted for in the cross-section.
• 0pair and ( 0 sin ) pair selections: A potential systematic uncertainty associated with the vertex-pointing requirement is assessed by varying the 0pair requirement from its default value of 0pair < 0.1 mm, to 0.08 and 0.14 mm, and the ( 0 sin ) pair requirement to ( 0 sin ) pair < 0.22 mm (from 0.2 mm) and including the resulting variation as an uncertainty. The uncertainty covers both the efficiency for signal muons to pass the selections and the effect of the vertex-pointing requirements on the HF background subtraction. The former is most relevant in the UPC and peripheral centrality intervals where the uncertainty is small or negligible. The latter is relevant in more central collisions. For the cross-section and normalized yields measurement, the largest systematic uncertainty, ≈2%, is obtained in the 0-5% centrality interval. Varying the 0pair and ( 0 sin ) pair requirements introduces no systematic variation of the dimuon yields with or ⊥ .
• Background shape parameterization: The corrected and ⊥ distributions are obtained by subtracting a background that is taken to be constant, consistent with fits to the distributions in Figure 6. The sensitivity of the results to the assumption of a constant background is evaluated by parameterizing the background by linear functions and performing the subtraction using the resulting function. The deviation from the default analysis result is included as a systematic uncertainty. This uncertainty applies only to the differential cross-sections. It is at most ≈0.2% and is negligible compared to the other uncertainties.
• Signal template variation for 0pair fits: The HF background levels estimated from the template fits depend on the shape of the signal template that, in the default analysis, is obtained from MC simulation. An alternative, data-driven approach uses for the signal template the 0pair distribution measured in the UPC sample, which is composed almost entirely of → + − pairs. An additional, centrality-dependent smearing determined from the MC simulation is applied to the UPC 0pair distribution to account for the poorer 33 0 resolution in the more central collisions. A ≈1% variation between the results obtained from the default analysis and from the alternative signal templates is included as a systematic uncertainty of the cross-sections and the normalized yields. For the differential cross-sections, the variation of the template produces an additive change in the background, and thus, in the background-subtracted results. The corresponding systematic uncertainty effectively accounts for the sensitivity of the signal fractions to the uncertainty in the shape of the signal template.
• Excess background: As discussed in Section 4.5, the difference between the excess observed at large or ⊥ after HF background subtraction and the P +P 8 DY estimate is taken as a systematic uncertainty to cover possible unidentified background sources. It is evaluated separately in each kinematic bin in the measurement, except for the measurements as a function of | | and |Δ |, where the low pair yields at largest | | and |Δ | values preclude a statistically viable determination. For those measurements, the excess fraction is evaluated averaged over | | and |Δ | and then applied in each bin. This uncertainty is treated one-sided (downwards only), and is typically within 4% for the cross-section and 2% for the normalized yields.
The final uncertainties in the cross-sections are obtained by summing the individual uncertainties listed above in quadrature. The resulting uncertainty in the cross-sections vary slightly with centrality but are typically ≈8%.

Theoretical calculations
The results of the measurements presented here are compared with three sets of theoretical calculations: STAR , QED calculations based on a generalization of the EPA method, and calculations using a fully quantum-mechanical treatment of the photon transverse momentum distribution through Wigner distributions. For all three, the calculations rely on distributions of impact parameters populated by Pb+Pb collisions in each centrality interval that are obtained from a Glauber MC simulation [43] that also provides the basis for the ATLAS centrality calibration.
While STAR was originally designed for UPC processes, it was recently updated to allow the calculation of → + − processes over fixed impact parameter ranges, although there is no treatment of the impact parameter dependence of the photon transverse momenta. Since it has already been established that the non-UPC and ⊥ distributions disagree with STAR , only the predictions for the total dimuon production cross-sections are compared with data. ⊥ distributions. These results are referred to as "PWF" in the rest of this paper. The top panel of Figure 15 shows the measured → + − cross-sections as a function of centrality for four different¯T intervals. Also shown for comparison are calculations from the STAR event generator. The cross-sections decrease slightly from central to peripheral events. This decreasing trend is qualitatively reproduced by STAR , but the predicted cross-sections from STAR are considerably different from the measurements. It has been observed [4] that STAR predicts too low a → + − yield in UPC collisions. One possible explanation [58] for the STAR deficit in UPC dimuon production is its model's truncation of the single-nucleus photon flux at transverse distances smaller than the nuclear radius. That truncation may have greater impact on → + − processes in non-UPC nuclear collisions, where the neglected regions have larger geometric overlap. In the bottom panel of Figure 15, the measured cross-sections in three different centrality intervals show a rapid decrease with¯T.

Cross-sections and relative yield measurements
The total cross-sections for → + − production in the Fid-and Fid-⊥ fiducial regions but with no centrality or ZDC requirement after subtraction of the HF and DY backgrounds are Fid- The cross-sections in the two fiducial regions differ by only 1% because they are dominated by the ultraperipheral contribution, for which the and ⊥ restrictions remove only a small fraction of the pairs. These cross-sections are larger than that reported in the recent UPC measurement [4], primarily due to the fact that the present measurement does not impose a minimum mass requirement of 10 GeV and also because it includes → + − pairs in non-UPC events. The total cross-section obtained using the same fiducial selections as the UPC measurement is consistent with the result reported in Ref. [4].   Figure 16 shows the normalized yields of → + − pairs, , calculated according to Eq. (6), in various centrality intervals and for the four¯T ranges presented in Figure 15. For a given¯T interval, it is observed that the values increase from peripheral to central events, i.e. the fraction of → + − pairs produced in central events per centrality percentile is larger than in peripheral events. It is also observed that at higher¯T, a larger fraction of → + − pairs is produced in more central events. This is compensated for by a reduction of the normalized yields in the UPC interval (not shown in Figure 16), since by construction the normalized yields add up to one (see Eq. (6)). This result indicates that the dimuon transverse momentum spectra become harder with decreasing Pb+Pb collision impact parameter. Figure 17 shows the | | dependence of the normalized yields for two centrality and three¯T ranges. The normalized yields decrease with rapidity, approaching zero at | | = 2.4. This behavior primarily results from the acceptance of the two muons used to form the pair: | | < 2.4, which restricts the pair rapidity to | | < 2.4. However, the | | dependent trends in the normalized yields can be compared between the different¯T intervals, since the acceptance effects are similar for the different¯T intervals. For | |<1.2 a stronger¯T dependence is observed, with the normalized yields increasing with¯T. While at forward rapidities, the normalized yields are consistent, within uncertainties, for the three¯T intervals shown in Figure 17. Figures 18 and 19 show two of the primary results of this measurement: the differential cross-sections as a function of and ⊥ , respectively, in different Pb+Pb centrality intervals, including the ET0-ET3 and UPC intervals. The main observation in the previous ATLAS measurement of non-UPC → + − production [7], namely a centrality-dependent broadening of the distributions, is confirmed. A similar broadening is observed in the ⊥ distributions. It is interesting to note that the and ⊥ distributions in the most peripheral ET0-ET3 intervals are visibly broadened compared to the UPC interval.

and ⊥ distributions
In addition to the broadening, a centrality-dependent suppression of the dimuon yield at small and ⊥ values is seen in the data. This suppression was present, but not statistically significant, in the results of the previous measurement [7]. The suppression is greater in the ⊥ distributions for reasons that are made clear below. The PWF and QED calculations reproduce the measured distributions in the UPC interval reasonably well. The calculations also qualitatively reproduce the centrality-dependent broadening observed in the data. In addition, both calculations predict a centrality-dependent depletion at small (or ⊥ ). However, the depletion observed in the data is better reproduced by the QED calculations, which describe well the observed suppression in the distributions in the most central collisions; agreement in the middle centrality intervals is not as good. In Ref. [13], the depletion at small values of acoplanarity is attributed to quantum interference effects. It is not, as yet, understood why the PWF calculations show less suppression.
Assuming that the broadening and depletion in the or ⊥ distributions results from a physical process having an associated transverse momentum scale, T , then, as argued in Ref.
[7], the modifications of the distributions, naively, should scale with the ratio of that scale to the muon momentum, T /¯T. In contrast, the ⊥ distributions, which effectively probe the component of the dimuon ì T perpendicular to the dimuon axis, should only depend on the scale, T . To test this hypothesis, Figure 20 shows distributions of (left panels) and ⊥ (right panels) for three different T intervals and three different centrality intervals chosen to improve the statistical significance of the data. The distributions show a clear dependence on¯T such that at higher¯T the distributions become narrower and show less suppression at = 0. Thus, integrating over¯T will smooth out some of the centrality-dependent features in the data. In contrast, the ⊥ distributions show no significant dependence on¯T. Thus, an integration over¯T will have little or no impact on the shape of the ⊥ distributions. For this reason, the ⊥ variable should be preferred for future measurements and theoretical calculations.

Characterizing the centrality dependence of ⊥ distributions
In Ref. [7], the broadening of the acoplanarity distributions was characterized by a transverse momentum scale, obtained from the RMS of the distribution. The observation that the and the ⊥ distributions have the most-probable value shifted away from zero means that the RMS does not fully capture the modifications observed in the data. Nonetheless, the moments provide a model-independent means to quantify modifications of the and ⊥ distributions. Figure 21 presents moments of the ⊥ distribution as a function of centrality for the¯T > 4 GeV selection. The moments are calculated from the ⊥ distributions prior to the background subtraction and then corrected, analytically, to remove the contribution from the constant background. Results are shown for the mean, the RMS, and the standard deviation of the ⊥ distributions.
A significant increase in the mean is observed between the UPC and the four ET0-ET3 intervals, and then there is a further steady increase in the more central Pb+Pb collisions. Similar behavior is observed in the RMS values, but the standard deviation shows a much smaller increase between UPC collisions and the most central collisions. The measured moments are compared in the figure with those obtained from the PWF calculation. The PWF predictions reproduce many of the trends seen in the data, but the mean and RMS values systematically lie below the data.   Figure 22 shows moments of the distributions compared with both the PWF and QED calculations. The moments calculated for the QED and PWF predictions are in excellent agreement except for the most central collisions where the QED results are slightly higher. Both calculations show a systematic difference from the data in more peripheral collisions, especially the 60-70% and 70-80% intervals. However, in those intervals, both calculations show weaker suppression near = 0 than is observed in the data.

Rapidity gap and event-plane dependence of ⊥
A specific prediction from Ref. [9] is that if the centrality-dependent modifications of the distributions result from the muons being deflected in magnetic fields generated during the Pb+Pb collision, then the broadening should vary as the hyperbolic tangent of the rapidity difference between the two muons, |Δ | ≡ | 1 − 2 |, where 1 and 2 are the rapidities of the two muons in the pair. Figure 23 shows the ⊥ distributions for two centrality intervals: 20-40% and 40-80%, and for three different |Δ | ranges: |Δ | ≤ 1, 1 < |Δ | ≤ 2, and |Δ | > 2. The results are presented in self-normalized form to allow them to be compared directly, even though the yields vary with |Δ |. At most, a weak variation with |Δ | is observed, but with a dependence opposite to what would be 0 50 100 150 [  Figure 23: Dimuon ⊥ distributions in three intervals of |Δ | in the 20-40% (left) and 40-80% (right) centrality intervals, respectively. The distributions are self-normalized to allow the distributions to be directly compared. The error bars show combined statistical and systematic uncertainties. expected from magnetic field effects. Namely, the suppression near ⊥ = 0 is greater for smaller |Δ |, while a tanh |Δ | dependence would have the opposite behavior and should vary more rapidly with |Δ |.
A second way to observe potential effects of magnetic fields on the → + − pairs is to study the dependence of the ⊥ distributions on the orientation of the dimuons in the transverse plane relative to the direction of the second-order event plane. In Pb+Pb collisions of intermediate centrality, the second-order event plane angle, Ψ 2 , is understood to be well-correlated with the direction of the impact parameter vector, which is, in turn, correlated with the direction of any magnetic fields generated during a Pb+Pb collision.
Dimuon yields were measured as a function of the quantity |2Δ | ≡ 2 − Ψ 2 , where represents the azimuthal orientation of the dimuon, and Ψ 2 is the second-order event plane angle obtained from the FCal using methods applied in previous ATLAS flow analyses [59,60]. The angle, , is calculated by rotating one of the muon azimuthal angles by and then averaging the result with the other muon, = 1 2 ( 1 + + 2 ) . Figure 24 shows differential cross-sections as a function of ⊥ in two intervals of |2Δ |, |2Δ | < /2 and /2 ≤ |2Δ | < , for two centrality intervals. A significant variation with |2Δ | is seen in the depletion near ⊥ = 0 for the 20-40% centrality interval. Otherwise, the distributions are indistinguishable at larger ⊥ values, which means that the overall broadening of the dimuon ⊥ distributions is the same in the two |2Δ | intervals.
To further explicate the results presented in Figures 23 and 24, Figure 25 shows results for the mean and standard deviation of the dimuon ⊥ distributions as a function of |Δ | (top) and |2Δ | (bottom) using finer bins in both variables. The average ⊥ values have, at most, a weak dependence on |Δ | or |2Δ |, while the standard deviations of the ⊥ distributions are constant within uncertainties. Thus, the mechanism responsible for the broadening of the ⊥ distributions does not appear to depend on either |Δ | or the direction of the muons relative to the impact parameter vector in the transverse plane. The absence of variation with |Δ | by values > 4 MeV of either the mean or standard deviation of the ⊥ distributions rules out magnetic broadening as a significant contribution to the observed modifications of the ⊥ distributions. However, the mechanism responsible for the suppression at ⊥ = 0 may vary with |Δ | or |2Δ |. Indeed, the ⊥ distributions in Figures 23  and 24 indicate such dependence, although the suppression near ⊥ = 0 is not easily seen in the calculated moments.   Figure 24: Differential cross-sections for → + − production as a function of ⊥ for two intervals of |2Δ |: 0 < |2Δ | < /2 and /2 < |2Δ | < in the 20-40% (top) and 40-80% (bottom) centrality intervals, respectively. The error bars show combined statistical and systematic uncertainties.

Conclusion
ATLAS has measured dimuon production via scattering processes in non-ultraperipheral Pb+Pb collisions at √ NN = 5.02 TeV. The measurements use data from the 2015 and 2018 Pb+Pb runs at the LHC corresponding to an integrated luminosity of 1.93 nb −1 . Backgrounds, dominated by heavy-flavor decays, are evaluated using template fits to the distribution of muon-pair 0 values. A much smaller background from DY processes, estimated using P +P 8 calculations implemented with nNNPDF2.0 nuclear PDFs, is subtracted from the data. Cross-sections and normalized yields of → + − pairs are measured as a function of pair rapidity,¯T, and centrality. The cross-sections vary weakly with centrality, decreasing from central to peripheral collisions. The STAR model, which was recently augmented to allow evaluation of cross-sections for (e.g.) → + − production within restricted impact parameter intervals, substantially underestimates the measured cross-sections.
Measurements of the and ⊥ distributions show a significant centrality dependence consistent with the results in Ref. [7]. However, with the improved statistical precision of this measurement, an additional depletion is observed in the acoplanarity and ⊥ distributions near zero values of the corresponding quantities.
The distributions are observed to vary strongly with¯T, while distributions of ⊥ are observed to be approximately independent of¯T, as would be expected if the broadening is determined by a physics effect with an intrinsic transverse momentum scale. The ⊥ distributions are, thus, better suited for studying the modifications of the dimuon alignment if the study integrates over any finite interval of¯T.
Comparisons of the measured and ⊥ distributions with two theoretical calculations are presented. Both the QED and the photon Wigner function-based calculations quantitatively reproduce the centrality-dependent broadening observed in the data. While both calculations also predict a suppression near = 0 and/or ⊥ = 0, the QED calculations describe the data more accurately; the PWF results show significantly less depletion than observed in the data.
Measurements of the ⊥ distributions as a function of |Δ | and |2Δ | suggest that the broadening of the and ⊥ distributions does not have significant contribution from interactions of the muons with magnetic fields generated in the quark-gluon plasma. Specifically, the data show no dependence of the broadening on |Δ |, while a tanh |Δ | dependence is predicted. Similarly, no dependence of the standard deviation of the ⊥ distributions on the orientation of the muon pair relative to the second-order event plane is observed. This result would be inconsistent with magnetic broadening in a field whose direction in the transverse plane is correlated with the impact parameter vector. While the overall broadening of the distributions does not depend on |Δ | or |2Δ |, the magnitude of the depletion at small ⊥ appears to depend on both quantities.      [46] S. Baker and R.