Observation of B$^0_s$ mesons and measurement of the B$^0_s$ / B$^+$ yield ratio in PbPb collisions at $\sqrt{s_\mathrm{NN}} =$ 5.02 TeV

The B$^0_s$ and B$^+$ production yields are measured in PbPb collisions at a center-of-mass energy per nucleon pair of 5.02 TeV. The data sample, collected with the CMS detector at the LHC, corresponds to an integrated luminosity of 1.7 nb$^{-1}$. The mesons are reconstructed in the exclusive decay channels B$^0_s$ $\to$ J/$\psi(\mu^+\mu^-)\phi($K$^+$K$^-)$ and B$^+$ $\to$ J/$\psi(\mu^+\mu^-)$K$^+$, in the transverse momentum range 7-50 GeV/c and absolute rapidity 0-2.4. The B$^0_s$ meson is observed with a statistical significance in excess of five standard deviations for the first time in nucleus-nucleus collisions. The measurements are performed as functions of the transverse momentum of the B mesons and of the PbPb collision centrality. The ratio of production yields of B$^0_s$ and B$^+$ is measured and compared to theoretical models that include quark recombination effects.


Introduction
Relativistic heavy ion collisions allow the study of quantum chromodynamics (QCD) at high energy densities and temperatures. Under such conditions, a state in which quarks and gluons are the relevant degrees of freedom, the quark-gluon plasma (QGP) [1,2], is formed [3] as predicted by lattice QCD calculations [4]. Multiple probes are necessary for characterizing the properties of the QGP medium. Among these, heavy quarks, which are abundantly produced at the CERN LHC, have the potential of providing novel insights into QCD calculations, serving as probes of the QGP [5,6]. As they traverse the QGP, these hard-scattered partons lose energy by means of elastic collisions and medium-induced gluon radiation [7][8][9][10]. The study of parton energy loss can provide insights into the energy density and diffusion properties of the QGP. The full reconstruction of beauty and charm gives access to their four-momenta and allows the study of the flavor and mass dependences of such processes. Although heavy flavor data became available at BNL RHIC and at the LHC [5], there are still large theoretical uncertainties on the hadronization of heavy quarks in the presence of the QGP [11]. To extract the transport properties of the QGP using heavy quarks, an understanding of the hadronization mechanism is necessary.
An enhancement in the strangeness content of the medium formed in heavy ion collisions [12][13][14][15][16][17][18] is expected when the medium temperature lies above the strange quark mass [19]. Through a quark-recombination mechanism [20-23], a corresponding yield enhancement of hadrons containing strangeness is expected relative to corresponding hadrons that do not contain strange quarks. Hints of the presence of such recombination processes for heavy mesons containing strangeness were recently reported for open-charm [24] and open-beauty [25] mesons in leadlead (PbPb) collisions at a nucleon-nucleon center-of-mass energy of √ s NN = 5.02 TeV at the LHC.
The production of B mesons employing fully reconstructed decays was studied at the LHC in proton-proton (pp) collisions at center-of-mass energies, √ s, of 7, 8 and 13 TeV [26-40] over wide transverse momentum (p T ) and rapidity (y) intervals, and in proton-lead ( In this Letter, we extend these results by employing the PbPb data set collected at the end of 2018, a sample larger by about a factor of three as compared to the first measurement with 2015 data, resulting in increased statistical precision and significance. The B + (ub) and B 0 s (sb) mesons are reconstructed via the exclusive decay channels B 0 s → J/ψ φ(1020) and B + → J/ψ K + , with J/ψ → µ + µ − and φ(1020) → K + K − . Throughout this Letter, unless otherwise specified, the y and p T variables given are those of the B mesons, and charge-conjugate states are implied. The production yields of the B 0 s and B + mesons in PbPb collisions at √ s NN = 5.02 TeV and their ratios are reported, using a data set corresponding to an integrated luminosity of 1.7 nb −1 collected by the CMS experiment. The measurements are performed as functions of the meson p T and of the PbPb collision centrality, where the latter characterizes the degree of overlap of the two lead nuclei and is related to the temperature of the created medium.

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 (45-150) µm in the transverse (longitudinal) impact parameter [44]. 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 and better than 6% in the endcaps. 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 to select PbPb collision events and to define their centrality class. Centrality, defined as the fraction of the total inelastic hadronic cross section with 0% representing collisions with the largest overlap of the two nuclei, is determined experimentally using the total energy in both HF calorimeters [45]. Events are filtered using a two-tiered trigger system [46]. The first level (L1), composed of custom hardware processors, uses information from the calorimeters and muon detectors to select events within a fixed time interval of less than 4 µs. The second level, known as the high-level trigger (HLT), consists of a farm of processors running a version of the full event reconstruction software optimized for fast processing. A detailed description of the CMS experiment and coordinate system can be found in Ref. [47].
Several samples of Monte Carlo (MC) simulated events are used to evaluate background components, signal efficiencies and detector acceptance corrections. The simulations include signal samples containing only the B meson decay channels being measured, and background samples of other b hadron decays also involving J/ψ mesons in the final state. Simulated protonproton collisions are generated with PYTHIA 8 v212 [48] tune CUETP8M1 [49] and propagated through the CMS detector model using the GEANT4 package [50]. The decay of the B mesons is modeled with EVTGEN 1.3.0 [51], and final-state photon radiation is simulated with PHOTOS 2.0 [52]. Each PYTHIA8 event is embedded into a simulated PbPb collision event generated with HYDJET 1.8 [53], which is tuned to reproduce global event properties, such as the charged hadron p T spectrum and particle multiplicity [54].
The PbPb collision events were collected with an L1 hardware trigger requiring the presence of two muon candidates, with no explicit momentum threshold, in coincidence with lead bunches crossing at the interaction point. At the HLT one of the muon candidates is reconstructed using information both from the muon detectors and the inner tracker, while only information from the muon detectors is required for the other muon candidate. The two muon candidates were required to have an invariant mass within 1 < m µµ < 5 GeV/c 2 . For the offline analysis, events have to pass a set of selection criteria designed to reject events from background processes (beam-gas collisions), as described in Ref. [54]. 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. The shapes of the clusters in the pixel detector must be compatible with those expected from particles produced by a PbPb collision [55]. In order to select inelastic hadronic collisions, the PbPb events are also required to have at least two towers in each of the HF detectors with energy deposits of more than 4 GeV per tower. Only events with centrality 0-90% are selected. Collision centrality bins are given in percentage ranges, with the 0-90% bin corresponding to the 90% of the collisions having the largest overlap of the two nuclei. The PbPb sample corresponds to an integrated luminosity of approximately 1.7 nb −1 . This value is indicative only, as the PbPb yield is normalized by the total number of minimum bias events sampled, N MB = 11.8 × 10 9 .

Candidate selection and signal extraction
Muon candidates are selected within kinematic constraints that ensure that their reconstruction efficiency stays above 10%. These limits are p µ T > 3.5 GeV/c for |η µ | < 1.2, p µ T > 1.5 GeV/c for 2.1 < |η µ | < 2.4, and p µ T > (5.47 − 1.89|η µ |) GeV/c in the 1.2 < |η µ | < 2.1 region. The muons are also required to match the trigger-level muon candidates, and soft muon selection criteria are applied to global muons (i.e. muons reconstructed using the combined information of the tracker and muon detectors), as defined in Ref. [56]. Two muons of opposite charge, with an invariant mass within ±150 MeV/c 2 of the world-average J/ψ meson mass [57], are selected to reconstruct a J/ψ candidate, with a mass resolution of typically 15-30 MeV/c 2 , depending on the dimuon rapidity and p T . The 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 φ(1020) meson candidates are formed with a common vertex constraint between two oppositely charged particle tracks with p T > 300 MeV/c, both required to pass standard selection criteria described in Ref. [54]. The invariant mass, assuming the world-average charged kaon mass [57] for the two tracks, with a resolution of ∼3.9 MeV/c 2 , is required to be within 15 MeV/c 2 of the world-average φ(1020) meson mass [57].
The B 0 s (B + ) meson candidates are constructed by combining a J/ψ candidate with a φ(1020) (single track) candidate, and requiring that they originate from a common vertex. In the kinematic vertex fit to the dimuon plus two-track (single-track) system, the tracks are assigned the mass of the charged kaon, and the invariant mass of the muon pair is constrained to the nominal J/ψ meson mass [57]. The B candidate selection is performed via multivariate discriminators, based on a boosted decision tree (BDT) method [58]. The selection is optimized separately for each meson, as well as each individual bin of p T . The discriminating variables employed include: the χ 2 probability of the decay vertex (the probability for the muon tracks from the J/ψ meson decay and the other charged particle track(s) to originate from a common vertex), the decay length (normalized by its uncertainty), the pointing angle (the angle between the line segment connecting the primary and decay vertices and the momentum vector of the B meson), the cosine value of the angle between the B meson displacement and momentum in the transverse plane, the two-dimensional (2D) distance between the primary and decay vertices of their daughter tracks (normalized by their uncertainties), and the p T of the daughter charged particle track(s). For the B 0 s meson, two additional selection variables are used: the distance along the z-direction between the primary vertex and the decay vertex (normalized by its uncertainty), as well as the absolute difference between the reconstructed φ(1020) meson invariant mass and its nominal value.
The BDT training is performed employing simulated B signal samples, and background samples taken from data sidebands (candidates with the invariant mass 0.15-0.35 GeV/c 2 away from the B meson nominal mass [57]). The signal samples are scaled to the number of B candidates predicted by fixed-order plus next-to-leading order logarithmic (FONLL) perturbative QCD calculations [59-61] for the integrated luminosity of the analyzed data sample. The background rejection achieved with the BDT selection is larger than 97%.
The raw B meson signal yields are extracted using an extended unbinned maximum likelihood fit to the invariant mass spectra. The signal shape is modeled by a sum of two Gaussian functions, with parameters determined from MC simulation, except for the common mean and the overall signal yield, which are free parameters of the fit. The combinatorial background, from uncorrelated combinations of J/ψ candidates with extra particles, gives rise to a falling contribution in the invariant mass spectrum that is modeled by an exponential function. An additional background source can arise from possible contamination from other b hadron decays. For the B + spectrum, partially reconstructed B decays, such as B 0 → J/ψ K * 0 → µ + µ − K + π − , lead to a heightened background in the invariant mass region below 5.2 GeV/c 2 . The decay B + → J/ψ π + , where the pion is misidentified as a kaon, results in a small peaking structure in the signal region. Such partially reconstructed and misreconstructed b hadron components are modeled from simulation, via an error function and a Gaussian function, respectively; both shapes and the normalization of the Gaussian function (relative to the signal) are fixed in the fit to the data. For the B 0 s meson case, such background contributions are found to be negligible, as a consequence of the tight selection on the mass of the φ(1020) candidate. The results of the fits to the invariant mass distributions are shown in Fig. 1.
The statistical significance of the B meson signals is estimated from the ratio of likelihoods obtained by fitting the data with the full model and the background-only model. Using the standard asymptotic approximation for the likelihood [62], the estimate obtained for the B 0 s meson corresponds to approximately 16 standard deviations. The estimate is obtained with the nominal signal and background model described above, and verified with the model variations presented in Section 5. This is the first observation of the B 0 s signal in heavy ion collisions.

Production yields
For each B meson, the production yield scaled by the overlap nuclear function, T AA [63], and the number of minimum bias events, N MB , is computed in each p T interval, according to where N obs denotes the raw signal yield, extracted in each p T interval of width ∆p T . The factor 1/2 accounts for the fact that the raw yield is measured for particles and antiparticles added together while the production yield is given for one species only, and B stands for the branching fraction of the corresponding decay chain. The factor T AA is equal to the number of nucleonnucleon (NN) binary collisions divided by the NN total inelastic cross section. The value of T AA can be interpreted as the NN-equivalent integrated luminosity per heavy ion collision. The T AA values for √ s NN = 5.02 TeV PbPb collisions for the centrality intervals 0-30%, 30-90%, and 0-90% are, respectively, (15.
The measurement is performed in a similar fashion as a function of the PbPb collision centrality. The centrality class is represented as the average number of participating nucleons, N part . This number, determined using the above Glauber model, is highly correlated with the impact parameter of the collision.
The acceptance and the efficiency factors, the last term in Eq. (1), are obtained employing 2D fine-grained maps of the acceptance α(p T , y) and the efficiency (p T , y). These maps are determined, for each centrality region, from MC-simulated samples of B meson signal events, generated within the analysis fiducial region given by: |y| < 1.5 and p T > 10 GeV/c; 1.5 < |y| < 2.4 and 7 < p T < 50 GeV/c. The acceptance corresponds to the fraction of generated events passing the muon and kaon selection thresholds specified in Section 3, while the efficiency is determined as the fraction of accepted events that pass the full analysis selection criteria, including the trigger. The maps are used to determine the 1/(α ) value for each B candidate in data, based on the kinematics (p T , y) of the candidates, and the corresponding average 1/(α ) obtained for the candidates within a 160 MeV/c 2 window centered on the B meson's worldaverage mass. The efficiency maps derived from simulation are corrected by data/MC scale factors for the muon reconstruction and trigger efficiencies, obtained by applying the tag-andprobe method using the J/ψ resonance [65].

Systematic uncertainties
The T AA -scaled yield measurements are affected by several sources of systematic uncertainties arising from the signal extraction, acceptance and efficiency, B, N MB , and T AA determinations. The uncertainty from the signal and background modeling is evaluated by considering the following fit variations: (i) using low-order polynomials for describing the combinatorial background, (ii) using a sum of three Gaussian functions with a common mean for alternatively describing the signal, (iii) varying the widths of the signal double-Gaussian function by 10% (to account for possible mismatches in the mass resolution between simulation and data), and (iv) fixing the double-Gaussian mean to its MC simulation value. The maximum of the signal variations and the maximum of the background variations are added in quadrature as the systematic uncertainty from fit modeling.
The systematic uncertainties associated with the limited size of the MC simulation samples are obtained by recomputing the 2D fine-grained α(p T , y) (p T , y) maps and the corresponding 1/(α ) averages over the data sample. The statistical uncertainties from the MC-derived product of acceptance and efficiency are propagated by generating 10000 2D maps (per bin in the differential measurement), where each entry is obtained by randomly sampling the α values from Gaussian distributions with widths given by the statistical uncertainty. The resulting systematic uncertainties are obtained as the width of the 1/(α ) distribution from the alternative 2D maps.
The efficiency and acceptance determinations, based on MC simulation, are further affected by potential disagreements between data and MC. The signal MC simulation is validated against data by inspecting the distributions of the variables employed in the selection. The signal distributions are extracted from the data employing the sPlot method [66], using the mass of the B meson candidates as discriminating variable, and are also cross-checked with a simple sideband-subtraction method. In particular, the sPlot-derived signal distributions for the BDT score are retrieved from the data, and the corresponding data/MC ratios formed. These ratios are in turn used to re-weight the MC simulation, and the resulting deviation in the 1/(α ) factors are assigned as systematic uncertainties. This procedure is employed using the higheryield channel, namely the B + . For the B 0 s meson, the limited size of the data sample yields results that are compatible between data and simulation, within the statistical uncertainties. As such, the B + channel is used as calibration mode for the B 0 s channel, and a procedure identical to that described above is applied to the B + sample but using only tracking related variables. The uncertainty in the efficiency of the muon trigger, reconstruction, and identification is evaluated using a tag-and-probe technique [65]. The derived data versus MC scale factors are employed in the determination of the nominal α 2D maps, while the associated asymmetric uncertainties are propagated as the systematic uncertainty in the 1/(α ) factors. The difference in the track reconstruction efficiency in data and simulation is estimated by comparing 3-prong and 5-prong D * decays, D * → Kππ(ππ) [54]. This results in 5 (10)% uncertainty in the efficiency determination for the B + (B 0 s ) decay for involving one (two) kaon(s). The determination of the 1/(α ) factors, for each analysis bin, involves averaging over the data. In order to reflect this statistical effect, the overall statistical uncertainty from the PbPb data is calculated using a bootstrap method [67], where data events are resampled 1000 times with replacement. The corrected yield is re-determined for each bootstrap sample, and the spread of its distribution is assigned as the statistical uncertainty.
Various consistency checks of the analysis procedure have been performed, which do not result in additional systematic uncertainties. The consistency of the fitting procedure was verified by generating and fitting pseudo-experiments using the likelihood model and parameters employed for fitting the data. The distributions of the pull (here defined as the difference between the generated and fitted parameter values divided by the uncertainty returned by the fit) show random behavior consistent with Gaussians of zero mean and unit width, which indicate the absence of bias in the fitting procedure. A consistency check was performed using MC simulation, where the signal yields are extracted with the fitting procedure and the acceptance and efficiency corrections applied, thus retrieving the generated yields. Potential biases in the acceptance and efficiency calculations associated with the shape of the B mesons kinematic distributions from simulation were probed by re-weighting the PYTHIA simulation according to varied p T shapes and recomputing the factors 1/(α ) . Such deviations were found to be negligible. This confirmed independence from the MC simulation is expected in view of the fine-grained 2D α MC maps and the 1/(α ) averaging based on data, and further attests the robustness of the analysis procedure. The effect of background contamination in the averaging procedure was studied employing the sPlot method and was found to be negligible.
The systematic uncertainties in the Glauber model normalization factor (T AA ) are derived from propagating the uncertainties in the event selection efficiency, and in the nuclear radius, skin depth, and minimum distance between nucleons in the Pb nucleus parameters of the Glauber model [54]. The uncertainties associated with the branching fraction of the decay chains B are obtained from Ref. [57]. The B factor is common to all bins in the production yield measurements, as is the case also of T AA and N MB in the p T -differential results; the corresponding uncertainties are denoted global uncertainties.
The total systematic uncertainties in the T AA -scaled yield measurement are computed as the sum in quadrature of the different contributions mentioned above, which are summarized in Tables 1 and 2.

Results
The differential T AA -scaled yields for B 0 s and B + mesons in PbPb collisions, Eq. (1), as a function of the meson p T are presented in Fig. 2 (left). On the horizontal scale, the abscissa of each data point is set to the mean value of the p T distribution, after background subtraction via sPlot, while the bin ranges are represented by the boxes. In Fig. 2 (right), the T AA -scaled yields are presented for different PbPb event centrality ranges.
The ratio of production yields of the B 0 s and B + mesons is formed as The terms T AA and N MB , that are common to the two meson samples, cancel in the ratio R. Correlated systematic uncertainties associated with the efficiency determination are partially reduced. For the tracking efficiency, a 5% uncertainty is assigned because of the different number of hadron tracks in the two decays. Muon efficiency uncertainties from the tag-and-probe method are determined by varying coherently the efficiency factors in the numerator and denominator within the corresponding uncertainties. For the determination of the uncertainties related to data/MC comparisons, an identical procedure to that used for the production yield is employed, but utilizing only tracking related variables, using the B + channel. The dominant uncertainty arises from data/MC agreement and the limited MC sample size. The R ratio results are shown in Fig. 3 as functions of the mesons' p T (left) and event centrality (right).
The p T dependence of the ratio R is compared to a transport model (TAMU) based on a Langevin equation that includes collisional energy loss and heavy quark diffusion in the medium [23] and an advanced Langevin hydrodynamic model with both heavy quark elastic and inelastic energy loss [68]. Both predictions account for recombination contributions to the B 0 s meson yields, which are expected to be more significant at low p T . Predictions form a quark combination model with the equal-velocity combination (EVC) approximation [69] are also shown for both PbPb and pp collisions. The ratio of the B 0 s and B + yields in pp collisions has been determined by LHCb [70], at different collision energies, for B mesons within 0.5 < p T < 40 GeV/c and 2 < η < 5. While no rapidity dependence of the ratio has been found, a decrease with p T and an increase with collision energy have been reported [70,71]. While the B 0 s to B + meson yield ratios obtained in charmonium decay channels cannot be directly compared with those measured in e + e − collisions at LEP [72], earlier LHCb measurements [32, 73,74] in pp collisions in other hadronic decay channels are compatible with the LEP results. The lack of sizable strangeness enhancement in the B 0 s meson hadronization in pp collisions, evident from the above agreement, makes pp collisions a robust baseline for comparison with our PbPb results. In Fig. 3, the pp results at 7 TeV by LHCb [70] are scaled by the same ratio of branching fractions [57] as used in Eq. (2). The results are compatible with both the model predictions for PbPb collisions and the pp reference results. However, no significant p T dependence can be established with the precision allowed by the current data. The results, while compatible, tend to lie systematically above the pp reference data. While the central value of R in the more central bin lies above that of the more peripheral bin, the results are compatible within the combined statistical and systematic uncertainty. These results are compatible with the ratio of nuclear modification factors, R AA , reported earlier by the CMS Collaboration [25], albeit with reduced statistical uncertainties.

Summary
The B 0 s and B + mesons are studied with the CMS detector at the LHC via the reconstruction of the exclusive hadronic decay channels B 0 s → J/ψ φ(1020) → µ + µ − K + K − and B + → J/ψ K + → µ + µ − K + . The measurements are performed within the B mesons' fiducial region given by transverse momentum p T > 10 GeV/c for rapidity |y| < 1.5 and 7 < p T < 50 GeV/c for 1.5 < |y| < 2.4. The first observation of the B 0 s meson in nucleus-nucleus collisions, with a statistical significance well surpassing five standard deviations, is attained. The production yields of B 0 s and B + mesons, scaled by the nuclear overlap function T AA and the number of minimum bias events N MB , in lead-lead (PbPb) collisions at a center-of-mass energy of 5.02 TeV per nucleon pair are presented as functions of the meson p T and, for the first time, of the event centrality. These results extend, and are compatible with, those previously reported by the CMS Collaboration [25,43], and are based on a three-fold larger PbPb data sample. The ratio of production yields of the two mesons in PbPb collisions is determined and it is found to be statistically compatible with the corresponding ratio in proton-proton (pp) collisions. The further investigation of possible hints of an enhancement of the ratio in PbPb, relative to pp collisions, will benefit from more precise PbPb and pp reference data taken at the same collision energy per nucleon. The larger PbPb data sets that should be accumulated in upcoming high-luminosity LHC heavy ion runs will provide greater precision and could help to further characterize the mechanisms of beauty hadronization in heavy ion collisions.

Acknowledgments
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 centres and personnel of the Worldwide LHC Computing Grid and other centres for delivering so effectively the computing infrastructure essential to our analyses. Finally, we acknowledge the enduring support for the construction and operation of the LHC, the CMS detector, and the supporting computing infrastructure provided by the following funding agencies:   [33] LHCb Collaboration, "Measurements of B + c production and mass with the B + c → J/ψπ + decay", Phys. Rev. Lett. 109 (2012) 232001, doi:10.1103/PhysRevLett.109.232001, arXiv:1209.5634.