Neutral pion and $\eta$ meson production in p–Pb collisions at $\sqrt{s_{\mathrm{NN}}}$  = 5.02 TeV

Abstract

Neutral pion and $\eta$ meson invariant differential yields were measured in non-single diffractive p–Pb collisions at $\sqrt{s_{\mathrm{NN}}}$  = 5.02 TeV with the ALICE experiment at the CERN LHC. The analysis combines results from three complementary photon measurements, utilizing the PHOS and EMCal calorimeters and the Photon Conversion Method. The invariant differential yields of $\pi ^{0}$ and $\eta$ meson inclusive production are measured near mid-rapidity in a broad transverse momentum range of $0.3< p_{\mathrm{T}} < 20 \hbox { GeV}/c$ and $0.7< p_{\mathrm{T}} < 20\hbox { GeV}/c$, respectively. The measured $\eta /\pi ^{0}$ ratio increases with $p_{\mathrm{T}}$ and saturates for $p_{\mathrm{T}}$ > 4 $\hbox {GeV}/c$ at $0.483 \pm 0.015_{\mathrm{stat}} \pm 0.015_{\mathrm{sys}}$. A deviation from $m_{\mathrm{T}}$ scaling is observed for $p_{\mathrm{T}}<$ 2 $\hbox {GeV}/c$. The measured $\eta /\pi ^{0}$ ratio is consistent with previous measurements from proton-nucleus and pp collisions over the full $p_{\mathrm{T}}$ range. The measured $\eta /\pi ^{0}$ ratio at high $p_{\mathrm{T}}$ also agrees within uncertainties with measurements from nucleus–nucleus collisions. The $\pi ^{0}$ and $\eta$ yields in p–Pb relative to the scaled pp interpolated reference, $R_{{\mathrm{pPb}}}$, are presented for $0.3 < p_{\mathrm{T}}<$ 20 $\hbox {GeV}/c$ and $0.7 < p_{\mathrm{T}}<$ 20 $\hbox {GeV}/c$, respectively. The results are compared with theoretical model calculations. The values of $R_{{\mathrm{pPb}}}$ are consistent with unity for transverse momenta above 2 $\hbox {GeV}/c$. These results support the interpretation that the suppressed yield of neutral mesons measured in Pb–Pb collisions at LHC energies is due to parton energy loss in the hot QCD medium.

1 Introduction

Proton-nucleus (pA) collisions are an important tool for the study of strongly interacting matter and the Quark–Gluon Plasma (QGP), complementing and extending measurements carried out with high energy collisions of heavy nuclei [1]. By using a proton instead of a heavy nucleus as one of the projectiles, measurements of pA collisions have unique sensitivity to the initial-state nuclear wave function, and can elucidate the effects of cold nuclear matter on a wide range of observables of the QGP [2, 3].

Measurements of inclusive distributions of hadrons at mid-rapidity at the LHC probe parton fractional momentum x in the range $10^{-4}<x<10^{-2}$, where nuclear modification to hadronic structure is expected to be sizable [2]. This range extends an order of magnitude smaller in x with respect to other colliders. Inclusive hadron measurements are also essential to constrain theoretical models of particle production ([4] and references therein).

Within the framework of collinearly-factorized perturbative QCD (pQCD), effects of the nuclear environment are parameterized using nuclear-modified parton distribution functions (nPDFs) [5,6,7,8,9,10], which have been determined from global fits at next-to-leading order (NLO) to data from deep inelastic scattering (DIS), Drell–Yan, and $\pi ^{0}$ production. Inclusive hadron measurements at the LHC provide new constraints on gluon nPDFs [5, 9, 11], and the flavor dependence of sea-quark nPDFs [12]. Hadron production measurements at the LHC are likewise needed to improve constraints on fragmentation functions (FFs) [13,14,15].

An alternative approach to the theoretical description of hadronic structure is the Color Glass Condensate (CGC) [16], an effective theory for the nuclear environment at low x where the gluon density is high and non-linear processes are expected to play a significant role. CGC-based calculations successfully describe measurements of particle multiplicities and inclusive hadron production at high $p_{\mathrm{T}}$ in pp, d−Au and p–Pb collisions at RHIC and at the LHC [17,18,19]. CGC calculations, with parameters fixed by fitting to DIS data, have been compared to particle distributions at hadron colliders, thereby testing the universality of the CGC description [19]. Additional measurements of inclusive hadron production at the LHC will provide new constraints on CGC calculations, and help to refine this theoretical approach.

Recent measurements of p–Pb collisions at the LHC indicate the presence of collective effects in such systems, which influence inclusive hadron distributions [3, 20,21,22,23,24,25]. Detailed study of identified particle spectra over a broad $p_{\mathrm{T}}$ range can constrain theoretical models incorporating such effects. For example, the EPOS3 model [26] requires the inclusion of collective radial flow in p–Pb collisions to successfully describe the $p_{\mathrm{T}}$ spectrum of charged pions, kaons, protons, $\Lambda$ and $\Xi$ baryons [27, 28]. Tests of this model with neutral pions and $\eta$ mesons will provide additional constraints to this approach.

The shape of the invariant production cross section of various hadron species in pp collisions can be approximated by a universal function of $m_{\mathrm{T}} =\sqrt{p_{\mathrm{T}} ^2+M^2}$ (“$m_{\mathrm{T}}$ scaling”) [29] where M is the hadron mass. This scaling has been tested with many different collision energies and systems [30,31,32], and is commonly utilized to calculate hadronic distributions in the absence of measurements. Violation of $m_{\mathrm{T}}$ scaling at low $p_{\mathrm{T}}$ in pp collisions at the LHC has been observed for $\pi ^{0}$ and $\eta$ mesons at $\sqrt{s}$  = 7 TeV [33], and at $\sqrt{s}$  = 8 TeV  [34]; this may arise from collective radial flow that is indicated in pp collisions for $\sqrt{s}>0.9~\hbox {TeV}$ [35]. However, a deviation from $m_{\mathrm{T}}$ scaling at very low $p_{\mathrm{T}}$ has also been observed in pA collisions at $\sqrt{s_{\mathrm{NN}}}$  = 29.1 GeV [36], where it was attributed to enhanced low $p_{\mathrm{T}}$ pion production from resonance decays. The simultaneous measurement of $\pi ^{0}$ and $\eta$ mesons over a broad $p_{\mathrm{T}}$ range is therefore important to explore the validity of $m_{\mathrm{T}}$ scaling in pA collisions. Precise measurements of $\pi ^{0}$ and $\eta$ mesons at low $p_{\mathrm{T}}$ also provide an experimental determination of the background for measurements of dilepton and direct photon production [37, 38].

Strong suppression of inclusive hadron yields at high $p_{\mathrm{T}}$ has been observed in heavy-ion collisions at RHIC [39,40,41,42,43,44] and the LHC [45,46,47,48,49]. This suppression arises from partonic energy loss in the QGP [50,51,52,53]. Measurements of p–Pb collisions, in which the generation of a QGP over a large volume is not expected, provide an important reference to help disentangle initial and final-state effects for such observables [3, 54, 55]. Suppression of inclusive hadron production is quantified by measuring $R_{\mathrm{pA}}$, the relative rate of inclusive production in pA compared to pp, scaled to account for nuclear geometry. Measurements at RHIC and at the LHC report $R_{\mathrm{pA}}$ consistent with unity for $p_{\mathrm{T}}$ > 2 $\hbox {GeV}/c$ [27, 28, 56,57,58,59,60,61]. Additional, precise measurements of the inclusive hadron production in p–Pb collisions will provide a new test of this picture.

This paper presents the measurement of $\pi ^0$ and $\eta$ $p_{\mathrm{T}}$ differential invariant yields, together with the $\eta$/$\pi ^{0}$ ratio in non-single diffractive (NSD) p–Pb collisions at $\sqrt{s_{\mathrm{NN}}}$  = 5.02 TeV. The measurement covers a range of $|y_\mathrm{lab}|<0.8$, where $y_\mathrm{lab}$ is the rapidity in the laboratory reference frame. The measured $\pi ^{0}$ spectrum is corrected for secondary neutral pions from weak decays. The inclusive $\pi ^{0}$ and $\eta$ yield suppression ($R_{{\mathrm{pPb}}}$) is determined using a pp reference that was obtained by interpolating previous measurements by the ALICE experiment of $\pi ^{0}$ and $\eta$ meson production in pp collisions at $\sqrt{s}$ = 2.76 TeV [47, 62], at 7 TeV [33], and at 8 TeV [34]. The results are compared to theoretical models incorporating different approaches, including viscous hydrodynamics, pQCD at NLO with nuclear-modified PDFs, and a color glass condensate model, as well as commonly used heavy-ion event generators.

The paper is organized as follows: the detectors relevant for this analysis are described in Sect. 2; details of the event selection are given in Sect. 3; photon and neutral meson reconstruction, the systematic uncertainties as well as the calculation of the pp reference for the nuclear modification factor are explained in Sect. 4; the results and comparisons to the theoretical models are given in Sect. 5 followed by the conclusions in Sect. 6.

2 Detector description

A comprehensive description of the ALICE experiment and its performance is provided in Refs. [63, 64]. The $\pi ^{0}$ and $\eta$ mesons were measured via their two-photon decay channels $\pi ^{0}\rightarrow \gamma \gamma$ and $\eta \rightarrow \gamma \gamma$ (branching ratio $\hbox {BR} = 98.823 \pm 0.034\%$ and $39.41 \pm 0.20\%$, respectively), and in case of the $\pi ^{0}$ also via the Dalitz decay channel $\pi ^{0} \rightarrow \gamma ^{*}\gamma \rightarrow e^{+}e^{-}\gamma$ ($\hbox {BR} = 1.174 \pm 0.035\%$) including a virtual photon $\gamma ^{*}$ [65]. Photon reconstruction was performed in three different ways, using the electromagnetic calorimeters, the Photon Spectrometer (PHOS) [66] and the Electromagnetic Calorimeter (EMCal) [67], and the photon conversion method (PCM). The PCM used converted $e^{+}e^{-}$ pairs reconstructed using charged tracks measured in the Inner Tracking System (ITS) [68] and the Time Projection Chamber (TPC) [69]. Each method of photon and neutral meson reconstruction has its own advantages, specifically the wide acceptance and good momentum resolution of PCM at low $p_{\mathrm{T}}$, and the higher $p_{\mathrm{T}}$ reach of the calorimeters [33, 47, 62, 70]. The combination of the different analysis methods provides independent cross-checks of the results, a broader $p_{\mathrm{T}}$ range of the measurement, and reduced systematic and statistical uncertainties.

The PHOS [66] is a fine-granularity lead tungstate electromagnetic calorimeter that covers $|\eta _{{\mathrm{lab}}}| < 0.12$ in the lab-frame pseudorapidity and ${260^\circ< \varphi < 320^\circ }$ in azimuth angle. During the LHC Run 1 it consisted of three modules at a radial distance of $4.6~\hbox {m}$ from the ALICE interaction point. The PHOS modules are rectangular matrices segmented into $64 \times 56$ square cells of $2.2 \times 2.2~\hbox {cm}^2$ transverse size. The energy resolution of the PHOS is $\sigma _E /E = 1.8\%/E \oplus 3.3\%/ \sqrt{E} \oplus 1.1\%$, with E in units of GeV. The EMCal [67] is a lead-scintillator sampling electromagnetic calorimeter. During the period in which the analyzed dataset was collected, the EMCal consisted of 10 modules installed at a radial distance of $4.28~\hbox {m}$ with an aperture of $|\eta _{{\mathrm{lab}}}|<0.7$ and $80^\circ< \varphi < 180^\circ$. The energy resolution of the EMCal is $\sigma _E /E = 4.8\%/E \oplus 11.3\%/ \sqrt{E} \oplus 1.7\%$ with energy E in units of GeV. The EMCal modules are subdivided into $24 \times 48$ cells of $6 \times 6~\text{ cm }^2$ transverse size. The material budget of the active volumes of both calorimeters is about 20 radiation lenghts ($X_{}0$). The amount of material of the inner detectors between the interaction point and the calorimeters is about $0.2\,X_0$ for PHOS and ranges between $0.55\,X_0$ to $0.8\,X_0$ for EMCal, depending on the module. The relative cell energy calibration of both calorimeters was obtained by equalization of the $\pi ^{0}$ peak position reconstructed in each cell with high-luminosity pp collisions.

The Inner Tracking System (ITS) consists of six layers of silicon detectors and is located directly around the interaction point, covering full azimuth. The two innermost layers consist of Silicon Pixel Detectors (SPD) positioned at radial distances of 3.9 and $7.6~\hbox {cm}$, followed by two layers of Silicon Drift Detectors (SDD) at 15.0 and $23.9~\hbox {cm}$, and two layers of Silicon Strip Detectors (SSD) at 38.0 and $43.0~\hbox {cm}$. While the two SPD layers cover $|\eta _{{\mathrm{lab}}}|<$ 2 and $|\eta _{{\mathrm{lab}}}|<$ 1.4, respectively, the SDD and the SSD subtend $|\eta _{{\mathrm{lab}}}|<$ 0.9 and $|\eta _{{\mathrm{lab}}}|<$ 1.0, respectively. The Time Projection Chamber (TPC) is a large ($\approx 85~\hbox {m}^3$) cylindrical drift detector filled with a $\hbox {Ne/CO}_2$ (90/10%) gas mixture. It covers $|\eta _{{\mathrm{lab}}}|<$ 0.9 over the full azimuth angle, with a maximum of 159 reconstructed space points along the track path. The TPC provides particle identification via the measurement of the specific energy loss ($\text{ d }E/\text{ d }x$) with a resolution of 5.5%. The material thickness in the range R < $180~\hbox {cm}$ and $|\eta _{{\mathrm{lab}}}|<$ 0.9 amounts to (11.4 ± 0.5)% of $X_{}0$, corresponding to a conversion probability of $(8.6 \pm 0.4)$% for high photon energies [64]. Two arrays of 32-plastic scintillators, located at $2.8< \eta _{{\mathrm{lab}}} < 5.1$ (V0A) and $-3.7< \eta _{\mathrm{lab}} < -1.7$ (V0C), are used for triggering [71].

3 Event selection

The results reported here use data recorded in 2013 during the LHC p–Pb run at $\sqrt{s_{\mathrm{NN}}}$  = 5.02 TeV. Due to the 2-in-1 magnet design of the LHC [72], which requires the same magnetic rigidity for both colliding beams, the nucleon-nucleon center-of-mass system was moving with $y_{\mathrm{NN}}=0.465$ in the direction of proton beam. About $10^8$ p–Pb collisions were recorded using a minimum-bias (MB) trigger, which corresponds to an integrated luminosity of $50~\mu \text{ b }^{-1}$. The ALICE MB trigger required a coincident signal in both the V0A and the V0C detectors to reduce the contamination from single diffractive and electromagnetic interactions [73].

The primary vertex of the collision was determined using tracks reconstructed in the TPC and ITS as described in detail in Ref. [64]. From the triggered events, only events with a reconstructed vertex ($\sim$98.5%) were considered for the analyses. Additionally, the z-position of the vertex was required to be within $\pm 10~\hbox {cm}$ with respect to the nominal interaction point. The event sample selected by the above-mentioned criteria mainly consisted of non-single diffractive (NSD) collisions. The neutral meson yields were normalized per NSD collision, which was determined from the number of MB events divided by the correction factor $96.4\% \pm 3.1\%$ to account for the trigger and vertex reconstruction efficiency [61, 73]. This correction factor was determined using a combination of different event generators and taking into account the type of collisions used in the analyses. This correction is based on the assumption that non-triggered events contain no neutral mesons at mid-rapidity; see Ref. [73] for details.

Pile-up events from the triggered bunch crossing, which have more than one p–Pb interaction in the triggered events, were rejected by identifying multiple collision vertices reconstructed by the SPD detector. The fraction of such pile-up events in the analyzed data sample was at the level of 0.3%.

4 Data analysis

4.1 Photon and primary electron reconstruction

Photons and electrons hitting the PHOS or the EMCal produce electromagnetic showers which deposit energy in multiple cells. Adjacent fired cells with energies above $E^{\mathrm{min}}_{\mathrm{cell}}$ were grouped together into clusters. Noisy and dead channels were removed from the analysis prior to clusterization. The clusterization process started from cells with an energy exceeding $E_{\mathrm{seed}}$. The choice of the values of $E_{\mathrm{seed}}$ and $E^{\mathrm{min}}_{\mathrm{cell}}$ was driven by the energy deposited by a minimum ionizing particle, the energy resolution, noise of the electronics, and optimizing the signal to background ratio of meson candidates. For PHOS, $E_{\mathrm{seed}} = 50~\hbox {MeV}$ and $E^{\mathrm{min}}_{\mathrm{cell}} = 15~\hbox {MeV}$ were chosen. The corresponding thresholds for EMCal were $E_{\mathrm{seed}} = 500~\hbox {MeV}$ and $E^{\mathrm{min}}_{\mathrm{cell}} = 100~\hbox {MeV}$. The photon reconstruction algorithm in PHOS separates the clusters produced by overlapping showers from close particle hits, via a cluster unfolding procedure. Due to a low hit occupancy in the calorimeters in p–Pb collisions, relatively loose selection criteria were applied for clusters to maximize the neutral meson reconstruction efficiency and minimize systematic uncertainties from photon identification criteria. The minimum number of cells in a cluster was set to three and two for PHOS and EMCal, respectively, to reduce contributions of non-photonic clusters and noise. Consequently, the energy threshold for PHOS and EMCal clusters was set to 0.3 and $0.7~\hbox {GeV}$, respectively.

Apart from the cluster selection criteria described above, additional detector-specific criteria were applied in the PHOS and EMCal analyses to increase the purity and signal to background ratio of the photon sample. The EMCal clusters were selected in $|\eta _{{\mathrm{lab}}}| <0.67$ and $80^\circ<\varphi <180^\circ$, which is the full EMCal acceptance during the LHC Run 1 p–Pb run. In the EMCal analysis, the purity of the photon sample was enhanced by rejecting charged tracks reconstructed in the TPC that are matched to a cluster in the EMCal. The matching criteria, based on the distance between the track and the cluster in $\eta$ and $\varphi$, depend on the track $p_{\mathrm{T}}$ to maximize purity at low $p_{\mathrm{T}}$ and statistics at high $p_{\mathrm{T}}$. The purity is further enhanced by requirements on the squared major axis of the cluster shape $\sigma _\mathrm{long}^{2}$ calculated as the principle eigenvalue of the cluster covariance matrix $s_{ij}$ via $\sigma _\mathrm{long}^{2} =(s_{\eta \eta }+s_{\varphi \varphi })/2 +\sqrt{(s_{\eta \eta }-s_{\varphi \varphi })^{2}/4+s^{2}_{\eta \varphi }}$ where $s_{ij} = \langle ij \rangle - \langle i \rangle \langle j \rangle$ are the covariance matrix elements, i, j are cell indices in $\eta$ or $\varphi$ axes respectively, $\langle ij \rangle$ and $\langle i \rangle$, $\langle j \rangle$ are the second and the first moments of the cluster cells weighted with the cell energy logarithm [62, 74]. Photon clusters in EMCal and PHOS were defined by the condition $0.1< \sigma _\mathrm{long}^2 < 0.5$ and $\sigma _\mathrm{long}^2 > 0.2$, respectively, which selected clusters with axial symmetry.

In addition to these requirements, a selection criterion on cluster timing was applied in order to exclude clusters from other bunch crossings. Since the minimum interval between colliding bunches was $200~\hbox {ns}$, $|t|<100~\hbox {ns}$ had to be fulfilled for PHOS. For EMCal the cell time of the leading cell of the cluster was required to be within $|t|<50~\hbox {ns}$ of the time of the triggered bunch crossing.

Photons converted into $e^{+}e^{-}$ pairs were reconstructed with a secondary-vertex algorithm that searches for oppositely-charged track pairs originating from a common vertex, referred to as $\hbox {V}^0$ [64]. Three different types of selection criteria were applied for the photon reconstruction: requirements on the charged track quality, particle identification criteria for electron selection and pion rejection, and requirements on the $\hbox {V}^{0}$ sample that exploit the specific topology of a photon conversion. Details of the PCM analysis and the selection criteria are described in Refs. [33, 47]. Electron identification and pion rejection were performed by using the specific energy loss $\text{ d }E/\text{ d }x$ in the TPC. Detailed requirements are listed in Table 1, where $n\sigma _{e}$ and $n \sigma _{\pi }$ are deviations of $\text{ d }E/\text{ d }x$ from the electron and pion expectation expressed in units of the standard deviation $\sigma _{e}$ and $\sigma _{\pi }$, respectively. In comparison to the previous analyses of the $\gamma \gamma$ decay channel ($\hbox {PCM}-\gamma \gamma$) [33, 47], the converted photon topology selection criteria were slightly modified to further increase the purity of the photon sample. The constant selection criterion on the $e^\pm$ transverse momentum with respect to the $\hbox {V}^0$ momentum, $q_{\mathrm{T}}$, was replaced by a two-dimensional selection in the ($\alpha$,$q_{\mathrm{T}}$) distribution, known as the Armenteros–Podolanski plot [75], where $\alpha$ is the longitudinal momentum asymmetry of positive and negative tracks, defined as $\alpha = (p_L^+ - p_L^-)/(p_L^+ + p_L^-)$. The fixed selection criterion on the reduced $\chi ^{2}$ of the converted photon fit to the reconstructed $\hbox {V}^0$ was changed to the $\psi _\mathrm{pair}$-dependent $\chi ^{2}$ selection, where $\psi _\mathrm{pair}$ is the angle between the plane that is perpendicular to the magnetic field (x–y plane) and the plane defined by the opening angle of the pair [76]. It is defined as $\psi _{{\mathrm{pair}}}= \text{ arcsin } \left( {{\Delta \theta }\over {\xi _{{\mathrm{pair}}}}}\right)$, where $\Delta \theta$ is the polar angle difference between electron and positron tracks, $\Delta \theta = \theta (\mathrm{e}^{+})-\theta (\mathrm{e}^{-})$, and $\xi _{{\mathrm{pair}}}$ is the total opening angle between them. For converted photons with vanishing opening angle between the $e^{+}e^{-}$ pair the $\psi _{{\mathrm{pair}}}$ distribution is peaked at zero, while it has larger or random values for virtual photons of the Dalitz decay or combinatorial background, respectively. The applied selection criteria on the converted photon for the $\hbox {PCM-}\gamma \gamma$ and $\hbox {PCM-}\gamma ^*\gamma$ decay channels are summarized in Table 1.

Table 1 Selection criteria of the converted photon reconstruction with PCM for the two-photon ($\hbox {PCM}-\gamma \gamma$) and the Dalitz decay channel ($\hbox {PCM}-\gamma ^{*} \gamma$)

Virtual photons ($\gamma ^{*}$) of the Dalitz decays were reconstructed from primary electrons and positrons with the ITS and the TPC for transverse momenta $p_{\mathrm{T}} >0.125$ $\hbox {GeV}/c$. Tracks were required to cross at least 70 TPC pad rows, with the number of TPC clusters to be at least 80% of the number expected from the geometry of the track’s trajectory in the detector. Track selection was based on $\chi ^{2}$ of the ITS and TPC clusters fit to the track. To ensure that the selected tracks came from the primary vertex, their distance of closest approach to the primary vertex in the longitudinal direction ($\hbox {DCA}_{z}$) was required to be smaller than $2~\hbox {cm}$ and $\hbox {DCA}_{xy} <0.0182~\hbox {cm} +0.0350~\hbox {cm}/p_{\mathrm{T}} ^{1.01}$ in the transverse plane with $p_{\mathrm{T}}$  given in $\hbox {GeV}/c$ which correspond to a 7 $\sigma$ selection [64]. In addition, in order to minimize the contribution from photon conversions in the beam pipe and part of the SPD, only tracks with at least one hit in any layer of the SPD were accepted. Electrons were identified by the TPC $\text{ d }E/\text{ d }x$ by requiring that tracks fall within $-4< n \sigma _{e} <5$ of the electron hypothesis. For the pion rejection at intermediate $p_{\mathrm{T}}$ the same $n \sigma _{\pi }$ selection as described for the conversion electron tracks was used while at high $p_{\mathrm{T}}$ the selection was not applied, to increase the efficiency.

For the neutral meson reconstruction via the Dalitz decay channel a $\gamma ^{*}$ is constructed from the primary $e^{+}e^{-}$ pairs and is treated as real $\gamma$ in the analysis, except with non-zero mass. The pion contamination in the primary electron sample was reduced by constraints on the $\gamma ^{*}$ invariant mass ($M_{\gamma ^{*}} < 0.015~\hbox {GeV}/c^{2}$ at $p_{\mathrm{T}} < 1~\hbox {GeV}/c$ and $M_{\gamma ^{*}} < 0.035~\hbox {GeV}/c^{2}$ at $p_{\mathrm{T}} > 1~\hbox {GeV}/c$) exploiting that most of the $\gamma ^{*}$ from $\pi ^{0}$ Dalitz decays have a very small invariant mass, as given by the Kroll–Wada formula [77]. Contamination of the $\gamma ^{*}$ sample by $\gamma$ conversions was suppressed by requiring the primary $e^{+}e^{-}$ pairs to satisfy $|\psi _\mathrm{pair}|<0.6 - 5\Delta \varphi$ and $0<\Delta \varphi <0.12$, where $\Delta \varphi =\varphi (e^{+})-\varphi (e^{-})$ is the difference between electron and positron azimuth angles.

4.2 Meson reconstruction

The $\pi ^{0}$ and $\eta$ meson reconstruction was done by pairing $\gamma \gamma$ or $\gamma ^{*}\gamma$ candidates and calculating their invariant mass in transverse momentum intervals. For simplicity, the notation PCM-EMC will stand for the method with one photon reconstructed via PCM and the second photon reconstructed in EMCal. PCM, EMC and PHOS refer to the methods with both photons reconstructed by the same methods. $\hbox {PCM-}\gamma ^{*}\gamma$ is the method of meson reconstruction via the Dalitz decay channel. In total, five different measurements (PCM, $\hbox {PCM-}\gamma ^{*}\gamma$, EMC, PCM-EMC and PHOS) were done for the $\pi ^{0}$ meson and three different ones (PCM, EMC and PCM-EMC) for the $\eta$ meson. The reconstruction of $\eta$ mesons is not accessible by PHOS due to the limited detector acceptance and the wider opening angle of the decay photons compared to the $\pi ^{0}$.

Examples of invariant mass distributions are shown in Figs. 1 and  2 for selected $p_{\mathrm{T}}$ intervals for $\pi ^{0}$ and $\eta$ mesons, respectively. The combinatorial background, estimated using the event mixing technique [78], was scaled to match the background outside the signal region and subtracted from the total signal. The shape of the combinatorial background was optimized by mixing events within classes of similar primary vertex position and for all methods except PHOS also similar photon multiplicity. In case of the EMC analysis a minimal opening angle selection between the two photons of 17 mrad between the cluster seed cells was applied, which corresponds to 1 cell diagonal at mid rapidity, in order to provide a good event mixed background description. For PCM and PCM-EMC an opening angle selection of 5 mrad was applied. The background-subtracted signal was fitted to reconstruct the mass position $(M_{\pi ^{0},\eta })$ and width of the $\pi ^{0}$ and $\eta$ mesons. In case of the PCM, $\hbox {PCM-}\gamma ^{*}\gamma$, EMC, and PCM-EMC analyses, the fit function consisted of a Gaussian function convolved with an exponential low-energy tail to account for electron bremsstrahlung [79] and an additional linear function to take into account any residual background. For the PHOS analysis a Gaussian function was used.

Fig. 1

The diphoton invariant mass distributions around the $\pi ^{0}$ mass for selected intervals in $p_{\mathrm{T}}$, without and with combinatorial background for each of the five measurements: PHOS, EMC, PCM, $\hbox {PCM-}\gamma ^{*}\gamma$, and PCM-EMC. The vertical lines correspond to the limits of the region used to compute the integration of the meson signal

Fig. 2

The diphoton invariant mass distribution around the $\eta$ mass for selected intervals in $p_{\mathrm{T}}$, without and with combinatorial background for each of the three measurements: EMC, PCM, and PCM-EMC. The vertical lines correspond to the limits of the region used to compute the integration of the meson signal

The reconstructed $\pi ^{0}$ and $\eta$ meson peak position and width versus $p_{\mathrm{T}}$ compared to GEANT3 [80] simulations are shown in Figs. 3 and 4, respectively. The reconstructed meson mass peak position and width for each method are in good agreement for data and MC. The $\pi ^{0}$ and $\eta$ meson peak position for EMC and PCM-EMC was not calibrated to the absolute meson mass, but the cluster energy in MC was corrected by a $p_{\mathrm{T}}$ dependent correction factor such that the $\pi ^{0}$ mass peak positions in data and MC match within 0.4% for EMC and 0.5% for PCM-EMC. The cluster energy correction factor was calculated with $\pi ^{0}$ mesons reconstructed with the PCM-EMC method where the energy resolution of converted photons is much better than the one of real photons detected in EMC. Deviations of the MC $\pi ^{0}$ peak position with respect to the measured one in data were fully assigned to the EMC cluster energy. The $\pi ^{0}$ mass peak positions in PHOS were also tuned in MC to achieve a good agreement with data, which was done with a cluster energy correction.

The $\pi ^{0}$ and $\eta$ raw yields were obtained by integrating the background-subtracted $\gamma \gamma$ or $\gamma ^{*}\gamma$ invariant mass distribution. The integration window around the reconstructed peak of the meson mass was determined by the fit function. The integration ranges, as shown in Table 2, were selected according to the resolution of respective methods.

Fig. 3

Reconstructed $\pi ^{0}$ mass width (top) and position (bottom) versus $p_{\mathrm{T}}$ for data and Monte Carlo simulation for all five methods

Fig. 4

Reconstructed $\eta$ mass width (top) and position (bottom) versus $p_{\mathrm{T}}$ for data and Monte Carlo simulation for PCM, EMC and PCM-EMC

Table 2 Integration windows for the $\pi ^{0}$ and $\eta$ meson invariant mass distributions, where $M_{\pi ^{0}}$ and $M_{\eta }$ are the reconstructed mass positions from the fit, and M is the nominal mass of the respective meson

The raw $\pi ^{0}$ and $\eta$ meson yields were corrected for secondary $\pi ^{0}$ mesons, reconstruction efficiency, and acceptance, to obtain the invariant differential yield [33, 47, 62]. The secondary $\pi ^{0}$ mesons from weak decays or hadronic interactions in the ALICE detector were subtracted by estimating the contribution in a cocktail simulation, using measured spectra of relevant particles as input. The $\hbox {K}^{0}_{S}$ meson is the largest source of secondary $\pi ^{0}$ mesons, followed by hadronic interactions. The contamination from secondaries is largest for low $p_{\mathrm{T}}$ and then steeply decreases with $p_{\mathrm{T}}$. This correction is of the order of 8.5, 4.4, 2.8, 7% at the corresponding lowest $p_{\mathrm{T}}$ and 1.4, 2.4, $<~1$, $<1~\%$ at high $p_{\mathrm{T}}$, for PHOS, EMC, PCM-EMC and PCM, respectively, and negligible for $\hbox {PCM}-\gamma ^{*}\gamma$. The PCM analysis is affected by events from bunch crossings other than the triggered one, referred to as out-of-bunch pile-up. In the PCM analysis a correction was applied, as described in Ref. [47], that is of the order of 10% for the lowest $p_{\mathrm{T}}$ bin and sharply declines to about 2% for high $p_{\mathrm{T}}$. The out-of-bunch pile-up contribution in PHOS, EMC and PCM-EMC is removed by time cuts. The $\hbox {PCM}-\gamma ^{*}\gamma$ analysis used Monte Carlo simulations to apply an additional correction for the remaining contamination ($\sim 2.5\%$) of the $\pi ^{0}$ $\rightarrow \gamma \gamma$ in the $\pi ^{0}$ $\rightarrow \gamma ^{*}\gamma$ decay channel. Furthermore, raw $\pi ^{0}$ and $\eta$ meson yield were corrected for acceptance and reconstruction efficiency using GEANT3 simulations with HIJING [81] (PCM and $\hbox {PCM}-\gamma ^{*}\gamma$) or DPMJET [82] (PHOS, EMC, PCM and PCM-EMC) as Monte Carlo event generators. The product of acceptance and efficiency was calculated in each $p_{\mathrm{T}}$ bin and normalized to unit rapidity and full azimuth angle $\Delta \varphi =2\pi$. A typical value of the acceptance $\times$ efficiency varies from a few percent to few tens of percent, depending on $p_{\mathrm{T}}$ and on the reconstruction method.

4.3 Systematic uncertainties

The systematic uncertainties of the $\pi ^{0}$ and $\eta$ invariant differential yields were evaluated as a function of $p_{\mathrm{T}}$ by repeating the analysis for variations on the selection criteria. The magnitude of the systematic uncertainty for each set of variations is quantified by the average of the largest significant positive and negative deviations, and is parametrized by a low order polynomial function to remove the statistical fluctuations. Tables 3 and  4 show all the sources of systematic uncertainties and their magnitude in two representative $p_{\mathrm{T}}$ bins for $\pi ^{0}$ and $\eta$ mesons, respectively. All contributions to the total systematic uncertainties within a given reconstruction method are considered to be independent and were added in quadrature. The systematic uncertainties of the $\eta$/$\pi ^{0}$ ratio were evaluated independently such that correlated uncertainties cancel out. All the sources to the total systematic uncertainty are briefly discussed in the following.

For each reconstruction method the material budget is a major source of systematic uncertainty. For the calorimeters the uncertainty comes from material in front of the PHOS and EMCal, resulting in $3.5\%$ for PHOS and $4.2\%$ for EMC. For the other methods, the material budget reflects the uncertainty in the conversion probability of photons [64], adding $4.5\%$ uncertainty for a reconstructed conversion photon.

The yield extraction uncertainty is due to the choice of integration window of the invariant mass distributions. The integration window is varied to smaller and larger widths to estimate the error. The yield extraction uncertainty for the $\pi ^{0}$ meson for the different methods is $\sim 2\%$, while for the $\eta$ meson it increases to $\sim 5\%$. The yield extraction uncertainty for PHOS is estimated by using the Crystal Ball function instead of a Gaussian to extract the yields, resulting in a contribution to the total systematic uncertainty of $2.2\%$ for low $p_{\mathrm{T}}$ and $2.5\%$ for higher $p_{\mathrm{T}}$.

The PCM $\gamma$ reconstruction uncertainty is estimated by varying the photon quality and Armenteros–Podolanski selection criteria. For PCM it is $0.9\%$ at low $p_{\mathrm{T}}$ and increases to $3\%$ for high $p_{\mathrm{T}}$. The uncertainty on the identification of conversion daughters in PCM is done by varying the TPC PID selection criteria. For PCM it is $0.8\%$ at low $p_{\mathrm{T}}$ and increases to $2.4\%$ for high $p_{\mathrm{T}}$, and for $\hbox {PCM}-\gamma ^{*} \gamma$ it is $2.7\%$ at low $p_{\mathrm{T}}$ and decreases to $2.3\%$ for high $p_{\mathrm{T}}$. The track reconstruction uncertainty is estimated by varying the TPC track selection criteria. This uncertainty slightly increases with increasing $p_{\mathrm{T}}$ and is $\sim$1%. The secondary $e^+/e^-$ rejection uncertainty reflects the uncertainty of the real conversion rejection from the $\gamma ^{*}$ sample and is only present in $\hbox {PCM}-\gamma ^{*} \gamma$. It is obtained varying the selection on $\psi _{\mathrm{pair}}\hbox {-}\Delta \varphi$ or requiring a hit in the second ITS pixel layer. The Dalitz branching ratio uncertainty ($3.0\%$) is taken from the PDG [65].

Table 3 Relative systematic uncertainties $(\%)$ of the $\pi ^{0}$ spectrum for the different reconstruction methods

The uncertainty on the cluster energy calibration is estimated from the relative difference between data and simulation of the $\pi ^{0}$ mass peak position and also includes the uncertainty from the cluster energy corrections for both calorimeters. In the PHOS analysis, the energy calibration is also verified by the energy-to-momentum E / p ratio of electron tracks reconstructed in the central tracking system. The residual deviation of $\pi ^{0}$ mass and E / p ratio of electrons is attributed to the systematic uncertainty of the energy calibration which contributes $4.9\%$ at low $p_{\mathrm{T}}$ and increases to $6.2\%$ for high $p_{\mathrm{T}}$. The uncertainty of the neutral meson spectra caused by the energy calibration uncertainty in EMC is estimated as $1.7\%$ at low $p_{\mathrm{T}}$ and increases to $2.5\%$ for high $p_{\mathrm{T}}$. The uncertainty on the cluster selection was estimated by varying the minimum energy, minimum number of cells and time of the EMCal clusterization process. For the EMC the $\sigma _{\mathrm{long}}$ selection and track matching criteria are varied to estimate the contribution to the cluster selection uncertainty. This uncertainty accounts for $4.6\%$ at low $p_{\mathrm{T}}$ and increases to $5.1\%$ for higher $p_{\mathrm{T}}$.

The $\pi ^{0}$($\eta$) reconstruction uncertainty is due to the meson selection criteria and was estimated by varying the rapidity window of the meson and the opening angle between the two photons. It is a minor contribution to the total error with a magnitude of $\sim 1\%$. A $p_{\mathrm{T}}$ dependent uncertainty from 2% at $1~\hbox {GeV}/c$ to smaller than 0.5% for $p_{\mathrm{T}}$ larger than $2~\hbox {GeV}/c$ is assigned for PHOS to the secondary $\pi ^{0}$ correction, and the other methods were not significantly affected by this contribution. The generator efficiency uncertainty quantifies the difference between different Monte Carlo generators that are used to calculate the reconstruction efficiency of the $\pi ^{0}$ and $\eta$ meson and affects photon reconstruction with the EMCal. It contributes $2.0\%$ to the $\pi ^{0}$ meson systematic uncertainty and $4.0\%$ to the $\eta$ meson systematic uncertainty. The uncertainty on the acceptance correction for PHOS is estimated to be $2.2\%$ and includes the uncertainty introduced by the bad channel map. For EMC this uncertainty is included in the generator efficiency correction.

For PCM and $\hbox {PCM}-\gamma ^{*} \gamma$, the uncertainty on the background estimation is evaluated by changing the event mixing criteria of the photons from using the $\hbox {V}^{0}$ multiplicity to using the charged track multiplicity. For PCM this contributes $0.1\%$ ($0.3\%$) for the $\pi ^{0}$ ($\eta$) meson and for $\hbox {PCM}-\gamma ^{*}\gamma$ it contributes $1.8\%$ at low $p_{\mathrm{T}}$ and increases to $2.0\%$ for high $p_{\mathrm{T}}$. For PHOS, the uncertainty of the background is estimated using different polynomial functions to scale the mixed event background. The contribution is of the order of $4.6\%$, increasing slightly towards high $p_{\mathrm{T}}$. The systematic uncertainty due to the out-of-bunch pile-up subtraction is $1.0\%$ for PHOS and it varies from $3.0\%$ at $0.35~\hbox {GeV}/c$ to $0.3\%$ at high $p_{\mathrm{T}}$ for PCM.

Table 4 Relative systematic uncertainties $(\%)$ of the $\eta$ spectrum for the different reconstruction methods

4.4 pp reference

In order to quantify cold nuclear matter effects in p–Pb collisions, we require inclusive $\pi ^{0}$ and $\eta$ distributions in pp collisions at the same collision energy. However, such distributions are not available at present for pp collisions at $\sqrt{s}$ = 5.02 TeV. Therefore, the pp reference was calculated by interpolating between the measured spectra at midrapidity at $\sqrt{s}$  = 2.76 TeV [47, 62], at $\sqrt{s}$  = 7 TeV [33] and at $\sqrt{s}$  = 8 TeV [34] assuming a power-law behavior for the evolution of the cross section in each $p_{\mathrm{T}}$ bin as a function of $\sqrt{s}$ given by $\text{ d }^{2}\sigma (\sqrt{s})/\text{ d }y\text{ d }p_{\mathrm{T}} \propto \sqrt{s}^{\alpha (p_{\mathrm{T}})}$ [61], where the fit parameter $\alpha (p_{\mathrm{T}})$ increases with $p_{\mathrm{T}}$ which reflects the hardening of hadron spectra with collision energy. The method was cross-checked using events simulated by PYTHIA 8.21 [83], where the difference between the interpolated and the simulated reference was found to be negligible.

The invariant differential spectra [33, 34, 62] were fitted either with a Tsallis function [33, 84]:

$$\begin{aligned}&\frac{1}{2\pi N_\mathrm{ev}}\,\frac{\mathrm{d}^2 N}{p_{\mathrm{T}}\mathrm{d} p_{\mathrm{T}}\mathrm{d}y } = \frac{A}{2\pi }\nonumber \\&\quad \cdot \frac{(n-1)(n-2)}{n T(n T+M(n-2))} \left( 1+\frac{m_{T}-M}{nT}\right) ^{-n}, \end{aligned}$$

(1)

where M is the particle mass, $m_{\mathrm{T}} =\sqrt{M^2 + p_{\mathrm{T}} ^2}$, and A, n and T are fitting parameters; or with a two component model (TCM) as proposed in Ref. [85]:

$$\begin{aligned} \frac{1}{2\pi N_\mathrm{ev}}\,\frac{\mathrm{d}^2 N}{p_{\mathrm{T}}\mathrm{d} p_{\mathrm{T}}\mathrm{d}y } = {A}_{\mathrm{e}} \exp {(-{E}_{\mathrm{T, kin}}/{T}_{\mathrm{e}})}+ {A} \left( 1 + \frac{p_{\mathrm{T}} ^{2}}{{T}^{2} n} \right) ^{-n} \end{aligned}$$

(2)

where ${E}_{\mathrm{T, kin}}=\sqrt{p_{\mathrm{T}} ^2+M^2}-M$ is the transverse kinematic energy of the meson, with M the particle mass, ${A}_{\mathrm{e}}$ and A are normalization factors, ${T}_{\mathrm{e}}$, T and n are free parameters. The parametrizations of the $\pi ^{0}$ and $\eta$ spectra at the different collision energies using the Tsallis or TCM fits were needed due to the different $p_{\mathrm{T}}$ binning of the various pp and p–Pb spectra. The fits were then evaluated in the used p–Pb binning. The systematic uncertainty for each bin was calculated as average uncertainty of adjacent bins in the original binning. The statistical uncertainties of the parametrized spectra were computed from the fits to the measured spectra with only statistical errors.

The PHOS, PCM, EMC and PCM-EMC pp references are based solely on their contribution to the published spectra [33, 34, 47, 62] in order to cancel part of the systematic uncertainties in the calculation of $R_{{\mathrm{pPb}}}$. The $\hbox {PCM-}\gamma ^{*}\gamma$ method used the same pp reference as the PCM. The PCM $\pi ^{0}$ measurement at $\sqrt{s}$  = 2.76 TeV was extrapolated for $p_{\mathrm{T}}$ > 10 $\hbox {GeV}/c$ using the published fit. The PCM $\eta$ measurements were also extrapolated for $p_{\mathrm{T}}$ > 6–8 $\hbox {GeV}/c$ using the published fits. The difference between the $\pi ^{0}$ spectrum at $y=0$ and at $y=-0.465$ has been evaluated with PYTHIA 8.21 to be 1% for $p_{\mathrm{T}}$  $>~2~\hbox {GeV}/c$ and 0.5% at $0.5~\hbox {GeV}/c$. This correction was applied to the pp reference spectrum. In each $p_{\mathrm{T}}$ bin, the systematic uncertainty of the interpolated spectrum was estimated by the largest uncertainty among the input spectra used for the interpolation process. The statistical error is obtained from the power-law fit.

5 Results

5.1 Invariant yields of $\pi ^{0}$ and $\eta$ mesons

The ALICE $\pi ^{0}$ and $\eta$ meson invariant differential yields were determined by combining the individual meson measurements via a weighted average as described in Refs. [86, 87]. The correlations among the measurements for PCM, PCM-EMC, EMC, and $\hbox {PCM}-\gamma ^{*}\gamma$ were taken into account using the Best Linear Unbiased Estimate (BLUE) method [88, 89]. The PCM, PHOS and EMC measurements are completely independent and are treated as uncorrelated. Due to different $p_{\mathrm{T}}$ reach, statistics, and acceptance, the binning is not the same for the various methods. For the combined result, the finest possible binning was chosen. Thus, yields were combined bin by bin and methods that did not provide the yield for the specific bin were not taken into account.

The invariant differential meson yields were normalized per NSD event, with the normalization uncertainty added in quadrature to the combined systematic uncertainties.

Fig. 5

Left: invariant differential $\pi ^{0}$ and $\eta$ yields produced in NSD p–Pb collisions at -$1.365~< y_{\mathrm {cms}}<~0.435$ at $\sqrt{s_{\mathrm{NN}}}$  = 5.02 TeV. The statistical uncertainties are represented as vertical error bars whereas the systematic uncertainties are shown as boxes. In addition, Tsallis fits to the measured yields are shown. The TCM fit to the $\langle T_\mathrm {pPb} \rangle$ scaled $\pi ^{0}$ and $\eta$ pp reference spectra (see Sect. 5.3 for details) is shown. Right: Ratios of the measured data to their corresponding Tsallis fits

Table 5 Fit parameters and $\chi ^2/$NDF of the Tsallis fits to the combined $\pi ^{0}$ and combined $\eta$ meson invariant differential yields

Fig. 6

Ratio between individual $\pi ^{0}$ (left) and $\eta$ (right) invariant differential yield measurements, and Tsallis fit to the combined meson yield. The statistical uncertainties are represented as vertical error bars whereas the systematic uncertainties are shown as boxes

The invariant differential $\pi ^{0}$ and $\eta$ yields measured in NSD p–Pb collisions at $\sqrt{s_{\mathrm{NN}}}$  = 5.02 TeV are shown in Fig. 5. The horizontal location of the data points is shifted towards lower $p_{\mathrm{T}}$ from the bin center by a few MeV and illustrates the $p_{\mathrm{T}}$ value where the differential cross section is equal to the measured integral of the cross section over the corresponding bin [90]. For the $\eta$/$\pi ^{0}$ ratio and $R_{{\mathrm{pPb}}}$ the bin-shift correction is done in y-coordinates. Fits with a Tsallis function (Eq. 1) to the combined NSD $\pi ^{0}$ and $\eta$ spectra with statistical and systematic uncertainties added in quadrature are also shown in Fig. 5. In each case the Tsallis fit leads to a good description of the meson yield. The resulting fit parameters and the $\chi ^2/\hbox {NDF}$ are listed in Table 5 for the $\pi ^{0}$ and $\eta$ meson. The small values of $\chi ^2/\hbox {NDF}$ arise from the correlation of systematic uncertainties. The ratios between the meson yields obtained in the various reconstruction methods and the Tsallis fit to the combined spectrum for $\pi ^{0}$ and $\eta$ are presented in Fig. 6. All measurements are consistent within uncertainties over the entire $p_{\mathrm{T}}$ range. The invariant differential yield of neutral pions is consistent with that of charged pions [61] over the entire $p_{\mathrm{T}}$ range.

5.2 $\eta$/$\pi ^{0}$ ratio and $m_{\mathrm{T}}$ scaling

A combined $\eta$/$\pi ^{0}$ ratio was calculated and is presented in Fig. 7. For this purpose, the $\pi ^{0}$ was measured with the same binning as the $\eta$ meson with the PCM, EMC and PCM-EMC methods. The $\eta$/$\pi ^{0}$ ratio was determined for each method separately to cancel out the common systematic uncertainties and then combined taking into account the correlations among the measurements using the BLUE method. The $\eta$/$\pi ^{0}$ ratio increases with $p_{\mathrm{T}}$ and reaches a plateau of $0.483 \pm 0.015_{\mathrm{stat}} \pm 0.015_{\mathrm{sys}}$ for $p_{\mathrm{T}}$ $> 4~\hbox {GeV}/c$. This value agrees with the $\eta$/$\pi ^{0}$ ratio of $0.48 \pm 0.03$ ($0.47 \pm 0.03$) for $p_{\mathrm{T}}$ $> 2~\hbox {GeV}/c$ measured by PHENIX [30] in pp (d-Au) collisions at $\sqrt{s_{\mathrm{NN}}}$  = 200 GeV and with results from pA collisions at fixed-target experiments E515 [91] (p–Pb at $\sqrt{s}=23.8~\hbox {GeV}$, $\eta /\pi ^{0} =0.47\pm 0.03$) and E706 [92] (p-Be at $\sqrt{s}=31.6~\hbox {GeV}$, $\eta /\pi ^{0} =0.45\pm 0.01$ and at $\sqrt{s}=38.8~\hbox {GeV}$, $\eta /\pi ^{0} =0.42\pm 0.01$). A comprehensive compilation of all measured $\eta$/$\pi ^{0}$ ratios [30] shows that this ratio reaches an asymptotic value of $R_{\eta /\pi ^{0}}\sim 0.4-0.5$ at high $p_{\mathrm{T}}$ in hadronic collisions. Figure 7 shows a good agreement between the $\eta$/$\pi ^{0}$ ratio measured in p–Pb and pp collisions at $\sqrt{s_{\mathrm{NN}}}$  = 5.02 TeV and $\sqrt{s}$  = 7 TeV with ALICE [33], respectively. To illustrate universality of the $\eta$/$\pi ^{0}$ ratio and its independence of the collision system or energy, Fig. 7 also shows the $\eta$/$\pi ^{0}$ ratio measured in d–Au collisions at $\sqrt{s_{\mathrm{NN}}}$  = 200 GeV with PHENIX [30] and in fixed-target p–Be and p–Au collisions at $\sqrt{s_{\mathrm{NN}}}$  = 29.1 GeV by the joint TAPS/CERES collaboration [36] in their corresponding $p_{\mathrm{T}}$ coverage.

Fig. 7

Left: $\eta$/$\pi ^{0}$ ratio as function of $p_{\mathrm{T}}$ measured in NSD p–Pb collisions at -$1.365< y_{\mathrm {cms}}< 0.435$ at $\sqrt{s_{\mathrm{NN}}}$  = 5.02 TeV. The statistical uncertainties are shown as vertical error bars. The systematic uncertainties are represented as boxes. For comparison, also the $\eta$/$\pi ^{0}$ ratios measured in $7~\hbox {TeV}$ pp collisions with ALICE [33], in d-Au collisions at $\sqrt{s_{\mathrm{NN}}}$  = 200 GeV with PHENIX [30], and in p–Au and p–Be collisions at $\sqrt{s_{\mathrm{NN}}}$  = 29.1 GeV with TAPS/CERES [36] are shown, as well as the ratio where the $\eta$ yield is obtained via $m_{\mathrm{T}}$ scaling from the measured p–Pb $\pi ^{0}$ yield. Right: Ratio of the measured $\eta$/$\pi ^{0}$ ratio to the one obtained via $m_{\mathrm{T}}$ scaling

To test the validity of $m_{\mathrm{T}}$ scaling, a comparison of the measured ratio to the ratio obtained via $m_{\mathrm{T}}$ scaling is shown in Fig. 7. For this purpose, the $\eta$ yield was calculated from the Tsallis parametrization to the combined $\pi ^{0}$ yield, $P_{\pi ^{0}}$, assuming $m_{\mathrm{T}}$ scaling $E\,\text{ d }^3 N^{\eta }/\text{ d }p^{3} = C_{m}\cdot P_{\pi ^{0}}\left( \sqrt{p^2_{\mathrm{T}}+m^2_{\eta }}\right)$, with $C_{m} = 0.483 \pm 0.015_{\mathrm{stat}} \pm 0.015_{\mathrm{sys}}$. The ratio of the $m_{\mathrm{T}}$-scaled $\eta$ yield to the $\pi ^{0}$ Tsallis fit is shown in Fig. 7 as a red curve.

Above $p_{\mathrm{T}} \sim 4~\hbox {GeV}/c$ the measured ratio agrees with the $m_{\mathrm{T}}$-scaled distribution. At lower $p_{\mathrm{T}}$ the measured ratio deviates from the $m_{\mathrm{T}}$ scaling prediction, reaching a 40% difference at $p_{\mathrm{T}}$  = 1 $\hbox {GeV}/c$. The TAPS/CERES data also supports a deviation from $m_{\mathrm{T}}$ scaling at low $m_{\mathrm{T}}$ while the PHENIX data were found to be consistent with $m_{\mathrm{T}}$ scaling, although this measurement starts only at $p_{\mathrm{T}}$ $\sim$ 2 $\hbox {GeV}/c$. The $m_{\mathrm{T}}$ scaling is often utilized in measurements of electromagnetic probes [38, 93] to describe decay photon spectra from heavier neutral mesons. The measurement reported here demonstrates that $m_{\mathrm{T}}$ scaling is not valid for the $\eta$ meson at low $p_{\mathrm{T}}$. Therefore, a measured $\eta$ yield, especially at low $p_{\mathrm{T}}$, is crucial for the study of direct photons and dileptons in pA collisions, since $m_{\mathrm{T}}$ scaling from the measured $\pi ^{0}$ yield overestimates the $\eta$ yield at low $p_{\mathrm{T}}$ considerably [94]. Measurements of heavier neutral mesons such as $\omega$ in a wide $p_{\mathrm{T}}$ range are thus also desirable.

5.3 Nuclear modification factor $R_{{\mathrm{pPb}}}$

The ratio of the yield of $\pi ^{0}$ or $\eta$ in pA collisions relative to that in pp collisions, also known as nuclear modification factors ($R_{\mathrm{pA}}$), are calculated using

$$\begin{aligned} R_{\mathrm{pPb}} (p_{\mathrm{T}})=\frac{\mathrm {d}^{2}N^{\text {pPb}}_{\pi ^{0},\eta }/ \mathrm {d}y \mathrm {d}p_{\mathrm{T}}}{ \langle T_{\text {pPb}} \rangle \cdot \mathrm {d}^{2}\sigma ^{\text {pp}}_{\pi ^{0},\eta }/\mathrm {d}y \mathrm {d}p_{\mathrm{T}}}, \end{aligned}$$

(3)

where $\mathrm {d}^{2}N^{\text {pPb}}_{\pi ^{0},\eta }/\mathrm {d}y \mathrm {d}p_{\mathrm{T}}$ are the $\pi ^{0}$ and $\eta$ invariant yields measured in p–Pb  collisions and $d^{2}\sigma ^{\text {pp}}_{\pi ^{0},\eta }/\mathrm {d}y \mathrm {d}p_{\mathrm{T}}$ are the interpolated invariant $\pi ^{0}$ and $\eta$ meson cross sections in pp collisions at $\sqrt{s_{\mathrm{NN}}}$  = 5.02 TeV, as described in Sect. 4.4. $\langle T_{\text {pPb}}\rangle$ is the average nuclear overlap function, $\langle T_\mathrm {pPb} \rangle = 0.0983\pm 0.0035~\text{ mb }^{-1}$ [58, 73].

Fig. 8

$\pi ^{0}$ (left) and $\eta$ (right) nuclear modification factors $R_{{\mathrm{pPb}}}$ measured in NSD p–Pb collisions at -$1.365< y_{\mathrm {cms}}< 0.435$ at $\sqrt{s_{\mathrm{NN}}}$  = 5.02 TeV compared to the nuclear modification factors of charged pions and charged kaons, respectively. The statistical uncertainties are shown as vertical error bars and the systematic uncertainties are represented as boxes. The overall normalization uncertainty is given as the solid black box around unity

In the absence of nuclear effects, $R_{{\mathrm{pPb}}}$ is unity in the $p_{\mathrm{T}}$ region where hard processes dominate particle production. The values of $R_{{\mathrm{pPb}}}$ were calculated for each individual method to cancel out the common systematic uncertainties and then combined using the BLUE method (Fig. 8). For the Dalitz $R_{{\mathrm{pPb}}}$ the PCM pp reference is used. This induces a correlation of the Dalitz $R_{{\mathrm{pPb}}}$ with the $R_{{\mathrm{pPb}}}$ from PCM. The NSD normalization uncertainty is added in quadrature to the overall normalization uncertainty together with the uncertainties of the $T_\mathrm {pPb}$ and of the inelastic pp cross sections. The fit to the reference $\pi ^{0}$ and $\eta$ spectra in pp collisions at $\sqrt{s}=5.02~\hbox {TeV}$ scaled by $\langle T_\mathrm {pPb} \rangle$ are also displayed in Fig. 5. The fit parameters are given in Table 6.

The values of $R_{{\mathrm{pPb}}}$ are consistent with unity for transverse momenta above 2 $\hbox {GeV}/c$ for the $\pi ^{0}$ and $\eta$ mesons. The $R_{{\mathrm{pPb}}}$ measurements for neutral and charged pions as well as the $R_{{\mathrm{pPb}}}$ measurements for $\eta$ mesons and charged kaons [61] agree with each other within uncertainties over the complete $p_{\mathrm{T}}$ range as shown in Fig. 8.

5.4 Comparisons to theoretical models

Comparisons of the $\pi ^{0}$ and $\eta$ meson transverse momentum spectra to several theoretical calculations are shown in Fig. 9. In the following, we discuss each model individually, compared with the experimental data.

Fig. 9

Comparison of several theoretical calculations to the invariant differential $\pi ^{0}$ and $\eta$ yields produced in NSD p–Pb collisions at -$1.365< y_{\mathrm {cms}}< 0.435$ at $\sqrt{s_{\mathrm{NN}}}$  = 5.02 TeV from Fig. 5. Theoretical calculations are shown for the EPOS3 model [26], CGC model [19], pQCD calculations at NLO [6, 13, 95] using EPPS16 nPDF [96] or using the nCTEQ nPDF [10] and DSS14 FF [15] for the $\pi ^{0}$ and using nCTEQ nPDF [10] and AESSS FF [97] for the $\eta$ meson, hydrodynamic framework (labeled as VISHNU) [98] using the iEBE-VISHNU code [99], DPMJET model [82], and HIJING model [81]. The blue band on the EPOS3 calculation shows the statistical errors of the prediction. The gray band on the pQCD calculation includes the uncertainties on the factorization, renormalization and fragmentation scales, as well as on the nPDF and FF. The ratios of the measured data and several theoretical calculations to the data Tsallis fits (Fig. 5) are shown in the right panel

pQCD calculations at NLO [6, 13, 95] using the EPPS16 nPDF [96] with the CT14 PDF [100] or using the nCTEQ nPDF [10] and DSS14 FF [15] reproduce the $\pi ^{0}$ spectrum in Fig. 9, within the uncertainties due to the nPDF, the FF and variation of the factorization, renormalization and fragmentation scales. The largest contribution to the systematic uncertainty is due to the uncertainties in the choice of scales. Note that the EPPS16 nPDF has larger uncertainties than EPS09 nPDFs. pQCD calculations at NLO [95] using the nCTEQ nPDF [10] and AESSS FF [97] reproduce the $\eta$ meson spectrum at intermediate $p_{\mathrm{T}}$ while it overestimates the spectrum up to a factor two at high $p_{\mathrm{T}}$. Inclusive $\eta$ meson production has been measured in pp collisions at different LHC energies [33, 34, 62], which could be used to improve the $\eta$ meson FF [97] utilizing global fits, similar to a recent calculation for pions and kaons [15, 101].

The HIJING model [81] combines a pQCD-based calculation for multiple jet production with low $p_{\mathrm{T}}$ multi-string phenomenology. The model includes multiple minijet production, nuclear shadowing of parton distribution functions, and a schematic mechanism of jet interactions in dense matter. The Glauber model for multiple collisions is used to calculate pA and AA collisions. Figure 9 shows that the central value of the model calculation for inclusive $\pi ^{0}$ is about 20% smaller than the measured value at intermediate $p_{\mathrm{T}}$, between 1 and 4 $\hbox {GeV}/c$, while it agrees with the $\eta$ meson production in p–Pb collisions. At lower and higher $p_{\mathrm{T}}$ values the calculation overestimates the $\pi ^{0}$ and $\eta$ yields by up to 60–80%.

The DPMJET event generator [82] based on the Gribov-Glauber approach is an implementation of the two-component Dual Parton Model. This model treats the soft and the hard scattering processes in an unified way, using Reggeon theory for soft processes and lowest order pQCD for the hard processes. DPMJET was tuned to reproduce RHIC measurements of hadron production at low and moderate $p_{\mathrm{T}}$ by introducing a new mechanism of percolation and chain fusion, though it overestimates inclusive hadron yields at high $p_{\mathrm{T}}$ at RHIC energies [102]. Comparison of the $\pi ^{0}$ and $\eta$ meson measurements with DPMJET calculations in Fig. 9 shows that the model reproduces the distributions for $p_{\mathrm{T}}$ < 1 $\hbox {GeV}/c$, but underestimates the yields by 40% at higher $p_{\mathrm{T}}$. This suggests that the model parameters may need to be adjusted for the new energy domain. Comparison of DPMJET model predictions to particle production measurements in pp collisions at LHC energies also shows that the energy dependence of hadron production predicted by the model does not agree with data [103].

The $\pi ^{0}$ invariant differential yield computed with the CGC model [19] with $\hbox {MV}^\gamma$ [104] as the initial condition agrees with the measurements in Fig. 9 for $p_{\mathrm{T}}$ < 5 $\hbox {GeV}/c$. The deviation seen at high $p_{\mathrm{T}}$ is similar to that observed for inclusive $\pi ^{0}$ production in pp collisions at LHC.

The iEBE-VISHNU package [99] consists of a 3+1 viscous hydrodynamical model coupled to a hadronic cascade model [98]. Fluctuating initial conditions in the transverse plane are generated using a Monte–Carlo Glauber model. Figure 9 shows that this model reproduces the $\pi ^{0}$ and $\eta$ meson inclusive spectra for 0.7 < $p_{\mathrm{T}}$ < 1.5 $\hbox {GeV}/c$. For lower momenta ($p_{\mathrm{T}}$ < 0.7 $\hbox {GeV}/c$) the model prediction is lower than the measured $\pi ^{0}$ yield by up to a factor of two at 0.35 $\hbox {GeV}/c$. For $p_{\mathrm{T}}$ > 1.5 $\hbox {GeV}/c$ the model predictions underestimate the $\pi ^{0}$ and $\eta$ meson yields by a factor 5 at 3.5 $\hbox {GeV}/c$. This comparison shows that additional mechanisms not included in the model, in particular jet production, are important to describe particle production in p–Pb collisions in this region.

The EPOS3 [26] event generator is based on 3D+1 viscous hydrodynamics, with flux tube initial conditions that are generated in the Gribov–Regge multiple scattering framework. The reaction volume is divided into a core and a corona part. The core is evolved using viscous hydrodynamics. The corona is composed of hadrons from string decays. Figure 9 shows that this model reproduces the $\pi ^{0}$ inclusive $p_{\mathrm{T}}$ spectrum well over the full measured range. The model also reproduces the charged pion and kaon inclusive spectra in pA collisions [26]. However, the $\eta$ meson spectrum is well-reproduced only for $p_{\mathrm{T}}$ < 4 $\hbox {GeV}/c$, while at $p_{\mathrm{T}}$  > 4 $\hbox {GeV}/c$ the calculations lie above the data, with the disagreement reaching a factor of two at 10 $\hbox {GeV}/c$. Note that the VISHNU theoretical predictions [98] and EPOS3 are within 10-20 and 30–40% for the $\pi ^{0}$ and $\eta$ mesons, respectively, for $p_{\mathrm{T}}$ < 1.5 $\hbox {GeV}/c$. The comparisons to VISHNU and EPOS3 shows that a picture incorporating viscous hydrodynamic flow is consistent with measured particle yields at low $p_{\mathrm{T}}$ in p–Pb collisions.

Comparison of the measured high-precision $\pi ^{0}$ and $\eta$ meson spectra with theoretical models in Fig. 9 clearly shows that different underlying pictures can describe the data qualitatively. However, systematic uncertainties of the theoretical models are not provided, or are sizable. Hydrodynamic models agree with the data at low $p_{\mathrm{T}}$, while jet production appears to be needed for a good description at $p_{\mathrm{T}} >4$ $\hbox {GeV}/c$. While the high $p_{\mathrm{T}}$ part of the spectra can be described by NLO pQCD calculations, the precise data presented here will help to reduce their uncertainties significantly, for instance providing additional constraints on identified-particle FFs. Improved theoretical uncertainties are needed in order to discriminate among the models.

The comparison of the $\eta$/$\pi ^{0}$ ratio to different theoretical predictions is shown in Fig. 10. The DPMJET and HIJING model calculations are very close to the $m_{\mathrm{T}}$ scaling prediction, i.e. they lie above the measured ratio for $p_{\mathrm{T}}$ < 4 $\hbox {GeV}/c$, and agree with it at larger $p_{\mathrm{T}}$. On the other hand, the EPOS3 model calculation is closer to the data at low $p_{\mathrm{T}}$ than the $m_{\mathrm{T}}$ scaling prediction, while for $p_{\mathrm{T}}$  > 4 $\hbox {GeV}/c$ it continues to increase instead of reaching the plateau observed in data. The prediction from the VISHNU hydrodynamical calculation is in agreement with the measured data and very close to the EPOS3 prediction. However, this comparison may only be relevant up to $p_{\mathrm{T}}$ of 1.5 $\hbox {GeV}/c$ where the calculation was able to reproduce the measured neutral meson spectra. This behavior highlights once more the importance of hydrodynamical flow in p–Pb collisions at the LHC.

Fig. 10

Comparison of different theoretical calculations to the $\eta$/$\pi ^{0}$ ratio measured in NSD p–Pb collisions at -$1.365< y_{\mathrm {cms}}< 0.435$ at $\sqrt{s_{\mathrm{NN}}}$  = 5.02 TeV from Fig. 7. Theoretical calculations are shown for the EPOS3 model [26] with statistical errors shown as a band, hydrodynamic framework (VISHNU) [98] using the iEBE-VISHNU code [99], DPMJET model [82] and HIJING model [81]

A comparison of the measured $\pi ^{0}$ and $\eta$ $R_{{\mathrm{pPb}}}$ to different model predictions is shown in Fig. 11. The NLO pQCD calculations for the $\pi ^{0}$ [6, 13, 95] utilize the EPPS16 nuclear PDF [96] or the nCTEQ nPDF [10], and DSS14 FF [15], and for the $\eta$ meson [95] the nCTEQ nPDF [10] and AESSS FF [97] are used. The central values of the NLO predictions for $\pi ^{0}$ and $\eta$ lie below the data for $p_{\mathrm{T}}$ $< 6~\hbox {GeV}/c$. While the uncertainties of $\pi ^{0}$ calculations using nCTEQ are small and show sizable difference, the uncertainties for $\pi ^{0}$ calculations using EPPS16 are large and in agreement with the data. The CGC prediction from Ref. [19] uses the $k_{\mathrm{T}}$ factorization method and is able to reproduce the measured $R_{{\mathrm{pPb}}}$.

Fig. 11

Comparison of different theoretical model calculations to the $\pi ^{0}$ (left) and $\eta$ (right) nuclear modification factors $R_{\mathrm{pA}}$ measured in NSD p–Pb collisions at -$1.365< y_{\mathrm {cms}}< 0.435$ at $\sqrt{s_{\mathrm{NN}}}$  = 5.02 TeV. The grey band shows a pQCD calculation at NLO using the EPPS16 nPDF [96], the CT14 PDF [100] and the DSS14 FF [15] including systematic uncertainties. Color Glass Condensate predictions using the $k_{\mathrm{T}}$ factorization method are also shown. NLO calculations using nCTEQ nPDF [10], and DSS14 FF ($\pi ^{0}$) [15] or AESSS FF ($\eta$) [97] are also shown

6 Conclusions

The $p_{\mathrm{T}}$ differential invariant yields of $\pi ^{0}$ and $\eta$ mesons were measured in NSD p–Pb collisions at $\sqrt{s_{\mathrm{NN}}}$  = 5.02 TeV in the transverse momentum range $0.3< p_{\mathrm{T}} < 20$ and $0.7< p_{\mathrm{T}} <20$ $\hbox {GeV}/c$, respectively. State-of-the-art pQCD calculations at NLO are able to describe the $\pi ^{0}$ spectrum within the uncertainties of the nPDF and the pQCD scale, whereas they describe the $\eta$ spectrum at intermediate $p_{\mathrm{T}}$ and overestimate it up to a factor of two at high $p_{\mathrm{T}}$. As the wealth of the $\eta$ measurements is already sizable at the LHC, it will be important to include them in global fits to reach a similar theoretical progress in the pQCD calculations of the $\eta$ meson.

The $\eta$/$\pi ^{0}$ ratio is constant with a value of $0.483 \pm 0.015_{\mathrm{stat}} \pm 0.015_{\mathrm{sys}}$ at $p_{\mathrm{T}} > 4$ $\hbox {GeV}/c$ which is consistent with the $\eta$/$\pi ^{0}$ measurements at lower-energy pp, pA and AA collisions. Universality of the $\eta$/$\pi ^{0}$ behavior at high $p_{\mathrm{T}}$ suggests that the fragmentation into light mesons is the same in all collisions systems. At $p_{\mathrm{T}} <2$ $\hbox {GeV}/c$, the $\eta$/$\pi ^{0}$ ratio shows a clear pattern of deviation from the ratio predicted by the $m_{\mathrm{T}}$ scaling, confirming a $m_{\mathrm{T}}$ scaling violation observed earlier in pA collisions at $\sqrt{s_{\mathrm{NN}}}$  = 29.1 GeV and in pp collisions at $\sqrt{s}$  = 7 TeV and $\sqrt{s}$  = 8 TeV. The presence of radial flow effects in small systems and contributions from heavier-mesons decays to the $\eta$ and $\pi ^{0}$ production spectra are among possible interpretations of the $m_{\mathrm{T}}$ scaling violation. The comparison to different model calculations suggests that hydrodynamical flow may help to describe the measured spectra at low $p_{\mathrm{T}}$. Theoretical calculations using DPMJET and HIJING are very close to the $m_{\mathrm{T}}$ scaling prediction and therefore overestimate the measured ratio. The $\eta$/$\pi ^{0}$ ratio is reproduced in the complete $p_{\mathrm{T}}$ range by the VISHNU calculations although any conclusions above 1.5 $\hbox {GeV}/c$ are difficult to extract as the spectra were underestimated by large factors. For $p_{\mathrm{T}}$ < 3 $\hbox {GeV}/c$, the $\eta$/$\pi ^{0}$ ratio calculated by EPOS3 is closer to the measured data than the $m_{\mathrm{T}}$ scaling prediction, and it agrees with the data in the intermediate $p_{\mathrm{T}}$ range 2 < $p_{\mathrm{T}}$ < 5 $\hbox {GeV}/c$. These model comparisons support the interpretation that radial flow plays a role in collisions of small systems at the LHC.

The measured nuclear modification factors $R_{{\mathrm{pPb}}}$ for the $\pi ^{0}$ and $\eta$ meson are consistent with unity at $p_{\mathrm{T}} >2$ $\hbox {GeV}/c$ which confirms previously reported measurements at RHIC [56, 57] and LHC [27, 28, 58,59,60]. Theoretical calculations based on the latest nPDFs and a model based on the CGC framework are able to describe $R_{{\mathrm{pPb}}}$ well. These results support the interpretation that the neutral pion suppression in central Pb–Pb collisions is due to parton energy loss in the hot QCD medium.

These data are an important input for theoretical models aiming at the description of particle production in small systems at LHC energies and provide additional constraints on nPDFs and identified FFs.

A Parameters of TCM fits

The parameters of the two-component model fits to the reference $\pi ^{0}$ and $\eta$ meson spectra in pp collisions at $\sqrt{s}=5.02~\hbox {TeV}$ shown in Fig. 5 are given in the Table 6. The $\pi ^{0}$ and $\eta$ meson references were calculated using Eq. 3, from the combined spectra in p–Pb collisions (Fig. 5), and combined $R_{{\mathrm{pPb}}}$ (Fig. 8).

Table 6 Fit parameters of the TCM fits to the reference $\pi ^{0}$ and $\eta$ spectra in pp collisions at $\sqrt{s}=5.02~\hbox {TeV}$