Transverse momentum distribution and nuclear modification factor of forward neutral pion in proton--lead collisions at $\sqrt{s_{NN}} = 5.02$TeV

The transverse momentum ($p_\text{T}$) distribution for inclusive neutral pions in the very forward rapidity region has been measured, with the Large Hadron Collider forward detector (LHCf), in proton--lead collisions at nucleon-nucleon center-of-mass energies of $\sqrt{s_{NN}} = 5.02$TeV at the LHC. The $p_\text{T}$ spectra obtained in the rapidity range $-11.0<y_\text{lab}<-8.9$ and $0<p_\text{T}<0.6$GeV (in the detector reference frame) show a strong suppression of the production of neutral pions after taking into account ultra-peripheral collisions. This leads to a nuclear modification factor value, relative to the interpolated $p_\text{T}$ spectra in proton-proton collisions at $\sqrt{s} = 5.02$TeV, of about 0.1--0.4. This value is compared with the predictions of several hadronic interaction Monte Carlo simulations.


I. INTRODUCTION
Measurements of particle production in hadronic interactions at high energies play an unique role in the study of strong interactions described by Quantum Chromodynamics (QCD). For example, as first discovered in measurements at HERA [1,2], it is still not well understood how a parton (dominantly gluon) density increases or even saturates when the momentum fraction that the parton itself carries is small (denoted as Bjorken-x). This even though the understanding of the behaviour of hadron constituents (partons) has been improving both theoretically and experimentally in the past few decades.
Such kind of phenomena at small Bjorken-x are known to be more visible in events at large rapidities. Since at large rapidities partons in projectile and target hadrons generally have large and small momentum fractions respectively, Bjorken-x in the target should be smaller than in mid rapidity events. Moreover, in case of nuclear target interactions, the parton density in the target nucleus is expected to be larger by ∼ A 1/3 . In these interactions partons in the projectile hadron would lose their energy while traveling in the dense QCD matter. Particle production mechanisms change accordingly when compared to those in nucleon-nucleon interactions.
These phenomena have been observed by several experiments at the Super Proton Synchrotron (SPS, CERN), at the Relativistic Heavy Ion Collider (RHIC, BNL), and at the Large Hadron Collider (LHC, CERN). The BRAHMS and STAR experiments at RHIC showed the modification of particle production spectra at forward rapidity [3,4] by comparing the p T spectra in deuteriumgold (d-Au) collisions at nucleon-nucleon center-of-mass energies √ s N N = 200 GeV with those in proton-proton (p-p) collisions at center-of-mass energies √ s = 200 GeV. Especially, a comparison of the experimental results between different rapidity regions by BRAHMS (pseudorapidity η = 2.2 and 3.2) and STAR (η = 4.0) indicates that particle production is strongly suppressed with increasing pseudorapidity. Also at the LHC the same suppression of hadron production has been found by the AL-ICE and LHCb experiments, both at mid and forward rapidities [5,6].
Here a naive question arises; in what way does particle production proceed within an extremely dense QCD matter at very forward rapidity regions which are thus characterised by p T 1 GeV ? Clearly all of the saturation based models cannot quantitatively predict the particle productions in this region which is characterised by a non-perturbative regime. Nevertheless, some mod-els exists that try to make quantitative predictions also near these regions. For example the black body hypothesis model predicts that meson production is strongly suppressed as a result of limiting fragmentation and that the p T distribution broadens [7]. Similarly hadronic interaction Monte Carlo (MC) simulations covering soft-QCD also predict a suppression of particle production. Experimental data that confirm the theoretical and phenomenological predictions and possibly constrain the remaining degrees of freedom in such models would thus be very welcome.
Understanding of particle production in very forward rapidity regions in nuclear target interactions is also of importance for ultrahigh-energy cosmic ray interactions, where parton density is expected to be much higher than at LHC energies due to the dependance of Bjorkenx ∝ 1/ √ s. In high energy cosmic-ray-observation, energy and chemical composition of primary cosmic rays are measured by analysing the cascade showers produced by the cosmic rays interacting with the nuclei in the earth atmosphere [8]. Secondary particles produced in the atmospheric interaction are, of course, identical to the forward emitted particles from the hadronic interactions at equivalent collision energy. In fact current modeling of particle production in nuclear interactions is limited by the available energy at the accelerators and is the cause of large systematic uncertainties in high energy cosmic ray physics.
The Large Hadron Collider forward (LHCf) experiment [9] is designed to measure the hadronic production cross sections of neutral particles emitted in very forward angles in p-p and proton-lead (p-Pb) collisions at the LHC, including zero degree. The LHCf detectors (see Sec. II) cover a pseudorapidity range larger than 8.4 and are capable of precise measurements of the forward high-energy inclusive-particle-production cross sections of photons, neutrons, and other neutral mesons and baryons. Therefore the LHCf experiment provides a unique opportunity to investigate the effects of high parton density which is the case in the small Bjorken-x region and in p-Pb collisions at high energies.
In the analysis presented in this paper, the focus is placed on the neutral pions (π 0 s) emitted into the direction of the proton beam (proton remnant side), the most sensitive probe to the details of the p-Pb interactions. From the LHCf measurements, the inclusive production rate and the nuclear modification factor for π 0 s in the rapidity range of −11.0 < y lab < −8.9 in the detector reference frame are then derived as a function of the π 0 transverse momentum.
The paper is organised as follows. Sec. II gives a brief description of the LHCf detectors. Section III and Sec. IV summarize the data taking conditions and the MC simulation methodology, respectively. In Sec. V the analysis framework is described, while the factors that contribute to the systematic uncertainties are explained in Sec. VI. Finally the analysis results are presented and discussed in Sec. VII. Concluding remarks are found in Sec. VIII.

II. THE LHCF DETECTOR
Two independent detectors called the LHCf Arm1 and Arm2 were assembled originally to study p-p collisions at the LHC [10,11]. In p-Pb collisions at √ s N N = 5.02 TeV, only the LHCf Arm2 detector was used to measure the secondary particles emitted into the proton remnant side. Hereafter the LHCf Arm2 detector is denoted as the LHCf detector for brevity. The LHCf detector has two sampling and imaging calorimeters composed of 44 radiation lengths of tungsten and 16 sampling layers of 3 mm thick plastic scintillators. The transverse sizes of the calorimeters are 25×25 mm 2 and 32×32 mm 2 . Four X-Y layers of silicon microstrip sensors are interleaved with the layers of tungsten and scintillator in order to provide the transverse profiles of the showers. Readout pitches of the silicon microstrip sensors are 0.16 mm [11].
The LHCf detector was installed in the instrumentation slot of the target neutral absorber (TAN) [12] located 140 m in the direction of the ALICE interaction point (IP2) from the ATLAS interaction point (IP1) and at a zero-degree collision angle. The trajectories of charged particles produced at IP1 and directed towards the TAN are bent by the inner beam separation dipole magnet D1 before reaching the TAN itself. Consequently, only neutral particles produced at IP1 enter the LHCf detector. The vertical position of the LHCf detector in the TAN is manipulated so that the LHCf detector covers the pseudorapidity range from 8.4 to infinity for a beam crossing angle of 145 µrad, especially the smaller calorimeter covers the zero-degree collision angle. After the operations in p-Pb collisions, the LHCf detector was uninstalled from the instrumentation slot of the TAN on April, 2013.
More details on the scientific goals of the experiment are given in Ref. [9]. The construction of the LHCf detectors (Arm1 and Arm2) is reported in Refs. [10,11] and the performance of the detectors has been studied in the previous reports [13,14].

III. EXPERIMENTAL DATA TAKING CONDITIONS
The experimental data used for the analysis in this paper were obtained in p-Pb collisions at √ s N N = 5.02 TeV with a 145 µrad beam crossing angle. The beam energies were 4 TeV for protons and 1.58 TeV per nucleon for Pb nuclei. Since the beam energies were asymmetric the nucleon-nucleon center-of-mass in p-Pb collisions shifted to rapidity = −0.465, with the proton beam traveling at θ = π and the Pb beam at θ = 0.
Data used in this analysis were taken in three different runs. The first (LHC Fill 3478) was taken on January 21, 2013 from 02:14 to 03:53. The second and third runs (LHC Fill 3481) were taken on January 21, 2013 from 21:03 to 23:36 and January 22, 2013 from 03:47 to 04:48, respectively. The integrated luminosity of the data was 0.63 nb −1 after correcting for the live time of data acquisition systems [15]. The average live time percentages for LHC Fill 3478 and 3481 were 12.1 % and 6.3 % respectively, the smaller live time percentage in Fill 3481 relative to Fill 3478 being due to a difference in the instantaneous luminosities. These three runs were taken with the same data acquisition configuration. In all runs the trigger scheme was essentially identical to that used in p-p collisions at √ s = 7 TeV. The trigger efficiency was greater than 99 % for photons with energies E > 100 GeV [14].
The multihit events that have more than one shower event in a single calorimeter may appear due to pileup interactions in the same bunch crossing and then could potentially cause a bias in the p T spectra. However, considering the acceptance of the LHCf detector for inelastic collisions ∼ 0.035, the multihit probability due to the effects of pileup is estimated only 0.4 % and is therefore producing a negligible effect. Detailed discussions of pileup effects and background events from collisions between the beam and residual gas molecules in the beam tube can be found in previous reports [14,16].
Also beam divergence can cause a smeared beam spot at the TAN leading to a bias in the measured p T spectra. The beam divergence at IP1 was ε/β * = 32 µrad [17] for the three fills mentioned, corresponding to a beam spot size at the TAN of roughly σ TAN = 4.5 mm. The effect of a non-zero beam spot size at the TAN is evaluated by comparing two p T spectra predicted by toy MC simulations; one assuming a beam spot size of zero and another assuming that the beam axis positions fluctuates following a Gaussian distribution with σ TAN = 4.5 mm. The π 0 yield at p T = 0.6 GeV is found to increase by a factor 1.8 at most. This effect is taken into account in the final results to the p T spectra.

IV. MONTE CARLO SIMULATIONS METHODOLOGY
MC simulations were performed in two steps: (I) p-Pb interaction event generation at IP1 explained in Sec. IV A and (II) particle transport from IP1 to the LHCf detector and consequent simulation of the response of the LHCf detector (Sec. IV B).
MC simulations which were then used for the validation of reconstruction algorithms and cut criteria, and the estimation of systematic uncertainties follow steps (I) and (II). These MC simulations are denoted as reference MC simulations. On the other hand, MC simulations used only for comparisons with measurement results in Sec. VII are limited to step (I) only, since the final p T spectra in Sec. VII are already corrected for detector responses and eventual reconstruction bias.

A. Signal modeling
Whenever the impact parameter of proton and Pb is smaller than the radius of each particle, soft-QCD induced events are produced. These p-Pb interactions at √ s N N = 5.02 TeV and the resulting flux of secondary particles emitted into forward rapidity region with their kinematics are simulated using various hadronic interaction models (dpmjet 3.04 [18], qgsjet II-03 [19], and epos 1.99 [20]). Dpmjet 3.04 also takes into account the Fermi motion of the nucleons in Pb nucleus and the Cronin effect [21]. Fermi motion enhances the π 0 yields at most by 5 % in the LHCf p T covered range p T < 0.6 GeV, while the Cronin effect is not significant in this p T range.
On the other hand, when the impact parameter is larger than the overlapping radii of each particle, socalled ultra-peripheral collisions (UPCs) can occur. In UPCs virtual photons are emitted by the relativistic Pb nucleus which can then collide with the proton beam [22]. The energy spectrum of these virtual photons follows the Weizsäcker-Williams approximation [23]. The sophia [24] MC generator is used to simulate the photon-proton interaction in the rest frame of the proton and then the secondary particles generated by sophia are boosted along the proton beam. For heavy nuclei with the radius R A , the virtuality of the photon |q 2 | < (ℏc/R A ) 2 can be neglected and the photons are regarded as real photon in the simulation in this analysis.
In these MC simulations, π 0 s from short-lived particles that decay within 1 m from IP1, mostly η mesons decaying into 3π 0 ( 10 % relative to all π 0 s), are also accounted for consistently with the treatment of the experimental data. B. Simulation of particle transport from IP1 to the LHCf detector and of the detector response The generated secondary particles are transported in the beam pipe from IP1 to the TAN, taking into account the bending of charged particles' trajectory by the Q1 quadrupole and the D1 beam separation dipole, particle decays, and particle interactions with the beam pipe wall and the Y-shaped beam-vacuum-transition chamber made of copper.
Finally the showers produced in the LHCf detector by the particles arriving at the TAN and the detector response are simulated with the cosmos and epics libraries [25]. The detector position survey data and random fluctuations due to electrical noise are taken into account. The simulations of the LHCf detector are tuned to the test beam data taken at SPS, CERN in 2007 and 2012 [13,26]. Since π 0 s decay into two photons very close to their point of creation at IP1, each photon's direction is geometrically calculated using the impact coordinates at the LHCf detector and the distance between IP1 and the detector itself. Photon four-momentum is then derived by combining the photon's energy as measured by the calorimeter with the previously obtained angle of emission. Candidate π 0 s are selected from events where the invariant mass of the two photons detected is within a narrow window around the π 0 rest mass.
The π 0 events are then classified in two categories: single-hit π 0 and multihit π 0 events. The former is defined as having a single photon in each of the two calorimeter towers only, while a multihit π 0 event is defined as a single π 0 accompanied by at least one additional background particle (usually a photon or a neutron) in one of the two calorimeters. In the analysis presented here, events having two particles within the same calorimeter tower (multihit events) are rejected when the energy deposit of the background particle is above a certain threshold [27]. Mostly then, only single-hit π 0 events are considered in this analysis. The final inclusive production rates reported at the end are corrected for this cut as described in Sec. V C.
Given the geometrical acceptance of the LHCf detector and to ensure a good event reconstruction efficiency, the rapidity and p T range of the π 0 s are limited to −11.0 < y lab < −8.9 and p T < 0.6 GeV, respectively. The reconstructed invariant mass of the reference MC simulations peaks at 134.8 ± 0.2 MeV and reproduces well the measured data which peaks at 134.7 ± 0.1 MeV, reproducing the π 0 mass. The uncertainties given for the mass peaks are statistical only.
Standard reconstruction algorithms used for this analysis are described in Ref. [14,27] and the π 0 event selection criteria that are applied prior to the reconstruction of the π 0 kinematics are summarized in Tab. I. Systematic uncertainties are discussed in Sec. VI.

B. Background subtraction
Background contamination of the π 0 events from hadronic events and the causal coincidence of two pho-tons not originated from the decay of a single π 0 are taken into account by subtracting the relevant contribution using a sideband method [27]. Figure 1 shows the reconstructed two-photon invariant mass distribution of the experimental data in the rapidity range −9.4 < y lab < −9.2. The sharp peak around 135 MeV is owing to π 0 events. The solid curve indicates the best fit of a composite physics model to the invariant mass distribution of the data; an asymmetric Gaussian distribution for the signal component and a third order Chebyshev polynomial function for the background component (dashed curve). The signal window is defined as the invariant mass region within the two solid arrows shown in Fig. 1, where the lower and upper limits are given bym − 3σ l andm + 3σ u , respectively.m denotes the expected mean, and σ l and σ u denote 1 σ deviations for lower and upper side of the signal component, respectively. The signal-rich rapidity-p T distributions are obtained by the events contained inside of the signal window. The remaining contribution of background events in the signal window is eliminated using the rapidity-p T distributions obtained from the background window, constructed from the two sideband regions, [m−6σ l ,m−3σ l ] and [m + 3σ u ,m + 6σ u ], that are defined as the invariant mass regions within the dashed arrows in Fig. 1. The raw rapidity-p T distributions must be corrected for (1) the reconstruction inefficiency and the smearing caused by finite position and energy resolutions, (2) geometrical acceptance and branching ratio of π 0 decay, and (3) the loss due to multihit π 0 .
First, an iterative Bayesian method [28] is used to si-multaneously correct for both the reconstruction inefficiency and the smearing. The unfolding procedure for the data is found in paper [27]. Next, a correction for the limited aperture of the LHCf calorimeters must be applied. The correction is derived from the rapidity-p T phase space. The determination of the correction coefficients follow the same method used in the π 0 event analysis in p-p collisions at √ s = 7 TeV [27]. Finally, the loss of multihit π 0 events, briefly mentioned in Sec. V A, is corrected for using the MC event generators. A range of ratios of multihit plus single-hit to single-hit π 0 events is estimated using three different hadronic interaction models (qgsjet II-03, dpmjet 3.04, and epos 1.99) for each rapidity and p T range. The observed p T spectra are then multiplied by the average of these ratios and the contribution to the systematic uncertainty is derived from the observed variations among the interaction models. The estimated range of the flux of multihit π 0 events lies within a range 0 %-10 % of the flux of single-hit π 0 events. The single-hit π 0 spectra are then corrected to represent the inclusive π 0 production spectra.
All the procedures were verified using the reference MC simulations mentioned in Sec. IV.

VI. SYSTEMATIC UNCERTAINTIES
We follow the same approach to estimate the systematic uncertainties as in Ref. [27]. Since systematic uncertainties on particle identification, single-hit selection, and position-dependent corrections for both shower leakage and light yield of the calorimeter are independent of the beam energy and type, the systematic errors are taken directly from Ref. [27].
Other terms deriving from the energy scale, beam axis offset, and luminosity are updated consistently to the current understanding of the LHCf detector and the beam configuration in p-Pb collisions. These terms are discussed in the following subsections. Table II summarises the systematic uncertainties for the analysis.

A. Energy scale
The uncertainty on the measured photon energy was investigated in the beam test at SPS [13] and also by performing a calibration with a radioactive source [10]. The estimated uncertainty on the photon energy from these tests is valued at 3.5 %. This uncertainty is dominated by the conversion factors that relate measured charge to deposited energy [13] and in fact, the reconstructed invariant mass of two photons reproduces the π 0 rest mass within the uncertainty of 3.5 % as shown in Fig. 1.
The systematic shift of bin contents due to energy scale uncertainties is estimated using two different p T spectra in which the photon energy is artificially scaled to the two extremes (±3.5 %). The ratios of the two extreme spectra to the non-scaled spectrum are assigned as systematic shifts of bin contents for each bin.

B. Beam axis offset
The projected position of the proton beam axis on the LHCf detector (beam centre) varies from Fill to Fill owing to the beam configuration, beam transverse position and crossing angles, at IP1. The beam centre on the LHCf detector can be determined by two methods; first we use the distribution of incident particle positions as measured by the LHCf detector and second we also use the information from the beam position monitors (BPMSW) installed ±21 m from IP1 [29].
From analysis results in p-p collisions at √ s = 7 TeV, the beam centre positions obtained by the two methods applied to LHC Fills 1089-1134 were found to be consistent within 1 mm. The systematic shifts to the p T spectra are then evaluated by taking the ratio of spectra with the beam centre displaced by ±1 mm to spectra with no displacement present. The fluctuations of the beam centre position modify the p T spectra by 5 %-20 % depending on the rapidity range.

C. Luminosity
The luminosity value used for the analysis is derived from on the online information provided by the ATLAS experiment. Since there is currently no robust estimation on the luminosity error by the ATLAS experiment, we assign a conservative ±20 % to the uncertainty. A more precise estimation of the luminosity will be reported in future by the ATLAS collaboration.

A. The QCD induced transverse momentum distribution
The p T spectra obtained from the data analysed are presented in Fig. 2 have all the corrections discussed in Sec. V C applied. The inclusive π 0 production rate is given as where σ pPb inel is the inelastic cross section for p-Pb collisions at √ s N N = 5.02 TeV and Ed 3 σ pPb /dp 3 is the inclusive cross section of π 0 production. The number of inelastic p-Pb collisions, N pPb inel , used for normalizing the production rates of Fig. 2 is calculated from N pPb inel = σ pPb inel Ldt, assuming an inelastic p-Pb cross section σ pPb inel = 2.11 b. The value for σ pPb inel is derived from the inelastic p-p cross section σ pp inel and the Glauber multiple collision model [30,31]. Using the integrated luminosities shown in Sec. III, N pPb inel is 9.33×10 7 . d 2 N pPb (p T , y) is the number of π 0 s detected in the transverse momentum interval (dp T ) and the rapidity interval (dy) with all corrections applied.
In Fig. 2, the filled circles represent the data from the LHCf experiment. The error bars and shaded rectangles indicate the 1 standard deviation statistical and total systematic uncertainties respectively. The total systematic uncertainties are given by adding all uncertainty terms except the one for luminosity in quadrature. The vertical dashed lines shown for the rapidity ranges greater than −9.2 indicate the p T threshold of the LHCf detector due to the photon energy threshold and the geometrical acceptance of the detector. The contribution from UPCs is presented as open squares (normalized to 1/2 for visibility). This UPC contribution is estimated with the MC simulations introduced in Sec. IV A using the UPC cross section from [32].
To obtain the soft-QCD component, the UPC contribution must be subtracted from the measured p T spectra. This is achieved by simply subtracting, point by point, the UPC induced p T spectra (open squares in Fig. 2) from the total p T spectra (filled circles in Fig. 2). Uncertainties in the subtracted results are obtained by adding the statistical and systematic errors in quadrature. The theoretical uncertainty on the UPC estimation ±5 % derives mainly from the knowledge of the virtual photon flux given by the relativistic Pb nucleus and of the virtual photon-proton cross section [32]. The inclusive production rates of π 0 s measured by LHCf after the subtraction of the UPC component are summarized in Appendix A. Figure 3 shows the LHCf data p T spectra after subtraction of the UPC component (filled circles). The size of the error bars correspond to 68 % confidence intervals (including both statistical and systematic uncertainties). The p T spectra in p-Pb collisions at 5.02 TeV predicted by the hadronic interaction models, dpmjet 3.04 (solid line, red), qgsjet II-03 (dashed line, blue), and epos 1.99 (dotted line, magenta), are also shown in the same figure for comparison. Predictions by the three hadronic interaction models do not include the UPC component. The experimental p T spectra are corrected for detector response, event selection and geometrical acceptance efficiencies, so that the p T spectra of the interaction models can be compared directly to the experimental spectra.
In Fig. 3, among the predictions given by the hadronic interaction models tested here, dpmjet 3.04 and epos 1.99 show a good overall agreement with the LHCf measurements. However qgsjet II-03 predicts softer p T spectra than the LHCf measurements and the other two hadronic interaction models. Similar features of these hadronic interaction models are also seen in p-p collisions at √ s = 7 TeV [27]. In Fig. 3, the p T spectra in p-p collisions at √ s = 5.02 TeV are also added and will be useful for the derivation of the nuclear modification factor described later in Sec. VII C. These spectra are multiplied by a factor 5 for visibility. The derivation of the p T spectra in p-p collisions at √ s = 5.02 TeV is explained in detail in Appendix B.

B. Average transverse momentum.
The average transverse momentum, p T , can be obtained by fitting an empirical function to the p T spectra in Fig. 3 for each rapidity range. Among several proposed approaches, we first choose an exponential to parametrize the p T spectra. The exponential function derived from a thermodynamical model by [33] has the form The parameter A [GeV −2 ] is a normalization factor and T [GeV] is the temperature of the π 0 s with a given transverse momentum p T . Using Eq. (2), p T is derived as a function of T : where K α (m π 0 /T ) is the modified Bessel function. Similarly we can also assume that the π 0 p T distributions follow the Gaussian distribution where the parameter A gives the absolute normalization and the Gaussian width σ Gauss [GeV] determines the average p T of the p T spectra. Thus p T is derived as a function of σ Gauss according to where f (p T ) is the right hand side of Eq. (4). For example, the left panel in Fig. 4 shows the experimental p T spectra (filled circles) and the best fit with the exponential (Eq. (2), dashed curve, blue) and with the On the other hand, p T can be simply estimated by numerically integrating the p T spectra. With this approach p T is given by where f (p T ) is the measured spectrum given in Fig. 3. In this approach, p T is obtained over the rapidity range −9.2 > y lab > −10.0 where the p T spectra are available down to 0 GeV. The data in the rapidity range −10.0 > y lab > −11.0 is not used here, since the bin content in 0 GeV < p T < 0.05 GeV is negative due to the subtraction of UPCs. Although the interval for the numerical integration is bounded from above by p upper T , the high p T tail contribution to p T is negligibly small.
The values of p T obtained by the three methods are in good agreement within the uncertainties. The specific values of p T for this paper, p T LHCf , are defined as follows. For the rapidity range −9.2 > y lab > −10.0, the values of p T LHCf are taken from the weighted average of p T from the exponential fit, the Gaussian fit, and the numerical integration. The uncertainty related to a possible bias of the p T extraction methods is derived from the largest value among the three methods. For the other rapidity ranges to where the numerical integration is not applicable, the weighted mean and the uncertainty are obtained following the same method above but using only the exponential and the Gaussian fit. Best-fit results for the three approaches mentioned above and the values of p T LHCf are summarised in Table III. Figure 5 shows the p T LHCf and the predictions by hadronic interaction models as a function of the rapidity y lab . The average p T of the hadronic interaction models is calculated by numerical integration. dpmjet 3.04 reproduces quite well the LHCf measurements p T LHCf , while epos 1.99 is slightly softer than both dpmjet 3.04 and the LHCf measurements. qgsjet II-03 shows the smallest p T among the three models and the LHCf measurements. These tendencies are also found in Fig. 3.

C. Nuclear modification factor
Finally the nuclear modification factor R pPb is derived. This factor quantifies the p T spectra modification caused by nuclear effects in p-Pb collisions. The nuclear modi-  fication factor is defined as Ed 3 σ pPb /dp 3 Ed 3 σ pp /dp 3 , where Ed 3 σ pPb /dp 3 and Ed 3 σ pp /dp 3 are the inclusive cross sections of π 0 production in p-Pb and p-p collisions at 5.02 TeV respectively. The UPC component is already subtracted in σ pPb . The average number of binary nucleon-nucleon collisions in a p-Pb collision is denoted as N coll , which is estimated as 6.9 from the MC simulations with the Glauber model [31]. Since there is no data at √ s = 5.02 TeV for p-p collisions Ed 3 σ pp /dp 3 is derived by scaling the p T spectra taken in p-p collisions at 7 TeV and 2.76 TeV. The details of the procedure are discussed in Appendix B. Figure 6 shows the nuclear modification factors R pPb from the LHCf measurements and the predictions by hadronic interaction models dpmjet 3.04 (red solid line), qgsjet II-03 (blue dashed line), and epos 1.99 (magenta dotted line). The LHCf measurements, although with a large uncertainty which increases with p T (mainly due to systematic uncertainties in p-Pb collisions at 5.02 TeV), show a strong suppression with R pPb equal 0.1 at p T ≈ 0.1 GeV rising to 0.3 at p T ≈ 0.6 GeV. All hadronic interaction models predict small values of R pPb ≈ 0.1, and they show an overall good agreement with the LHCf measurements within the uncertainty. Clearly other analyses which are more sensitive to exclusive π 0 signals are needed, for example diffractive dissociation, to investigate the reason for this strong suppression. However the measured R pPb dependency on p T and rapidity may hint to an understanding of the break down of the π 0 production mechanism.

VIII. CONCLUSIONS
The inclusive production of neutral pions in the rapidity range −8.9 > y lab > −11.0 has been measured by the LHCf experiment in p-Pb collisions at √ s N N = 5.02 TeV at the LHC in 2013. Transverse momentum spectra of neutral pions measured by the LHCf detectors have been compared with the predictions of several hadronic interaction models. Among the hadronic interaction models tested in this paper, dpmjet 3.04 and epos 1.99 show the best overall agreement with the LHCf data in the rapidity range −8.9 > y lab > −11.0, while qgsjet II-03 shows softer p T spectra relative to the the LHCf data and the other two hadronic interaction models. These tendencies are also recognized in the comparison on the average p T distribution as a function of rapidity y lab .
The nuclear modification factor, R pPb , derived from the LHCf measurements indicates a strong suppression of the π 0 production in the nuclear target relative to those in the nucleon target. All hadronic interaction models present an overall good agreement with the LHCf measurements within the uncertainty.
As a future prospect, additional analyses which are sensitive to exclusive π 0 spectra, are needed to reach a better understanding of this strong suppression and its break down. sured in p-Pb collisions at a given collision energy to the reference p T spectra in p-p collisions at the same collision energy. In this analysis, since a measurement in p-p collisions at √ s = 5.02 TeV is not available, the reference p T spectra are made by scaling the p T spectra measured in the p-p collisions at √ s = 7 TeV and 2.76 TeV. First the average p T at 5.02 TeV is estimated by scaling the average p T obtained at 7 TeV and 2.76 TeV. Figure 7 shows the average p T in p-p collisions at √ s = 7 TeV (filled circles) and 2.76 TeV (open circles) as a function of rapidity loss ∆y ≡ y beam − y cms , where y beam is beam rapidity and y cms is the rapidity of the center-of-mass frame. For the proton beam with E = 3.5 TeV and 1.38 TeV, y beam gives 8.917 and 7.987, respectively. In the following we assume y cms is positive. According to the scaling law proposed by several authors [34][35][36] (Feynman scaling), the average p T as a function of ∆y should be independent of the center-ofmass energy in the projectile fragmentation region. Thus the average p T can be directly compared among different collision energies. The values of the average p T at 7 TeV are taken from measurements by LHCf [27] in which the associated ∆y points are modified to take into account event population for each rapidity bin. These weighted bin centers are estimated using the MC simulation by epos 1.99. The values of the average p T at 2.76 TeV are obtained by a similar analysis on the data that was taken in p-p collisions at √ s = 2.76 TeV on February 13, 2013.
These data were taken with essentially the same data acquisition configuration as at 5.02 TeV. Although the two measurements in Fig. 7 have limited overlap on the ∆y range owing to the smaller collision energy at 2.76 TeV, the p T spectra at 7 TeV and 2.76 TeV follow mostly a common line. A linear function fit is then made to these measurements. The solid line and shaded area in Fig. 7 show the best-fit linear function and the 1 standard deviation uncertainty obtained by a chi-square fit to the data points p T best-fit (∆y) = 216.3 + 116.0∆y.
The minimum chi-square value is 11.1 with a number of degrees of freedom equal to 9. With the fitted result in Eq. 8 and y beam for the proton beam with E = 2.51 TeV, 8.585, the average p T at a given rapidity y cms at 5.02 TeV can be evaluated as p T (y cms )| 5.02 TeV = 216.3 + 116.0(8.585 − y cms ), (9) where we assume the proton beam travels to the positive rapidity direction. Note that the rapidity range of the reference p T spectra at 5.02 TeV is enclosed by the data points taken at 7 TeV and 2.76 TeV.
The absolute normalization scaling among three collisions is then estimated for 2.76 TeV, 5.02 TeV, and 7 TeV energies. Since the systematic uncertainty of the LHCf measurements on the luminosity is ±20 % and ±6.1 % at 2.76 TeV and 7 TeV respectively, the predictions by MC simulations are used instead for the estimation. According to dpmjet 3.04, qgsjet II-03 and epos 1.99, the relative normalization at 5.02 TeV to at 7 TeV or 2.76 TeV, defined as R norm ≡ dp 3 1 σ inel E d 3 σ dp 3 √ s=5.02 TeV (10) dp 3 1 σ inel E d 3 σ dp 3 √ s=7 TeV or 2.76 TeV , is mostly unity in the rapidity and p T ranges covered by LHCf. Therefore we apply the measured absolute normalization at 7 TeV to the reference p T spectra at 5.02 TeV without scaling. The uncertainty on the normalization is taken from the luminosity error ±6.1 %.
Accordingly the average p T and normalization of p T spectra at 5.02 TeV can be scaled from 7 TeV and 2.76 TeV. With these two values, the p T spectra in pp collisions at 5.02 TeV can be effectively derived. In the analysis of this paper, the expected p T spectra are presumed to follow a Gaussian distribution with the width of the distribution σ Gauss equal to 2 p T / √ π. The expected p T spectra in p-p collisions take into account the rapidity shift −0.465 explained in Sec. III for consistently with the asymmetric beam energies in p-Pb collisions.