Nuclear modification factor of neutral pions in the forward and backward regions in $p$-Pb collisions

The nuclear modification factor of neutral pions is measured in proton-lead collisions collected at a center-of-mass energy per nucleon of $8.16~{\rm TeV}$ with the LHCb detector. The $\pi^0$ production cross section is measured differentially in transverse momentum ($p_{\rm T}$) for $1.5<p_{\rm T}<10.0~{\rm GeV}$ and in center-of-mass pseudorapidity ($\eta_{\rm c.m.}$) regions $2.5<\eta_{\rm c.m.}<3.5$ (forward) and $-4.0<\eta_{\rm c.m.}<-3.0$ (backward) defined relative to the proton beam direction. The forward measurement shows a sizable suppression of $\pi^0$ production, while the backward measurement shows the first evidence of $\pi^0$ enhancement in proton-lead collisions at the LHC. Together, these measurements provide precise constraints on models of nuclear structure and particle production in high-energy nuclear collisions.

Neutral pion production is an important probe of nuclear effects in heavy ion collisions. In proton-lead (p-Pb) collisions at the Large Hadron Collider (LHC), π 0 production is particularly sensitive to cold nuclear matter (CNM) effects on the initial state of the bound nucleons in the colliding nucleus. Pion production in p-Pb collisions can be described using the collinear factorization framework where CNM effects are encoded into nuclear parton distribution functions (NPDFs) [1][2][3][4] and the parton-parton collision dynamics are described by perturbative quantum chromodynamics (PQCD) [5]. These NPDFs are determined using fits to data and are poorly constrained for partons with momentum fraction x smaller than about 10 −4 . At low x and low momentum transfer Q, the parton number density may become so large that parton saturation effects become sizable. Parton saturation would lead to nonlinear evolution of parton densities and is often described using the color glass condensate (CGC) effective field theory [6]. Measurements of π 0 production in p-Pb collisions at forward and backward rapidities with the LHCb detector can provide constraints on NPDFs for x between 10 −6 and 10 −1 , potentially helping identify the onset of parton saturation effects [7].
Measurements of angular correlations in small-collision systems at the LHC [8-12] and the Relativistic Heavy Ion Collider (RHIC) [13,14] show evidence of collective flow, which cannot be described using only collinear factorization and NPDFs. In addition, the LHCb experiment recently measured the nuclear modification factor of inclusive charged particles in p-Pb collisions at a center-of-mass energy per nucleon of √ s NN = 5 TeV, observing a much larger enhancement at backward pseudorapidities than predicted by NPDF calculations [15]. Similar enhancements have been observed in charged-particle production at RHIC [16][17][18][19] and have been attributed to the Cronin effect [20]. The Cronin effect is often described as initial-state multiple scattering within the nucleus [21], although radial collective flow and final-state recombination could produce similar enhancements [22,23]. However, measurements of inclusive charged-particle and π 0 production at central rapidity at the LHC show no evidence of a Cronin-like enhancement [24][25][26][27], potentially pointing to interplay between enhancing effects and low-x effects such as parton saturation [28,29]. Untangling the causes of the effects observed at RHIC and the LHCb experiment requires measurements of the nuclear modification factor of identified particles, such as neutral pions, over a broad rapidity range. Multiple scattering in the nucleus would result in similar enhancements for all particle species, while radial flow is expected to have a larger transverse momentum (p T ) hardening effect on higher mass particles included in charged-particle measurements [30]. Alternatively, recombination is expected to produce a particularly large baryon enhancement [23]. The PHENIX collaboration recently reported enhancements of π 0 production in small-collision systems, observing a system size dependence qualitatively consistent with radial flow, although still quantitatively consistent with NPDF predictions [31]. A study of π 0 production with the LHCb detector will help differentiate between contributions from nuclear parton density effects, initial state multiple scattering, and final-state collective flow.
This Letter presents a measurement of the nuclear modification factor of π 0 meson production in p-Pb collisions at √ s NN = 8. 16 TeV with the LHCb detector. The nuclear modification factor is given by where A = 208 is the atomic mass number of the lead nucleus and dσ p-Pb /dp T and dσ pp /dp T are the π 0 differential production cross sections in p-Pb and proton-proton (pp) collisions, respectively. The π 0 cross section is measured in p-Pb collisions at √ s NN = 8.16 TeV at both forward and backward rapidities relative to the proton beam. The LHCb detector collected only a small sample of unbiased pp collisions at √ s = 8 TeV, hence a pp reference cross section is constructed by interpolating between measurements performed using pp collisions at √ s = 5 and 13 TeV. The nuclear modification factor is measured for 1.5 < p T < 10.0 GeV in the pseudorapidity regions 2.5 < η c.m. < 3.5 and −4.0 < η c.m. < −3.0, where η c.m. is the pseudorapidity in the center-of-mass frame related to the laboratory frame pseudorapidity η lab by η c.m. = η lab − 0.465 in p-Pb collisions and η c.m. = η lab in pp collisions (natural units are used throughout this Letter). The LHCb detector is a single-arm forward spectrometer covering the pseudorapidity range 2 < η lab < 5, described in detail in Refs. [32,33]. The detector includes a high-precision tracking system consisting of a silicon-strip vertex detector (VELO) surrounding the interaction region, a large-area silicon-strip detector located upstream of a dipole magnet, and three stations of silicon-strip detectors and straw drift tubes placed downstream of the magnet. Photons and electrons are identified by a calorimeter system consisting of scintillating-pad and preshower detectors, an electromagnetic (ECAL) and a hadronic calorimeter. Particularly important for this analysis is the ECAL, which consists of alternating layers of lead and scintillator and has an energy (E) resolution of 13.5%/ E/GeV ⊕ 5.2% ⊕ (0.32 GeV)/E [34]. Simulated data samples are used to model the detector response to neutral pion reconstruction. In the simulation, p-Pb collisions are generated using EPOS-LHC [35], while pp collisions are generated using Pythia [36] with a specific LHCb configuration [37]. Decays of unstable particles are described by EvtGen [38] and final-state radiation is generated using Photos [39]. The interaction of the generated particles with the detector, and its response, are implemented using the Geant4 toolkit [40, 41] as described in Ref. [42].
Proton-lead collisions are recorded in the forward configuration in which the proton beam travels from the interaction region toward the spectrometer, and in the backward configuration in which the lead ion beam travels toward the spectrometer. The forward (backward) datasets correspond to an integrated luminosity of 328 ± 9 µb −1 (267 ± 7 µb −1 ). The forward configuration data are used to perform measurements for positive η c.m. , while the backward sample is used for negative η c.m. . All p-Pb events must have at least one reconstructed track in the VELO, while the pp samples consist of unbiased events.
Neutral pions are reconstructed using their decays to pairs of photons. To ensure that both photons are in the detector acceptance, π 0 candidates must have 2.5 < η lab < 4.0. Photons may reach the ECAL, or they may convert to electron-positron pairs in the detector material. Photons that reach the ECAL or that convert in the material downstream of the dipole magnet are reconstructed as energy deposits (clusters) in the ECAL and are referred to as ECAL photons. Photons that convert upstream of the dipole magnet are referred to as converted photons. Only π 0 candidates reconstructed from combinations of one ECAL photon and one converted photon are considered. This reconstruction method provides better momentum resolution than ECAL photon pairs alone and does not suffer from systematic effects due to the small opening angle between π 0 decay photons. ECAL photon clusters must contain at least two ECAL cells and must be far from the extrapolated trajectories of all reconstructed tracks. Converted photons are reconstructed from pairs of good-quality tracks with p > 2 GeV. To minimize the effects of energy LHCb pPb, √ s NN = 8.16 TeV 2.5 < η CM < 3.5 2.0 < p T (π 0 ) < 2.2 GeV Data Fit Comb. Brem. loss via bremsstrahlung radiation, only photons that convert downstream of the VELO are considered. Converted photons are required to have a small invariant mass, with a maximum that varies based on the position of the conversion vertex along the beam axis. The reconstructed converted photon momentum is also required to point to a reconstructed primary vertex. The ECAL photon must have p T > 400 MeV and the converted photon must have p T > 500 MeV. An example π 0 candidate invariant mass [M (γγ)] distribution in one interval of p T and η c.m. is shown in Fig. 1. The π 0 yield is determined using a binned maximum likelihood fit to the M (γγ) distribution. The π 0 signal is modeled by a two-sided Crystal Ball function [43]. The parameters describing the tails of the signal distribution are determined from simulation, while the Gaussian mean and width are left to vary in the fit to the data. The combinatorial background is modeled using charged particles as proxies for neutral pions. Charged particles reconstructed in the tracking system are given the π 0 mass, and their decays to two photons are simulated. The simulated photons are combined with reconstructed ECAL photons to form background candidates. The mass distribution of the proxy background candidates accurately describes combinatorial backgrounds in simulation and is used as a background model in the fit. An additional background component arises when a converted photon is combined with its own bremsstrahlung radiation to form a π 0 candidate, producing a peak at low mass. This background is modeled by convolving the reconstructed converted photon mass distribution with the sum of two positive half-normal distributions. The width of the narrower half-normal distribution is left to vary in the fit, while the larger width is fixed using simulation.
The π 0 yields are corrected for the effects of the detector response using simulation and an iterative unfolding procedure. First, correction factors are calculated for each p T interval from simulation and are used to correct the measured π 0 yields. The corrected π 0 spectrum in data and the true π 0 spectrum in simulation are then fit to Hagedorn functions [44]. The ratio of these fits is used to weight the true π 0 p T spectrum in simulation. The procedure is then repeated with the weighted simulated data sample. The procedure consistently converges after three iterations. Differences in reconstruction efficiency between data and simulation are measured and corrected for [45]. The ECAL photon reconstruction efficiency is the product of the cluster reconstruction efficiency and the photon cluster identification efficiency. The former is estimated with a tag-and-probe method using converted photons. The tag electron must be matched to an ECAL cluster and be identified as an electron. The cluster efficiency is then the fraction of probe electron tracks with a matching ECAL cluster. The photon cluster identification efficiency is evaluated by repeating the photon reconstruction without any photon cluster identification requirements in a subset of events in each data sample and measuring the fraction of π 0 mesons that pass the requirements. The converted-photon efficiency is measured using the ratio of the yields of π 0 mesons reconstructed as an ECAL photon and a converted photon to those reconstructed as two ECAL photons. This ratio is measured for 1.0 < p T < 1.5 GeV where ECAL photons from π 0 decays are always reconstructed as two separate clusters. The result is then combined with the measured ECAL photon efficiency correction factors to calculate a converted-photon efficiency correction. The corrected efficiencies range from about 0.5% at 1.5 GeV to about 5% at 10 GeV.
The π 0 cross section in pp collisions at √ s = 8. 16 TeV is estimated by interpolating between the measured pp cross sections at √ s = 5 and √ s = 13 TeV. The interpolation is performed independently in each p T interval. The cross section is interpolated using the functional form σ(s) = as b . A linear interpolation in √ s is also considered, as is a relative placement interpolation method [46]. In the latter, placement factors are calculated using simulation assuming a linear or power-law dependence of the cross section on √ s. The placement factors are then applied to data to determine the interpolated cross section.
The systematic uncertainties are summarized in Table 1. The largest sources of systematic uncertainty come from the interpolation between pp cross sections and the π 0 fit model. The maximum variation of interpolation results using the different pp interpolation methods is taken as the interpolation uncertainty. The fit model uncertainty is estimated by varying the default fit parameters by their uncertainties from the simulation. The variance of the resulting yields is taken as a systematic uncertainty. The fit model uncertainty also includes a contribution from the background model, which is estimated by repeating the fit using a polynomial background and taking the difference relative to the default result as a systematic uncertainty. An additional contribution to the fit model uncertainty is assigned to the yield extraction method itself by using an alternative method, where a background-only fit is performed in the mass region [0, 50] ∪ [250, 400] MeV. The background is then subtracted, and the background-subtracted distribution is integrated over [50,250] MeV to estimate the π 0 yield. The difference between the alternative and default yield extraction methods is taken as a systematic uncertainty.
Smaller systematic uncertainties come from unfolding, the luminosity estimate, and the efficiency correction factors. The unfolding uncertainty is estimated by using a closure test in simulation. The unfolded yields agree with the true yields in simulation to within about 1% in most p T intervals. An additional unfolding uncertainty arises from Table 1: Relative systematic and statistical uncertainties in dσ/dp T and R p-Pb in percent. The ranges correspond to the minimum and maximum values of the associated uncertainties across all p T intervals and both η c.m. regions. The dσ/dp T ranges cover the uncertainties for each of the p-Pb and pp samples. All sources of systematic uncertainty are fully correlated across p T intervals.

Source
dσ/dp T [%] differences in π 0 p T resolution in data and simulation. This difference is estimated to be less than 10% by comparing the fitted widths of the π 0 peaks in data and simulation. The resolution is varied in the unfolding by ±10%, resulting in a systematic uncertainty of less than 1% in every p T interval. The efficiency correction uncertainty arises from the finite size of the simulated data samples and results in a global uncertainty of about 1%-2%. The luminosity has been measured in pp collisions with a precision of 2% and in p-Pb collisions with a precision of 2.6% in the forward configuration and 2.5% in the backward configuration. The luminosity uncertainty is 50% correlated between datasets. The differential cross sections have an additional 4% uncertainty due to uncertainties in the detector material budget. This uncertainty is fully correlated between datasets and cancels in the nuclear modification factor. The fully corrected π 0 differential cross sections and nuclear modification factor are shown in Figs. 2 and 3, respectively. The nuclear modification factor shows a Cronin-like enhancement at backward pseudorapidity and a strong suppression at forward pseudorapidity. These measurements are compared to next-to-leading order PQCD calculations [47] using the EPPS16 [2] and nCTEQ15 [3] NPDF sets and the DSS14 π 0 fragmentation functions [48]. The NPDFs have been reweighted to incorporate LHCb D 0 production data [49-51], resulting in considerably smaller uncertainties than calculations using the default EPPS16 and nCTEQ15 NPDFs (see Supplemental Material [51]). The measurement uncertainties in the forward region are much smaller than the NPDF uncertainties, indicating that this measurement can provide powerful constraints on NPDFs at low x. In addition, the forward results present tension with the CGC calculation [52]. The enhancement in the backward direction between 2 and 4 GeV is larger than predicted by the PQCD calculation, suggesting that effects not described by NPDFs contribute to the enhancement.
The π 0 nuclear modification factor is also compared to the charged-particle nuclear modification factor measured by the LHCb experiment in p-Pb collisions at √ s NN = 5 TeV [15].
The forward π 0 measurement agrees with the charged-particle measurement, and the enhancement at backward η c.m. is smaller than that seen for charged particles at the LHCb experiment. Because the charged-particle measurement includes heavier mesons  In conclusion, the π 0 nuclear modification factor is measured at forward and backward rapidities in p-Pb collisions at √ s NN = 8. 16 TeV. The measured nuclear modification factor has a total uncertainty of less than 6% in most p T intervals, which will provide strong constraints on models of nuclear structure and particle production in heavy ion collisions. In particular, the forward measurement is sensitive to NPDFs for x as low as 10 −6 and can provide useful constraints in future NPDF analyses. Furthermore, the backward measurement shows the first evidence of enhanced π 0 production in proton-ion collisions at the LHC. These measurements will help constrain nuclear parton densities and models of collectivity in small-collision systems, as well as explain the origin of Cronin enhancement at the LHC.  [7] A. Morreale and F. Salazar, Mining for gluon saturation at colliders, Universe 7 (2021) 312, arXiv:2108.08254. [21] A. Accardi, Cronin effect in proton nucleus collisions: A Survey of theoretical models, arXiv:hep-ph/0212148.   [51] See Supplemental Material at the end of this Letter for additional plots.

Supplemental Material for LHCb-PAPER-2021-053
Nuclear modification factor of neutral pions in the forward and backward regions in pPb collisions Additional Comparisons to Theory Fig. S1 shows comparisons of the measured π 0 nuclear modification factor to NLO PQCD predictions calculated using the nominal EPPS16 [2] and nCTEQ15 [3] NPDF sets. These NPDF analyses do not incorporate LHCb D 0 data [50], which provide strong constraints at low x. As a result, calculations using these NPDF sets have much larger uncertainties than those incorporating LHCb measurements.