Measurements of longitudinal and transverse momentum distributions for neutral pions in the forward-rapidity region with the LHCf detector

The differential cross sections for inclusive neutral pions as a function of transverse and longitudinal momentum in the very forward rapidity region have been measured at the Large Hadron Collider (LHC) with the Large Hadron Collider forward detector (LHCf) in proton-proton collisions at $\sqrt{s}=$ 2.76 and 7 TeV and in proton-lead collisions at nucleon-nucleon center-of-mass energies of $\sqrt{s_\text{NN}}=$ 5.02 TeV. Such differential cross sections in proton-proton collisions are compatible with the hypotheses of limiting fragmentation and Feynman scaling. Comparing proton-proton with proton-lead collisions, we find a sizable suppression of the production of neutral pions in the differential cross sections after subtraction of ultra-peripheral proton-lead collisions. This suppression corresponds to the nuclear modification factor value of about 0.1-0.3. The experimental measurements presented in this paper provide a benchmark for the hadronic interaction Monte Carlo simulation codes that are used for the simulation of cosmic ray air showers.

sponds to the nuclear modification factor value of about 0.1-0.3. The experimental measurements presented in this paper provide a benchmark for the hadronic interaction Monte Carlo simulation codes that are used for the simulation of cosmic ray air showers.

I. INTRODUCTION
Observations of high-energy cosmic rays with energy above 10 14 eV provide key information for a yet unestablished origin(s) and acceleration mechanism(s) for cosmic rays. The compilation of current observations reveals kinks in the energy spectrum that agree with the turning points in the mass composition [1] at ∼ 3 × 10 15 eV (the so called 'knee') and provide a consistent description of the transition from Galactic to extragalactic cosmic rays at ∼ 5 × 10 18 eV (the so called 'ankle'). In particular, a cut-off feature of ultrahighenergy cosmic rays (UHECRs) at ∼ 5 × 10 19 eV is supposed to existence of Greisen-Zatsepin-Kuzmin [2,3] cutoff, while the source and propagation of the UHECRs is still a mystery [4]. In order to grasp the experimental signature of the source of UHECRs, and to understand consistent picture of transition from galactic component around 10 14 eV, many extensive air-shower experiments, including on-going UHECR observatories (i.e. Auger [5] and Telescope Array [6]) have collected the data on energy spectrum, mass composition, and arrival direction of UHECRs high energy cosmic rays over the past few decades [7][8][9].
It is important to note that critical parts of the analysis still depend on Monte Carlo (MC) simulations of air shower development that are sensitive to the choice of hadronic interaction models. Therefore different hadronic interaction models, which simultaneously predict the soft and hard QCD interactions, provide different viewpoints even using exactly the same data compilation [1,10]. Currently the lack of knowledge about forward particle production in hadronic collisions at high energy hinders the interpretation of observations of highenergy cosmic rays [10].
Here it should be remarked that the LHC at CERN has so far reached 13 TeV centre-of-mass energy in protonproton (p + p) collisions. This energy corresponds to the cosmic ray energy 9.0 × 10 16 eV in the target rest frame which is well above the first turning point in the mass composition of primary cosmic rays from proton dominated to light nuclei dominated, namely the knee at approximately 3×10 15 eV [11]. The data provided by the LHC in the forward region, defined as the fragmentation region of a projectile particle, should thus provide a useful benchmark for the MC simulation codes that are used for the simulation of air showers.
The energy in the laboratory frame converted from the collision energy in p + p collisions at √ s = 7 TeV (E lab = 2.6 × 10 16 eV) is two orders of magnitude lower than the ankle region where a transition from galactic to extragalactic cosmic rays may occur. However, extrapolation from the LHC energy range to a higher energy range can be achieved by using a scaling law in the forward rapidity region. One possibility for such a scaling law is the hypothesis of limiting fragmentation [12][13][14], which specifies that the secondary particles will approach a limiting distribution of rapidity in the rest frame of the target hadron. In this case the fragmentation of a colliding hadron would occur independently of the centerof-mass energy and then the differential cross sections as a function of rapidity (hereafter rapidity distributions) in the fragmentation region, namely the forward rapidity region, would form a limiting distribution.
Understanding particle production in nucleon-nucleus or nucleus-nucleus interactions is also of importance for ultrahigh-energy cosmic ray interactions, where parton density in nuclei is expected to be enhanced by ∝ A 1/3 . The presence of a high gluon density in the nucleus is known to greatly modify the absolute yield and the momentum distribution of the particles that are produced [15].
The LHCf experiment [16] is designed to measure the hadronic production cross sections of neutral particles at very forward angles in p + p and proton-lead (p + Pb) collisions. The LHCf experiment also provides a unique opportunity to investigate all the effects mentioned in the previous paragraph, namely, the limiting fragmentation, the Feynman scaling [17], and the high parton density in nuclear target. In a previous publication [18] we presented the π 0 production cross sections as a function of the transverse momentum (hereafter p T distributions) in p + p collisions at √ s = 7 TeV. However tests of the limiting fragmentation and the Feynman scaling predictions were not performed. Conversely, in the analysis of this paper, the comparison of the LHCf data taken in p + p collisions at √ s = 2.76 and 7 TeV respectively makes it possible to perform these tests. In addition the analysis presented in this paper has updates that lead to a deeper understanding of forward π 0 production compared to our previous publications [18,19]: the upper range for p T analysis is extended to 1.0 GeV and differential cross sections as a function of longitudinal momentum (hereafter p z distributions) as well as p T distributions are presented.
The paper is organized as follows. In Sec. II, the LHCf detectors are described. Sections III and IV summarize the conditions for taking data and the MC simulation methodology, respectively. In Sec. V, the analysis framework and the factors that contribute to the systematic uncertainty of the results are explained. In Sec. VI the analysis results are presented and compared with the predictions of several hadronic interaction models. In Sec. VII the analysis results for p+p and p+Pb collisions are described. Finally, concluding remarks are found in Sec. VIII.

II. THE LHCF DETECTOR
Two independent detectors called LHCf Arm1 and LHCf Arm2 were assembled to study p+p and p+Pb collisions at the LHC [20]. In p + p collisions at √ s = 7 TeV, both LHCf Arm1 and LHCf Arm2 detectors were operated to measure the neutral secondary particles emitted into the positive and negative large rapidity regions, respectively. In p+p collisions at √ s = 2.76 TeV and p+Pb collisions at √ s NN = 5.02 TeV, only the LHCf Arm2 detector was used to measure the neutral secondary particles emitted into the negative rapidity region (the proton remnant side in p + Pb collisions). Here the rapidity y is defined as y = tanh −1 (p z /E) [21]. The LHCf detectors each consist of 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 20×20 mm 2 and 40×40 mm 2 for Arm1, and 25×25 mm 2 and 32×32 mm 2 for Arm2. The smaller and larger calorimeters are hereafter called the Small Calorimeter and the Large Calorimeter, respectively. Four X-Y layers of position-sensitive detectors are interleaved with the layers of tungsten and scintillator in order to provide the transverse profiles of the showers. Scintillating fiber (SciFi) belts [22] are used for Arm1 and silicon microstrip sensors [23] are used for Arm2. Readout pitches are 1 mm and 0.16 mm for Arm1 and Arm2, respectively. The Front Counters, additional components of the LHCf detectors, are simple thin plastic scintillators (80×80 mm 2 ) and are installed in front of the LHCf calorimeters. They act as monitors for beam-beam collision rates with a higher detection efficiency than the LHCf calorimeters.
The LHCf detectors were installed in the instrumentation slots of the target neutral absorbers (TANs) [24] located ±140 m from the ATLAS interaction point (IP1) in the direction of the LHCb interaction point for Arm1 and in the direction of the ALICE interaction point for Arm2, and at a zero-degree collision angle. The trajectories of charged particles produced at IP1 and directed towards the TANs are deflected by the inner beam separation dipole magnets D1 before reaching the TANs themselves. Consequently, only neutral particles produced at IP1 enter the LHCf detectors. The vertical positions of the LHCf detectors in the TANs are manipulated so that the LHCf detectors cover the pseudorapidity range from 8.4 to infinity for a beam crossing half angle of 145 µrad.
The Small Calorimeter effectively covers the zero-degree collision angle. Following p + Pb collision operation, the LHCf detectors were removed from the TAN instrumentation slots in April, 2013 in order to protect them from radiation damage when the LHC is operated at high luminosity.
LHCf triggers are generated at three levels [25]. The first level trigger is generated from beam pickup signals when a bunch passes IP1. A shower trigger is then generated when signals from any successive three scintillation layers in any calorimeter exceeded a predefined threshold. The shower trigger threshold is chosen to detect photons greater than 100 GeV with an efficiency of > 99 %. A second level trigger is generated when a shower trigger has occurred and the data acquisition system is activated. The highest level trigger, or third level trigger, is generated when a specified combination of shower triggers, front counter triggers and data acquisition trigger has occurred. The live time efficiency of the data acquisition systems is defined as the ratio of the number of second level triggers to the number of shower triggers. The efficiency depends on the luminosity during the data taking period and is always less than unity due to pileup. The final results shown are corrected for the live time efficiency.
More details on the scientific goals of the experiment are given in Ref. [16]. The performance of the LHCf detectors has been studied in previous reports [25,26].

III. EXPERIMENTAL DATA TAKING CONDITIONS
The experimental data used for the analysis in this paper were obtained at three different collision energies and colliding particle configurations. Data taking conditions are explained in the subsections below, ordered according to the dates of the operation periods with the earliest first.
A. p + p collisions at √ s = 7 TeV The data in p + p collisions at √ s = 7 TeV with a zerodegree beam crossing angle were obtained from May 15 to 22, 2010(LHC Fills 1104, 1107, 1112, and 1117. The events that were recorded during a luminosity optimization scan and a calibration run were removed from the data sets for this analysis. The integrated luminosities for the data analysis reported in this paper were derived from the counting rate of the Front Counters [27] and were 2.67 nb −1 (Arm1) and 2.10 nb −1 (Arm2) after taking the live time efficiencies into account.
Pileup interactions in the same bunch crossing may increase the multi-hit events that have more than one shower event in a single calorimeter, leading to a potential bias in the momentum distributions of π 0 s. The contamination of multi-hit events due to pileup interactions is estimated to be only 0.2 % and therefore produces a negligible effect [18]. Detailed discussions of background events from collisions between the beam and residual gas molecules in the beam tube can be found in a previous report [25].
The data in p + Pb collisions were obtained at √ s NN = 5.02 TeV with 145 µrad beam crossing half angle and with only the Arm2 detector recording data on the proton remnant side. The beam energies were 4 TeV for protons and 1.58 TeV per nucleon for Pb nuclei. Because of the asymmetric beam energies where the proton beam travels at θ = π and the Pb beam at θ = 0, the nucleon-nucleon center-of-mass in p + Pb collisions is shifted to rapidity −0.465 where A and Z are the mass and atomic numbers, respectively [28]). Data used in this analysis were taken in two different fills; during LHC Fill 3478 on January 21, 2013 and during LHC Fill 3481 on January 21 and 22. The integrated luminosity of the data was 0.63 nb −1 after correcting for the live time efficiencies of the data acquisition systems [29]. The trigger scheme was essentially identical to that used in p + p collisions at √ s = 7 TeV. The bunch spacing in p + Pb collisions (200 ns), which was smaller than the gate width for analog to digital conversion in the LHCf data acquisition system (500 ns) created the possibility of integrating two or at most three signal pulses from the pileup of successive p + Pb collisions. The probability for this to occur was estimated from the timing distribution for shower triggers and was less than 5 %. Contamination by successive collisions is not corrected for in this study, while it is considered in the beam-related systematic uncertainty. The contamination of multi-hit events due to pileup interactions is negligible (0.4 %).
It should be remarked that beam divergence causes a smeared beam spot at the TAN, leading to a bias in the measured momentum distributions. The effect of a non-zero beam spot size at the TAN was evaluated with MC simulations (see Ref. [19]). This effect is taken into account in the final results reported for the p T and p z distributions.
The data in p + p collisions at √ s = 2.76 TeV were obtained with a 145 µrad beam crossing half angle and beam energy 1.38 TeV for each proton. Data used in this analysis were taken during LHC Fill 3563 on February 13, 2013. The integrated luminosity for this data was 2.36 nb −1 after correcting for the live time efficiencies of the data acquisition system [30]. The trigger scheme, trigger efficiency, and contamination of multi-hit events were mostly the same as the p + Pb collision data at √ s NN = 5.02 TeV. The effects of beam divergence were dealt with in the same way as was described for p + Pb collisions at √ s NN = 5.02 TeV (Sec. III B).

IV. MONTE CARLO SIMULATIONS METHODOLOGY
MC simulations have been performed in two steps: (I) Event generation in p + p and p + Pb collisions at IP1 (Sec. IV A) and (II) particle transport from IP1 to the LHCf detectors and consequent simulation of the response of the LHCf detectors (Sec. IV B).
MC simulation events are generated following steps (I) and (II) and are used for the validation of reconstruction algorithms, determination of cut criteria, and determination of the response matrix for momentum distributions unfolding. Conversely, MC simulations that are used only for comparison with the measurement results in Sec. VI are limited to step (I) only, since the final p T and p z distributions in Sec. VI are already corrected for detector response and eventual reconstruction bias. The statistical uncertainties of the MC simulations used in this paper are negligibly small compared to the statistical uncertainties of the LHCf data.

A. Collision event modeling
Collision event modeling of p + p hadronic interactions at √ s = 2.76 and 7 TeV are simulated and the resulting fluxes of secondary particles are generated with several event generators: dpmjet 3.06 [31], qgsjet II-04 [32], sibyll 2.1 [33], epos lhc [34] and pythia 8.185 [35,36]. Hereafter the version number for these event generators is omitted for simplicity, unless otherwise noted.
In the analysis of this paper, we use the integrated interface crmc 1.5.3 [37] for executing the first four event generators, whereas the fifth event generator, pythia, serves as its own front end for the generation of protonproton hadronic interaction events.
Events in p + Pb collisions are divided into two categories according to the value of the impact parameter: (1) general hadronic interactions and (2) ultraperipheral collisions (UPCs). Category (1) occurs when the impact parameter between p and Pb is smaller than the sum of their radii. These inelastic p + Pb interactions at √ s NN = 5.02 TeV are simulated using the hadronic interaction models dpmjet, qgsjet, and epos with the crmc interface. sibyll was not used because it only supports nuclei lighter than Fe. pythia also does not support heavy ion collisions and thus was also not used for p + Pb collisions. Category (2) p + Pb UPCs occur when the impact parameter is larger than the sum of p and Pb radii. The UPC events are simulated by the combination of starlight [38] for the virtual photon flux, sophia 2.1 [39] for low-energy photon-proton interactions, and either dpmjet 3.05 [31] or pythia 6.428 [35] for highenergy photon-proton interactions. The UPC simulation distributions used in this analysis are taken from the average of two UPC simulations; one using dpmjet 3.05 and the second using pythia 6.428 for the highenergy photon-proton interaction. Differences between these two UPC simulations are taken into account as a systematic uncertainty in the UPC simulation. See Ref. [40] for more details.
In both p + p and p + Pb collisions, the MC events used for the determination of the response matrix for unfolding the momentum distributions (Sec. V B) are simulated by pythia at the requisite beam energies. A single π 0 with energy larger than 100 GeV and possible associated background particles are selected from the secondary particles produced. There is no significant dependence of the unfolding performance on the choice of event generator for the MC simulation events that are used for the response matrix. This was verified by repeating event simulations with three of the event generators; dpmjet, pythia and epos.
In all of the MC simulations, the π 0 s from short-lived particles that decay within 1 m of IP1, e.g. η, ρ, ω, etc. ( 10 % for each relative to all π 0 s), are accounted for consistently in the treatment of LHCf data. The 145 µrad beam crossing half angle is also taken into account for p + p collisions at √ s = 2.76 TeV and for p + Pb collisions at √ s NN = 5.02 TeV. B. Simulation of particle transport from IP1 to the LHCf detector and of the detector response Transport of secondary particles inside the beam pipe from IP1 to the TAN, the electromagnetic and hadronic 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 [41].
Secondary particles produced by the interaction between IP1 collision products and the beam pipe are also taken into account in this step. The secondary particles from beam pipe interaction events generally have energy well below 100 GeV and thus provide no bias to the momentum distributions of collision events that focus only on energies above 100 GeV. The survey data for detector position and random fluctuations due to electrical noise are also taken into account in this step. See Ref. [18] for more details.

V. ANALYSIS FRAMEWORK
A. π 0 event reconstruction and selection The standard reconstruction algorithms consist of four steps: hit position reconstruction, energy reconstruction, particle identification, and π 0 event selection.

Position reconstruction
Hit position reconstruction starts with a search for multi-hit and single hit events. A multi-hit event is defined to have more than one photon registered in a single calorimeter. A single-hit event is defined to have a single hit in each of the two calorimeters in a given detector, Arm1 or Arm2.
Therefore multi-hit event candidates should have two or more distinct peaks in the lateral-shower-impactdistribution of a given calorimeter and are then identified using the TSpectrum algorithm [42] implemented in root [43]. TSpectrum provided the basic functionality for peak-finding in a spectrum with a continuous background and statistical fluctuations.
The MC simulation estimated efficiencies for identifying multi-hit events are larger than 70 % and 90 % for Arm1 and Arm2, respectively [25]. Given the list of shower peak position candidates that have been obtained above, the lateral distributions are fit to a Lorenzian function [44] to obtain more precise estimates of the shower peak positions, heights, and widths. In the case of multi-hit events, two peaks are fit using superimposed Lorenzian functions. Multi-hit events with three or more peaks are rejected from the analysis. Conversely, singlehit events, not having two or more identifiable peaks in a single calorimeter but having a single hit in each calorimeter are correctly selected with an efficiency better than 98 % for true single-photon events with energy greater than 100 GeV for both Arm1 and Arm2.

Energy reconstruction
The photon energy is reconstructed using the measured energy deposited in the LHCf calorimeters. The charge information in each scintillation layer is first converted to a deposited energy by using the calibration coefficients obtained from the SPS electron test beam data taken below 200 GeV [26]. The sum of the energy deposited in the 2nd to 13th scintillation layers is then converted to the primary photon energy using an empirical function. The coefficients of the function are determined from the response of the calorimeters to single photons using MC simulations. Corrections for shower leakage effects and the light-yield collection efficiency of the scintillation layers are carried out during the energy reconstruction process [20]. In the case of multi-hit events, the reconstructed energy based on the measured energy deposited is split into two energies, primary and secondary. Fractions of the energy for the primary and secondary hits are determined according to the peak height and width of the corresponding distinct peaks in the lateral-showerimpact-distribution.

Particle identification
Particle identification (PID) is applied in order to efficiently select pure electromagnetic showers and to reduce hadron (predominantly neutron) contamination. PID in the study of this paper depends only on the parameter L 90% . L 90% is defined as the longitudinal distance, in units of radiation length (X 0 ), measured from the 1st tungsten layer of the calorimeter to the position where the energy deposition integral reaches 90 % of the total shower energy deposition. Events with an electromagnetic shower generally have a L 90% value smaller than 20 X 0 , while events with a hadronic shower generally have L 90% larger than 20 X 0 . The threshold L 90% value as a function of the photon energy is defined in order to keep the π 0 selection efficiency at 90 % over the entire energy range of the individual photons. PID criteria are determined by MC simulations for each calorimeter.

π 0 event selection
The π 0 are then identified by their decay into two photons, leading to the distinct peak in the invariant mass distribution around the π 0 rest mass. The invariant mass of the two photons is calculated using the reconstructed photon energies and incident positions. The π 0 events used in the analysis of this paper are classified into two categories: Type-I π 0 and Type-II π 0 events. A Type-I event is defined as having a single photon in each of the two calorimeters of Arm1 or Arm2 (the left panel of Fig. 1). A Type-II event is defined as having two photons in the same calorimeter (the right panel of Fig. 1). Note that Type-II events were not used in the previous analyses [18,19], and thus are taken into account for the first time in this paper. As detailed in Sec. V B, the phase spaces covered by Type-I and Type-II events are complementary. In particular, the inclusion of Type-II events extends the p T upper limit for analysis from 0.6 GeV in the previous analyses to 1.0 GeV. Figure 2 shows the reconstructed two-photon invariant mass (M γγ ) distributions of LHCf data in the rapidity range 8.8 < y < 10.8. The left and right panels of Fig. 2 show the distributions for Type-II events in the Arm2 small calorimeter and Arm2 large calorimeter respectively. The sharp peaks around 135 MeV are due to π 0 events. The distributions in Fig. 2 are based only on data from p + p collisions at √ s = 7 TeV during LHC Fill 1104. Similar invariant mass distributions are obtained from other fills and from Arm1. Kinematic quantities of the π 0 s (four-momenta, p T , p z and rapidity) are reconstructed by using the photon energies and incident positions measured by the LHCf calorimeters, and are used for producing the p T and p z distributions. The projected position of the proton beam axis on the LHCf detector (beam center) is used in order to derive the correct p T and p z values of each event. The beam center position is obtained from the LHCf position-sensitive detectors of   The π 0 event selection criteria that are applied prior to the reconstruction of the π 0 kinematics are summarized in Table. I. Type-I events accompanied by at least one additional background particle in one of the two calorimeters (usually a photon or a neutron) and not originating in a π 0 decay are denoted as multi-hit π 0 events and are rejected as background events. Similarly, Type-II events accompanied by at least one additional background particle in the calorimeter used for π 0 identification are rejected. Figure 3 shows diagrams of all types of multi-hit events that are rejected. Panels (a) and (b) show the multi-hit Type-I π 0 events and panels (c) and (d) show the multi-hit Type-II π 0 events. Red and green arrows indicate a background particle not originating in a π 0 decay and two photons originating in a π 0 decay, respectively. The final inclusive production rates reported in this paper are corrected for these cut efficiencies and will be discussed in Sec. V B. show the multi-hit Type-I π 0 events and panels (c) and (d) show the multi-hit Type-II π 0 events. Red and green arrows indicate a background particle not originating in a π 0 decay and two photons originating in a π 0 decay, respectively. The raw p T and p z distributions of π 0 s are corrected for: (1) contamination by background events, (2) reconstruction inefficiency and the smearing caused by finite position and energy resolutions, (3) geometrical acceptance and the branching ratio of π 0 decay, and (4) the efficiency of the multi-hit π 0 cut. We now discuss each of these corrections in some detail.

Background contamination
First, the background contamination of the π 0 events from hadronic events, and from the coincidence of two photons not originating from the decay of a single π 0 are estimated using a sideband method [18]. As shown in Fig. 2 for instance, the reconstructed two-photon invariant mass distributions of LHCf data are fit to a composite physics model (solid blue curve). The model consists of an asymmetric Gaussian distribution for the π 0 signal component and a third order Chebyshev polynomial function for the background component. The fit is performed over the two photon invariant mass range 0.08 < M γγ < 0.18 GeV. The π 0 signal window is defined by the two dashed vertical lines in Fig. 2 that are placed ±3σ from the mean value. Here the mean value and the standard deviation are obtained from the best-fit asymmetric Gaussian distribution. The background window is defined as the region within ±6σ distance from the peak value and excluding the π 0 signal window. The fraction of the background component included in the π 0 signal window can be estimated using the ratio of the integral of the best-fit third order Chebyshev function over the π 0 signal window divided by the integral over the π 0 signal and background windows. The width of the asymmetric Gaussian function comes from the detector response, predominantly from shower leakage near the edges of the calorimeters. The reconstructed energy is corrected for shower leakage.

Reconstruction inefficiency and smearing in position and energy resolution
Second, a spectrum unfolding is performed to simultaneously correct for both the reconstruction inefficiency and the smearing caused by finite position and energy resolution. The contamination by background events that has been estimated by the sideband method is taken into account in the unfolding process. We follow basically same unfolding procedure as in the previous analyses [18,19], although the unfolding algorithm is based on a fully Bayesian unfolding method [45] instead of an iterative Bayesian unfolding method [46]. The calculation of the "a posteriori" probability in multi-dimensional space (the measured spectrum multiplied by the true spectrum) is achieved using a Markov Chain Monte Carlo simulation [47]. The convergence of the Markov Chain Monte Carlo simulation is ensured by the Gelman-Rubin test [48]. Production of the MC events used for the calculation of the response matrix for the unfolding is explained in Sec. IV A.

Geometric acceptance and branching ratio corrections
Thirdly, the limiting aperture of the LHCf calorimeters is estimated by using MC simulations. The procedure for performing MC simulations is given in Ref. [18]. Figure 4 shows the acceptance efficiency as a function of the π 0 p z and p T . The acceptance efficiency has been obtained by taking the ratio of the p z -p T distribution of π 0 s that are within the aperture of the LHCf calorimeters divided by the distribution of all simulated π 0 s. The fiducial cuts [18] and reconstructed energy cut (both of the π 0 decay photons must have E > 100 GeV) are also applied to the accepted π 0 events. Dashed curves in Fig. 4 indicate lines of constant π 0 rapidity. The acceptance efficiencies in Fig. 4 are purely kinematic and do not depend upon a particular hadronic interaction model. The aper-ture correction is achieved by dividing, point by point, the distributions before the acceptance correction by the acceptance efficiency. The branching ratio inefficiency is due to π 0 decay into channels other than two photons. The branching ratio for π 0 decay into two photons is 98.8 % and is taken into account by increasing the π 0 acceptance efficiency by 1.2 %. 4. Loss of events due to the multi-hit π 0 cut Fourth, in order that the reported π 0 distributions represent inclusive cross sections it is necessary to correct the data for the loss of events due to the multihit cut (Sec. V A 4). The correction factor is defined are the number of expected multihit and single-hit π 0 events in the i-th bin respectively. The factors f multihit i are estimated using hadronic interaction models introduced in Sec. IV A and are in the range 1.0 < f multihit i < 1.1 over all the p T and p z bins. LHCf p T and p z distributions are then multiplied by the average of these factors for the various interaction models and their contribution to the systematic uncertainty is derived from the observed variations amongst the interaction models. Consequently, the single-hit π 0 distributions are corrected to represent inclusive π 0 production distributions. All the procedures just described have been verified using the MC simulations introduced in Sec. IV A.

C. Systematic uncertainties
Systematic uncertainties are determined by three factors: (1) possible biases in event reconstruction, (2) uncertainty of the LHC machine conditions, and (3) an interaction model dependence.

Systematic uncertainties in event reconstructions and unfolding of distributions
Uncertainties related to biases in event reconstruction are mainly due to five causes: (1) single-hit/multi-hit separation, (2) PID, (3) energy scale uncertainty, (4) position-dependent corrections for both shower leakage and the light yield of the calorimeters, and (5) unfolding of distributions. For the first four terms, we follow the same approaches to estimate the systematic uncertainties as we used in the previous study [18].
Concerning the unfolding process, the uncertainty is estimated by adding the following three components in quadrature. First, the uncertainty due to a possible dependence of the unfolding procedure on the shape of the p T or p z distributions to be unfolded is estimated from MC simulations; we estimate the variation of the ratios of the unfolded distributions to the true distributions among the three true distributions predictions by dpmjet, qgsjet, and epos. The second component is a dependence of the unfolding procedure on the event generator used in the generation of the response matrix for unfolding, which is negligible as we mentioned in Sec. IV A. Finally, the third component is the systematic uncertainty in the unfolding algorithm itself. This is evaluated by comparing two unfolded distributions, one obtained by a fully Bayesian unfolding method and the second obtained by the iterative Bayesian unfolding method. The uncertainty in the first component is 10 % over the all p T and p z bins, and the uncertainties in the other two components make no significant contribution. Thus we assign 10 % for the systematic uncertainty in the unfolding of p T and p z distributions.

Systematic uncertainties in the LHC machine conditions
The LHC machine conditions introduce systematic uncertainties in beam position and luminosity. The beam position at the LHCf detectors varies from fill to fill owing to variations of the beam transverse position and the crossing angles at IP1. The beam center positions at the LHCf detectors obtained for LHC Fills 1089 to 1134 by the LHCf position-sensitive detectors and by the beam position monitors (BPMSW) installed ±21 m from IP1 [49] are consistent with each other within ±1 mm.
The systematic shifts to the p T and p z distributions are then evaluated by taking the ratios of distributions with the beam center displaced by ±1 mm to distributions with no displacement present. The evaluated systematic shifts to the p T and p z distributions are 5-20 % depending on the p T and p z values.
The uncertainty in the luminosity depends on the collision configuration. For the data in p + p collisions at √ s = 7 TeV, the luminosity value used for the analysis is derived from the counting rate of the Front Counters. Considering the uncertainties in both the calibration of the Front Counters ±3.4 % and in the beam intensity measurement ±5.0 % during the Van der Meer scans, we estimate an uncertainty of ±6.1 % in the luminosity for p + p collisions at √ s = 7 TeV [27]. For the p + p collision data at √ s = 2.76 TeV and p + Pb collision data at √ s NN = 5.02 TeV, LHCf data were taken simultaneously with data taken by the ATLAS experiment. The luminosity values used for this data analysis were then provided by the LHCf Front Counters and also by the ATLAS collaboration. The luminosity uncertainties in p+p collisions at √ s = 2.76 TeV and in p + Pb collisions at √ s NN = 5.02 TeV are estimated to be ±3.1 % [30] and ±20 % [29], respectively. Pileup of successive p + Pb collisions due to the small bunch spacing (200 ns) relative to the data acquisition time (500 ns) amounts to < 5 % systematic uncertainty of p T and p z distributions (see Sec. III B), and may provide a slight shift of the absolute normalization for the p T and p z distributions. This effect is not corrected for in this study, but is taken into account as uncertainty related to the LHC machine condition.

Systematic uncertainties depending on the interaction models used in the MC simulations
The analysis in this paper unavoidably relies on the predictions given by MC simulations. First, we correct LHCf data for the loss of multi-hit π 0 events (Sec. V B 4).
The correction factors f multihit show a systematic uncertainty of less than 10 % among the hadronic interaction models. Second, for p + Pb collisions only, the contamination from UPC induced π 0 events in LHCf data is derived from MC simulations (Sec. IV A). The comparison of the predicted p T and p z distributions of π 0 s between two UPC MC simulations, one using dpmjet 3.05 and the other one using pythia 6.428 for highenergy photon-proton interaction, show a systematic uncertainty of roughly 3-20 %.
In summary, there are 10 systematic uncertainties. The first four (1) Single/multihit selection, (2) PID, (3) energy scale and (4) position-dependent correction are explained in Ref. [18] and we follow the same approaches as we used in Ref. [18]. The remaining six systematic uncertainties and the text containing their explanations are: (5) Unfolding uncertainty is explained and evaluated in Sec. V C 1, (6) Offset of beam axis is explained in TABLE II. Summary of the systematic uncertainties. Numerical values indicate the maximum variation of bin contents in the pT and pz distributions due to systematic uncertainties. Note that the uncertainty in contamination of successive p + Pb collisions and in UPC π 0 simulation pertain only to p + Pb collisions.
the 1st paragraph of Sec. V C 2, 5-20 % shifts in p T or p z distributions are obtained, (7) Luminosity uncertainty is explained in the 2nd paragraph of Sec. V C 2. (8) Contamination of successive p + Pb collisions is explained in the 3rd paragraph of Sec. V C 2, (This uncertainty is due to contamination and thus only a positive error is quoted.) (9) The uncertainty in multihit π 0 events ±10 %, and (10) the uncertainty in UPC ±(3-20 %) are found in Sec. V C 3. Table II summarizes the systematic uncertainties of the π 0 p T and p z distributions.

VI. ANALYSIS RESULTS
A. Results in p + p collisions at √ s = 7 TeV The inclusive production rate of neutral pions as a function of p T and p z is given by the expression [21] 1 2πp T dp T dp z .
(1) σ inel is the inelastic cross section for p + p collisions at √ s = 7 TeV. Ed 3 σ/dp 3 is the inclusive cross section for π 0 production. The number of inelastic collisions, N inel , used for the production rate normalization is calculated from N inel = σ inel Ldt, taking the inelastic cross section σ inel = 73.6 mb [18]. The uncertainty in σ inel is estimated to be ±3.0 mb by comparing the values of σ inel reported in Refs. [50][51][52][53].
Using the integrated luminosities Ldt, reported in Sec. III A, N inel is (2.67 ± 0.11) × 10 8 for Arm1 and (2.10 ± 0.09) × 10 8 for Arm2. d 2 N (p T , y) is the number of π 0 s produced within the transverse momentum interval dp T and the rapidity interval dy. Similarly d 2 N (p T , p z ) is the number of π 0 s produced within dp T and the longitudinal momentum interval dp z .
Experimental p T and p z distributions measured inde-pendently with the Arm1 and Arm2 detectors are combined following a pull method [54] and the final p T and p z distributions are then obtained by minimizing the value of the chi-square function defined by where the index i represents the p T or p z bin number running from 1 to n (the total number of p T or p z bins), and the index a indicates the type of distributions: a = 1 Arm1 Type-I events, a = 2 Arm1 Type-II events with the Large Calorimeter, a = 3 Arm2 Type-I events, a = 4 Arm2 Type-II events with the Small Calorimeter, and a = 5 Arm2 Type-II events with the Large Calorimeter. Note that Arm1 Type-II events with the Small Calorimeter are not used for this analysis since the energy reconstruction accuracy for these events is still being investigated. R obs a,i is the inclusive production rate in the ith bin of the ath distribution, which corresponds to the second and third terms Eq. (1). R comb i is the inclusive production rate in the ith bin obtained by combining all R obs a,i 's for a =1-5. σ a,i is the uncertainty of R obs a,i . The σ a,i are calculated by quadratically adding the statistical uncertainty and the systematic uncertainty in the energy scale. The energy scale uncertainty has been estimated with test beam data taken at the SPS and is uncorrelated bin-by-bin unlike the other systematic uncertainties [18]. The systematic correction S a,i modifies the number of events in the ith bin of the ath distribution: The coefficient f j a,i is the systematic shift of the ith bin content of the ath distribution due to the jth systematic uncertainty term. The systematic uncertainty consists of seven uncertainties related to the single-hit/multihit separation, the PID, the energy scale (owing to the invariant mass shift of the measured π 0 events), the positiondependent correction, the unfolding procedure, the beam center position, and the loss of multihit π 0 events. These uncertainties are assumed to be fully uncorrelated between the Arm1 and Arm2 detectors, while correlations between Type-I and Type-II events and bin-bin correlations have been accounted for. The coefficients ε j a , which should follow a Gaussian distribution, can be varied, within the constraints of the penalty term given by to achieve the minimum χ 2 value for each chi-square test.
Note that the uncertainty in the luminosity determination, ±3.1 %-±20 %, not included in Eq. (3) and Eq. (4), can cause independent shifts of all the p T and p z distributions.
The LHCf p T distributions (filled circles) are obtained from the best-fit R comb and are shown in Fig. 5. The 68 % confidence intervals incorporating the statistical and systematic uncertainties, except for the luminosity uncertainty, are indicated by the error bars. LHCf p T distributions are corrected for the influences of the detector response, event selection efficiencies and geometrical acceptance efficiencies, and thus LHCf p T distributions can be compared directly to the predicted p T distributions from hadronic interaction models. For comparison, the predictions from various hadronic interaction models are also shown in Fig. 5: dpmjet (solid red line), qgsjet (dashed blue line), sibyll (dotted green line), epos (dashed-dotted magenta line), and pythia (default parameter set, dashed-double-dotted brown line). For these hadronic interaction models, the inelastic cross section used for the production rate normalization is taken from the predefined value in each model. Figure 6 presents the ratios of the inclusive production rates predicted by the hadronic interaction models listed above to those obtained by LHCf data. Shaded areas have been taken from the statistical and systematic uncertainties. In Fig. 6, the denominator and the numerators, namely the inclusive production rate for LHCf data and for the hadronic interaction models, respectively, are properly normalized by the inelastic cross section for each, thus we do not apply any other normalization to the ratios. The inclusive production rates of π 0 s measured by LHCf and the ratios of π 0 production rate of MC simulation to data are summarized in Appendix.
In the comparisons in Fig. 5 and 6, qgsjet has good overall agreement with LHCf data, while epos produces a slightly harder distribution than the LHCf data for p T > 0.5 GeV. These two models are based on the parton-based Gribov-Regge approach [55,56] and are tuned by using the present LHC data (ALICE, ATLAS, CMS, and TOTEM) [32,34]. The prediction of sibyll agrees well with the LHCf data for 8.8 < y < 9.2 and p T < 0.4 GeV, while the absolute yield of sibyll is about half that of the LHCf data for y > 9.2. The predictions of dpmjet and pythia are compatible with LHCf data for 9.0 < y < 9.8 and p T < 0.2 GeV, while for p T > 0.2 GeV they become significantly harder than both LHCf data and the other model predictions. Generally the harder distributions appearing in sibyll, dpmjet, and pythia can be attributed to the baryon/meson production mechanism that is used by these models. For example the popcorn approach [57,58] implemented in the Lund model is known to produce hard distributions of forward mesons [59]. Indeed, by only changing the tuning parameters of the popcorn approach in dpmjet one obtains softer meson distributions and consequently p T distributions that are compatible with LHCf data. However such a crude tune may bring disagreements between the model predictions and other experimental results, e.g. forward neutron p z and p T distributions.
The LHCf p z distributions are shown in Fig. 7. The models are also shown in Fig. 7. Figure 8 presents the ratios of p z distributions predicted by the hadronic interaction models to the LHCf p z distributions. Shaded areas have been taken from the statistical and systematic uncertainties. The same conclusions for the comparisons are obtained as those found for Fig. 5 and 6. There is again an overall agreement between LHCf data and the qgsjet prediction, especially for 0.0 < p T < 0.2 GeV. The epos prediction is compatible with LHCf data for p T < 2 TeV, while showing a hard slope for p T > 2 TeV in all p T regions. The predictions by dpmjet and pythia agree with LHCf data for p T < 0.2 GeV and p z < 1.6 TeV, while showing a harder distribution for the higher p z regions. sibyll predicts a smaller production of π 0 s for p T < 0.2 GeV and becomes similar with dpmjet and pythia with increasing p T .
The inclusive production rates of π 0 s as a function of p T and p z are given by Eq. (1). Using the inelastic cross section σ inel = (62.5 ± 5.0) mb [21] and the integrated luminosities reported in Sec. III C, N inel is calculated as (1.60 ± 0.13) × 10 8 . The uncertainty on σ inel is estimated by comparing the σ inel value with the present experimental result [60]. Note that only the LHCf Arm2 detector was operated in p+p collisions at √ s = 2.76 TeV and that only Type-I events are used for the analysis since Type-II event kinematics are outside the calorimeter acceptance for √ s = 2.76 TeV. LHCf p T distributions are shown in Fig. 9. The p T distributions predictions for the hadronic interaction models are also shown in Fig. 9 for comparison.  presents the ratios of p T distributions predicted by the hadronic interaction models to the LHCf p T distributions. qgsjet provides the best agreement with LHCf data, although it is slightly softer than the LHCf data for y > 9.2. The prediction of epos shows a harder behavior than both qgsjet and LHCf data. sibyll tends to have generally a smaller π 0 yield and a harder distribution compared to qgsjet and epos, leading to the smaller and larger yields with respect to LHCf data in the p T regions below and above 0.1 GeV. dpmjet and pythia predict larger π 0 yields than both LHCf data and other models over the entire rapidity range. The same discussion on the popcorn model in the previous Section VI A can be applied to the predictions of sibyll, dpmjet, and pythia.
LHCf p z distributions are shown in Fig. 11. Figure 12 presents the ratios of p z distributions predicted by the hadronic interaction models to LHCf p z distributions. The same tendencies found in Fig. 7 are present here also qgsjet gives the best agreement for 0.0 < p T < 0.4 GeV and epos has a harder behavior especially for 0.2 < p T < 0.4 GeV. The predictions of dpmjet and pythia are significantly harder than LHCf data for p T < 0.4 GeV and show poor overall agreement with LHCf data. This can be explained by the popcorn model in a way similar to the harder p T distributions of the sibyll, dpmjet and pythia models found in Fig. 7   The inclusive π 0 production rate in p + Pb collisions is given as where σ pPb inel is the inelastic cross section, Ed 3 σ pPb /dp 3 is the inclusive cross section of π 0 production in p + Pb collisions at √ s NN = 5.02 TeV, and y lab is the rapidity in the detector reference frame. The number of inelastic p + Pb collisions, N pPb inel , used for normalizing the production rates is calculated from N pPb inel = σ pPb inel Ldt, assuming the inelastic p + Pb cross section σ pPb inel = (2.11 ± 0.11) b [61]. The value for σ pPb inel is derived from the inelastic p + p cross section σ pp inel and the Glauber multiple collision model [61,62]. The uncertainty on σ pPb inel is estimated by comparing the σ pPb inel value with other calculations and experimental results presented in Refs. [63,64]. Using the integrated luminosities described in Sec. III, N pPb inel is (9.33 ± 0.47) × 10 7 . Note that only the LHCf Arm2 detector (proton remnant side) was operated in p + Pb collisions at √ s NN = 5.02 TeV. Figure 13 shows LHCf p T distributions with both statistical and systematic errors (filled circles and error bars). The p T distributions in p + Pb collisions at √ s NN = 5.02 TeV predicted by the hadronic interaction models, dpmjet (solid red line), qgsjet (dashed blue line), and epos (dotted magenta line), are also shown in the same figure for comparison. The expected UPC contribution discussed in Sec. IV A is added to the hadronic interaction model predictions for consistency with the treatment of LHCf data, and the UPC p T distribution is shown for reference (dashed-double-dotted green line).
In Fig. 13, dpmjet shows good agreement with LHCf data at −8. in UPCs. In fact the UPC simulation reproduces such a bump. Figure 14 presents the ratios of LHCf p T distributions to the p T distributions predicted by hadronic interaction models taking the UPC contribution into account in the p T distributions.
The p z distributions are shown in Fig. 15. Figure 16 presents the ratios of LHCf p z distributions to the p z distributions predicted by the hadronic interaction models. A similar tendency to that found in p+p collisions at √ s = 7 TeV is found for LHCf data relative to model predictions. Concerning the comparison of hadronic interaction models with LHCf data, qgsjet shows a very good agreement at p T < 0.2 GeV. However at p T > 0.2 GeV, there are no models giving a consistent description of LHCf data within uncertainty over all p z bins, although epos shows a certain compatibility with LHCf data for p T > 0.4 GeV and for p z < 3 TeV. The dpmjet predictions agree with LHCf data at p T < 0.6 GeV and p z < 2 TeV, while showing a harder distribution at higher p z similar to p + p collisions. Again note the characteristic bump found in the LHCf data at p z ∼ 1.2 TeV and p T < 0.4 GeV, originating from the channel γ+p → π 0 +p via baryon resonances in UPCs.

A. Average transverse momentum
According to the scaling law proposed in Ref. [65], the average transverse momentum, denoted p T , as a function of rapidity should be independent of the center-ofmass energy in the projectile fragmentation region. Here we obtain and compare the p T distributions as func- tions of rapidity for p + p and p + Pb collisions. In the study of this paper, p T is obtained by three methods discussed below. The first two methods use analytic distributions that are fit to the LHCf data and the third uses numerical integration of the LHCf data.
The first method uses the fit of an empirical Gaussian distribution to the LHCf p T distributions for each rapidity range in Fig. 5, 9, and 13. The second method uses a Hagedorn function. Here we pay attention to the fact that soft scattering dominates the measured π 0 events for p T <∼ 1 GeV thus excluding from the analysis a powerlaw distribution that is used predominantly for hard scattering at p T >∼ 1 GeV. These methods do not necessarily require that the measured p T distribution be available down to 0.0 GeV, although the best-fit distribution may then include a systematic uncertainty in its fit [66]. Detailed descriptions of the parametrization and derivation of p T by using the best-fit Gaussian distribution can be found in Ref. [18]. In a Hagedorn function [66], the invariant cross section of identified hadrons, namely π 0 s in this paper, with a given mass m [GeV] and temperature T [GeV] can be written as where A [GeV −3 ] is a normalization factor and K 1 is the modified Bessel function. Approximately half of the π 0 measured with the LHCf detector are daughters from the decay of parent baryons and mesons and are not directly produced. Thus the measured p T distribution is no longer a thermal distribution of prompt π 0 s and so we set m as a free parameter as well as A and T in the fit of a Hagedorn function to the p T distribution. Equation (6) converges by n ≈ 5 and the computation is in fact stopped at n = 10. The p T value is calculated by using the modified Bessel functions K 5/2 and K 2 as functions of the ratio of the best-fit m and T values [66], For reference, Fig. 17 shows LHCf p T distributions (filled black circles) and the best fits of the Gaussian distributions and the Hagedorn functions. The left panel of Fig. 17 shows the results for 9.2 < y < 9.     for LHCf data is plotted after subtraction of the UPC component where the systematic uncertainty in the simulation of UPC events has been taken into account. The best-fit Gaussian distribution and the Hagedorn function reproduce the LHCf p T distributions within the total uncertainties and are also compatible with each other.
Finally, for the third method, p T is obtained by numerically integrating the p T distributions in Fig. 5, 9, and 13. The LHCf p T distributions in p + Pb collisions have already had the UPC component subtracted. In this approach, p T is calculated only in the rapidity range where the p T distributions are available down to 0.0 GeV. The high-p T tail that extends beyond the data (p T p T ) has a negligible contribution to p T . The final p T values obtained in this analysis, denoted p T LHCf , have been determined by averaging the p T values calculated with the three above described independent approaches: Gaussian, Hagedorn and numerical integration. The uncertainty of p T LHCf for each rapidity bin is assigned to fully cover the minimum and maximum p T values obtained by the three approaches. The p T LHCf values are summarized in Table. III.
In Fig. 18, p T in p + p collisions at √ s = 2.76 and 7 TeV, and in p + Pb collisions at √ s NN = 5.02 TeV are presented as a function of rapidity loss ∆y ≡ y beam − y, where y beam is the beam rapidity for each collision energy. The shift in rapidity by y beam allows a direct comparison to be made between the p T results at different collision energies. We see that for ∆y > −1.3, p T at √ s = 2.76 TeV (open red circles) has slightly smaller values than at √ s = 7 TeV (filled black circles), although the two sets of data are mostly compatible at the ±10 % level. For reference, the SppS UA7 results for p +p collisions at √ s = 630 GeV [67] (open magenta squares) show a rapid roll off of p T at low ∆y compared to LHCf data. The LHCf and UA7 results are particularly incompatible for −0.3 < ∆y < 0.3. The comparison of the LHCf data with the UA7 results indicates that p T may depend on the center-of-mass energy. However, in order to firmly confirm a center-of-mass energy dependence of p T , we need to have experimental data at a lower collision energy, e.g., √ s < 1 TeV and with a wider range of rapidity. Approved plans are underway to obtain this data by installing the LHCf detector at the RHIC ZDC position [68]. The p T values obtained from p + Pb collisions at √ s NN = 5.02 TeV (filled blue triangles) are consistent with those from p + p collisions at √ s = 7 TeV within the systematic uncertainties present. The predictions from dpmjet (thick solid red line) and qgsjet (thin solid blue line) in p+p collisions at √ s = 7 TeV and p + Pb collisions at √ s NN = 5.02 TeV have been added to Fig. 18 for reference. The predictions at √ s = 2.76 TeV are excluded in Fig. 18, since these curves mostly overlap with those at 7 TeV. LHCf data in p + p collisions at √ s = 7 TeV are close to the predictions made by dpmjet at large ∆y (small y) and become close to those made by qgsjet at small ∆y (large y). These relations between LHCf data and the model predictions are consistent with the p T distributions shown in Fig. 5 and 9. prediction obtained from qgsjet (thin dashed blue line) is smaller than LHCf data for ∆y > −0.5 and approaches the LHCf results on decreasing ∆y. This tendency was already found in Fig. 13; the prediction from qgsjet shows an overall agreement with LHCf p T distributions at y lab < −9.8.

B. Limiting fragmentation
The hypothesis of limiting fragmentation [12][13][14] asserts that the number of fragments of a colliding hadron will follow a limiting rapidity distribution in the rest frame of the target hadron. In this case the rapidity distribution of the secondary particles in the forward rapidity region would be independent of the center-of-mass energy. In this paper, a test of the limiting fragmenta-tion hypothesis is performed by using LHCf data in p + p collisions at √ s = 2.76 and 7 TeV. The normalized rapidity distribution of π 0 s, (1/σ inel )(dσ/dy), in this analysis can be obtained by using very similar methods that were used for the derivation of the average p T in Sec. VII A.
The first method uses the fit of an empirical distribution to the LHCf p T distributions in Fig. 5 and 9 in each rapidity range. As we discussed in Sec. VII A, two distributions are chosen to parametrize the p T distributions: a Gaussian distribution and a Hagedorn function. The rapidity distribution is derived by integrating the best-fit Gaussian distribution and Hagedorn function along the p T axis from 0.0 GeV to infinity.
The rapidity distribution can also be obtained by numerically integrating the p T distributions in Fig. 5 and 9. In this approach, the derivation of the (1/σ inel )(dσ/dy) value is possible only in the rapidity range where the p T distributions are available down to 0.0 GeV. Again, the final rapidity distribution is derived by averaging the rapidity distributions obtained by the above three methods. The estimated uncertainty is obtained from the minimum and maximum values for each rapidity bin. Figure 19 shows the rapidity distributions as functions of the rapidity loss ∆y (i.e., y beam −y) in p+p collisions at √ s = 2.76 TeV (open red circles) and 7 TeV (filled black circles). The rapidity distributions for both collision energies mostly appear to lie along a common curve in the rapidity range −1.8 < ∆y < −0.8. LHCf data are consistent at the ±15 % level with the hypothesis of limiting fragmentation in the very forward region.
For comparison the experimental results from the UA7 experiment [67] are also shown in Fig. 19. The extrap-olated LHCf curve at 7 TeV to higher ∆y (i.e., lower y) could be compatible with the UA7 results, at least for ∆y 0.5.
The predictions of dpmjet (thick red curve) and qgsjet (thin blue curve) at √ s = 7 TeV have been added to Fig. 19 for reference. The predictions at √ s = 2.76 TeV have been omitted, since these curves mostly overlap with those at 7 TeV since limiting fragmentation holds in dpmjet and qgsjet. The best agreement with LHCf data at √ s = 2.76 and 7 TeV is obtained by the qgsjet model. The dpmjet predictions generally give a larger π 0 yield and a harder p T distribution especially for y > 9.8 at √ s = 7 TeV and for y > 9.4 at 2.76 TeV.

C. Feynman scaling
In Ref. [17] Feynman proposed that the production cross sections of secondary particles as functions of the Feynman-x variable (defined by x F ≡ 2p z / √ s) were independent of the incident energy in the forward region. If the so-called Feynman scaling holds, the differential cross section as a function of x F (hereafter x F distribution) (x F /σ inel )(dσ/dx F ) should be independent of the center-of-mass energy for x F 0.2. Here the rapidity distribution introduced in Sec. VII B can be rewritten as where x E ≡ 2E/ √ s and dy = dp z /E are used for the second form. Considering p z ≈ E in the forward region, x E can be considered as x F and thus the right hand side of Eq. (8) becomes approximately (x F /σ inel )(dσ/dx F ). Consequently, the limiting fragmentation hypothesis that states (1/σ inel )(dσ/dy) is independent of the center-ofmass energy in each rapidity bin can be rewritten as Feynman scaling which states (x F /σ inel )(dσ/dx F ) is independent of the center-of-mass energy in each x F bin. In this paper, we test the Feynman scaling hypothesis by comparing LHCf data in p + p collisions at √ s = 2.76 and 7 TeV.
In Fig. 20, we compare the x F distributions in the p T range 0.0 < p T < 0.4 GeV. Other p T ranges are excluded from the comparison, since LHCf data at √ s = 2.76 TeV are unavailable outside this range. The x F distributions at √ s = 2.76 and 7 TeV are compatible with each other at the ±20 % level. In Fig. 21, we further compare the x F distributions for the reduced p T ranges: 0.0 < p T < 0.2 GeV and 0.  that Feynman scaling holds at the ±20 % level at these center-of-mass energies in the very forward region. Besides a test of the Feynman scaling, we find in Fig. 21 that the yield of π 0 s at √ s = 2.76 TeV relative to 7 TeV is slightly larger for 0.0 < p T < 0.2 GeV and slightly smaller for 0.2 < p T < 0.4 GeV. This tendency means that the p T distributions at √ s = 2.76 TeV are softer than those at 7 TeV, leading to the small p T values at 2.76 TeV relative to those at 7 TeV as already found in Fig. 18.

D. p T dependence of the xF distributions
In hadronic interactions at large rapidities, partons from the projectile and target hadrons generally have large and small momentum fractions respectively, since the momentum fraction that the parton itself carries relative to the parent projectile and target hadrons, i.e., the Bjorken-x variable or x Bj , is proportional to e ±y (+y for projectile and −y for target). Here we note that a parton (dominantly gluon) density, rapidly increases with decreasing x Bj when x Bj < 0.01 with the target approaching the blackbody limit where the gluon density is saturated. In the blackbody regime, the partons cannot go through the target nuclear medium without interaction and suffer transverse momenta transfers proportional to the saturation momentum scale Q s . The Q s values in the very forward region are ∼ 1 GeV in p + p collisions and ∼ 10 GeV in p + Pb collisions, although the calculation of Q s itself suffers from both theoretical and experimental uncertainties and is also dependent on the impact parameter of the colliding hadrons [15].
In the p T region below Q s , the x F distribution in the forward region can be asymptotically written [69] as where α is the leading exponent. In the blackbody regime, the x F distribution of the leading hadron is strongly suppressed and thus α increases relative to the value found for a dilute target. Conversely, α decreases with increasing p T when p T approaches or exceeds the saturation momentum scale Q s . Figure 22 shows the best-fit leading exponent α in each p T range in p + p and p + Pb collisions. The leading exponent in p + p collisions at √ s = 7 TeV (filled black circles) is α ≈ 3.7 at p T < 0.6 GeV and decreases to α ≈ 3.0 at 0.6 < p T < 1.0 GeV. The reduction of α with increasing p T can be understood as much of the target staying in the blackbody regime for p T < 0.6 GeV and then gradually escaping from the blackbody regime for p T > 0.6 GeV. The leading exponent in p + p collisions at √ s = 2.76 TeV (open red circles) is slightly smaller than that at 7 TeV though with large uncertainty. The comparison between √ s = 2.76 TeV and 7 TeV may indicate that the upper p T limit of the measurement at 2.76 TeV is near the saturation momentum at 2.76 TeV and that the suppression due to the gluon density is weaker than at 7 TeV, although the calculated Q s at 2.76 TeV is only slightly different from the Q s at 7 TeV. The leading exponents in p + Pb collisions at √ s NN = 5.02 TeV (filled blue triangles) are rather flat along the p T axis, within uncertainties that are generally larger than those in p + p collisions. This may indicate that the saturation momentum in p + Pb collisions is well above the measured p T range and also that the x F distributions in p + Pb collisions are suppressed relative to those for p + p collisions.

E. Nuclear modification factor
The effects of high gluon density in the target are inferred from the comparison of the leading exponent α between in p + p and p + Pb collisions (see the preceding Sec. VII D). Here we introduce the nuclear modification factor that quantifies the p T spectra modification caused by nuclear effects in p + Pb collisions with respect to p + p collisions. The nuclear modification factor R pPb is defined as where Ed 3 σ pPb /dp 3 and Ed 3 σ pp /dp 3 are the inclusive cross sections of π 0 production in p + Pb collisions at √ s NN = 5.02 TeV and in p+p collisions at √ s = 5.02 TeV, respectively. These cross sections are obtained from Eq. (5) and Eq. (1), with the subtraction of the expected UPC contribution applied to the cross section for p + Pb collisions. The uncertainty in the inelastic cross section σ pPb inel is estimated to be ±5 % [19]. The average number of binary nucleon-nucleon collisions in a p + Pb collision, N coll = 6.9, is obtained from MC simulations using the Glauber model [61]. The uncertainty of σ pp inel / N coll is estimated by varying the parameters in the calculation with the Glauber model and is of the order of ±3.5 % [19]. Finally the quadratic sum of the uncertainties in σ pPb inel and σ pp inel / N coll is added to R pPb . Since there is no LHCf data for p + p collisions at exactly √ s = 5.02 TeV, Ed 3 σ pp /dp 3 is derived by scaling the p T distributions taken in p + p collisions to other collision energies. The derivation follows three steps. First, the p T at √ s = 5.02 TeV is estimated by interpolating the measured p T values at 7 TeV. The uncertainty of the interpolated p T values is estimated to be ±10 % according to the differences between the measured p T values at √ s = 2.76 and 7 TeV for −1.7 < ∆y < −0.8 (see Fig. 18). Second, the absolute normalization of the p T distribution value in each rapidity range for √ s = 5.02 TeV, i.e., (1/σ inel )(dσ/dy), is determined by interpolating the rapidity distribution at 7 TeV (see Fig. 19). The uncertainty of the absolute normalization is estimated to be ±15 % according to the discussion in Sec. VII B and is taken into account in the interpolated normalization. Finally, the p T distributions at √ s = 5.02 TeV are produced by assuming that the p T distribution follows either a Gaussian distribution or a Hagedorn function and by using the p T values obtained in the first step and the normalization obtained in the second step. The difference of the p T distribution produced using a Gaussian distribution or a Hagedorn function gives the systematic uncertainty. Note that the rapidity shift −0.465 explained in Sec. III B is also taken into account for the p T distribution in p + p collisions at √ s = 5.02 TeV. Figure 23 shows the nuclear modification factors R pPb obtained from LHCf data and the predictions from the hadronic interaction models, dpmjet (solid red curve), qgsjet (dashed blue curve), and epos (dotted magenta curve). LHCf data show a strong suppression with R pPb equal to 0.3 at y lab ∼ −8.8 down to < 0.1 at y lab ∼ −10.8, although a large uncertainty is present due to systematic uncertainties in the estimation of the p T values in p + p collisions at √ s = 5.02 TeV. All hadronic interaction models, which employ different approaches for the nuclear effects, predict small values of R pPb 0.15. Within the uncertainties the hadronic interaction models show an overall good agreement with R pPb estimated from LHCf data.

VIII. CONCLUSIONS
The inclusive production of π 0 s was measured with the LHCf detector in p + p collisions at √ s = 2.76 and 7 TeV and in p + Pb collisions at √ s NN = 5.02 TeV. In p + p collisions at √ s = 7 TeV, differential cross sections as a function of p T and p z for the π 0 s were measured by two independent LHCf detectors, Arm1 and Arm2, with consistent results. Conversely, only the LHCf Arm2 detector was used in p + p collisions at √ s = 2.76 TeV and in p + Pb collisions.
In p + p collisions, qgsjet II-04 shows an overall agreement with LHCf data, while epos lhc distributions have a slightly harder behavior than LHCf data for p T > 0.5 GeV. dpmjet 3.06 and pythia 8.185 have in general shown a harder momentum distributions and a poor agreement with LHCf data. In p + Pb collisions, dpmjet 3.06 showed good agreement with LHCf data for −8.8 > y lab > −10.0 and p T < 0.3 GeV, while qgsjet II-04 and epos lhc did better reproducing the LHCf data for p T > 0.4 GeV than dpmjet 3.06.
The average values of p T , denoted p T , at y > 8.8 in p + p collisions and at y lab > 8.8 in p + Pb collisions were calculated using the LHCf p T distributions. The p T values obtained have been shown to be independent of the center-of-mass energy at the 10 % level. Tests of limiting fragmentation and Feynman scaling hypotheses using LHCf data in p + p collisions show that both hypotheses hold in the forward region at the 15 %-20 % level. The leading exponent α and the nuclear modification factor R pPb derived from LHCf data indicate a strong suppression of π 0 production from the nuclear target relative to that from the nucleon target. Within the uncertainties all of the hadronic interaction models presented gave an overall good agreement with R pPb estimated by LHCf data. According to the analysis in this paper, we expect that the number of particles leading to an electromagnetic component in air showers would follow the limiting rapidity distribution and Feynman scaling hypotheses. Combining the results for forward π 0 s in this paper with the recent results for forward neutrons in Ref. [70] strongly constrain models for air shower production at the TeV scale.
As a future prospect, additional analyses using correlations between forward π 0 s and other particles (e.g., two-particle angular correlations) are needed to reach a better understanding of the forward meson production mechanism and the strong suppression of π 0 production in p + Pb collisions compared to p + p collisions. The AT-LAS and LHCf Collaborations have taken p + p data at √ s = 13 TeV and p + Pb data at √ s NN = 5.02 TeV with common triggers. This data could provide the possibility for performing analyses of two particle correlations.

ACKNOWLEDGMENTS
The LHCf Collaboration acknowledges CERN staff and the ATLAS Collaboration for their essential contributions to the successful operation of LHCf. We thank S. Ostapchenko and T. Pierog for numerous discussions and for confirming the results of the MC simulations. We are grateful to C. Baus, T. Pierog and R. Ulrich for providing the crmc program codes and useful comments. This work has been partly supported by a Grant-in-Aid for Scientific research by MEXT of Japan, a Grant-in-Aid for a JSPS Postdoctoral Fellowship for Research Abroad, a Grant-in-Aid for Nagoya University GCOE "QFPU" from MEXT, and Istituto Nazionale di Fisica Nucleare (INFN) in Italy. Part of this work was performed using the computer resources provided by the Institute for the Cosmic-Ray Research (University of Tokyo), CERN, and CNAF (INFN).

APPENDIX: DATA TABLES
The inclusive production rates of π 0 s measured by LHCf with all corrections applied are summarized in Tables IV-XVIII. The ratios of π 0 production rate of MC simulation to data are summarized in Tables XIX-LV Production rate for the π 0 production in the rapidity range 9.0 < y < 9.2 in p + p collisions and in the rapidity range −9.0 > y lab > −9.2 in p + Pb collisions. The rate and corresponding total uncertainty are in units of [  VI. Production rate for the π 0 production in the rapidity range 9.2 < y < 9.4 in p + p collisions and in the rapidity range −9.2 > y lab > −9.4 in p + Pb collisions.
TABLE VIII. Production rate for the π 0 production in the rapidity range 9.6 < y < 9.8 in p + p collisions and in the rapidity range −9.6 > y lab > −9.8 in p + Pb collisions.