Measurement of B$^0_\mathrm{s}$ meson production in pp and PbPb collisions at $\sqrt{s_\mathrm{NN}} =$ 5.02 TeV

The production cross sections of B$^0_\mathrm{s}$ mesons and charge conjugates are measured in proton-proton (pp) and PbPb collisions via the exclusive decay channel B$^0_\mathrm{s}\to\mu^+\mu^-$K$^+$K$^-$ at a centre-of-mass energy of 5.02 TeV per nucleon pair and within the rapidity range $|y|\lt$ 2.4 using the CMS detector at the LHC. The pp measurement is performed as a function of transverse momentum (p$_\mathrm{T}$) of the B$^0_\mathrm{s}$ meson in the range of 7 to 50 GeV/$c$ and is compared to the predictions of perturbative QCD calculations. The B$^0_\mathrm{s}$ production yield in PbPb collisions is measured in two p$_\mathrm{T}$ intervals, 7 to 15 and 15 to 50 GeV/$c$, and compared to the yield in pp collisions in the same kinematic region. The nuclear modification factor ($R_\mathrm{AA}$) is found to be 1.5 $\pm$ 0.6 (stat) $\pm$ 0.5 (syst) for 7-15 GeV/$c$, and 0.87 $\pm$ 0.30 (stat) $\pm$ 0.17 (syst) for 15-50 GeV/$c$, respectively. Within current uncertainties, the results are consistent with models of strangeness enhancement and a suppression as observed for the B$^+$ mesons.


Introduction
Relativistic heavy ion collisions allow the study of quantum chromodynamics (QCD) at high energy density and temperature. Under such extreme conditions, a state consisting of deconfined quarks and gluons, the quark-gluon plasma (QGP) [1,2], is predicted by lattice QCD calculations [3]. The study of the phenomenon in which the outgoing partons interact strongly with the QGP and lose energy by means of elastic collisions and mediuminduced gluon radiation [4][5][6][7][8] can provide insights into the energy density and diffusion properties of the QGP. Heavy quarks are effective probes to study these properties of the medium. Charm and beauty quarks that are primarily produced in hard scatterings at the early stages of the collision are expected to carry the full evolution history of the QGP formation [8]. On the other hand it is expected [9] that, via the process gg → ss, an enhancement of strangeness in a thermally and chemically equilibrated QGP should occur if its temperature is above the strange quark mass. Measurements at the BNL RHIC of the production of strange baryons and mesons, using different collision systems and beam energies, provide systematic support for this expectation [10][11][12][13][14]. Because of the interplay between the predicted enhancement of strange quark production and the quenching mechanism of beauty quarks, the measurement of strange beauty particles is important for studying the mechanisms of beauty hadronization in heavy ion colli-E-mail address: cms -publication -committee -chair @cern .ch. sions. In the presence of a medium with increased strangeness content [15,16], the relative yield of B 0 s mesons with respect to nonstrange beauty mesons at transverse momentum (p T ) below ∼15 GeV/c [8,17] can be enhanced in nucleus-nucleus collisions compared to proton-proton (pp) interactions. This can happen if recombination is a significant factor of beauty hadronization in the QGP [18][19][20]. The recombination processes, which are considered markers for the presence of a deconfined medium, were most recently tested in the open charm sector by the ALICE Collaboration [21]. A possible hint for an enhancement in the relative yield of D + s mesons with respect to nonstrange charmed mesons for p T < 8 GeV/c in central PbPb collisions at a center-of-mass energy of √ s NN = 5.02 TeV per nucleon pair was observed.
The production of B 0 s mesons was previously measured at the CERN LHC by the CMS Collaboration in pp collisions at a center-ofmass energy of √ s = 7 TeV [22] and in proton-lead (pPb) collisions at √ s NN = 5.02 TeV [23]. In this letter, we report the first measurement of exclusive B 0 s meson decays ever performed in nucleusnucleus collisions and in pp collisions at 5.02 TeV. The pp measurement is performed as a function of p T and compared to the predictions of fixed-order plus next-to-leading order logarithmic (FONLL) perturbative QCD calculations [24][25][26]. The nuclear modification factor (R AA ) of B 0 s mesons, which is defined as the ratio of the yield in PbPb collisions with respect to that in pp collisions scaled by the corresponding number of binary nucleon-nucleon (NN) collisions, is shown. The comparison between the R AA of B 0 s mesons and that of B + mesons measured by CMS at the same energy [27] is also presented. The B 0 s meson and its charge conjugate are measured in the rapidity range |y| < 2.4 via the reconstruction of the decay channel B 0 s → J/ψφ → μ + μ − K + K − , which has the branching fraction B = (3.12 ± 0.24) × 10 −5 [28]. The pp measurement is performed as a function of the B 0 s p T in three intervals, 7-15, 15-20, and 20-50 GeV/c. The PbPb production yield and the R AA measurement are performed in two p T intervals, 7-15 and 15-50 GeV/c, inclusively for all events (i.e., 0-100% centrality, the degree of overlap of the two colliding nuclei). Throughout the letter, unless otherwise specified, the y and p T variables given are those of the B 0 s mesons. This analysis does not distinguish between the charge conjugates.

Experimental apparatus and data sample
The central feature of the CMS detector is a superconducting solenoid, which provides a magnetic field of 3.8 T. Within the solenoid volume are a silicon tracker that measures charged particles in the pseudorapidity range |η| < 2.5, a lead tungstate crystal electromagnetic calorimeter, and a brass and scintillator hadron calorimeter. For charged particles of 1 < p T < 10 GeV/c and |η| < 1.4, the track resolutions are typically 1.5% in p T and 25-90  μm in the transverse (longitudinal) impact parameter [29]. Muons are measured in the range |η| < 2.4, with detection planes made using three technologies: drift tubes, cathode strip chambers, and resistive-plate chambers. The muon reconstruction algorithm starts by finding tracks in the muon detectors, which are then fitted together with tracks reconstructed in the silicon tracker to form "global muons". Matching muons to tracks measured in the silicon tracker results in a relative p T resolution for muons with 20 < p T < 100 GeV/c of 1.3-2.0% in the barrel (|η| < 1.2) and better than 6% in the endcaps (1.6 < |η| < 2.4). For muons with higher p T up to 1 TeV/c, the p T resolution in the barrel is better than 10% [30]. The hadron forward (HF) calorimeter uses steel as an absorber and quartz fibers as the sensitive material. The two halves of the HF are located 11.2 m away from the interaction point, one on each end, providing together coverage in the range 3.0 < |η| < 5.2. In this analysis, the HF information is used for performing an offline event selection. A detailed description of the CMS experiment and coordinate system can be found in Ref. [31].
Several Monte Carlo (MC) simulated event samples are used to evaluate background components, signal efficiencies, and detector acceptance corrections. The simulations include samples containing only the B 0 s meson decay channels being measured, and samples with inclusive (prompt and nonprompt) J/ψ mesons. Proton-proton collisions are generated with pythia8 v212 [32] tune CUETP8M1 [33] and propagated through the CMS detector using the Geant4 package [34]. The decay of the B 0 s mesons is modeled with evtgen 1.3.0 [35], and final-state photon radiation in the B 0 s decays is simulated with photos 2.0 [36]. For the PbPb MC samples, each pythia8 event is embedded into a PbPb collision event generated with hydjet 1.8 [37], which is tuned to reproduce global event properties, such as the charged-hadron p T spectrum and particle multiplicity. For both samples, the signal p T shape is reweighted to match the one from FONLL. For both p p and PbPb data and MC samples, the dimuon and ditrack mass distributions/resolutions are consistent.
Events were collected with the same trigger during the pp and PbPb data acquisition, requiring the presence of two muon candidates (with no explicit momentum threshold) in coincidence with a bunch crossing. For the offline analysis, events have to pass a set of selection criteria designed to reject events from background processes (beam-gas collisions and beam scraping events) as described in Ref. [38]. Events are required to have at least one reconstructed primary interaction vertex, formed by two or more tracks, with a distance from the center of the nominal interaction region of less than 15 cm along the beam axis. In PbPb collisions, the shapes of the clusters in the pixel detector have to be compatible with those expected from particles produced by a PbPb collision [39]. In order to select hadronic collisions, the PbPb events are also required to have at least three towers in each of the HF detectors with energy deposits of more than 3 GeV per tower. The combined efficiency for this event selection, including the remaining non-hadronic contamination, is (99 ± 2)%. Values higher than 100% are possible, reflecting the potential presence of ultra-peripheral (i.e., non-hadronic) collisions in the selected event sample. The PbPb sample corresponds to an integrated luminosity of approximately 351 μb −1 . This value is indicative only, as the PbPb yield is normalized by the total number of minimum bias events sampled, N MB [38]. The pp data set corresponds to an integrated luminosity of 28.0 pb −1 , which is known to an accuracy of ±2.3% from the uncertainty in the calibration based on a van der Meer scan [40]. The average number of additional collisions per bunch crossing is approximately 0.9 for pp and less than 0.01 for PbPb data. The presence of multiple collisions is found to have a negligible effect on the measurement.

Signal extraction
The analysis procedure is common for pp and PbPb data. Kinematic limits are imposed on the single muons so that their reconstruction efficiency stays above 10%. These limits are p μ T > 3.5 GeV/c for |η μ | < 1.2, p μ T > 1.8 GeV/c for 2.1 ≤ |η μ | < 2.4, and linearly interpolated in the 1.2 < |η μ | < 2.1 region. The muons are also required to match the muons that triggered the event online, and to pass selection criteria optimized for low p T (the so-called soft selection [30]). Two muons of opposite sign (OS), with an invariant mass within ±150 MeV/c 2 of the world-average J/ψ meson mass [28] are selected to reconstruct a J/ψ candidate, with a mass resolution of typically 18-55 MeV/c 2 , depending on the dimuon rapidity and p T . The OS muon pairs are fitted with a common vertex constraint and are kept if the p-value of the χ 2 of the fit is greater than 1%, thus lowering the background from charm and beauty hadron semileptonic decays. Similarly, the φ meson candidates are formed with a common vertex constraint between two OS charged-particle tracks with p T > 300(150) MeV/c for PbPb (p p) sample, both required to pass standard selections [38]. The invariant mass, with a resolution of ∼3.9 (3.4) MeV/c 2 for PbPb (p p) data, is required to be within 15 MeV/c 2 of the world-average φ meson mass [28]. The B 0 s meson candidates are constructed by combining the J/ψ and φ candidates and requiring that they originate from a common vertex. Without using particle identification, assumptions need to be made about the masses of the charged particles. The difference between the natural width (according to PDG [41]) and the measured width (reflecting detector resolution) of the peaks is much bigger for the J/ψ meson than for the φ meson. Therefore, in calculating the mass of the B 0 s candidates, the two charged particles are always assumed to have the mass of charged kaons, and the muon pair is assumed to have the mass of a J/ψ meson.
The B 0 s candidates are selected according to their daughter charged particle track kinematics, the χ 2 probability of their decay vertex (the probability for the muon tracks from the J/ψ meson decay and the other charged particle tracks to originate from a common vertex), the distance between the primary and decay vertices (normalized by its uncertainty), and the pointing angle (the angle between the line segment connecting the primary and decay vertices and the momentum vector of the B 0 s meson). The selection is optimized separately for pp and PbPb results as well as each individual p T bin, using a multivariate technique that employs the boosted decision tree (BDT) algorithm [42], in order to , which is about 5σ away from the PDG B 0 s mass value. The optimal selection criterion is the working point with the highest signal significance are the expected signal (background) candidate yields from the simulated signal (data sidebands) within the mass range |M μμKK − M B 0 s ,PDG | < 0.08 GeV/c 2 .
The raw yields of B 0 s mesons in pp and PbPb collisions are extracted using an extended unbinned maximum likelihood fit to the invariant mass distribution of the B 0 s candidates in the mass range 5-6 GeV/c 2 . The estimation of the statistical uncertainties of the fitted raw yields is based on the second derivatives of the negative log-likelihood function. Examples of fits to the invariant mass distributions in pp and PbPb collisions are shown in Figs. 1 and 2 for the p T regions 7-15 and 15-50 GeV/c, respectively. The signal shape is modeled by two Gaussian functions with a common mean (which is a free parameter together with the amplitude), and different widths individually determined from MC simulations for the pp and PbPb results. The relative contribution of the two Gaussian functions to the signal yield is also fixed at the value given by the MC sample. The background is dominated by random combinations of prompt and nonprompt J/ψ candidates with extra particles and it is modeled by a first-order polynomial, as determined by studies of the inclusive J/ψ MC sample. Peaking structures that could arise from the background contamination of other B meson decays (e.g., B 0 → J/ψ K * 0 ) were found to be negligible as a consequence of the tight selection on the mass of the φ candidate.
The differential cross section for B 0 s production in |y| < 2.4 is computed in each p T interval according to dσ B 0 s dp T |y|<2.4 for pp data, and for PbPb data according to The N pp,PbPb is the raw signal yield extracted in each p T interval of width p T , (α, ) pp,PbPb represents the corresponding acceptance times efficiency, and B is the branching fraction of the decay chain. For the pp cross section, L represents the integrated luminosity, and for the PbPb cross section, N MB is the number of minimum bias events and T AA is the nuclear overlap function [43]. The T AA is equal to the number of NN binary collisions divided by the NN total inelastic cross section, and it can be interpreted as the NN-equivalent integrated luminosity per heavy ion collision. The T AA value for inclusive PbPb collisions at √ s NN = 5.02 TeV is (5.6 ± 0.2) mb −1 as estimated from an MC Glauber model [38,43].
Assuming that, in the kinematic region accessible by the present measurement, the B 0 s and B 0 s production cross sections are equal, the factor 1/2 accounts for the fact that the yields are measured for particles and antiparticles added together, but the cross section is given for one species only.

Systematic uncertainties
The cross section measurements are affected by several sources of systematic uncertainties arising from the signal extraction, corrections, B, L, N MB , and T AA determination. Unless mentioned otherwise, the same procedures were used to estimate the uncertainties for the p p and PbPb results. The uncertainty of the signal modeling is evaluated by considering four fit variations: (i) increasing/decreasing the width parameters determined from simulation by 4% (the maximum relative statistical uncertainty of the fitted width parameter among all p T bins from p p and PbPb data); (ii) using a single Gaussian function; (iii) using a sum of three Gaussian functions with a common mean, and, (iv) fixing the mean of the Gaussian function to the value determined from simulation. The uncertainty in the modeling of the background shapes is also evaluated by varying the probability distribution functions used to describe the background to a higher-order polynomial and exponential function. The maximum of the signal variations and the maximum of all the background variations are propagated as systematic uncertainties. For the p p results, the systematic uncertainty due to the selection of the B 0 s meson candidates is estimated by comparing the BDT-obtained nominal result with the results using a cut-based method (a rectangular cut) that uses the Genetic Algorithm to determine the best cut value for each parameter [42]. The same signal and background shape parametrization are used, and the same analysis parameters are optimized as in the BDT nominal method. The significance is similar for the two methods (∼8) for the p p bins. This provides an estimate of the potential difference between different selection criteria. The full difference between the two methods is propagated as a systematic uncertainty. For the PbPb results, because of the small signal in data, in order to minimize the impact of statistical fluctuations, a different approach was taken. In this case, the B 0 s selection uncertainty was estimated using the pp data sample, as the full difference in the yield between the p p results with the BDT trained on the p p sample (the nominal result) and the results with the BDT trained on the PbPb sample (the selection used for the PbPb results).
The bin-by-bin systematic uncertainties associated with the acceptance correction are estimated by varying the shape of the generated B 0 s meson p T and y spectra. For the purpose of the systematic studies only, both data and MC are split into four p T and y bins. The ratio between data and simulated p T spectra (including their statistical uncertainties) is used to generate pseudoexperiments ('toys'). Each toy is fit with a polynomial, which is A new acceptance value is calculated for each modified shape, for each kinematic bin. The root mean square (RMS) of all acceptances determined via toys is propagated as the systematic uncertainty by choosing the maximum RMS value emerging from the p T and y shape variations. Because of the small signal available, for the PbPb results the p p ratio is used to generate the toys. There is also an uncertainty assigned to account for potential bias in the efficiency calculations from the FONLL simulations of the B 0 s meson p T shape. This uncertainty is calculated as the difference between the nominal results and those obtained by generating the pythia p T shape. An additional uncertainty comes from the finite size of the MC samples. This is determined by the statistical uncertainty of the simulated signal, after applying all selection criteria.
The uncertainty in the efficiency of the muon trigger, reconstruction, and identification is evaluated bin-by-bin using control samples in data [44]. A relative systematic uncertainty of 4% per hadron track in pp collisions [29] and 6% in PbPb collisions [38] is also considered, to account for the uncertainty in the track reconstruction efficiency. This uncertainty propagates to 8% and 12% for the B 0 s measurement in pp and PbPb, respectively. The systematic uncertainty in the cross section measurement is computed as the sum in quadrature of the different contributions mentioned above.  Table 1.

Results
In Fig. 3   spectrum has a smaller uncertainty than that of the FONLL calculation. The nuclear modification factor R AA , shown in Fig. 4, is computed as: The B 0 s meson R AA is 1.5 ± 0.6(stat) ± 0.5(syst) for 7-15 GeV/c, and 0.87 ± 0.30(stat) ± 0.17(syst) for 15-50 GeV/c, respectively. In the bottom panel of Fig. 4, the R AA of B + mesons from a previous measurement [27] is also shown. Compared to the B + mesons, there is an indication of an enhancement for B 0 s mesons, which would be the expectation in the presence of a contribution from beauty recombination with strange quarks in heavy ion collisions. However, the B 0 s R AA values are compatible with unity and their large uncertainties do not exclude a significant suppression. The p T dependence of R AA is compared to the B 0 s prediction of a perturbative QCD based model that includes both collisional and radiative energy loss, (CUJET3.0) [45][46][47], and a transport model based on a Langevin equation that includes collisional energy loss and heavy quark diffusion in the medium, (TAMU) [17,48]. The difference between the two models below p T ∼ 15 GeV reflects the contribution from recombination processes, which are included in the TAMU but not in the CUJET3.0 model. The results measured for p T > 7 GeV/c have the power to disentangle the two models, albeit after an increase in precision, which can be achieved with a bigger data sample.
To further quantify the significance of a possible enhancement of the B 0 s /B + ratio in PbPb with respect to pp collisions, the ratio between the B 0 s and the B + R AA is also calculated, canceling the systematic uncertainty sources that are common to both measurements (acceptance, tracking efficiency, and muon-related). The B + R AA with a wider p T binning (15-50 GeV/c) is obtained by a B + yield weighted average of the results from three p T bins (15-20, 20-30 and 30-50 GeV/c) presented in previous work [27]. The result is shown in Fig. 5. The ratio is 4.0 ± 1.8(stat) ± 1.3(syst) for 7-15 GeV/c, and 1.8 ± 0.7(stat) ± 0.3(syst) for 15-50 GeV/c, respectively. Assuming a Gaussian distribution with mean and width equal to that of the R AA ratio and its uncertainty (including statistical and systematic components added in quadrature), the hypoth- The global systematic uncertainty, represented by the grey box at R AA = 1, comprises the uncertainties in the integrated luminosity measurement and T AA value. Two B 0 s theoretical calculations are also shown for comparison: TAMU [17,48] and CUJET3.0 [45][46][47]. The line width of the theoretical calculation from Refs. [17,48] represents the size of its statistical uncertainty. esis of the ratio values being consistent with unity (no enhancement) is tested with a χ 2 test. The resulting p-values are 18% and 28% for 7-15 and 15-50 GeV/c, respectively. This shows that, with a p-value cutoff of 5%, the scenario of no enhancement cannot be rejected. This analysis demonstrates the capability of performing a fully reconstructed B 0 s measurement in PbPb collisions with the CMS detector.

Summary
The first measurement of the differential production cross section of  7 to 50 GeV/c, inclusively for 0-100% event centrality. A hint of an enhancement of the B 0 s /B + ratio in PbPb with respect to pp collisions is seen. More precise measurements of the B 0 s and B ± mesons R AA with the upcoming high-luminosity LHC heavy ion runs could provide further constraints on the relevance of recombination, a marker of deconfined matter, for beauty hadron production, and unambiguous information about the mechanisms of beauty hadronization in heavy ion collisions.

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