Measurement of energy flow, cross section and average inelasticity of forward neutrons produced in $\mathrm{\sqrt{s} = 13 TeV}$ proton-proton collisions with the LHCf Arm2 detector

In this paper, we report the measurement of the energy flow, the cross section and the average inelasticity of forward neutrons (+ antineutrons) produced in $\sqrt{s} = 13$ TeV proton-proton collisions. These quantities are obtained from the inclusive differential production cross section, measured using the LHCf Arm2 detector at the CERN Large Hadron Collider. The measurements are performed in six pseudorapidity regions: three of them ($\eta>10.75$, $8.99<\eta<9.21$ and $8.80<\eta<8.99$), albeit with smaller acceptance and larger uncertainties, were already published in a previous work, whereas the remaining three ($10.06<\eta<10.75$, $9.65<\eta<10.06$ and $8.65<\eta<8.80$) are presented here for the first time. The analysis was carried out using a data set acquired in June 2015 with a corresponding integrated luminosity of $\mathrm{0.194~nb^{-1}}$. Comparing the experimental measurements with the expectations of several hadronic interaction models used to simulate cosmic ray air showers, none of these generators resulted to have a satisfactory agreement in all the phase space selected for the analysis. The inclusive differential production cross section for $\eta>10.75$ is not reproduced by any model, whereas the results still indicate a significant but less serious deviation at lower pseudorapidities. Depending on the pseudorapidity region, the generators showing the best overall agreement with data are either SIBYLL 2.3 or EPOS-LHC. Furthermore, apart from the most forward region, the derived energy flow and cross section distributions are best reproduced by EPOS-LHC. Finally, even if none of the models describe the elasticity distribution in a satisfactory way, the extracted average inelasticity is consistent with the QGSJET II-04 value, while most of the other generators give values that lie just outside the experimental uncertainties.


Introduction
In the last decade, gigantic cosmic ray air shower experiments, like the Pierre Auger Observatory [1] and the Telescope Array [2], have contributed to fundamental progress in the knowledge of Ultra-High Energy Cosmic Rays (UHECRs), i.e. cosmic rays with energies above 10 18 eV [3]. Thanks to the large statistics and high quality results obtained by these experiments, various theories have been suggested to explain the mechanisms responsible for acceleration and propagation of UHECRs. However, despite this significant progress, the current measurements leave several open questions on the nature of UHECRs, mainly relative to their astrophysical sources, acceleration mechanisms and mass composition at high energies.
In gigantic cosmic ray air shower experiments, the properties of the primary cosmic ray particles are indirectly reconstructed from the detection of Extensive Air Showers (EASs), i.e. the showers of secondary particles that primary particles form when they interact in the atmosphere. In particular, mass composition is not a direct information, but can be derived from the measurement of the depth at which the number of particles in the shower reaches its maximum, or from the simultaneous measurements of the number of muons and electrons at the ground level. In order to extract this information, the experimental results must be compared to simulations performed under different assumptions on the relative mass abundance of the primary cosmic rays. Nevertheless, due to the lack of high energy calibration data, the predictions of hadronic interaction models used to simulate the EASs are subject to large systematic errors, which in turn constitute the major source of uncertainty on mass composition measurements [4]. To reduce this systematic contribution, accurate information on several parameters used to model EAS development are needed: inelastic cross section, multiplicity of secondaries and forward energy distributions, from which one can derive the average inelasticity [5] [6]. This information can be extracted from measurements carried out at hadron colliders and the CERN Large Hadron Collider (LHC) [7] is currently the best candidate to accomplish this task. During LHC Run II, the machine operated p-p collisions at √ s = 13 TeV, which is the highest energy ever achieved at a collider and in the range relevant for UHECRs (equivalent to the collision of a 0.9 × 10 16 eV proton with a nucleon of a nucleus in the atmosphere). At LHC, inelastic cross section and multiplicity of secondaries are mainly accessible to the four central detectors [8][9][10][11] and roman pot detectors like TOTEM [12] and ATLAS ALFA [13], while forward energy distributions and average inelasticity can be measured only with dedicated forward detectors. Some of the results obtained during LHC Run I from the first group of experiments (e.g. [14][15][16][17]) were used to tune the so-called post-LHC models, like QGSJET II-04 [18], EPOS-LHC [19] and SIBYLL 2.3 [20]. However, the level of agreement among these models is still far from satisfactory [21], thus highlighting the importance of the second group of experiments. Even if several LHC experiments have forward subdetectors (see [22] for a complete list of forward detectors operating during Run II), LHCf [23] is the only one that has been specifically designed to accurately measure the distributions of very forward neutral particles produced in proton-proton collisions.
Because of the key role that forward baryons have in the development of the air showers and in the abundance of the muonic component [24], a large activity of the LHCf collaboration is dedicated to neutron measurements. A first result relative to the inclusive differential production cross section of very forward neutrons+antineutrons (hereafter simply called neutrons) produced in p-p collisions at √ s = 13 TeV was already published [25]. Here this analysis is extended from three to six pseudorapidity regions, in order to have enough data points to derive three important quantities that are directly connected to EAS development: neutron energy flow, cross section and average inelasticity. It is worth mentioning that these are the first measurements of this type at such a high energy and that their direct comparison to generator expectations gives an important indication on the validity of these generators for air shower simulation.

The LHCf experiment
The LHCf experiment consists of two detectors [26], Arm1 and Arm2, placed in two regions on the opposite sides of LHC Interaction Point 1 (IP1). These regions, called Target Neutral Absorber (TAN), are located at a distance of 141.05 m from IP1, after the dipole magnets that bend the two proton beams. In this position, the two detectors are capable of detecting the neutral particles that are produced in hadronic collisions with a pseudorapidity η > 8.4.
The measurements reported in this paper are obtained using the LHCf Arm2 detector only, because, due to the different geometric acceptance of the two arms, it offers a better coverage of the pseudorapidity regions where the core of the energy flow is expected 1 . The Arm2 detector consists of two calorimetric towers with different transverse sizes: the small tower (25 mm × 25 mm) and the large tower (32 mm × 32 mm). Each tower is made out of 16 Gd 2 SiO 5 (GSO) scintillator layers (1 mm thick), interleaved with 22 tungsten (W) plates (7 mm thick), for a total length of about 21 cm, which is equivalent to 44 X 0 or 1.6 λ I . In addition, 4 xy imaging layers, made out of silicon microstrip detectors with a read-out pitch of 160 µm, are inserted at different depths. In this way, for all particles hitting the detector, it is possible to reconstruct the incident energy (from the scintillator layers) and the transverse position (from the imaging layers). Concerning the reconstruction performances for incident hadrons at the two different energies of 1 TeV and 6.5 TeV, the detection efficiency increases from 52% to 72%, the energy resolution worsens from 28% to 38%, while the position resolution improves from 300 µm to 100 µm [27].
Thanks to the detector design, photons and hadrons are easily separated using the different features of electromagnetic and hadronic showers. However, being unable to distinguish among hadrons, the experimental measurement contains not only neutrons and antineutrons, but also Λ 0 , K 0 L and other neutral and charged hadrons 2 . As described later, this residual background is subtracted at the end of the analysis by the usage of Monte Carlo simulations.

Experimental data set
The results discussed in this paper are based on an experimental data set, relative to p-p collisions at √ s = 13 TeV, which was acquired from 22:32 of June 12th to 1:30 of June 13th 2015 (CEST). This time period corresponds to LHC Fill 3855, a special low luminosity and high β * fill specifically provided for LHCf operations. Low luminosity (3 − 5 × 10 28 cm −2 s −1 ) was required in order to keep a small average number of collisions per bunch crossing (0.007 − 0.012). In this way, given the 15% acceptance of the calorimeter for inelastic collisions, event pile-up at the detector level remains below a negligible 1%. High β * (19 m) was required in order to keep the protons composing the beam almost parallel at the interaction point. In this way, the scattering angle (and therefore the pseudorapidity) of the detected secondary particles can be accurately reconstructed. In Fill 3855, 29 bunches were made to collide at IP1 with a downward half crossing angle of 145 µrad, a value chosen to increase the detector coverage down to a pseudorapidity of 8.4. In addition, there were 6 and 2 non-colliding bunches circulating in the clockwise and counter-clockwise beams, respectively. Using the instantaneous luminosity measured from the ATLAS experiment [28] and taking into account LHCf data acquisition live time, the recorded integrated luminosity was estimated to be 0.194 nb −1 .

Simulation data set
Monte Carlo simulations with the same experimental configuration used for LHC data taking are necessary for four different purposes: estimation of correction factors and systematic uncertainties, validation of the whole analysis procedure, energy spectra unfolding, and data-model comparison. A description of all the simulation data sets, detailing how they were generated, which models were used and which effects were considered, can be found in [25]. Here it is important to remind that, depending on the purpose, a given sample is generated taking into account one or more of the following steps: hadronic collision, transport of secondary products from IP1 to TAN, and detector interactions. The Cosmos 7.633 [29] and EPICS 9.15 [30] libraries were used for all these simulation steps, except for the samples needed for the final data-model comparison. In this last case, QGSJET II-04, EPOS-LHC, SIBYLL 2.3 and DPMJET 3.06 [31] data sets were simulated using CRMC [32], an interface tool for event generators, whereas the PYTHIA 8.212 [33] data set was simulated using its own dedicated interface.

Analysis
The analysis procedure presented in this section is divided in four main steps. First, after defining the event selection criteria, the events in the sample are reconstructed in order to derive the detector level energy distributions for each pseudorapidity region (section 5.1). Second, the correction factors that must be applied to the spectra before and after unfolding are estimated (section 5.2). Third, the effect of migration between energy bins is unfolded using a deconvolution procedure that takes into account the detector response (section 5.3). Fourth, all necessary systematic uncertainties are evaluated and associated to the unfolded spectra (section 5.4). Finally, these distributions are used to derive the experimental measurements, which are then compared to the various generator predictions. These results are presented and discussed in the last two sections, relative to the inclusive differential production cross section (section 6) and to several quantities connected to air shower development (section 7).

Event reconstruction
Apart from small refinements, event reconstruction algorithms and event selection criteria are similar to the ones described in [25]. The incident energy is reconstructed from the total energy deposit in the calorimeter, properly weighting the release in each scintillator layer for the energy lost in the corresponding front tungsten layer and applying position dependent correction factors for shower lateral leakage and light collection efficiency. The transverse position is determined by fitting the shower lateral profile reconstructed in the imaging layers, always looking for a singlehit event, i.e. a single peak in each tower. In case of a multihit event, i.e. more simultaneous hits in the same tower, the transverse position is obviously misreconstructed, an effect that is later corrected for in the analysis. An event is accepted if it satisfies an offline selection that mimics the hardware trigger, but with a slightly higher threshold respect to the one used by the discriminators during data taking. This offline selection requires a raw energy deposit above 850 MeV in at least three consecutive scintillator layers, ensuring a good selection efficiency for hadrons with incident   Table 1: Definition of the six pseudorapidity regions used in the analysis: the coverage is expressed in terms of pseudorapidity η, scattering angle θ and azimuthal angle φ.
energies above 500 GeV. Particle IDentification (PID) exploits the different development of electromagnetic and hadronic showers in the calorimeter to distinguish between photon and hadron candidates. This discrimination is based on the variable where L X% represents the longitudinal depth at which the fraction of the energy deposit with respect to the total release in the calorimeter is X%. Finally, after taking into account the beam center projection on the detector plane, the event is selected if it is within one of the six pseudorapidity regions defined for this analysis. Three of these pseudorapidity intervals (η > 10.75, 8.99 < η < 9.21 and 8.80 < η < 8.99) were already considered in [25], whereas the remaining three (10.06 < η < 10.75, 9.65 < η < 10.06 and 8.65 < η < 8.80) are new additions to this work. The six regions superimposed to the detector area are shown in figure 1 and more details on their definition are given in table 1. Due to the slightly different definition of the three regions already included in [25], and to the refined reconstruction algorithm and selection criteria, the correction factors and the systematic uncertainties discussed here may differ of a few percent from the values reported there.

Correction factors
Correction factors account for the same effects and approximately amount to the same values described in [25]. All contributions are summarized in table 2 and are shortly discussed in the following. Note that correction factors are divided into two groups, depending on whether they are applied before or after unfolding. The reason for this distinction is that, as explained in section 5.3, the final result is obtained making use of energy unfolding based on a single hadron response matrix. Thus, the measured distributions must be corrected before unfolding for all the effects that deviate from the single hadron case. The remaining corrections are applied after unfolding in order to simplify the analysis procedure. Four correction factors are applied before unfolding. The first correction is due to beam background, caused by two different processes: the interaction of primary protons with the residual gas in the beam line and the interaction of secondary particles from collisions with the beam pipe. The former, estimated using non-colliding bunches, leads to a correction of about −1%, whereas the latter, estimated using dedicated simulations, ranges from −10% to 0%, depending on the energy. The second correction accounts for the limited efficiency (hadron identification) and purity (photon contamination) of the PID selection. Correction factors, estimated via template fit of the L 2D distributions of simulated photons and hadrons to the experimental data, range between −15% and +10%. The uncertainty of this correction, obtained from the combination of a dedicated confidence interval fit and different template fit strategies, is generally below 20%. The third correction is due to multihit events, which cannot be treated as a simple background to be subtracted, because this choice would not lead to inclusive measurements. Instead, the corrections applied should restore the ideal case distributions where one is able to reconstruct the n simultaneous hits in the same tower as n independent events. Exploiting two simulation samples generated using QGSJET II-04 and EPOS-LHC, the correction factors, estimated from their average, are in a range between −20% and +10%, with an uncertainty, given by half of their deviation, mostly below 10%. The fourth correction accounts for fake events, i.e. the events that, due to detector misreconstructions, are incorrectly included in the measured distributions. Correction factors, estimated employing the same simulation samples and the same computation method just described, range from −10% to −1%, whereas the uncertainty is always smaller than 1%.
Two correction factors are applied after unfolding. The first correction accounts for missed events, i.e. the events that, due to detector misreconstructions, are incorrectly excluded from the measured distributions. Correction factors, estimated employing the same simulation samples and the same computation method used for fake events, range from +30% to more than +100%, whereas the uncertainty is always smaller than 1%. As expected, missed events are the dominant contribution to correction factors because of the limited hadron detection efficiency, which decreases even more at low energies. Two factors contribute to the uncertainty on fake and missed corrections. The first one is the dependence on the simulation model used to generate the collisions, which, as discussed in this section, is negligible respect to the other contributions. The second one is the dependence on the simulation model used to describe the interaction with the detector, which, as discussed in section 5.4, is not negligible. The last correction is due to hadron contamination, i.e. the fact that a non-negligible fraction of hadrons reaching the detector are not actually neutrons, but other particles that cannot be distinguished from them. Exploiting five simulation samples generated using all the models discussed in section 4, the correction factors, estimated from their average, are in a range between −60% and −5%, with an uncertainty, given by their maximum deviation, mostly below 20%.

Spectra unfolding
Given the limited hadron energy resolution, spectra unfolding must be applied to deconvolute the measurements from the detector response. The iterative Bayesian method [34], implemented in the RooUnfold package [35], was used for this purpose. The best estimation for the final unfolded distributions was obtained choosing a flat prior as initial hypothesis and constructing the response matrix from a single neutron flat energy sample simulated using the DPMJET 3.04 model. As discussed in section 5.4, in order to estimate the uncertainties due to the deconvolution, the unfolding procedure was repeated several times with different priors, response matrices, and test input spectra. As a general convergence criteria, the iterative method was stopped when ∆χ 2 , the χ 2 variation between the outputs of two consecutive iterations, was below 1. This threshold was chosen as a compromise between convergence requirements and the number of iterations needed to reach that convergence. As an example, in the case of the data points in the final distributions, this convergence criteria corresponds to 51, 35, 21, 13, 11 and 10 iterations for pseudorapidity regions A, B, C, D, E and F, respectively. Table 3: The (minimum, maximum) values of the relative systematic uncertainties for each contribution discussed here, separately reported for regions A, B, C, D, E and F. The top five and the bottom four uncertainties are respectively applied before and after unfolding. Note that the two energy independent systematic uncertainties, relative to integrated luminosity (1.9%) and detection efficiency (6%), applied after unfolding, are not reported in this table.

Systematic uncertainties
Systematic uncertainties are due to the same effects and approximately amount to the same values described in [25]. Considering all the uncertainty sources as independent, the final uncertainty is obtained from their quadrature sum. All contributions are summarized in table 3 and are shortly discussed in the following. Note that systematic uncertainties are divided into two groups, depending on whether they are applied before or after unfolding. The reason for this distinction is that some uncertainties, which are relative to experimental conditions or detector response, can be estimated only on the measured distributions, whereas other uncertainties, which are energy independent or unfolding related, can be directly applied on the unfolded distributions. Thus, while the latter are immediately added to the final results, the former must be fully propagated through the unfolding procedure. Three uncertainties are applied before unfolding. The first uncertainty is relative to the energy scale, which is known with an error of ±4.5%, determined from the detector calibration at SPS, the LHC data taking conditions, and the observed π 0 mass peak shift.
The contribution to the final uncertainty, estimated by shifting the energy scale for the corresponding error and comparing them with the nominal distribution, increases from a value of about 10% at low energy up to more than 100% at high energy, thus representing the dominant contribution to the systematic uncertainty. The second uncertainty is associated to the beam center projection on the detector plane, which was determined by fitting the position distribution of high energy hadrons with a two dimensional function, leading to an error of ±0.3 mm. The corresponding uncertainty, obtained by shifting the beam center for the corresponding error on the two axes and comparing them with the nominal distribution, is generally below 10%. The third uncertainty comes from the position resolution, which is not taken into account in the deconvolution in order to simplify the unfolding procedure. Exploiting two simulation samples generated with QGSJET II-04 and EPOS-LHC, the corresponding uncertainty, estimated from the ratio between the spectra obtained using true and reconstructed position, is mostly below 10%. Note that, for each event, these three uncertainties affect the reconstructed energy (energy scale and position resolution) and/or the reconstructed pseudorapidity (beam center and position resolution).
Three uncertainties are applied after unfolding. The first uncertainty is associated to the integrated luminosity, as derived from ATLAS measurements, and has an energy independent value of 1.9%. The second uncertainty comes from the detection efficiency, which is affected by the error on the hadron-detector interaction cross section at high energies. In order to estimate the impact of this parameter on the experimental measurements, its value was extracted both for the DPMJET 3.04 model in EPICS and for the QGSP_BERT 4.0 model in GEANT4 [36]. Then, assuming that the difference on the hadron-detector interaction cross section leads to a similar difference on the detection efficiency, the corresponding uncertainty is estimated from the maximum deviation below 6.5 TeV, amounting to an energy independent value of 6%. The third uncertainty accounts for three different contributions relative to spectrum unfolding: method uncertainty, due to the deviation between unfolded spectrum and true distribution, estimated using different generators; prior uncertainty, due to the dependence of the unfolded spectrum on the prior chosen as input of the algorithm, estimated using different priors; algorithm uncertainty, due to the dependence of the unfolded result on the deconvolution algorithm, estimated using different approaches (Bayesian, SVD [37] and TUnfold [38]). Each of these three uncertainties is mostly below 20%, except at very high energies, where they can increase to more than 100%, due to the small number of events.
Finally, additional sources of uncertainties are associated to the correction factors discussed in section 5.2: in this case, each of these contributions belongs to the group of uncertainties applied either before or after unfolding, depending on where the corresponding correction is considered in the analysis.

Results
In this section, the measurement of neutron production is presented and the comparison with model predictions is discussed. Distributions are expressed in terms of the inclusive differential production cross section dσ n /dE for six pseudorapidity intervals. In each interval, the distribution is corrected to take into account the limited coverage of the detector in terms of the azimuthal angle. The data are then compared with the generator predictions, using for each model its own inelastic cross section.
The dσ n /dE unfolded distributions are shown in figure 2. The measurements are limited to energies above 500 GeV because, as described in section 5.1, the hardware trigger was optimized to have a good detection efficiency above this value. The total uncertainty is given by the quadratic sum of all statistical and systematic sources: the systematic un-certainties are the dominant contribution everywhere, whereas the statistical uncertainties are generally smaller than 5%. The numerical values of the inclusive differential production cross section are also summarized in appendix A.
At a first glance, it is interesting to note that all experimental distributions exhibit a peak structure, whose position progressively moves to lower energy as pseudorapidity decreases, from about 5 TeV for η > 10.75 to about 1.5 TeV for 8.65 < η < 8.80. Furthermore, it is found that all these distributions, except the one relative to the most forward region, can be fitted by a very simple function like an asymmetric gaussian -a fact that could give a hint on the underlying physical process responsible for neutron production.
Considering the data-model comparison, the largest deviation between model distributions and experimental data is at η > 10.75, as already noted in [25]. In this region, LHCf measurements indicate a peak structure at around 5 TeV which is not predicted by any generator and, in addition, all models underestimate the total production cross section by at least 20% (QGSJET II-04). No model completely reproduces the experimental measurements even in the other five regions, although deviations are smaller than at η > 10.75. Depending on the pseudorapidity region, the generators showing the best overall agreement with data are either SIBYLL 2.3 or EPOS-LHC. At 10.06 < η < 10.75, EPOS-LHC reproduces well the peak position but it underestimates neutron production, whereas, at 9.65 < η < 10.06, SIBYLL 2.3 reproduces well the peak position but it overestimates neutron production. In all the remaining three regions, it is interesting to note that SIBYLL 2.3 is compatible with the experimental measurements in the energy range between 1.5 TeV and 2.5 TeV, where the production peak is located. However, SIBYLL 2.3 is the model showing the best overall agreement with data only at 8.99 < η < 9.21, whereas EPOS-LHC leads to better results at 8.80 < η < 8.99 and 8.65 < η < 8.80.

Discussion
The results presented in section 6 are interesting by themselves and useful for model tuning. However, in order to test the validity of generators in the context of cosmic ray measurements, it is better to express these results in terms of quantities that are directly connected to air shower physics. In this section, the distributions shown in figure 2 are therefore used to extract three important parameters: energy flow, cross section and average inelasticity of forward neutrons.
The differential energy flow dE n /dη and the differential cross section dσ n /dη are expressed as a function of pseudorapidity in the following manner. For each region, the corresponding mean pseudorapidity, η, and pseudorapidity interval, ∆η, are given by the average value and the distance of the two extremes, respectively 3 . In a similar way, for each energy bin i of the dσ n /dE distribution, the mean energy, E i , and the energy interval, dE i , are given by the average value and the distance of the two extremes, respectively. Thus, dE n /dη and dσ n /dη are given by where the TOTEM measurement of σ inel = (79.5 ± 1.8) mb was used for the total inelastic cross section (and the corresponding uncertainty) of p-p collisions at √ s = 13 TeV [39].
Note that these quantities must be corrected to take into account the contribution of neutrons below 500 GeV, which are not included in the dσ n /dE distributions. Correction factors were estimated from five simulation samples generated using all the models discussed in section 4, taking the average as best estimation and the maximum deviation as uncertainty. For dE n /dη, the corrections are at most 1.5% and absolute uncertainties below 1%, whereas, for dσ n /dη, the corrections are at most 8% and absolute uncertainties below 5%. Another important aspect is how the several sources of uncertainty acting on dσ n /dE contribute to the uncertainty of these two derived quantities. First, all contributions are assumed to be independent, so that they can be added in quadrature. Then, the uncertainties are divided in two groups: bin-by-bin independent (only statistical) and bin-by-bin fullycorrelated (all systematic). For each bin-by-bin fully-correlated source, its contribution to the final uncertainty is given by the maximum deviation between the quantity derived from the nominal distribution and the one derived from the distributions corresponding to all possible shifts induced by that uncertainty. Figure 3 shows the differential energy flow and the differential production cross section measured using the LHCf Arm2 detector. The numerical values of dσ n /dη and dE n /dη are also summarized in table 4. As expected, in both cases all models underestimate these quantities at η > 10.75, whereas for the other regions the disagreement between generators and data is less serious. In particular, EPOS-LHC is the model showing the best overall agreement with the experimental measurements, QGSJET II-04 is consistent with data at 9.65 < η < 10.75, but underestimates these quantities at 8.65 < η < 9.21, whereas SIBYLL 2.3 reasonably reproduces dσ n /dη, but strongly overestimates dE n /dη. Furthermore, even if there are no data points for η values between 9.21 and 9.65, experimental measurements seem to indicate that the peak in dE n /dη must be located here, if this distribution has such kind of structure as predicted by all models except EPOS-LHC.
The average inelasticity 1 − k can be extracted from the elasticity distribution, where the elasticity is defined as k = 2E/ √ s, and E is the energy of the leading particle. The leading particle is the particle (not necessarily a baryon) that, on one of the two sides of the collision (let us say η > 0), carries the largest fraction of the primary proton momentum. Note that, being sensitive to neutrons only, the LHCf measurements are restricted to the specific case where the leading particle is a neutron. To distinguish the quantities that are measured here from their general definition, in the following they are indicated as k n and 1 − k n , instead of k and 1 − k . The role of neutrons as forward leading particles was investigated with the use of generators. Despite the large variations among the different models, two important features were clearly highlighted. The first one is that, depending on the energy bin, the fraction of events where the leading particle is a neutron ranges from about 5 to 35%. The second one is that, for energies above half the beam energy, almost 100% of the neutrons produced from the collisions are leading particles. In order to obtain the elasticity distribution, the dσ n /dE contributions of all the six regions are summed in a single histogram. Then, the x axis is rescaled to the beam energy and the y axis is multiplied for the bin width, so that the distribution represents the total production    Table 4: Differential energy flow dE n /dη and differential cross section dσ n /dη of neutrons produced in p-p collisions at √ s = 13 TeV, measured using the LHCf Arm2 detector. Upper and lower uncertainties are also reported. The values are relative to the experimental measurements with (E > 0 GeV) and without (E > 500 GeV) the simulation-driven correction factors for the limited detection efficiency below 500 GeV. The last two columns correspond to the numbers used for the experimental points shown in figure 3.
cross section σ n as a function of elasticity k n . At this point, a correction must be applied to take into account two different effects: the first one is due to the fact that the detector has a limited pseudorapidity coverage; the second one is due to the fact that not all neutrons are leading particles. These two effects are considered together in a single correction factor that was obtained from five simulation samples generated using all the models discussed in section 4, taking the average as best estimate and the maximum deviation as uncertainty.
Corrections range between 1% and 70%, whereas absolute uncertainties go from 5% to 70%. The several sources of uncertainties acting on the dσ n /dE distributions contribute TeV. These quantities, measured using the LHCf Arm2 detector, are only relative to the events where the leading particle is a neutron. Black markers represent the experimental data with the quadratic sum of statistical and systematic uncertainties. Solid lines (left) and full circles (right) refer to model predictions at the generator level, obtained using only the events where the leading particle is a neutron. In order to compare this approach to the general case, 1 − k , the average inelasticity obtained using all the events independently of the nature of the leading particle, is also reported as open circles in the right figure. to the uncertainty on the elasticity distribution in a similar way to what was previously described, i.e. assuming that all contributions are independent and dividing them in bin-bybin independent (only statistical) and bin-by-bin fully-correlated (all systematic) sources. Note that, differently from the previous case, the term bin does not refer to the energy bin, but to the pseudorapidity bin, because the summation index is on pseudorapidity and not on energy. The elasticity distribution cannot directly be used to extract the error on the average inelasticity, because systematic uncertainties are correlated both on energy and on pseudorapidity. The entire procedure must therefore be repeated to extract the uncertainty, but an average value is computed instead of building a histogram, so that both sources of correlation are correctly considered in the estimation of the uncertainty. This value is then corrected to take into account the contribution of neutrons below 500 GeV, which are not included in the dσ n /dE distributions. The correction factor, estimated in a similar way to what was previously discussed, amounts to a value of (0.4 ± 0.4)%. Figure 4 shows the inclusive production cross section as a function of elasticity and the average inelasticity extracted from that distribution, measured using the LHCf Arm2 detector. In the left plot, the contribution to the error bars of σ n is dominated by the uncertainty on dσ n /dE for large values of k n and by the uncertainty on elasticity correction factors for small values of k n . As can be seen, no generator shows a satisfactory agreement with the experimental measurements for all the elasticity range. In particular, for 0.3 < k n < 0.6, the data lies between SIBYLL 2.3 and EPOS-LHC, QGSJET II-04 predictions. Outside this region, the model most consistent with the data is SIBYLL 2.3 for k n < 0.3, 0.6 < k n < 0.8 and EPOS-LHC for k n > 0.8. The average inelasticity extracted from the elasticity distribution with the method just discussed is 1 − k n = (0.536 +0.031 −0.037 ). In the right plot, this value is compared to model predictions for two different quantities: 1 − k n (experimental case) and 1−k (general case). As can be seen, all generators predict similar numbers for these two quantities, indicating that, even if limited only to the events where the leading particle is a neutron, this measurement provides a relevant information also for the general definition of inelasticity. With a 1 − k n value of 0.533, QGSJET II-04 is the only model consistent with the experimental measurement, even if all the remaining generators except PYTHIA 8.212 lie just outside the error bars. As a final remark, it is important to stress that, despite the relative agreement found on the average inelasticity, the elasticity distribution predicted by all models differs from the one observed in the data.

Summary
The LHCf experiment measured the inclusive differential production cross section of forward neutrons (+ antineutrons) produced in √ s = 13 TeV proton-proton collisions. , whereas at lower pseudorapidities the data-model agreement is better, but still not satisfactory. In this last case, depending on the pseudorapidity region, the generator showing the best overall agreement with data is either SIBYLL 2.3 or EPOS-LHC. From the inclusive differential production cross section, three quantities directly connected to air shower development were derived: neutron energy flow, cross section and average inelasticity. The first and the second quantities are well reproduced by EPOS-LHC, except for the most forward region, whereas QGSJET II-04 gives a consistent value with the data for the third quantity, with most generators lying just outside the experimental uncertainties. However, despite the relative agreement observed on the average inelasticity of leading neutrons, no generator satisfactorily reproduces the corresponding elasticity distribution.
The LHCf results are the first measurements of forward neutron production at such a high energy and, given the significant deviation between generators and data, they suggest that all models would benefit from a tuning that takes these results into account. A further development of this analysis is possible thanks to the fact that the LHCf and ATLAS experiments had common data acquisition during the operations where the data considered here was recorded. By exploiting the ATLAS information in the central region, it will be possible to investigate the different mechanisms responsible for neutron production in the forward region [40] and the amount of correlation between central and forward hadron production [41].
A  Table 5: Inclusive differential neutron production cross section for p-p collisions at √ s = 13 TeV, relative to three of the six regions used in the analysis. These regions are located on the small tower of the LHCf Arm2 detector. Upper and lower uncertainties, expressed as the quadratic sum of statistical and systematics contributions, are also reported. ( Table 6: Inclusive differential neutron production cross section for p-p collisions at √ s = 13 TeV, relative to three of the six regions used in the analysis. These regions are located on the large tower of the LHCf Arm2 detector. Upper and lower uncertainties, expressed as the quadratic sum of statistical and systematics contributions, are also reported.