Hadronic interaction studied by TA

Two studies by the Telescope Array group related to the hadronic interaction observed with Extensive Air Showers are reviewed. (1) Inelastic p-air cross section σinel p−air = 567.0±70.5 [stat] +29 −25 [sys] mb and total pp cross section σtot p−p = 170 +48 −44 [stat] +19 −17 [sys] mb were determined using the 5 years of TA hybrid data with one of the 3 FD stations. These results at the highest energy √ sNN=95 TeV showed good agreements with the extrapolation from the previous measurements and model predictions. (2) The signal sizes of SD were compared between data and MC using 7 years of TA SD data in the energy range from 1018.8 eV to 1019.2 eV. It was found that the data/MC ratios exceed unity and the deviation becomes larger when the expected fraction of muon signal, defined as muon purity P, is higher. The results support the muon excess (with respect to MC) problem reported by the previous observations.


Introduction
Observations of Ultra-High Energy Cosmic Rays (UHECR) through extensive air showers (EAS) are the only way to access high-energy particle physics beyond the LHC energy. So far the data from LHC are used to study hadronic interaction up to √ s=13 TeV, which corresponds to cosmic-ray equivalent energy E CR ∼10 17 eV in case of proton primaries, while UHECR observations can access up to √ s=440 TeV or E CR =10 20 eV. It is known that some EAS observables such as number of muons are not well reproduced even using the post-LHC interaction models in air shower simulation [1]. Investigation of high-energy interactions using EAS observations are thus an important theme for both cosmic-ray physics and high-energy particle physics.
The Telescope Array (TA), the largest EAS observatory in the Northern hemisphere, has been observing cosmic rays with E >10 18 eV since 2008 [2]. The surface detector (SD) array of TA covering 700 km 2 area near the Delta city in Utah, USA, observes EASs with more than 95% live time [3]. The sky over the SD array is monitored by 38 fluorescence telescope detectors (FD) installed at 3 stations [4]. FDs image the EAS developments and measure the three-dimensional EAS reconstruction although their observation time is limited to dark and clear nights. The depth of shower maximum development in the atmosphere, X max , used in this paper is measured using the FD data. In addition, analyzing the events triggered by both the SD and FD detectors, called hybrid events, the accuracy of geometrical reconstruction, and hence the X max resolution, is improved. and hence the X max resolution, is improved.
In this paper, two studies by TA related to the hadronic interaction are reviewed. One is the determination of * e-mail: sako@icrr.u-tokyo.ac.jp proton-air inelastic and proton-proton total cross sections (σ inel p−air and σ tot p−p ) at √ s NN = 95 TeV [5], where √ s NN designates the center-of-momentum energy per nucleon. The other study is data-Monte Carlo (MC) comparison of the SD signal size for various shower geometries, which is related to the muon density in EASs [6]. The two studies are reviewed in Sec.2 and Sec.3, and a summary is given in Sec.4.

Method using FDs
The proton-proton cross-section, σ tot p−p , is one of the most fundamental quantities in particle physics. However, except for the unitarity limit, the evolution of the crosssection for hadronic interactions is not well understood. Cross-section measurements using FDs was pioneered by the Fly's Eye experiment in 1984 [7], and was updated by later experiments [8] [9].
The basic idea is to determine σ inel p−air using FD data assuming the exponential slope of X max distribution, labeled as Λ m , deep in the atmosphere is proportional to the mean free path of proton-air interaction, λ p−air . By selecting deep penetrating showers within a certain energy range, the corresponding primary particles are dominated by protons. The proportional coefficient, K, that relates Λ m =Kλ p−air is determined using air shower MC simulations. Finally, σ inel p−air is calculated with a relation σ inel p−air = 14.45m p /λ p−air , where 14.45 and m p are the mean mass number of air and the mass of the proton (nucleon), respectively. In principle, the K factor is interaction model dependent, however, as seen in Sec.2.2, this dependency is small. Furthermore, once Λ m is determined by an experimental group, λ p−air can be updated using K with latest knowledge on the hadronic interaction. Historical improvement in the measurement of K is shown in Fig.1. It is proposed in [5] to scale σ inel p−air =530±66 mb reported by Fly's Eye group [7] down to 392±49 mb according to the K factor update from 1.6 to 1.2.
Derivation of σ tot p−p requires another theoretical step to correct for the nuclear effect. Glauber calculation gives σ inel p−air considering the multiple nucleon effect inside a nucleus with inputs of σ tot p−p , nucleon distribution function and the elastic slope parameter B [12] [13]. Fixing the nucleon distribution function and introducing another theoretical assumption of a σ tot p−p -B relation, σ tot p−p is determined using σ inel p−air obtained above. More details on the methodology are found in [10] [13].

TA result on σ inel p−air
In the TA data analysis, data obtained from the SD array and the Middle Drum station, one of the 3 FD stations, from May 2008 to May 2013 were used. Using the hybrid reconstruction, the estimated X max resolution is about 23 g/cm 2 . A total of 438 events in the energy range from 10 18.3 eV to 10 19.3 eV remained after event selection. The mean energy of all events,10 18.68 eV, corresponds to the center of momentum energy of 95 TeV that is the highest energy where the cross section is determined so far. Note that because the recent composition study of TA concluded a light mass composition at this energy [14], proton dominance of deep X max events is a reasonable assumption.
The energy and interaction model dependences of K factor were studied using the CONEX one-dimensional air shower simulation code [15]. Validity of using the 1-D simulation code was confirmed by comparing with calculations using fully three-dimensional CORSIKA code, TA detector response simulation and TA reconstruction code. Fig.2 shows the K factor determined by the CONEX simulation as a function of primary energy. QGSJET II-04 [16] was used as the high-energy hadronic interaction model. No energy dependence is found and K is treated as a constant in the energy range of interest. The K factors were derived using four different high-energy hadronic interaction models QGSJET-01 [17], SIBYLL [18] and EPOS-LHC [19] in addition to QGSJET II-04, to be 1.22, 1.18, 1.19 and 1.15, respectively. Because the uncertainty of each value is only ±0.01, model dependency dominates the uncertainty of K determination, but it is still only ±3% from the mean value.
The distribution of X max measured by TA is shown in Fig.3. Deep penetrating events in the range from 790 g/cm 2 and 1000 g/cm 2 were fit with exponential function. The best fit function is drawn in the figure and its slope Λ m was found to be (50.47±6.26 [stat]) g/cm 2 . Using the K factor and Λ m obtained above with a relation σ inel p−air = 14.45m p K/Λ m , was obtained for the proton-air cross section at √ s NN = 95 TeV. The dominant elements of the systematic uncertainty accounted here are the interaction model dependence (±3% corresponding to ±17 mb), 20% Helium contamination (-18 mb) and 1% photon contamination (+23 mb). Detail of these estimations is found in [5]. σ inel p−air result obtained from the TA data is plotted in Fig.4 with the previous results together with the prediction lines of 4 interaction models. 1

TA result on σ tot p−p
The function of σ inel p−air (σ tot p−p , B) based on the Glauber calculation is shown in Fig.5. Three contour levels of σ inel p−air , a solid curve corresponding to 567 mb obtained in Sec.2.2 and two dotted curves showing the upper and lower uncertainty bands, are drawn. The shaded area is not allowed due to the unitarity constraint. A dashed straight line shows the theoretical prediction of the σ tot p−p -B relation using the QCD-inspired Block, Halzen, and Stanev (BHS) fit to the pp andpp data from the Tevatron [21].  Taking the intersections between the BHS prediction line and three σ inel p−air curves and projecting them on the horizontal axis, σ tot p−p and its error are obtained to be The result is plotted together with the previous works in Fig.6.

Comparison to recent LHC results
This section discusses the relationship to the latest LHC results after the TA publication. The analysis in Sec.2.3 relied on the theoretical BHS fit to relate σ tot p−p and B based on the Tevatron data at √ s=1.8 TeV. The TOTEM collaboration at LHC measured B = (20.36±0.19) GeV −2 at √ s=13 TeV and this is clearly larger than a smooth extrapolation from the data below the Tevatron energy [22]. The deviation does not exist in the √ s=2.76 TeV TOTEM result, but becomes clear in the √ s=7 and 8 TeV included. This includes the statistical (outer/thinner error bar) and the systematic (inner/thicker error bar).

V. PROTON-PROTON CROSS SECTION
From the TA proton-air cross section result we can determine the total proton-proton cross section. The process of inferring σ p−p from σ inel p−air is described in details in [35], and [36].
The σ p−p is calculated from the measured cross section, also known as the inelastic cross section σ inel p−air , using both Glauber Formalism [37] and the relation: Where σ total p−air is the total cross section, σ el p−air is the elastic cross section and σ qel p−air is the quasielastic cross section. The quasielastic cross section corresponds to scattering processes in which nuclear excitation occurs without particle production.
The relation between the σ inel p−air and the σ p−p is highly dependent on the forward scattering elastic slope B.
This is shown in the B, σ total p−p plane in Fig. 8. Here the solid and dotted curves represent a constant value of σ inel p−air that reflects the Telescope Array measured value and the statistical fluctuations.
There have been many theories predicting the relationship between B and σ p−p . However many of these models either failed to describe the elastic scattering data, or the elastic slope energy dependence from the Tevatron ([35,38,39]). A more updated theory using the single pomeron exchange model while describing the Tevatron data correctly is not consistent with the Unitarity constraint ([35,40]). Here the unitarity constraint is shown by solid grey shaded area in Fig. 8. A more recent prediction is the BHS fit [5]. It is consistent with unitarity while using a QCD inspired fit to the pp andpp data from the Tevatron. The dashed line in Fig. 8 shows the BHS prediction. Here     √ s=14 TeV in their paper before LHC [21] that is closer to the TOTEM result than the smooth extrapolation. A BHS prediction of σ tot p−p = (107.9±1.2) mb is also in a good agreement with the TOTEM measurement σ tot p−p = (110.6±3.4) mb. These good agreements with the latest (even unexpected) LHC results validate the usage of the BHS fit.

Muon measurements 3.1 Muon measurement by SD
The number density of muons in an EAS carries important information about mass of the primary particle. However there is clear evidence that the number of muons observed in the high-energy EAS is higher than the prediction of MC simulations. An exhaustive summary of this muon problem is found in the report of Hadron Interaction Working Group from this conference [1]. So far the latest hadronic interaction models reasonably agree with the accelerator data and no clear explanation is found to explain the muon problem. The collection of observational information over a wide range of parameter space is needed. Because the TA SD is composed of thin plastic scintillators [3], the signal size from muon passage is essentially identical to that of electromagnetic particles, and the latter dominate the signal due to a higher flux at the ground level. Because of the same (or similar) response to different type of charged particles, TA SDs do not have high sensitivity of particle identification at the detector level. The fraction of muon density in EAS is a function of geometrical condition because the electromagnetic components attenuate faster than muons after shower maximum and at distance far from the core. In this study, the TA collaboration performed a data-MC comparison for different geometrical condition of each SD as shown in Fig.7. The zenith angle (θ), distance from the shower axis (R) and azimuthal angle of the SD with respect to the primary arrival direction (ϕ) are used as geometrical parameters. When the signal size is calculated using MC simulation, the fraction of muon signal after detector simulation, muon purity P, is also recorded. Higher P is expected for larger R, larger θ and larger ϕ. The data-MC comparison was finally performed as a function of P.

TA dataset and MC
In this analysis, TA SD data obtained in 7 years from May 2008 to May 2015 were used. The energy range was limited between 10 18.8 eV and 10 19.2 eV to avoid possible energy dependence. The energy scale used in the data analysis is the one calibrated using the FD designated as E FD , while there is another energy scale E S D determined by SD alone. A relation E FD = E S D /1.27 between two scales is used in the TA analyses.
MC simulations were performed using CORSIKA 6.960 with QGSJET II-03 and FLUKA2008.3C for highand low-energy hadronic interactions, respectively, and EGS4 for electromagnetic interaction. A technique of thinning [23] and dethinning [24] was used to accelerate calculation speed. The detector response was calculated using GEANT4 and the final output was formatted the same as the experimental data so that the output was processed in the identical pipeline as the experimental data. To extract P, the detector simulation was individually performed for different type of particles at the detector. It should be noted that in the MC data analysis, energy scaling from E S D to E FD was NOT applied because this correction factor may also arise from the muon problem.
The primary protons were distributed between 10 16.55 eV and 10 20.55 eV following the energy spectrum by the HiRes experiment [25]. An isotropic distribution up to θ=60 • was also assumed. In addition, 0.05 muons/SD/±32µs were randomly added to account for the accidental background. Simulations of iron primary and with different hadronic interaction models were also performed and are discussed in Sec.3.3.
The signal size in units of vertical equivalent muons (VEM) in the SD area of 3 m 2 after analysis process as a function of R is shown in Fig.8 (top). Events with reconstructed geometries 30 • < θ <45 • , 150 • < |ϕ| <180 • and 10 18.8 eV <E S D <10 19.2 eV are selected. Results are plotted for the different type of particles together with the sum of all particles. The muon purity P, fraction of muon signal size to total, defined from this result is shown in Fig.8 (bottom). An increase of P as a function of R is found as expected and it rises up to 0.7 when R=3 km.

Results of TA data analysis
A data and MC comparison performed under the same conditions as Fig.8 is shown in Fig.9. The systematic uncertainty of the experimental data, ±(22-24)%, is dominated by the energy determination uncertainty using the FD. A clear excess of the experimental data to the MC prediction is found. The N data /N MC ratio shown in Fig.9 (bottom) shows an increasing trend with R. This suggests the cause of excess is due to muons because P also increases at larger R. A drop of the ratio at the furtherest bin, >4 km, is because at this distance the signal is dominated by background.
The same analysis was repeated using different hadronic interaction models and iron primaries in the MC. The model dependent ratios are shown in Fig.10. Increasing ratio with R is found for all hadronic interaction models. A comparison between proton and iron primaries using QGSJET II-03 is shown in Fig.11. Because the heavier nuclear showers contain more muons than lighter ones, the MC simulation with iron primary results a smaller ratio than the proton primary. However, the ratio is still >1.
Finally, the N data /N MC ratios were calculated for various geometrical conditions in θ and ϕ and plotted as a function of muon purity P in Fig.12. In this analysis R was fixed in the range between 2 km and 4 km where on average high P is expected as shown in Fig.8 (bottom). A positive correlation is found between the ratio and R and this also suggests that the excess is due to muons.
Summarizing all results it was concluded that part of the discrepancy in signal sizes between the data and the MC is due to a muon excess in the data. This result is qualitatively consistent with the reports by the Auger group  ( [26] for example) although the threshold energy of muon detection and R range used in the analysis are different. The muon excess found in the TA data will also explain the cause of the discrepancy between E S D and E FD . Though the muon purity, P, at the core distance of 800 m, where the energy scale parameter is defined, is 20%, the muon excess may partly cause misunderstanding of the conversion relation.

Summary
Two studies by TA related to the hadronic interactions are reviewed in this paper. The proton-air inelastic cross section σ inel p−air and the proton-proton total cross section σ tot p−p at √ s NN = 95 TeV were determined using the 5 years of TA hybrid data with one of the 3 FD stations. Our results are consistent with the extrapolation of the previous results and various theoretical predictions, but constrain them at the highest energy point. The assumed relation between σ tot p−p and the elastic slope B is investigated with the recent results at LHC. It is found that the BHS theory used in our analysis successfully predicted the LHC results of σ tot p−p and B. This reinforces the results of the TA analysis. Now TA has 10 years data and using 3 FD stations the statistics can be increased by a factor of 5.7. Analysis using this full dataset is ongoing and a significant improvement in the statistical uncertainty, that dominates the current results, is assured. Further increase of statistics and hence extension to higher energy are also anticipated using the coming TA×4 data [27].
A new analysis technique to investigate the muon contents in EASs was developed. Instead of extracting number of muons in the observed data, signal sizes of SD were compared between data and MC. According to the geometrical conditions of SD, the expected muon purity P was defined according to MC simulation. The TA SD data of 7 years were analyzed in a limited energy range from 10 18.8 eV to 10 19.2 eV. When the signal was compared as a function of R, an excess was observed in the experimental data and it increased with R following the increase of P. The same trend was observed using different hadronic interaction models and assuming iron primary. Though the amount of excess decreases with iron primary, it is still larger than unity. Results are divided in different θ and ϕ ranges with a large R range of 2-4 km. The data/MC ratio shows a correlation with P. All these measurements support an excess of muons in the experimental data compared to the MC predictions. This fact is also anticipated to explain the possible SD, FD discrepancy observed by TA.