Production of $\varUpsilon(\textrm{nS})$ mesons in Pb+Pb and $pp$ collisions at 5.02 TeV

A measurement of the production of vector bottomonium states, $\varUpsilon(\textrm{1S})$, $\varUpsilon(\textrm{2S})$, and $\varUpsilon(\textrm{3S})$, in Pb+Pb and $pp$ collisions at a center-of-mass energy per nucleon pair of 5.02 TeV is presented. The data correspond to integrated luminosities of 1.38~$\textrm{nb}^{-1}$ of Pb+Pb data collected in 2018, 0.44 $\textrm{nb}^{-1}$ of Pb+Pb data collected in 2015, and 0.26 $\textrm{fb}^{-1}$ of $pp$ data collected in 2017 by the ATLAS detector at the Large Hadron Collider. The measurements are performed in the dimuon decay channel for transverse momentum $p_\textrm{T}^{\mu\mu}<30$~ GeV, absolute rapidity $|y^{\mu\mu}|<1.5$, and Pb+Pb event centrality 0-80%. The production rates of the three bottomonium states in Pb+Pb collisions are compared with those in $pp$ collisions to extract the nuclear modification factors as functions of event centrality, $p_\textrm{T}^{\mu\mu}$, and $|y^{\mu\mu}|$. In addition, the suppression of the excited states relative to the ground state is studied. The results are compared with theoretical model calculations.


Introduction
Quantum chromodynamics (QCD) predicts that at high temperatures and energy densities, hadronic matter undergoes a phase transition and turns into a state of deconfined quarks and gluons known as quark-gluon plasma (QGP). This state of matter is typically thought to be created in the collisions of two heavy nuclei at ultra-relativistic energies. In such collisions, heavy-flavor quarks, especially charm and bottom, are produced at an early stage in hard scattering processes and hence can probe QGP over its full evolution.
Formation of the QGP and the consequent modification to the heavy-quark potential is expected to lead to different quarkonium states dissolving at different temperatures of the medium [1]. This effect is known as sequential suppression [2]. While the excited states are dissociated just above the transition temperature

Data selection and simulation samples
The results presented in this paper were obtained using data recorded in 2017 at a center-of-mass energy of 5.02 TeV as well as Pb+Pb data collected in 2015 and 2018 at 5.02 TeV per nucleon-nucleon pair.
The integrated luminosity of the analyzed collision samples is 0.26 fb −1 . The events were collected using a dimuon trigger which requires at least two spatially separated muon candidates at L1, while both satisfy the criterion of T > 4 GeV in the HLT. For the Pb+Pb analysis, the integrated luminosity is 0.44 nb −1 for the 2015 data sample and 1.38 nb −1 for the 2018 data sample. Pb+Pb events were collected using triggers which require at least one muon with T > 4 GeV at both L1 and the HLT, and at least one additional muon satisfying T > 4 GeV in the HLT, without requiring matching to L1. Muon pairs were required to fulfill the following criteria: at least one reconstructed muon matching the HLT's dimuon trigger, and both muons matching the HLT without an L1 trigger requirement; both muons satisfy the Medium identification criteria (described in Ref. [29]) without any requirements on TRT hits; both muons are associated with the primary vertex reconstructed using all of the tracks in each event; the muon pair has a dimuon mass 7.7 < < 12.3 GeV, dimuon transverse momentum | T | < 30 GeV and dimuon rapidity | | < 1.5; the selected muon pair is re-fitted to a common vertex, with the vertex fit quality satisfying 2 < 100 and the significance of the transverse displacement of the re-fitted vertex relative to the primary vertex satisfying | / ( )| < 3, where ( ) is the primary vertex resolution. The dimuon rapidity requirement was chosen to be | | < 1.5 because at more forward rapidities the mass resolution starts to deteriorate quickly.
Monte Carlo (MC) simulation is used to study the (nS) acceptance, the fit model used in signal extraction, and the closure of muon reconstruction and identification corrections (due to residual biases associated with the yield correction procedure). The unpolarized prompt ( ) in pp events were generated with the CTEQ6L1 [30] parton distribution function set in Pythia 8 [31] with subsequent decay in muon pairs. Pythia 8 implements prompt (nS) production sub-processes using the non-relativistic QCD Color Octet mechanism [32]. Prompt (nS) production includes prompt production from the hard interactions, as well as the radiative feed-down from → (nS) decays. The production of nonprompt / , used for determining per-muon corrections, was simulated in Pythia 8 by forcing¯production and and retaining only events consistent with / decays. The Pb+Pb MC sample was created by overlaying simulated Pythia 8 events with recorded minimum-bias Pb+Pb events, so that the "data overlay" simulation samples contain the same level of underlying-event activity as is present in the Pb+Pb data. The response of the ATLAS detector was simulated [33] using Geant4 [34]. The MC events are reconstructed with the same algorithms as used in data.

.1 Centrality definition in Pb+Pb
The transverse energy measured in the forward calorimeter, Σ FCal T , in minimum bias events, is used to estimate the degree of overlap between the two colliding Pb nuclei. Each centrality class corresponds to a fixed percentile in the Σ FCal T distribution of minimum-bias events using the procedure described in Ref. [21], where a table containing exact centrality bin definitions can be found. A Monte Carlo Glauber-based model [35] is used to calculate the mean number of participant nucleons, ⟨ part ⟩, and the mean nuclear overlap function, ⟨ AA ⟩, for each centrality class.

Corrections to raw invariant mass distributions
Before the dimuon invariant mass distributions are fit to extract Upsilon yields, each candidate dimuon pair is corrected with a weight that accounts for the (nS) dimuon acceptance, trigger efficiency, and reconstruction efficiency.
The kinematic acceptance A is defined as the probability that both muons from → + − decay pass the fiducial selection ( T > 4 GeV and | | < 2.4). The kinematic acceptance is calculated from a generator-level simulation separately for different (nS) states as described in Ref. [36]. The acceptance correction also accounts for final-state radiation from one or both of the decay muons. In principle, the acceptance could depend on the spin-alignment of the (nS). In this analysis, (nS) mesons are assumed to be produced unpolarized, following the previous measurements in collisions [37][38][39]. These measurements are consistent with no polarization, but have large uncertainties. No extra systematic uncertainty due to polarization is added in this study.
The dimuon reconstruction efficiency, reco ( 1 2 ), is determined as the product of two single-muon reconstruction efficiencies. The single-muon reconstruction efficiency is factorized into ID track reconstruction efficiency and MS reconstruction efficiency. For collisions, the values of the reconstruction efficiency are obtained from the / → + − Pythia 8 simulation, and additional data-to-MC efficiency scale factors are derived using a / → + − tag-and-probe method [29] using data to account for residual differences between data and simulation. The same / → + − tag-and-probe method is employed to measure the muon reconstruction efficiency in Pb+Pb collisions. The ID reconstruction efficiency is obtained directly from Pb+Pb data, with a requirement on the transverse displacement of the / vertex to suppress potential biases from displaced muons. The MS reconstruction efficiency is obtained from the / → + − Pythia 8 simulation overlaid with minimum-bias Pb+Pb data, and additional data-to-MC scale factors are determined in Pb+Pb data to account for the small difference between data and simulation. The ID reconstruction efficiency is found to be larger than 99% in both and Pb+Pb collisions, with a weak centrality dependence at low T in the latter case. The MS reconstruction efficiency at T = 4 GeV is about 65% in the barrel region (| | < 1.05) and 75% in the endcap region (1.05 < | | < 2.4), and the MS efficiency increases with T and saturates at 95% around T = 7 GeV in both the barrel and endcap regions. Due to the absorption of most hadronic activity in the ATLAS calorimeters, the MS reconstruction efficiency has no centrality dependence. For a given muon pair, the dimuon trigger efficiency in collisions is factorized as the product of two single-muon efficiencies. The dimuon trigger efficiency in Pb+Pb collisions is the combined efficiency of either of the two muons matching HLT trigger while the other one matching HLT trigger without L1 trigger requirement.
The single-muon trigger efficiency is determined in data and simulation using a / → + − tag-and-probe method similar to that used for the muon reconstruction efficiency. Two different muon-trigger logic schemes are used in this analysis: 1) a full-chain muon trigger, which requires the formation of an L1 muon candidate that is subsequently confirmed in the HLT, and 2) a full-scan muon trigger, which is only performed in the HLT via a muon candidate search of the full MS system without any requirement at L1. The values of the full-chain muon trigger efficiency are determined from MC simulation, with data-to-MC scale factors determined from the data to take into account the difference between data and simulation. The same values are used for and Pb+Pb collisions, except that an additional centrality-dependent correction is applied to Pb+Pb data. The centrality-dependent correction is determined as the ratio of the trigger efficiency measured in Pb+Pb collisions as a function of centrality to the trigger efficiency in data. The full-chain muon trigger efficiency plateau value is found to be 70% in the barrel and 90% in the endcaps. The centrality-dependent correction factor is about 90% (100%) for the 0-10% (60-80%) centrality range. The full-scan muon trigger is only used in Pb+Pb collisions and its efficiency is determined directly from Pb+Pb data, using the tag-and-probe method. The full-scan trigger plateau is found to be 90% in both the barrel and endcap regions. The dimuon trigger efficiency in collisions is factorized as the product of two full-chain muon trigger efficiencies, and in Pb+Pb collisions the factorization form consisting of full-chain and full-scan muon trigger efficiencies, as detailed in Ref.

Upsilon signal extraction
Upsilon states are reconstructed in the + − decay channel and their yields are determined via unbinned maximum-likelihood fits to the weighted dimuon invariant mass distributions, following the same procedure for both and Pb+Pb data. Each of the three (nS) state signal shapes is described by a sum of Crystal Ball [40] and Gaussian functions.
The probability distribution function for the fit is defined as a normalized sum of three signal components and a background component as where (nS) ( ) = G ( ; nS , nS ) + (1 − ) CB ( ; nS , 1.7 nS , , ); (nS) and bkg are the Upsilon and background yields, respectively; and G and CB are Gaussian and Crystal Ball functions, respectively, with representing relative weight of the Gaussian function. The quantities nS and nS are the Gaussian function's mean and width for each Upsilon state, and and are CB tail parameters. The line shape parameters are assumed to be the same for all three states since it is mostly determined by detector effects. The background shape, bkg ( ), is represented by a second-order polynomial at high T ( T > 6 GeV) and as a product of an error function and an exponential function at low T . The determination of the yield corrected for acceptance and efficiencies, corr (nS) , proceeds in several steps. First, a resonance-dependent weight, total ( (nS)), is determined for each selected dimuon candidate as where A ( (nS)) is the acceptance for (nS) → + − decay, reco is the muon reconstruction efficiency, trig is the trigger efficiency, and pvAsso is the efficiency related to the primary-vertex association. Next, an unbinned maximum-likelihood fit to the weighted dimuon invariant mass distribution ( ) is performed to extract the (nS) yields. Three fits with an acceptance value corresponding to each state are performed to extract the yields for the three states.
The mean and Gaussian width of the (1S) signal, 1 and 1 , are left unconstrained in each T , | | and centrality range, while the means and widths of (2S) and (3S) in the same range are fixed to (1S) parameters scaled by the respective PDG [41] mass ratios. The widths of the CB functions are set by scaling the corresponding Gaussian widths by a constant factor, 1.7, which was determined in a previous analysis [42] and validated with MC studies. The relative weight of the Gaussian and CB functions, , is not constrained, and has the same value for all three Upsilon states. For the collision analysis, the CB function parameter , which defines the point at which the low-mass tail transitions from a Gaussian shape to a power-law shape, was fixed to the value obtained from (1S) MC samples, while , which describes the shape of the tail, was a free parameter. For Pb+Pb collisions, both and were fixed to the values from the fit for collisions, for each kinematic selection. The nominal signal fit model described above was validated by fitting the signal MC samples in various T , | | and centrality intervals.
The background functional form varies with dimuon T . For T > 6 GeV the background is parameterized as a second-order polynomial. However, for lower T ( T < 6 GeV), and for the integrated over T fits, in order to describe the turn-on behavior of the dimuon acceptance caused by the single-muon transverse momentum requirement, an error function multiplied by an exponential function is used. The background model parameters are initialized using a background-enriched sample, which consists of events with same-sign muon pairs, as well as a control sample consisting primarily of muons from b-hadron decays. The control sample is produced by requiring at least one of the two muons to satisfy | 0 |/ 0 > 2 or | 0 sin( )| > 0.2 mm, where 0 and 0 are the distances of closest approach of the muon to the primary vertex in the plane perpendicular to the beam and in the beam direction, respectively. After using the background distributions to determine the relevant parametrizations, the full model including signal and background contributions is used to fit the data, with the slope of the exponential function allowed to float. Figure 1 shows an example of the fit to the dimuon mass plots for (left) and 0-80% centrality Pb+Pb (right) collisions for the inclusive T and | | selection. The lower panels show the pull distribution, which represents the distance between data points and fit function normalized by the data points statistical uncertainty: (Data − Fit)/ (Data). The goodness of the fit is assessed by calculating the reduced chisquare, 2 /NDF, summing the squared deviations of the data points from the fit (weighted by the inverse errors) and then dividing by the number of degrees of freedom (NDF) in the fit. Typically, 2 /NDF for the fits varies from ≈2.5 to ≈1 for and from ≈2 to ≈1 for Pb+Pb, indicating that the data and model agree within statistical uncertainties. The extracted values of 2 /NDF decrease with T , indicating that the fit quality improves with T . Some relatively large 2 /NDF values are also found and can be explained by deviations of the functions used for background description from the actual background at the edges of the fit range.

Systematic uncertainties
The main sources of point-by-point uncorreleated systematic uncertainties pertain to the corrections for muon reconstruction and trigger efficiencies, and the yield extraction. For the lowest T range, the final-state radiation correction is also a significant contribution. Other, subdominant sources of systematic uncertainty considered in this study are the primary-vertex association uncertainty and bin migration due to momentum resolution. Finally, the Upsilon states are assumed to be unpolarized and no extra systematic uncertainty is assigned to cover this assumption.
The systematic uncertainty in the MS reconstruction efficiency is dominated by the uncertainty in the data-to-MC scale-factor determination. The scale-factor uncertainty is evaluated following Ref. [29] by changing the tag-muon selection criteria and varying the line-shapes in the efficiency extraction fit procedure in the data. The uncertainty related to the ID reconstruction efficiency, which is close to 1, is estimated by comparing the results while varying this efficiency in both up and down directions by one standard deviation.
The systematic uncertainty in the muon trigger efficiency is also dominated by the tag-and-probe efficiency determination procedure. For Pb+Pb collisions, an additional systematic uncertainty associated with the centrality-dependent correction is included. This uncertainty is evaluated by comparing the centralitydependence-corrected Pb+Pb efficiency with the efficiency as a function of . Individual variations described above are added in quadrature to form the total systematic uncertainty of the efficiency corrections.
The sensitivity of the signal extraction to the choice of a particular fit model is evaluated by varying the line-shape of each fit component. The maximum variation between the recalculated values and the central value is used to estimate their uncertainty. Eight variations are considered, and these can be categorized into three groups: signal resolution (width of the peak), shape of the final-state radiation tail, and background shape. The final uncertainty from the fit model is obtained by calculating the difference between the maximum and the minimum yield from the eight line-shape variations and dividing by √ 12, assuming a flat distribution.
The uncertainty in the acceptance due to the final-state radiation correction was calculated by comparing the result of a fully simulated acceptance calculation with that obtained using an MC sample designed for a high-precision determination of the acceptance. The signal Monte Carlo samples were processed with a fast simulation [33] which relies on a parameterization of the calorimeter response [43]. This uncertainty is only important for the lowest T range. The global uncertainty of the integrated luminosity for the 2017 data is 1.6%, derived using methods described in Ref. [44]. Primary vertex association uncertainty, which results mostly from small discrepancies between data and MC, was studied by varying the primary vertex association requirements. Since the primary-vertex association affects all Upsilon state yields in the same way, it is treated as a global uncertainty together with those of the luminosity and AA . The combined systematic uncertainty for the luminosity and primary-vertex association in data is 2.6%. For Pb+Pb collisions, the global systematic uncertainty of ⟨ AA ⟩ is estimated by varying the Glauber model parameters as detailed in Ref. [45]. The combined systematic uncertainty for ⟨ AA ⟩ and primary-vertex association in Pb+Pb collisions is 3.7%.

Systematic uncertainties in
and Pb+Pb collisions are summarized in Table 1. While some systematic uncertainties for AA values and excited-state to ground-state double ratios are correlated (e.g. the acceptance) and cancel out in the ratios, most of the systematic uncertainties are not completely correlated and are estimated by directly studying their effects on the ratios.

Differential cross-section
Differential (nS) production cross-sections in collisions are measured according to the relation where B ( (nS) → + − ) is the dimuon decay branching fraction, corr (nS) is the (nS) yield corrected for acceptance and efficiencies, Δ T and Δ are the bin widths in T and , and ∫ L is the integrated luminosity.
The (nS) differential cross-sections in collisions at 5.02 TeV, multiplied by the respective dimuon branching fractions, are shown as a function of T in the left panel of Figure 2. The per-event yields of (nS) states in Pb+Pb collisions are defined by where evt is the total number of minimum-bias Pb+Pb collisions in each centrality class. In particular, this number is 1.02 × 10 9 for the 0-10% centrality interval. Per-event Upsilon yields in Pb+Pb collisions divided by ⟨ AA ⟩ are shown in the right panel of Figure 2. The results for (3S) mesons are not shown because their peaks are not statistically significant in Pb+Pb collisions.

Nuclear modification factor
The modifications of bottomonium production yields in Pb+Pb collisions relative to the system are quantified by the nuclear modification factor AA , which can be defined for each centrality interval as where AA is the observed per-event yield of bottomonium states in Pb+Pb collisions, and is the bottomonium production cross-section in collisions at the same collision energy. Figure 3 shows the AA values of (nS) as functions of ⟨ part ⟩ (top), dimuon T (bottom left), and | | (bottom right). The centrality-integrated results are also shown in the right panel of the top plot. In addition to the results for (1S) and (2S), only the combined result for the two excited states, (2S+3S), is presented because the (3S) peak is not statistically significant in the Pb+Pb data. The (nS) states are observed to be suppressed over the whole kinematic range investigated, and the AA values of (2S) and (2S+3S) are always lower than those of (1S). The AA value decreases with ⟨ part ⟩ for all three states. No strong T or | | dependence is observed. When no statistically significant non-zero yield was extracted for a particular kinematic selection, 95% confidence level upper limit was calculated.

Excited-state to ground-state double ratios
The suppression of different Upsilon states can be compared by constructing an excited-state to ground-state double ratio of nuclear modification factors. The advantage of measuring the double ratios is that the  acceptance and efficiency corrections partially cancel out, and the overall systematic uncertainty is reduced. Although defined in terms of the individual nuclear suppression factors, the double ratio can be understood as being defined as the ratio of the yields of excited states (2S), (3S) or of combined yield of the two excited states ( (2 + 3 )) to the yield of the ground state (1S) in Pb+Pb collisions, divided by the same ratio in collisions.       [8]. This model uses a framework with coupled transport equations for open heavy-flavor and quarkonium states in order to describe their transport inside the quark-gluon plasma, including regeneration. Cold nuclear matter effects are included by using nuclear parton distribution functions for the initial primordial heavy-flavor production. A calibrated (2 + 1)-dimensional viscous hydrodynamic model is used to describe the bulk QCD medium. The model depends on the choice of nucleus parton distribution function (nPDF) and two coupling constant parameters, s and pot s , which are varied by ±10% from their nominal values.

Theory comparisons and discussion
All three models are in agreement with the data within experimental and theoretical uncertainties. It is notable that all three theoretical results were calculated after the publication of a CMS measurement of suppression in Ref.
[11] at the same beam energy. Figure 8 shows the ATLAS results presented in this paper compared with CMS published values [12] for the (1S) (left) and (2S) (right) states. These measured nuclear modification factors, as well as the double ratios, are consistent between ATLAS and CMS within uncertainties. The kinematic selections are similar, though ATLAS uses a narrow rapidity window | | < 1.5 compared with | | < 2.4 in CMS.
Published results by the ALICE Collaboration [13,14] are at more forward rapidities, and thus no direct comparison is made. The sequential melting of Upsilon states is thus confirmed by the present measurement. Suppression of the (1S) observed at the level of approximately 0.25-0.30 in central Pb+Pb collisions, indicates suppression of the directly produced (1S) and not just suppression of higher states feeding down to the (1S) thus providing important information on the highest temperatures achieved in the QGP.
More recently, CMS presented a measurement which includes the first observation of the (3S) in Pb+Pb collisions at the LHC [48]. The (3S) is even more suppressed than the (1S) and (2S), further quantifying the sequential melting of quarkonium states in QGP. The combined nuclear modification factor for the two excited states (2S+3S) presented in this paper is consistent with the (3S) CMS measurement.

Conclusions
This paper presents a measurement of (nS) yields in and Pb+Pb data at 5.02 TeV per nucleon-nucleon collision. The measurement uses data from collisions collected in 2017 with a total integrated luminosity of 0.26 fb −1 and Pb+Pb collisions collected in 2015 and 2018 with total integrated luminosities of 0.44 nb −1 and 1.38 nb −1 respectively, recorded by the ATLAS experiment at the LHC. The and Pb+Pb measurements are used to obtain the nuclear modification factor for (1S), (2S), and (2S+3S), as well as excited-state to ground-state double ratios of AA for (2S) and (2S+3S) as functions of transverse momentum, rapidity, and centrality. Both the (1S) and (2S) yields are suppressed with increasing centrality in Pb+Pb compared to those in collisions, and the excited states are shown to be stronger suppressed than the ground state, resulting in double ratios smaller than one. The AA value for (2S+3S) appears to be systematically lower than that for (2S), and so are the corresponding double ratios, which indicates that the (3S) state is suppressed more than (2S).
The measured nuclear modification factors and double ratios are found to be consistent with the previous CMS measurement [12] at similar kinematical conditions. All theoretical calculations considered in this paper describe the data well and incorporate deconfinement as a key ingredient in the suppression of the Upsilon yields. Further discrimination between the different implementations of these deconfinement effects require additional precision data from upcoming runs at RHIC and the LHC.
The crucial computing support from all WLCG partners is acknowledged gratefully, in particular from CERN, the ATLAS Tier-1 facilities at TRIUMF (Canada)