Measurement of very forward neutron energy spectra for 7 TeV proton--proton collisions at the Large Hadron Collider

The Large Hadron Collider forward (LHCf) experiment is designed to use the LHC to verify the hadronic-interaction models used in cosmic-ray physics. Forward baryon production is one of the crucial points to understand the development of cosmic-ray showers. We report the neutron-energy spectra for LHC $\sqrt{s}$ = 7 TeV proton--proton collisions with the pseudo-rapidity $\eta$ ranging from 8.81 to 8.99, from 8.99 to 9.22, and from 10.76 to infinity. The measured energy spectra obtained from the two independent calorimeters of Arm1 and Arm2 show the same characteristic feature before unfolding the difference in the detector responses. We unfolded the measured spectra by using the multidimensional unfolding method based on Bayesian theory, and the unfolded spectra were compared with current hadronic-interaction models. The QGSJET II-03 model predicts a high neutron production rate at the highest pseudo-rapidity range similar to our results and the DPMJET 3.04 model describes our results well at the lower pseudo-rapidity ranges. However no model perfectly explains the experimental results in the whole pseudo-rapidity range. The experimental data indicate the most abundant neutron production rate relative to the photon production, which does not agree with predictions of the models.


Introduction
The forward particle production process induced by collisions of high-energy particles is a poorly understood phenomenon in high-energy physics. Though it is important to understand the development of cosmicray showers in the atmosphere, the validity of hadronicinteraction models have not been sufficiently verified at * corresponding author Email address: kawade@stelab.nagoya-u.ac.jp (K. Kawade) energies for ultra-high-energy cosmic rays (UHECRs, > 10 18 eV) because of lack of experimental data in this energy range. This lack of data results into a large uncertainty in the interpretation of the energy and chemical composition of UHECRs. Forward baryons play a very important role in the development of cosmic-ray showers. If forward baryons carry more collision energy, cosmic-ray showers develop much deeper in the atmosphere, and vice versa. However, in the energy range of UHECRs, the predictions by current models differ significantly among themselves.
The excess of muons at ground level is reported as one of the problems in the cosmic-ray shower observations. The number of muons observed by the surface detector array of the Pierre Auger Observatory (PAO) [1] is higher than the one expected based on the energy determined by the fluorescence detectors even if a heavy primary mass is assumed [2]. It is suggested that the number of (anti) baryons generated in the forward region is strongly related to the number of muons observed by PAO at the ground [3]. Therefore, baryon production at the very forward region is quite important to understand cosmic-ray showers.
In this paper, we report the results of analyzing the data of the Large Hadron Collider forward (LHCf) experiment for forward neutron spectra. Forward baryon spectra at the laboratory equivalent energy of 2.5×10 16 eV ( √ s=7 TeV) will be a crucial input to improve the hadronic-interaction models used in the air shower analyses.

LHCf experiment
The LHCf experiment is designed to use the LHC to verify the hadronic-interaction models used in cosmicray experiments [4,5]. Two independent detectors named Arm1 and Arm2 were installed in the detector installation slots of the TANs located 140 m away from the interaction point 1 (IP1). Because charged particles are swept away by the D1 bending magnets, LHCf can measure only neutral particles in the very forward region of the LHC (pseudo rapidity |η| > 8.4). Both detectors have two different sampling calorimeters with 44 radiation lengths (1.6 hadron-interaction lengths) of tungsten plates and 16 layers of sampling scintillators [5]. Four layers of the position sensors (SciFi in Arm1 and silicon micro-strip sensor in Arm2) can measure the hit position transverse to the beam direction. The transverse dimensions of the calorimeters are 20 mm × 20 mm and 40 mm × 40 mm in Arm1, and 25 mm × 25 mm and 32 mm × 32 mm in Arm2. The details of the detector performance during the 2009-2010 proton-proton collisions are reported in [6].
The performance of the LHCf detectors for hadron measurements was studied by Monte Carlo (MC) simulations and confirmed by using 350 GeV proton beams at CERN-SPS [7]. Depending on the incident-neutron energy, the energy resolution and position resolution are about 40% and 0.1-1.3 mm, respectively. The detection efficiency for neutrons was estimated to be 70%-80% for neutrons above 500 GeV.
In this paper we assume hadronic showers are produced by neutrons. According to the EPOS 1.99 gener-ator [8], about 10% of other hadrons, i.e., Λs and K 0 s, are also included in the data and this fraction and the energy dependence are model dependent.

Data used in analysis
The data used in this analysis were obtained on May 15, 2010 from proton-proton collisions at √ s = 7 TeV (LHC Fill # 1104). The typical luminosity corresponding to this fill derived from the counting rate of the LHCf front counters [9] was (6.3-6.5) ×10 28 cm −2 s −1 . The data set was the same as one used in the previously published photon analysis results, and additional details can be found in [10]. The trigger for LHCf events was generated when signals from any three successive scintillation layers in any calorimeter exceeded a predefined threshold (typically 130 minimum ionizing particles (MIPs)). Data acquisition (DAQ) was performed under 85.7% (Arm1) and 67.0% (Arm2) average livetimes.
Taking the DAQ live-times into account, the integrated luminosities of the data set were 0.68 nb −1 for Arm1 and 0.53 nb −1 for Arm2, each with ±6.1% uncertainty. The numbers of inelastic collisions were about 48M and 38M collisions for Arm1 and Arm2, respectively.
MC predictions were conducted with the generators DPMJET 3.04 [11], EPOS 1.99 [8], PYTHIA 8.145 [12], QGSJET II-03 [13], and SYBILL 2.1 [14] and compared with the experimental results. In the MC simulations, the COSMOS (v7.49) and EPICS (v8.81) [15] libraries that are used in air-shower and detector simulations were used to simulate the flight of particles from the IP1 to the detectors and the response of the detectors. About 10M inelastic collisions were simulated for each model.

Event reconstruction
Initially, the offline-event selection was applied when energy depositions equivalent to more than 200 MIPs were recorded for three successive layers in addition to the experimental trigger. The positions at which particles hit the detector were determined by using the position sensors. Because reconstruction of the events is difficult at the edge of the calorimeters due to large fluctuations in the energy deposition, events within 2 mm from the edge were discarded from the analysis. The lateral shower leakage caused by the limited lateral size of the detectors deteriorates the energy resolution. This position-dependent leakage effect was corrected by using the transverse-hit position measured by the position sensors.
Particle identification (PID) between neutrons and photons was based on the difference in the longitudinal shape of the shower development. Two simple parameters called L 20% and L 90% were introduced to characterize the shower shape. These parameters were defined as the depths containing 20% and 90%, respectively, of the total energy deposited within the layers. Considering the correlation between L 20% and L 90% an optimized parameter L 2D was defined as L 2D = L 90% − 1/4 × L 20% to improve the selection efficiency and purity than previous analyses. Figure 1 shows the L 2D distributions. Two distinct peaks are identified in the observed L 2D distribution indicated by the closed circles. The histograms correspond to the MC prediction of pure photons (red) and pure neutrons (blue) and they are called 'templates' hereafter. The templates were produced by accumulating the MC simulation as a mixture of five models, DP-MJET 3.04, EPOS 1.99, PYTHIA 8.145, QGSJET II-03, and SYBILL 2.1, with same statistics for each.
To obtain neutron spectra, only events with the parameter L 2D exceeding a certain threshold were identified as neutron-like events. The effects of the neutron selection efficiency and the photon contamination were corrected by using the efficiency ǫ and purity P (i.e., multiplying by P/ǫ ) determined by MC simulations and the template fitting method [6,10,16], respectively. To estimate the photon contamination, the templates for photons and hadrons were independently scaled to reproduce the experimental results (called method A). The open circles in Figure 1 represent the fitting result with method A. Neutron selection criteria, as indicated in Figure 1, were chosen to maximize ǫ × P. To cope with energy dependence, events were separated into eight different energy ranges according to reconstructed energy.
After correcting the PID efficiency and purity, we obtain the production rate of neutrons as a function of obtained energy that was determined from the total deposited energy in the calorimeter and from the identity of the particle. Details of the event reconstruction for neutrons are summarized in [7]. About 0.3 million neutron-like events passed the PID selection for each arm.
To combine the results of Arm1 and Arm2, we selected events that occurred within the common rapidity regions. Events within 6 mm (η > 10.76) from the beam center were selected for the small towers. The large towers of Arm1 and Arm2 were divided into two regions. The inner region "A" was defined by a radius of 28-35 mm (8.99 < η < 9.22) from the beam center, whereas the outer region "B" was defined by a radius of 35-42 mm (8.81 < η < 8.99). For the analysis, we used azimuthal-angle intervals dφ of 360 • for the small towers and 20 • for the large towers.

Energy scale
In order to determine energy scale and estimate its systematic uncertainty, we followed the previous analyses [10,17]. We observed invariant mass excesses of 8.1% (Arm1) and 3.7% (Arm2) compared with the π 0 mass reconstructed in the MC simulations. These excesses were corrected in this analysis. Based on this mass excess correction and known calibration uncertainty, totally ±5.6% (Arm1) and ±4.4% (Arm2) were assigned as systematic uncertainties with respect to the central value of the mass shift. In addition, according to the differences between the SPS beam test and MC simulation in the reconstructed energy of 350 GeV proton showers [7], +2.0% (Arm1) and -3.8% (Arm2) errors were quadratically added to the respective energy scale uncertainties.

PID
Method A described in section 3.2 did not perfectly -500 -1000 -1500 -2000 -2500 -3000 -3500 3500< Small tower -8.0% -3.4% -1.8% -1.1% 0.0% 0.0% 0.0% -0.7% Large tower A -12.8% -9.7% -4.2% -1.7% -0.7% -0.6% 0.1% -0.8% Large tower B -17.7% -11.5% -3.6% -1.2% -0.7% -0.6% 0.4% -1.3% reproduce the experimental results. To estimate systematic effects from these differences, we used a more artificial method (Method B) that allows longitudinal displacements and modifications in width of the distributions to fit the experimental results. The systematic uncertainty from the PID process was estimated by comparing the results using the Method A and Method B. The energy-dependent PID systematic uncertainty is mostly 1% above 1.5 TeV and 12% at most below this energy. The relative differences in the neutron production rate defined as are summarized in Table 1. Here, ǫ A (ǫ B ) and P A (P B ) are the efficiency and purity determined by Method A (Method B), respectively. Final results are given using the Method A, while the differences shown in Table 1 are taken as a part of the systematic uncertainties.

Multi-hit
When two or more particles enter one of the LHCf calorimeters, these events are called 'multihit' events. Because discriminating between single and multihit events for neutron-like events was difficult because of large fluctuation of hadronic showers, rejection of multihit events causes a large systematic uncertainty. Because the multihit event rate was predicted by MC simulations to be less than a few percent, all events used for this analysis were treated as single-hit events. The systematic effects on the energy spectra from the multihit events were studied by using MC simulations. We tested the difference of the spectra with current method and those with ideal multihit reconstruction by using the MC study. The difference was less than 6% above 500 GeV as summarized in Table 2 and was taken into account as a part of the systematic uncertainties.

Position resolution
Because the resolutions of transverse hit position for neutrons are 0.1-1.3 mm depending on the neutron energy [7], some events were reconstructed as hits at positions far from the true hit position. Thus, particles hitting outside of the fiducial area may be contaminated and vice versa. The effect of the position resolution was estimated by using MC simulations.
The difference between neutron-energy spectra selected by the reconstructed position and the true hit positions were calculated with EPOS 1.99, QGSJET II-03, SYBILL 2.1, DPMJET 3.04, and PYTHIA 8.145. Because no significant dependence on model was found, the average of the predictions by the five different models was assigned to the systematic uncertainty. The systematic error from these effects was less than 8.4% at energies above 500 GeV and are summarized in Table 3.

Other systematic errors
As similar as the previous study [10], the other systematic errors, such as the integrated luminosity (±6.1%) and position of beam center (typically ±3-10%) were taken into account. The systematic uncertainty from pile-up events (0.2%) was negligibly small. Figure 2 shows the energy spectra of forward neutrons measured by the LHCf Arm1 and LHCf Arm2 detectors. The energy scale correction described above was applied in these spectra (-8.1% for Arm1 and -3.7% for Arm2). The vertical axes were normalized to the number of inelastic collisions per GeV. The vertical bars represent the statistical (they are very small) and systematic uncertainties except the energy scale and luminosity uncertainties 1

. The closed and open circles
show the results of the Arm1 and Arm2 detectors, respectively. A small difference in the detection efficiencies of the Arm1 and Arm2 detectors were not corrected  2.0% Large tower A -54.2% -1.0% 0.9% -2.5% -0.6% -3.3% -3.6% -5.0% Large tower B -11.0% -0.0% 1.4% 0.6% -1.1% -1.0% -3.7% -5.7% here because it is treated in the unfolding process discussed in Section 4.2. In spite of the difference in detector response, data from both Arms show the same characteristic feature of the spectra. Figures 3 compares the energy spectra measured by the Arm1 detector with the MC predictions. Colored lines indicate MC predictions by EPOS 1.99 (magenta), QGSJET II-03 (blue), SYBILL 2.1 (green), DPMJET 3.04 (red), and PYTHIA 8.145 (yellow). The model spectra were obtained from full detector simulations and taking account of the same reconstruction process as the experimental data. None of the models perfectly matches the experimental data. Experimental result indicates the hardest spectra than any other MC predictions.

Spectra unfolding
To estimate the true energy distribution, we used the multidimensional-spectra unfolding method [18] with the variables energy and transeverse momentum (p T ). To create training samples for the unfolding process, we used MC simulations with neutrons having a flat energy spectrum from 50 to 3500 GeV and an uniform injection to the detector plane. The training samples were reconstructed by using the same method as the experimental data. Performance of the unfolding method was checked by applying the unfolding process to the MC spectra and comparing the results with the true spectra. The upper panel in Figure 4 shows the unfolded spec- tra for the DPMJET 3.04 and EPOS 1.99 models together with the true spectra at the small tower of Arm1.
Here "true spectra" means the true neutron energy distribution of the MC events after acceptance and trigger threshold were applied. The bottom panel shows the ratio of the unfolded spectra to the true spectra. The differences between the unfolded and true spectra were mostly within 20% except at the highest energy bins (50-100%). These systematic differences were due to the choice of the flat energy distribution as a training sample. We found that the difference did not strongly depend among five input models. Thus we applied another correction; dividing the unfolded spectra by the average of differences. The differences among the five models, typically ±10%, were considered as a part of the systematic uncertainties. Figure 5 shows the unfolded experimental spectra measured by the Arm1 and Arm2 detectors for each rapidity range. The hatched areas show the Arm1 systematic errors, and the bars represent the Arm2 systematic errors. The detection efficiency of neutrons and the correction in the PID efficiency and purity were also considered. The results below 500 GeV were not shown because of the large systematic errors. The unfolded spectra from both Arms show good agreement within systematic errors.
The experimental spectra were combined according to the previously used method [17]. It was assumed that the systematic uncertainties caused by the energy scale, PID correction, beam center position, multihit events, and position resolution have bin-to-bin correlations while the other elements are independent between bins. The systematic uncertainties except the unfolding processes were thought to be fully uncorrelated between Arm1 and Arm2. Because we treated the systematic uncertainties of the unfolding processes identical in Arm1 and Arm2, these errors were quadratically added after the combining process.
The differential neutron production cross sections dσ n /dE were calculated from the unfolded experimental spectra by using where dN(∆η, ∆E) means the number of neutrons observed in the each rapidity range, ∆η, and each energy bin, ∆E. L is the integrated luminosity corresponding to the data set. The cross sections are summarized in Table A.5. Figure 6 shows the combined Arm1 and Arm2 spectra together with the model predictions. The experimental results indicate that the highest neutron production rate compared with the MC models occurs at   the most forward rapidity. The QGSJET II-03 model predicts a neutron production rate similar to the experimental results at the largest rapidity range. However, the DPMJET 3.04 model predicts neutron production rates better at the smaller rapidity ranges.
The neutron-to-photon ratios (N n /N γ ) at three different rapidity regions are extracted after unfolding and summarized in Table 4. Here N n and N γ are the number of neutrons and number of photons, respectively, with energies greater than 100 GeV. The numbers of photons were obtained from the previous analysis [10] and the same analysis for the pseudo-rapidity range 8.99-9.22 defined in this study. The experimental data show the highest neutron ratio compared with the hadronicinteraction models for all rapidity regions.

Summary and discussions
An initial analysis of neutron spectra at the very forward region of the LHC is presented in this paper. The data were acquired in May 2010 at the LHC from √ s = 7 TeV proton-proton collisions with integrated luminosity of 0.68 nb −1 and 0.53 nb −1 for the LHCf Arm1 and Arm2 detectors, respectively.
The neutron energy spectra were analyzed in the three different rapidity regions. The results obtained from the two independent calorimeters of Arm1 and Arm2 are consistent each other. The measured spectra were combined and unfolded by using a two-dimensional unfolding method based on Bayesian theory.
The experimental results, both in folded and unfolded were compared with the MC predictions of QGSJET II-03, EPOS 1.99, DPMJET 3.04, PYTHIA 8.145, and SYBILL 2.1; however, no model perfectly reproduces the experimental results in the whole pseudorapidity range. Moreover, compared with the hadronicinteraction models, the experimental results show the most abundant neutron production relative to the photon production rate.
Because of the difference in the phase space, neutrons at lower rapidity carry more energy than those at higher rapidity. EPOS and QGSJET II 2 , current standard models in the air shower analysis, underestimate the neutron production to about 30% at the lower rapidity regions. It is interesting to investigate how these differences affect to the development of air showers.
The differential cross sections for neutron production at very forward rapidity were measured at the Intersecting Storage Ring for proton-proton collisions at √ s = 30.6-62.7 GeV [19,20] and by the PHENIX experiment at the Relativistic Heavy Ion Collider for proton-proton collisions at √ s = 200 GeV [21]. The results of previous experiments were consistent with Feynman x (x F ) scaling and do not depend on the collision energy. To accurately extrapolate the hadronic-interaction models to the high-energy range (independently of whether or not the x F scaling in the neutron production is still relevant at energies in the TeV range) is also an important issue. LHCf will extend the energy range to test the Feynman scaling with the proton-proton collision data obtained at √ s= 0.9, 2.76, 7 TeV and to be obtained at √ s=13 TeV in 2015.