Measurement of $B^+$, $B^0$ and $\Lambda_b^0$ production in $p\mkern 1mu\mathrm{Pb}$ collisions at $\sqrt{s_\mathrm{NN}}=8.16\,{\rm TeV}$

The production of $B^+$, $B^0$ and $\Lambda_b^0$ hadrons is studied in proton-lead collisions at a centre-of-mass energy per nucleon pair of $\sqrt{s_\mathrm{NN}}=8.16\,{\rm TeV}$ recorded with the LHCb detector at the LHC. The measurement uses a dataset corresponding to an integrated luminosity of $12.2\pm0.3\,\mathrm{nb}^{-1}$ for the case where the proton beam is projected into the LHCb detector (corresponding to measuring hadron production at positive rapidity) and $18.6\pm0.5\,\mathrm{nb}^{-1}$ for the lead beam projected into the LHCb detector (corresponding to measuring hadron production at negative rapidity). Nuclear effects are probed through double-differential cross-sections, forward-to-backward cross-section ratios and nuclear modification factors of the beauty hadrons. The double-differential cross-sections are measured as a function of the beauty-hadron transverse momentum and rapidity in the nucleon-nucleon centre-of-mass frame. Forward-to-backward cross-section ratios and nuclear modification factors indicate a significant nuclear suppression at positive rapidity. The ratio of $\Lambda_b^0$ over $B^0$ production cross-sections is reported and is consistent with the corresponding measurement in $pp$~collisions.


Introduction
Charm and beauty quarks provide a unique probe of nuclear matter in heavy-ion collisions [1]. They are produced at early times of the collisions and experience the whole evolution of the nuclear medium before hadronization [2]. Their kinematics and hadronization provide information on the extent of thermalization effects and on transport coefficients. The hard scale provided by the heavy-quark mass is larger than the Quantum Chromodynamics (QCD) scale, Λ QCD . Therefore, heavy-quark production can be addressed with perturbative QCD down to zero transverse momentum (p T ).
The characterization of the extended color-deconfined thermodynamic system, the quark-gluon plasma, using heavy-quark observables in heavy nucleus-nucleus collisions, requires an understanding of background effects. Therefore, it is mandatory to identify and constrain other QCD effects that may appear in nuclear collisions. Among these effects, the modification of the parton distribution functions [3][4][5][6][7] or, alternatively, the breakdown of collinear factorization in the gluon-dense nuclear wave function [8,9] are discussed most extensively. Besides the modification of the nuclear wave function compared to that of free nucleons, coherent gluon radiation at small angles may modify final-state heavy-quark kinematic distributions [10]. Furthermore, the nuclear effect that is responsible for the change of hadronization patterns as a function of final-state particle multiplicities in small collision systems (pp and proton-nucleus collisions), first observed for strange-hadron production [11], is not yet fully understood. Measurements sensitive to hadronization fractions in the heavy-flavor sector can contribute to a better understanding. Studies of hadronization in heavy nuclear collisions may help to explain the puzzle of heavyflavor hadron collective behaviour that was observed recently in pp and proton-lead collisions [12][13][14]. These measurements in small collision systems still require a common reconciliation with the global theoretical picture of heavy-ion collisions based on fluid dynamics, or might result in modifications to the fluid-based description.
Observables related to charm hadrons have been extensively studied at the high-energy frontier of heavy-ion collisions at RHIC and the LHC [1]. Recently, the first measurements of Λ + c baryon 1 production in proton-lead collisions have been performed at the LHC [15,16]. The measurements of charm-baryon production were the last important step towards the evaluation of the total charm production cross-section without relying on assumptions about charm fragmentation functions based on measurements made before the start of the LHC. 2 Beauty hadrons are not yet explored experimentally to the same extent in heavy-ion collisions due to lower production rates. Theoretically, computations of the production of beauty hadrons are more reliable than charm hadrons since the larger beauty-quark mass allows for a better separation of energy scales with respect to Λ QCD . The LHCb collaboration has recently studied the production of J/ψ mesons from beauty-hadron decays (nonprompt J/ψ ) in proton-lead collisions [17]. This measurement is sensitive to beauty-quark production down to vanishing transverse momentum with good precision.
This article presents measurements of the production cross-sections of fully reconstructed B + , B 0 and Λ 0 b hadrons in proton-lead collisions recorded by the LHCb experiment, as a function of the hadron kinematics down to p T = 2 GeV/c, which is lower than the hadron masses. The measurement of heavy-quark production at low p T helps to constrain the gluon wave function in the nucleus in the small Bjorken x region [18][19][20][21], where x is the fraction of the nucleon momentum carried by the interacting gluon. In addition, production measurements of fully reconstructed beauty hadrons in heavy-ion collisions can test whether the hadronization fractions in nuclear collisions are the same as those measured in pp collisions [22][23][24][25]. and rapidity, y. The rapidity is defined in the nucleon-nucleon centre-of-mass frame, using the proton beam direction as the direction of the z-axis of the coordinate system. Since the energy per nucleon in the proton beam is larger than in the lead beam, the nucleon-nucleon centre-of-mass system has a rapidity in the laboratory frame of 0.465. During the data taking in 2016, the LHC provided collisions with two configurations, inverting the direction of the proton and lead beams. The LHCb forward spectrometer covers the positive (negative) rapidity ranges when the proton (lead) beam direction is projected into the detector from the interaction region, denoted as "pPb" ("Pbp") configuration.
The dataset corresponds to an integrated luminosity of 12.2 ± 0.3 nb −1 for the pPb configuration and 18.6 ± 0.5 nb −1 for the Pbp configuration, calibrated using dedicated luminosity runs [42]. The double-differential cross-section of the production of a H b hadron as a function of p T and y is computed as where, for a given interval of p T and y, N (H b ) + N (H b ) is the sum of the observed signal yields in a particular decay mode and its charge-conjugated decay mode, B(H b ) is the product of the branching fractions for the beauty decay and the subsequent charm decay, L is the integrated luminosity, and is the total detection efficiency of the final state particles. The measurements are carried out in the kinematic range 2 < p T < 20 GeV/c and 1.5 < y < 3.5 for the pPb configuration, and in the range 2 < p T < 20 GeV/c and −4.5 < y < −2.5 for the Pbp configuration. The p T intervals used to study the efficiency and signal yield are 2-4 GeV/c, 4-7 GeV/c, 7-12 GeV/c and 12-20 GeV/c, and the rapidity regions are split into two equal size intervals, −4.5 < y < −3.5 and −3.5 < y < −2.5 for the pPb configuration, and 1.5 < y < 2.5 and 2.5 < y < 3.5 for the Pbp configuration. The range p T < 2 GeV/c is not considered due to the small signal yield with the current sample. This restriction is not related to any detector limitation specific to the collision system, but to the limited integrated luminosity and the small production cross-section. Nuclear effects are quantified by the nuclear modification factor, R pPb , where A Pb = 208 is the mass number of the lead ion, d 2 σ pPb (p T , y)/dp T dy the H b production cross-section in proton-lead collisions as defined in Eq. (1), and d 2 σ pp (p T , y)/dp T dy the H b reference production cross-section in pp collisions at the same nucleon-nucleon centre-of-mass energy. In the absence of nuclear effects, the nuclear modification factor is equal to unity. To quantify the relative forward-to-backward production rates, the forward-to-backward ratio, R FB , is measured, which is the ratio of cross-sections in the positive and negative y intervals corresponding to the same absolute value range, 3 Selections, signal yields and efficiency

Candidate reconstruction and selection
The B + cross-section is measured in the B + → J/ψ K + mode, with J/ψ → µ + µ − , and in the purely hadronic mode B + → D 0 π + , with D 0 → K + π − . The cross-sections of the B 0 and Λ 0 b hadrons are studied in the hadronic decays modes, the candidates are reconstructed from a sample selected by a hardware trigger requiring a minimum activity in the scintillating-pad detector. This hardware trigger selection has an efficiency of 100% for the signal. The intermediate charm-hadron candidates are reconstructed using tracks that are identified as pion, kaon and proton candidates by the LHCb particle identification system [27]. The tracks used to form the D 0 (D − and Λ + c ) candidates are required to have p T > 300 MeV/c, and at least one of them has to satisfy p T > 500 MeV/c (p T > 1000 MeV/c). They must also have momentum p > 3 GeV/c (p > 10 GeV/c for protons) and pseudorapidity in the range 2 < η < 5. In addition, they are required to be separated from any primary vertex by requiring χ 2 IP > 16, where χ 2 IP is the difference between the χ 2 values of a given PV reconstructed with and without the considered track. The tracks are required to form a vertex of good quality. Further requirements are imposed to ensure that this vertex is consistent with charm-hadron decays by requiring a minimum reconstructed decay time and a reconstructed mass within an interval centred on the known values of the hadron mass [43]: [1834. 8, 1894.8] MeV/c 2 , [1844.6, 1894.6] MeV/c 2 and [2268.5, 2304.5] MeV/c 2 for D 0 , D − and Λ + c candidates, respectively. Each mass interval corresponds to six times the experimental resolution on the reconstructed mass. A charm-hadron candidate, inconsistent with originating from the PVs as ensured by the requirement χ 2 IP > 4, is then combined with a positively identified pion of the appropriate charge to form a beauty hadron. This pion is required to have p T > 500 MeV/c and to be separated from any PV with the condition χ 2 IP > 16. Reconstructed beauty hadrons with a good-quality vertex and a significant displacement from any PV are selected and are further required to point back to a PV by imposing χ 2 IP < 16. The offline selected beauty-hadron candidates are also required to match an online vertex, reconstructed from two, three or four tracks, with a large sum of the transverse momenta of the tracks and a significant displacement from the PVs.
The B + candidates studied with the B + → J/ψ K + decay are obtained from a data sample that contains J/ψ candidates reconstructed by the online software trigger [34]. The muons used to reconstruct a J/ψ meson are identified by the muon detector and information from all subsystems combined by a neural network. The J/ψ candidate must have a well-reconstructed vertex, a mass in the range [3056.9, 3136.9] MeV/c 2 , and pass the hardware trigger that selects muons with p T > 500 MeV/c. The J/ψ candidate with a reconstructed decay vertex significantly separated from all PVs is combined with a kaon track to form a B + candidate. The K + candidate must be positively identified and is required to have a transverse momentum p T > 500 MeV/c and to be separated from all PVs with the requirement χ 2 IP > 16. The reconstructed B + candidate is required to have a good-quality vertex, be displaced from the PVs and point back to a PV by requiring χ 2 IP < 16.

Signal yield determination
The signal yields for each decay mode are determined from extended unbinned maximumlikelihood fits to their mass distributions. The fits are used to calculate per-candidate weights with the sPlot method [44]. The weights are then used to determine the signal yields in each p T and y bin. As a cross-check, fits are also performed in individual p T and y bins, and the results are consistent with those obtained using the sPlot method. The signal mass distribution is described by a Crystal Ball (CB) function [45] for the B + → J/ψ K + , B 0 → D − π + and Λ 0 b → Λ + c π − decays. For the B + → D 0 π + decay an additional Gaussian function, which shares the peak position with the CB function, is necessary to achieve a satisfactory fit quality. The tail parameters for the CB function and the fractions of the CB and the Gaussian components are fixed to values obtained from fits to mass spectra of simulated signal decays. The mean and width of the Gaussian core in the CB function, and the width of the separate Gaussian component are free parameters determined from data. The combinatorial background is described by an exponential function, with parameters allowed to vary in the fits.
The contribution of misidentified background from B where the K ± (π + ) meson is reconstructed as a π ± (K + ) candidate, is described by an empirical function obtained using simulation. Due to the small branching fraction of the misidentified background compared to the signal and the suppression from the PID requirement, the contribution relative to the signal mode in the selected sample is expected to be around or below 5% depending on the decay mode.
The shape for each component, except that of the combinatorial background, is constrained to be the same for the fits to pPb and Pbp data. The yields for each contribution in the fit model are free parameters determined from data with the constraint that the ratio of misidentified background to signal yield is the same in pPb and Pbp data. The signal yields for each decay model considered in this analysis are summarized in Table 1 for the kinematic range 2 < p T < 20 GeV/c and 1.5 < y < 3.5 (−4.5 < y < −2.5) in the pPb (Pbp) sample. The mass distributions and the fit projections are shown in Figs. 1 to 4 for the decays The higher combinatorial background level in the Pbp sample compared to the pPb sample is due to higher charged track multiplicities seen by the LHCb detector in the Pbp beam configuration.

Efficiency
The total efficiency is the product of the geometrical acceptance of the detector and the efficiencies of the reconstruction, the selection, the PID and the trigger requirements. Table 1: Signal yields in the range 2 < p T < 20 GeV/c and 1.5 < y < 3.5 (−4.5 < y < −2.5) for pPb (Pbp) collisions. Uncertainties are statistical only.

Decay
pPb Pbp It is about a few percent in the low-p T region, and 20% in the high-p T region. These efficiencies, except for the PID, are evaluated using samples of simulated signal decays, in bins of the beauty-hadron p T and y. The reconstruction efficiency obtained from simulated signals is corrected using a data-driven method which is detailed in the next paragraph. The occupancy distribution in the minimum bias simulation sample is weighted to reproduce that in data, in order to model correctly the PV reconstruction efficiency. For the decays B + → D 0 π + , B + → J/ψ K + and B 0 → D − π + and subsequent charm-hadron decays, the angular distributions of the final state particles are well described by EvtGen.
For the Λ 0 b → Λ + c π − decay, the Dalitz-plot distribution of the Λ + c → pK − π + decay in simulation is described by a mixture of uniform phase space and resonant contributions of ∆(1232) ++ → pπ + and K * (892) 0 → K − π + . The Λ + c Dalitz-plot distribution in the simulation is corrected to match that in the background subtracted data.
The track reconstruction efficiency from simulation is corrected using a tag-andprobe approach. For this method, J/ψ candidates in data are formed combining a fully reconstructed "tag" track with a "probe" track reconstructed using a subset of the tracking detectors [27,48]. The single-track reconstruction efficiency is obtained as the fraction of the probe tracks that are matched to fully reconstructed tracks, in bins of the track momentum and pseudorapidity. The ratio of the tag-and-probe efficiency between protonlead data and simulation is used to correct the simulation efficiencies. The correction factors are determined for the pPb and Pbp samples separately.
The PID efficiency for each track is determined with a tag-and-probe method [49,50] using calibration samples of proton-lead data. The track PID efficiency depends on the detector occupancy. Since the occupancy distribution is found to be consistent between the calibration samples and the beauty-signal events, the efficiency is parametrized as a function of track momentum and pseudorapidity. The pion and kaon PID efficiencies are calibrated using D 0 → K − π + decays, where the D 0 flavor is tagged by the charge of the pion in D * + → D 0 π + decays, the proton PID efficiency is studied using Λ → pπ − decays and the PID efficiency for muons is obtained using J/ψ → µ + µ − decays. For each beauty candidate, the product of the single-track PID efficiencies, measured as a function of the track momenta and pseudorapidity, gives the combined PID efficiency for all the tracks in the final state. The efficiency is then averaged over all beauty-hadron candidates for each bin of p T and y.

Systematic uncertainties
The various sources of systematic uncertainties, and their quadratic sum, on the crosssections for B + , B 0 and Λ 0 b hadrons are summarized in Tables 2 and 3 for the pPb and Pbp data samples, respectively. The ranges in the tables correspond to the minimum and maximum values over the p T and y bins of the measurement. The cross-section of the B + hadron is measured in the two decay modes, B + → J/ψ K + and B + → D 0 π + , which give consistent results within statistical uncertainties.
The uncertainty on the b-hadron signal yields is studied by using alternative fit models or different fitting ranges for the mass distributions. The nominal CB function for the signal mass distribution is replaced by a combination of a Gaussian function plus a CB function or vice-versa for the B + → D 0 π + decay, giving a relative change of 2% on the signal yields for all the decay modes. A second-order polynomial is employed to replace the exponential function for the combinatorial background, which results in a difference of 1% for the signal yields at the maximum. The effect of partially reconstructed background is studied by fitting the mass distribution in a smaller region where its contribution is reduced or absent. The signal yields change by at most 1% for all the decay channels. The effect of the misidentified background is studied by fixing its branching fraction relative to that of the signal [43], corrected by the PID selection efficiency. The change in signal yields amounts to 0.1%. The maximum value among all these effects, 2%, is quoted as the systematic uncertainty, and is considered as a global uncertainty for all decay modes    and all p T and y bins. The corrections to the track reconstruction efficiency are limited in precision by the size of the calibration data sample, which results in a systematic uncertainty dominating in most of the analysis bins. This effect is studied by generating sets of correction factors according to Gaussian distributions centered on their nominal values and with widths equal to the statistical uncertainties. The standard deviation of the variations of the corrected efficiency in simulation is assigned as uncertainty, labelled as "Tracking sample size" in the summary tables. It ranges from 2.0% to 9.5% for pPb and from 4.6% to 17.8% for Pbp, depending on the decay modes and the beauty-hadron p T and y bins. The larger uncertainty for the Pbp sample, where the LHCb detector accepts particles produced in the lead beam direction, is due to higher background that makes the signal yield determination in the calibration data sample more difficult. The tag-and-probe method used to calculate the tracking efficiency has an uncertainty estimated to be 0.8% per track [48], giving a total value of 2.4% (3.2%) for a three-(four-)track decay mode. Since the tracking efficiency is measured using muons, an additional uncertainty of 1.5% per track is introduced for hadrons, to account for the possible imperfect modeling of the amount of interactions with the detector material. Labelled as "Hadron tracking" in the summary tables, the result is equal to 1.5% for B + → J/ψ K + and to 4.5% (6%) for three-(four-)track hadronic decays. The uncertainties related to the track reconstruction efficiency method and to the hadron-detector interactions are fully correlated among different hadron species and between the pPb, Pbp and pp datasets.
Several sources of systematic uncertainties are associated with the PID efficiencies. The contribution due to the limited size of the data calibration samples is determined by varying the single-track PID efficiencies within their uncertainties for all momentum and pseudorapidity bins simultaneously, and calculating the resulting spread of the PID efficiencies on the b-hadron signal decays. Since large samples are available for the kaon, pion, and proton calibration, the resulting systematic uncertainties are found to be small and in the range of 0.2%-0.7% (0.1%-0.5%) for the B + → D 0 π + decay, B 0 → D − π + and Λ 0 b → Λ + c π − decays in pPb (Pbp) collisions. They are labelled as "PID sample size" in the summary tables. For B + → J/ψ K + decays, the smaller size of the muon calibration samples results in a systematic uncertainty between 1.4% and 2.7% for the pPb data and between 0.7% and 2.1% for the Pbp data. For each bin of track momentum and pseudorapidity, the possible difference in track kinematics between the PID sample and the b-hadron sample is counted as a second source of systematic uncertainty. The effect is studied by varying the default binning scheme using finer bins, and determining the changes of the PID efficiencies on the b-hadron signal decays. The result is labelled as "PID binning" in the summary tables and is found to be at most 1.4%. The systematic uncertainty related to a possible difference of detector occupancy between the PID samples and the b-hadron samples is studied by weighting the occupancy in the PID samples to match that of the signal beauty sample, and the resulting change of the efficiency is found to be negligible.
The imperfect modeling of b-hadron kinematic distributions and decay properties in the simulation introduces systematic uncertainties on the reconstruction and selection efficiencies. The two-body invariant mass distributions of the Λ + c decay products, or Dalitz-plot distribution, for the Λ 0 b → Λ + c π − mode in simulation is weighted to match data, and the uncertainty on the Dalitz-plot distribution is counted as a source of systematic uncertainty. Its magnitude is studied by pseudoexperiments. For each pseudoexperiment, a sample is constructed by randomly sampling Λ 0 b candidates from data allowing for repetition, and this sample is used to correct the Dalitz-plot distribution in the simulation. The root-mean-square value of the efficiencies corrected with multiple pseudoexperiments is quoted as the systematic uncertainty. It is found to be in the range 0.8%-3.1% for the different Λ 0 b p T and y bins and is labelled as "Dalitz structure" in the summary tables. The distributions of variables used to select candidates show good agreement between data and simulation. The effect of the residual differences is quantified by weighting the reconstructed b-hadron decay-time distribution in simulation to match that in data, and studying the corresponding variation of the selection efficiency. The result, labelled as "Selection" in the summary tables, amounts to 1% for the two B + decay modes, and to 3% and 2% for the B 0 and Λ 0 b decay modes. Simulation and data also show reasonable agreement in the beauty-hadron p T and y distribution, even if a modest discrepancy in the p T distribution is observed, especially for the Λ 0 b baryon. Due to the limited data sample size it is not possible to accurately determine the b-hadron p T and y distributions from data directly. However, as the crosssection is measured differentially in bins of p T and y, the small discrepancy on these kinematic distributions has a reduced impact. A systematic uncertainty is evaluated as the change in the reconstruction efficiency after reweighting the p T and y distributions in simulation to match data using a finer binning scheme. The result, labelled as "Kinematics" in the summary tables, ranges from a fraction of a percent to a few percent depending on the decay modes and the beauty-hadron p T and y bins.
The muon trigger efficiency is validated using a large sample of J/ψ → µ + µ − decays obtained with an unbiased trigger selection [32]. The result is compared with the trigger efficiency estimated in simulation, showing a difference of at most 1%, which is quoted as the systematic uncertainty due to the trigger selection for the B + → J/ψ K + decay. Thanks to the loose requirement applied by the online event selection, the overall trigger efficiency for the purely hadronic decay modes is found to be above 99% for the offline selected candidates. A systematic uncertainty of 1% is assigned.
The finite sizes of the simulated b-hadron signal samples introduce uncertainties on the efficiency, which are propagated to the cross-section. Labelled as "Simulation sample size", these uncertainties range from subpercent to a few percent depending on the decay modes and the p T and y bins. The uncertainties due to the integrated luminosity of the pPb and Pbp datasets are of 2.6% and 2.5%, respectively. The uncertainties on the branching fractions of the b-hadron decays and of the intermediate charm-hadron decays are also sources of systematic uncertainty, and are evaluated using the uncertainties on the measured values [43].
The dominant systematic effect is the uncertainty on the track reconstruction efficiency which, however, largely cancels in the cross-section ratios. For the Λ 0 b → Λ + c π − decay, the branching fraction is also a large source of systematic uncertainty, but cancels for the nuclear modification factor measurements. The systematic uncertainties are considered to be fully correlated among all kinematic bins for a particular decay mode, except that labelled as "Simulation sample size" which is uncorrelated. Table 4: Differential cross-sections of B + , B 0 and Λ 0 b production in bins of p T and y, d 2 σ dp T dy ( µb/[GeV/c]), and in bins of y integrated over 2 < p T < 20 GeV/c, dσ dy ( µb). The first uncertainty is statistical and the second systematic.

Cross-sections
The B + cross-sections measured in the J/ψ K + and D 0 π + decay modes are consistent and their weighted average is reported. The weights are calculated using the statistical uncertainties combined with the systematic uncertainty due to the limited sample size of the simulation samples. The systematic uncertainties due to luminosity, kinematics, track reconstruction efficiency and kaon PID efficiency are entirely or strongly correlated, while those due to simulation sample size, muon and pion PID efficiencies, trigger selection and branching fractions are uncorrelated between the two decay modes. The double-differential cross-section of the averaged B + production in four rapidity bins as a function of p T and integrated over p T as a function of rapidity are shown in Fig. 5 and reported in Table 4. The same quantities for B 0 production are displayed in Fig. 6 and listed in Table 4. The measured cross-sections increase towards central rapidity both at positive and at negative rapidity. A good precision is achieved in the B + sample due to the averaging over two decay channels, which allows for improved precision with respect to the measurement in each single B + decay mode. The double-differential cross-section of Λ 0 b production is shown in Fig. 7 in four rapidity bins as a function of p T and integrated over p T as a function of rapidity, and is listed in Table 4. The trend observed as a function of the two variables is similar to that of the B mesons.    In order to probe the hadronization in proton-lead collisions, ratios of B 0 over B + and Λ 0 b over B 0 production cross-sections are studied with results shown in Fig. 8. Both ratios show no significant rapidity dependence within experimental uncertainties. The ratio between meson species is consistent with being independent of y and p T of the beauty hadrons. Most interestingly, the baryon-to-meson ratio shows a p T dependence with a significantly lower value at the highest p T compared to the p T -integrated measurement. However, the current uncertainties do not allow to draw firm conclusions. The production ratio, averaged over the kinematic range in the analysis, is measured to be 0.41 ± 0.06 (0.39 ± 0.05) for the pPb (Pbp) sample. The value is consistent with that measured by the LHCb collaboration in pp collisions [22][23][24][25].
The cross-sections are used to calculate forward-backward ratios and nuclear modification factors. In the following, the experimental results on these nuclear modification observables are compared with calculations using the HELAC-onia generator [51][52][53] with two different nuclear parton distribution function (nPDF) sets, nCTEQ15 [6] and    Table 5: Forward-backward ratios, R FB , of B + , B 0 and Λ 0 b production in bins of p T and integrated over 2.5 < |y| < 3.5. The first uncertainty is statistical and the second systematic.

Forward-backward ratios
The forward-backward production ratio of B + mesons is shown in Fig. 9 as a function of p T and y, while the corresponding values are reported in Table 5. A significant suppression of the production in the pPb sample with respect to that in the Pbp data is measured at the level of 20% when integrating over p T . Within the experimental uncertainty, no dependence as a function of p T is observed. The HELAC-onia calculations using EPPS16 and nCTEQ15 are in agreement with the experimental data. The EPPS16 * set exhibits the smallest uncertainties and is also in agreement with data.
The R FB ratio as a function of p T for B 0 mesons and the p T -integrated value is shown in Fig. 10 and given in Table 5. A significant suppression is observed when integrating over the considered p T range, consistent with the value measured for B + mesons. No significant dependence on p T is seen within the current experimental uncertainties.
In Fig. 11, the forward-backward cross-section ratio, R FB , of Λ 0 b production is shown. The numerical values are summarized in Table 5. The observed central value of R FB for the Λ 0 b baryon is consistent with the measured value for the two b-meson species and with the no-suppression hypothesis. A significant suppression of Λ 0 b production in pPb data compared to Pbp data is observed for the most precisely measured bin, between 7 and 12 GeV/c. The R FB measurement of Λ 0 b baryons is consistent with the modifications observed for the beauty mesons within the uncertainties for all kinematic bins. In Fig. 12, the values of R FB as a function of p T and as a function of y for the three hadrons are compared directly.   Figure 9: Forward-backward ratio, R FB , for B + mesons as a function of (left) p T and (right) y in proton-lead collisions compared with HELAC-onia calculations using different nPDF sets. For the data points, the vertical bars (boxes) represent the statistical (total) uncertainties. The rapidity range 2.5 < |y| < 3.5 represents the common range between the pPb and Pbp data samples studied in this analysis.  Figure 10: Forward-backward ratio, R FB , of B 0 mesons as a function of (left) p T and as a function of (right) y in proton-lead collisions. The vertical bars (boxes) represent the statistical (total) uncertainties.

Nuclear modification factors
In order to gain insight into potential modifications of the b-quark hadronization in pPb and Pbp collisions with respect to pp collisions, the Λ 0 b /B 0 cross-section ratio shown in Fig. 8 is divided by the corresponding measurement in pp collisions at √ s = 7 TeV [23]. Neglecting the dependence on the collision energy of the hadronization with respect to the experimental uncertainties, the quantity corresponds to the ratios of nuclear modification factors  If the overall nuclear effects for B 0 mesons and Λ 0 b baryons are identical, R Λ 0 b /B 0 pPb is expected to be unity. This double ratio is presented as a function of p T and y in Fig. 13 and in Table 6. At positive rapidity, the value of the ratio in all kinematic bins is consistent with unity. At negative rapidity (Pbp), the lowest p T bin exhibits a value smaller than one by more than two standard deviations and the third bin exceeds one by about two standard deviations. The p T -integrated value in the rapidity range −3.5 < y < −2.5 is about two standard deviations away from unity. However, more data are required to test whether there are different nuclear effects in beauty mesons and baryons. It would be interesting to check from the theory side whether deviations from unity are expected from models of quark recombination effects in heavy-flavor production in heavy-ion collisions.
The R p Pb modification factor for B + production is shown in Fig. 14, with the numerical values given in Table 7. The values are reported integrated over the considered p T range for the two y intervals, −3.5 < y < −2.5 and 2.5 < y < 3.5. They are also given as a 1.11 ± 0.14 ± 0.03 0.97 ± 0.17 ± 0.05 ( 7,12) 0.91 ± 0.13 ± 0.03 1.44 ± 0.21 ± 0.07 (12,20) 0.81 ± 0.21 ± 0.03 0.89 ± 0.22 ± 0.07 ( 2,20) 0.92 ± 0.09 ± 0.03 0.78 ± 0.11 ± 0.04 function of p T for both pPb and Pbp collisions. For the pp reference cross-section, an interpolation between existing pp cross-section measurements by the LHCb collaboration at 7 TeV [54] and 13 TeV [55] is performed. A power-law function is used following the approach of Refs. [17,56,57], which yields a prediction of B + production at √ s = 8.16 TeV consistent with an extrapolation of the measured value at √ s = 7 TeV using a FONLL calculation [58,59]. The interpolation takes into account the correlations provided in Ref. [55]. The measurement of R p Pb for nonprompt J/ψ production at the same collision energy by the LHCb collaboration [17] is also shown.
At positive rapidity, a significant suppression by more than 20% is observed integrating over the whole p T range, whereas at negative rapidity, the result is consistent with unity. The measurement is also consistent with that of nonprompt J/ψ production obtained in a similar kinematic range. The p T -differential result at positive rapidity shows a significant suppression, at the level of 25% for the lowest p T bin. The ratio tends to increase for high p T , however, the current experimental uncertainties that grow also with p T do not allow to establish a significant p T dependence. At negative rapidity, all values are consistent with a nuclear modification factor of one. The experimental data points are in good agreement with the three considered nPDF sets. At positive rapidity, the experimental uncertainties are smaller than the nPDF ones for the integrated values as well as for the three lowest   Figure 14: Nuclear modification factor, R pPb , for B + mesons as function of (top) y and as a function of p T in (bottom left) pPb and (bottom right) Pbp compared with HELAC-onia calculations using different nPDF sets as well as with the measurement of R pPb for nonprompt J/ψ production. For the data points, the vertical bars (boxes) represent the statistical (total) uncertainties. Table 7: Nuclear modification factor, R pPb , of B + production in pPb and Pbp collisions, in bins of p T for the range 2.5 < |y| < 3.5. The first uncertainty is statistical and the second systematic.

Conclusions
The differential production cross-sections of B + , B 0 mesons and Λ 0 b baryons in proton-lead collisions at √ s NN = 8.16 TeV are measured in the range 2 < p T < 20 GeV/c within the rapidity ranges 1.5 < y < 3.5 and −4.5 < y < −2.5. The cross-sections and the derived nuclear modification factors and forward-backward ratios of b-hadron production are measured for the first time with exclusive decay modes at transverse momenta smaller than the mass of the hadrons. They represent the first measurement of beauty-hadron production with different exclusive decay channels in nuclear collisions in that kinematic regime. The results with fully reconstructed beauty hadrons confirm the significant nuclear suppression of beauty-hadron production at positive rapidity measured via nonprompt J/ψ mesons. The observed experimental uncertainties at positive rapidity are smaller than those achieved in a weighting of nPDFs with heavy-flavor data. Therefore, this measurement can serve as a valuable input for future fits of nPDF, assuming that modifications of nPDFs are the dominant source of nuclear effects in proton-lead collisions at the LHC. Finally, the unique measurement of Λ 0 b production constrains the fragmentation of the beauty quark in a nuclear environment. The baryon-to-meson cross-section ratio in proton-lead collisions is found to be compatible with the equivalent ratio measured in pp collisions, and more data will be needed to study whether nuclear effects modify beauty baryon and meson production differently. These findings are important steps towards a better understanding of heavy-flavor production in nuclear collision systems and will serve as an input for the characterization of the quark-gluon plasma with heavy-flavor observables.