Measurement of the Proton-Air Cross Section with Telescope Array's Black Rock Mesa and Long Ridge Fluorescence Detectors, and Surface Array in Hybrid Mode
R.U. Abbasi, M. Abe, T. Abu-Zayyad, M. Allen, R. Azuma, E. Barcikowski, J.W. Belz, D.R. Bergman, S.A. Blake, R. Cady, B.G. Cheon, J. Chiba, M. Chikawa, A. di Matteo, T. Fujii, K. Fujisue, K. Fujita, R. Fujiwara, M. Fukushima, G. Furlich, W. Hanlon, M. Hayashi, N. Hayashida, K. Hibino, R. Higuchi, K. Honda, D. Ikeda, T. Inadomi, N. Inoue, T. Ishii, R. Ishimori, H. Ito, D. Ivanov, H. Iwakura, H.M. Jeong, S. Jeong, C.C.H. Jui, K. Kadota, F. Kakimoto, O. Kalashev, K. Kasahara, S. Kasami, H. Kawai, S. Kawakami, S. Kawana, K. Kawata, E. Kido, H.B. Kim, J.H. Kim, J.H. Kim, M.H. Kim, S.W. Kim, S. Kishigami, V. Kuzmin, M. Kuznetsov, Y.J. Kwon, K.H. Lee, B. Lubsandorzhiev, J.P. Lundquist, K. Machida, H. Matsumiya, T. Matsuyama, J.N. Matthews, R. Mayta, M. Minamino, K. Mukai, I. Myers, S. Nagataki, K. Nakai, R. Nakamura, T. Nakamura, Y. Nakamura, Y. Nakamura, T. Nonaka, H. Oda, S. Ogio, M. Ohnishi, H. Ohoka, Y. Oku, T. Okuda, Y. Omura, M. Ono, R. Onogi, A. Oshima, S. Ozawa, I.H. Park, M.S. Pshirkov, J. Remington, D.C. Rodriguez, G. Rubtsov, D. Ryu, H. Sagawa, R. Sahara, Y. Saito, N. Sakaki, T. Sako, N. Sakurai, K. Sano, T. Seki, K. Sekino, et al. (45 additional authors not shown)
MMeasurement of the Proton-Air Cross Section with Telescope Array’s Black RockMesa and Long Ridge Fluorescence Detectors, and Surface Array in Hybrid Mode
R.U. Abbasi, ∗ M. Abe, T. Abu-Zayyad, M. Allen, R. Azuma, E. Barcikowski, J.W. Belz, D.R. Bergman, S.A. Blake, R. Cady, B.G. Cheon, J. Chiba, M. Chikawa, A. di Matteo, † T. Fujii, K. Fujisue, K. Fujita, R. Fujiwara, M. Fukushima,
7, 11
G. Furlich, W. Hanlon, ‡ M. Hayashi, N. Hayashida, K. Hibino, R.Higuchi, K. Honda, D. Ikeda, T. Inadomi, N. Inoue, T. Ishii, R. Ishimori, H. Ito, D. Ivanov, H.Iwakura, H.M. Jeong, S. Jeong, C.C.H. Jui, K. Kadota, F. Kakimoto, O. Kalashev, K. Kasahara, S.Kasami, H. Kawai, S. Kawakami, S. Kawana, K. Kawata, E. Kido, H.B. Kim, J.H. Kim, J.H. Kim, M.H. Kim, S.W. Kim, S. Kishigami, V. Kuzmin, § M. Kuznetsov,
20, 8
Y.J. Kwon, K.H. Lee, B.Lubsandorzhiev, J.P. Lundquist, K. Machida, H. Matsumiya, T. Matsuyama, J.N. Matthews, R. Mayta, M. Minamino, K. Mukai, I. Myers, S. Nagataki, K. Nakai, R. Nakamura, T. Nakamura, Y. Nakamura, Y. Nakamura, T. Nonaka, H. Oda, S. Ogio,
10, 26
M. Ohnishi, H. Ohoka, Y. Oku, T. Okuda, Y. Omura, M. Ono, R. Onogi, A. Oshima, S. Ozawa, I.H. Park, M.S. Pshirkov,
20, 29
J. Remington, D.C. Rodriguez, G. Rubtsov, D. Ryu, H. Sagawa, R. Sahara, Y. Saito, N. Sakaki, T. Sako, N. Sakurai, K. Sano, T.Seki, K. Sekino, P.D. Shah, F. Shibata, T. Shibata, H. Shimodaira, B.K. Shin, H.S. Shin, J.D.Smith, P. Sokolsky, N. Sone, B.T. Stokes, T.A. Stroman, T. Suzawa, Y. Takagi, Y. Takahashi, M.Takamura, M. Takeda, R. Takeishi, A. Taketa, M. Takita, Y. Tameda, H. Tanaka, K. Tanaka, M. Tanaka, Y. Tanoue, S.B. Thomas, G.B. Thomson, P. Tinyakov,
20, 8
I. Tkachev, H. Tokuno, T.Tomida, S. Troitsky, Y. Tsunesada,
10, 26
Y. Uchihori, S. Udo, T. Uehama, F. Urban, T. Wong, K.Yada, M. Yamamoto, K. Yamazaki, J. Yang, K. Yashiro, M. Yosei, Y. Zhezher,
7, 20 and Z. Zundel (The Telescope Array Collaboraton) Department of Physics, Loyola University Chicago, Chicago, Illinois, USA The Graduate School of Science and Engineering, Saitama University, Saitama, Saitama, Japan High Energy Astrophysics Institute and Department of Physics and Astronomy, University of Utah, Salt Lake City, Utah, USA Graduate School of Science and Engineering, Tokyo Institute of Technology, Meguro, Tokyo, Japan Department of Physics and The Research Institute of Natural Science,Hanyang University, Seongdong-gu, Seoul, Korea Department of Physics, Tokyo University of Science, Noda, Chiba, Japan Institute for Cosmic Ray Research, University of Tokyo, Kashiwa, Chiba, Japan Service de Physique Thorique, Universit Libre de Bruxelles, Brussels, Belgium The Hakubi Center for Advanced Research and Graduate School of Science,Kyoto University, Kitashirakawa-Oiwakecho, Sakyo-ku, Kyoto, Japan Graduate School of Science, Osaka City University, Osaka, Osaka, Japan Kavli Institute for the Physics and Mathematics of the Universe (WPI),Todai Institutes for Advanced Study, University of Tokyo, Kashiwa, Chiba, Japan Information Engineering Graduate School of Science and Technology, Shinshu University, Nagano, Nagano, Japan Faculty of Engineering, Kanagawa University, Yokohama, Kanagawa, Japan Interdisciplinary Graduate School of Medicine and Engineering,University of Yamanashi, Kofu, Yamanashi, Japan Earthquake Research Institute, University of Tokyo, Bunkyo-ku, Tokyo, Japan Academic Assembly School of Science and Technology Institute of Engineering, Shinshu University, Nagano, Nagano, Japan Astrophysical Big Bang Laboratory, RIKEN, Wako, Saitama, Japan Department of Physics, Sungkyunkwan University, Jang-an-gu, Suwon, Korea Department of Physics, Tokyo City University, Setagaya-ku, Tokyo, Japan Institute for Nuclear Research of the Russian Academy of Sciences, Moscow, Russia Faculty of Systems Engineering and Science, Shibaura Institute of Technology, Minato-ku, Tokyo, Japan Department of Engineering Science, Faculty of Engineering,Osaka Electro-Communication University, Neyagawa-shi, Osaka, Japan Department of Physics, Chiba University, Chiba, Chiba, Japan Department of Physics, Yonsei University, Seodaemun-gu, Seoul, Korea Faculty of Science, Kochi University, Kochi, Kochi, Japan Nambu Yoichiro Institute of Theoretical and Experimental Physics, Osaka City University, Osaka, Osaka, Japan Department of Physical Sciences, Ritsumeikan University, Kusatsu, Shiga, Japan Quantum ICT Advanced Development Center, National Institute forInformation and Communications Technology , Koganei, Tokyo, Japan Sternberg Astronomical Institute, Moscow M.V. Lomonosov State University, Moscow, Russia Department of Physics, School of Natural Sciences,Ulsan National Institute of Science and Technology, UNIST-gil, Ulsan, Korea a r X i v : . [ a s t r o - ph . H E ] J un Graduate School of Information Sciences, Hiroshima City University, Hiroshima, Hiroshima, Japan Institute of Particle and Nuclear Studies, KEK, Tsukuba, Ibaraki, Japan Department of Research Planning and Promotion, Quantum Medical Science Directorate,National Institutes for Quantum and Radiological Science and Technology, Chiba, Chiba, Japan CEICO, Institute of Physics, Czech Academy of Sciences, Prague, Czech Republic Department of Physics and Institute for the Early Universe,Ewha Womans University, Seodaaemun-gu, Seoul, Korea
Ultra high energy cosmic rays provide the highest known energy source in the universe to mea-sure proton cross sections. Though conditions for collecting such data are less controlled than anaccelerator environment, current generation cosmic ray observatories have large enough exposuresto collect significant statistics for a reliable measurement for energies above what can be attained inthe lab. Cosmic ray measurements of cross section use atmospheric calorimetry to measure depthof air shower maximum ( X max ), which is related to the primary particle’s energy and mass. Thetail of the X max distribution is assumed to be dominated by showers generated by protons, al-lowing measurement of the inelastic proton-air cross section. In this work the proton-air inelasticcross section measurement, σ inelp − air , using data observed by Telescope Array’s Black Rock Mesa andLong Ridge fluorescence detectors and surface detector array in hybrid mode is presented. σ inelp − air is observed to be 520 . ± . +25 . − [Sys.] mb at √ s = 73 TeV. The total proton-protoncross section is subsequently inferred from Glauber formalism and is found to be σ totpp = 139 . +23 . − . [Stat.] +15 . − . [Sys.] mb. I. INTRODUCTION
Ultra high energy cosmic rays (UHECRs) offer aunique opportunity as testing grounds for physics be-yond the standard model, as they represent a class ofparticles in the energy frontier beyond what can be gen-erated in human-made accelerators. In addition to ques-tions concerning their astrophysical nature, such as lo-cation of sources, composition, acceleration mechanisms,and propagation modes, fundamental aspects regardingthe nature of matter can be investigated as well. In par-ticular, UHECRs provide a way to measure the protoninteraction cross section at energies beyond what canbe achieved in the lab to test standard model predic-tions of how the cross section evolves with energy beyondwhat is measured in accelerators. Whereas acceleratorsare highly controlled environments, specially designed tomaximize integrated luminosity, cosmic ray experimentsmust rely on natural accelerators in the universe whichcan not be tuned to deliver a desired luminosity. The onlychoice for UHECR detectors is to increase their apertureto collect more events given a fixed interval of collectiontime.UHECR detectors do not directly observe the primaryparticle of interest due to the extremely low flux of thespectrum ( ∼ − eV − m − sr − s − ) [1]. Instead,the primary particle enters the Earth’s atmosphere andquickly interacts with an air molecule generating an ex-tended air shower which generates copious amounts offluorescence light with secondary particles reaching theground. Telescope Array collects ∼ ∗ [email protected] † Currently at INFN, sezione di Torino, Turin, Italy ‡ [email protected] § Deceased with energies > √ s >
43 TeV) with the surfacedetector array [2], which runs continuously day and night(100% duty cycle), and ∼
700 events per year per eachmonocular fluorescence detector station [3], which onlyrun on clear, moonless nights ( ∼
10% efficiency), in thesame energy range.In an accelerator experiment cross sections are mea-sured through careful design and control of the sourceand target, using either colliding beams or fixed targetsetup. Cross section, σ , in colliding beams is determinedby understanding the acceptance of the detector andmeasuring the event rate for a given beam luminosity, L = σ − dN/dt . A cosmic ray measurement of cross sec-tion is more akin to a fixed target calorimeter, with thecosmic ray flux acting as the beam and the atmospherethe target material. In the case of UHECRs, the mea-surement of cross section is made in the lab frame andthe atmosphere can be treated as a fixed target since anincoming proton has Lorentz factor γ in excess of 10 for E ≥ E > ∼ √ s > ∼
43 TeV). As accelerator de-signs are improved over time, human-made acceleratorsare closing the energy gap between center of mass ener-gies than can be achieved in the lab and what is providedby nature.Ultra High Energy Cosmic Ray detectors have beenreporting on the proton-air cross section measurementbeyond the capability of particle accelerators since1984 [15–22]. This work presents the second TelescopeArray report on the proton-air cross section. The firstresult was reported in 2015 using the Middle Drum (MD)fluorescence detector and the surface detector in hybrid
14 15 16 17 18 19 20 21(eV) log ) - s r - s - m . ( e V . E J Tibet IIICASA-BLANCAAkenoTunka-25YakutskIceTopTelescope ArrayTunka-133HiRes-MIAHiRes2Fly's EyeAugerHiRes1
Knee Ankle
12 12.5 13 13.5 14 14.5 15(eV)) (log pp sEquivalent CoM Energy T eva t r on L HC ( T e V ) L HC ( T e V ) FIG. 1. The cosmic ray spectrum starting at the knee, ob-served by recent experiments. Cosmic ray energies are mea-sured in the lab frame (i.e., as a fixed target measurement),while the highest energy accelerator based measurements aredone using colliding beams, reported in the center of momen-tum (CoM) frame. Recent accelerator experiments such as theTevatron and LHC are closing the energy gap between humanbuilt accelerators and astrophysical accelerators. Data from[4–14]. mode [23]. In this paper, we are reporting on the inelas-tic proton-air cross section, at √ s = 73 TeV, using nearlynine years of data observed by Black Rock Mesa (BRM)and Long Ridge (LR) fluorescence detectors (FDs) andthe surface detector (SD) array in hybrid mode. Notethat the BRM and LR detectors used in this analysis, arecloser in distance than MD to the surface detector arrayas shown in Figure 2. This enables us to study the inelas-tic proton-air cross section with higher statistical powerfor lower energy events. The technique used to analysethese events is similar to that used in the first proton-aircross section report [23]. The statistical power, on otherhand, increased by a factor of four. Note that all thesystematic sources are revisited and updated in additionto using the most recent hadronic high energy models.The proton-proton cross section is also calculated inthis work using Glauber formalism [24] and BHS fit [25].The inelastic proton-air and the total proton-proton crosssection are compared to previous cosmic ray experimentalresults and to predictions from models. II. DETECTOR DESCRIPTION
Telescope Array (TA) is a cosmic ray observatory thatdeploys multiple types of detectors to record the pas-sage of extended air showers caused by ultra high energycosmic ray primaries as they impact in the Earth’s atmo-sphere. The primary way TA observes air showers is byusing surface detectors which detect the energy depositedby high energy particles as they pass through them orusing fluorescence detectors which observe the UV lightgenerated in the atmosphere as the shower particles in-teract and exchange energy with air molecules. SDs do
FIG. 2. The Telescope Array detector configuration. Thefilled squares are the 507 SD scintillators on a 1.2 km grid.The SD scintillators are enclosed by three fluorescent detec-tors shown in filled triangles together with their field of view insolid lines. The northernmost fluorescence detector is calledMiddle Drum while the southern fluorescence detectors arereferred to as Black Rock Mesa and Long Ridge. The filledcircle in the middle equally spaced from the three fluorescencedetectors is the Central Laser Facility used for atmosphericmonitoring and detector calibration. not measure the development of the air shower in thesky, while FDs do. The shower size, number of chargedparticles at depth X ( N ( X )), of an air shower can beparameterized using the Gaisser-Hillas function [26] N ( X ) = N max (cid:18) X − X X max − X (cid:19) X max − X λ exp (cid:18) X max − Xλ (cid:19) (1)where N is the number of particles at slant depth X . Theparameters N max , X max , X , and λ describe the showershape. N max is the maximum number of shower particlesand the slant depth at which this occurs is denoted by X max . λ and X correspond to the proton-air interactionlength and slant depth of first interaction, respectively.In practice, when fitting real shower profiles using theGaisser-Hillas function, λ and X are fixed parameters,while X max is observed by the FDs. To get an accu-rate measure of X max , a monocular FD measurementis not sufficient. To improve X max resolution, simulta-neous observation of a shower by multiple FD stationsmust be employed or simultaneous observation by a FDstation and the SD ground array. In the case of mul-tiple FD stations, the independently measured shower-detector planes provide a strong constraint on the showertrack, leading to greatly improved geometrical resolu-tion. Similarly, a shower observed by a single FD sta-tion, along with the arrival time and core location on theground provided by the SD array delivers the same ben-efit. Resolution on X max improves dramatically from 84and 52g/cm g/cm for showers with energies of 1 EeVto 100 EeV respectively [27] to ∼
20 g/cm for E > . Thecenter of the SD array is located at 39 ◦ (cid:48) (cid:48)(cid:48) N 112 ◦ (cid:48) (cid:48)(cid:48) W, 1370 masl. Three FD stations overlooking theSD array are located outside the array boundary. All FDstations are located ∼
21 km away from the center of theSD array. Middle Drum station is located on the northend of the array, Black Rock Mesa is located at the south-east border and Long Ridge is located at the southwestborder. This work uses FD data from the Black RockMesa and Long Ridge detectors. While the general oper-ation of all FDs is similar, the design and location of theMiddle Drum detector relative to the SD array borderresults in different low energy acceptance than the BlackRock Mesa and Long Ridge detectors. The descriptionof FD equipment used in this analysis that follows is forBlack Rock Mesa and Long Ridge.Each FD station is comprised of twelve telescopes con-sisting of a multi-segmented 6.8 m mirror, a 16x16 pho-tomultiplier tube (PMT) array camera, electronics todigitize PMT signals at 40 MHz, trigger on air showertrack candidates, and readout and communications witha remote DAQ which controls event readout and storageamong all of the electronics racks. The telescopes arearranged in a two ring configuration providing zenith an-gle coverage in two bands. Six telescopes are assignedto ring 1, observing 3 ◦ − ◦ in elevation angle and sixare assigned to ring 2 observing 17 ◦ − ◦ . Azimuthalcoverage of 108 ◦ is the same for both rings. On clear,moonless nights the FDs scan the skies for potential airshower events. Because of this constraint, FD collectionefficiency is about 10%, whereas properly operating SDshave 100% operating efficiency since SDs can operate inall weather conditions 24 hours a day. Each FD elec-tronics rack has a track finder module which implementstemporal-spatial pattern recognition algorithms to deter-mine if a track has been observed. If a set of tube triggersmeets the criteria an event level trigger is generated andcommunicated to the remote DAQ which forces readoutof all mirrors for storage and offline analysis. BRM andLR electronics utilize FADC electronics which allows dig-itization of PMT signals at an equivalent 14-bit, 10 MHzsampling rate, allowing observation of the time develop-ment of an event with 100 ns time resolution. Above10 . eV, showers are seen with distance of closest ap-proach (impact parameter) >
25 km. Further detailsabout the construction and design of the BR and LRstations can be found in [28, 29].Each surface detector is composed of two layers of 3 m plastic scintillator, 1.2 cm thick. Wavelength shiftingfibers are embedded in grooves in the scintillator lay-ers and optically coupled to a PMT (one for each layer).PMT signals are digitized by a 12 bit FADC operatingat 50 MHz sampling rate. Onboard electronics deployedwith each SD scan for signals above threshold ( > > . . and becomes 100% efficientwith no zenith angle dependence. Refer to [30] for fur-ther information regarding the technical details of TA’ssurface detector array.The second additional cut added to this analysis isa zenith angle cut. Because the atmosphere acts as acalorimeter, fluorescence detectors are limited in theirability to reconstruct showers based upon the atmo-spheric overburden available for air showers to develop.Fluorescence detectors far above sea level in general willaccept fewer showers for reconstruction than those closerto sea level, because the higher detector has less atmo-sphere for the shower to reach maximum size before hit-ting the ground. Thus, showers with small zenith angle(close to vertical) traverse less atmosphere than thosewith large zenith angle, reducing the chance that theshower will achieve maximum size before penetrating theground or falling below the field of view of the detector.Figure 4 shows Telescope Array hybrid X max acceptanceof QGSJET II.4 protons. As seen in the figure, eventswith zenith angle less than 30 degrees show a break inacceptance roughly corresponding to the vertical depthof ground level, indicated by the dashed line showing thevertical depth of the Central Laser Facility (CLF) at thecenter of the SD array. Events that have zenith anglegreater than 30 degrees, are sufficiently inclined to pro-vide enough slant depth to reach shower maximum inthe atmosphere. Indeed, these events show no significantbreak in the X max acceptance.The proton-air cross section measurement uses infor-mation from the deep tail of the X max distribution underthe assumption that only light primaries such as protons,and possibly some helium contamination, populate thisregion of the distribution. To ensure the highest level ofproton purity in the tail of the X max distribution usedto determine σ inelp − air , we search for the minimum zenithangle cut which results in nearly flat X max acceptance forall X max in the energy range 18 . ≤ log ( E/ ev) < . X max acceptance shows a break in the deep X max re-gion for some range of zenith angles, those events mustbe removed because showers induced by proton primariesmay be lost in the X max distribution tail.Analysis, data, and Monte Carlo used for this work isidentical to that used in [31] except for energy binningand the additional zenith angle cut described above. Theresultant data set contains 1975 events with a resolutionin X max of ∼
20 g/cm and an average energy of 10 . eV. For further details concerning the hybrid analysisprocedure refer to [31]. III. DATA TRIGGER, RECONSTRUCTION,AND SELECTION.
The FD and SD data streams are collected indepen-dent of each other. To create a hybrid data stream thestreams are searched for coincident triggers that occurwithin 500 µ s. For this set of hybrid events, SD recon-struction proceeds as described in [30] to determine theshower core location and arrival time. FD reconstruc-tion is performed as described in [3] to determine theshower-detector plane for each FD station that observesa shower. This determines the shower-detector planeangle, ψ , impact parameter, core location, and arrivaltime. A hybrid reconstruction takes the additional stepof casting the individual SDs into “pixels” that observethe shower in a similar way FD PMTs do. This allowsus to use them in the shower-detector plane fit. Becauseof their accurate measure of the shower track positionand arrival time on the ground, these points provide anadditional constraint on the track geometry. Once thehybrid shower geometry is determined, the shower pro-file is measured by each observing FD station using thisimproved measure of the shower track. Shower profile re-construction determines the shower size measured by thenumber of charged particles as a function of atmosphericdepth (equation 1) and proceeds as described in [3]. Theshower profile is used to determine the primary particleenergy and X max , both of which are the essential inputsto the proton-air cross section measurement.The data used for this analysis was collected from 2008May 27 to 2016 November 29, nearly nine years, and isthe same data used the BR/LR hybrid X max measure-ment in [31]. That analysis examined X max for eventswith E > = 10 . eV and resulted in 3330 events afterapplying all quality cuts to the data. The present anal-ysis imposes two more cuts on the data required for agood quality cross section measurement: here we restrictanalysis to events with energy 18 . ≤ log ( E/ eV) < . > ◦ . The rational for these additionalcuts is described below.Abbasi et al. [31] demonstrated that below 10 . eV,the TA RMS of the X max distribution σ ( X max ) is con-sistent with light composition ranging between 52 and63 g/cm . Above this energy σ ( X max ) begins to de-crease. Due to changing zenith angle acceptance andfalling statistics it is premature to say if this narrowing of σ ( X max ) is astrophysical in nature or caused by selec-tion bias. We can compare TA’s observed mean (cid:104) X max (cid:105) and RMS σ ( X max )of the X max distribution to MonteCarlo predictions by randomly sampling the X max dis-tributions of individual elements such as proton, helium,nitrogen, and iron according to data statistics. To dothis the simulated X max distributions of each of thoseelements is randomly sampled N times, where N is thenumber of events observed in the data for the given en-ergy bin, a distribution of X max is therefore generated,and (cid:104) X max (cid:105) and σ ( X max ) of the distribution is recorded.Note the that X max distributions are fully simulated withall acceptance effects present in reconstructed data. Thisprocedure is repeated 5000 times. We can then measurethe 68%, 90%, and 95% confidence intervals of the jointexpectation of (cid:104) X max (cid:105) and σ ( X max ) for each element asshown. We can then compare the predictions of (cid:104) X max (cid:105) and σ ( X max ) to what is observed in the data, which isshown in figure 3 for two energy bins. The figure alsoshows (cid:104) X max (cid:105) and σ ( X max ) observed by TA, as well asthe systematic and statistical uncertainties.Figure 3a shows (cid:104) X max (cid:105) and σ ( X max ) observed for18 . ≤ log ( E/ eV) < . (cid:104) X max (cid:105) and σ ( X max ) of the data are closest to the prediction ofQGSJET II.4 protons. Additionally, the predictions fromMonte Carlo simulation have relatively small dispersionand are easily distinguishable because of the relativelylarge statistics in this energy bin. The TA hybrid datais tested against these single element models and, in thisenergy bin, it is easy to see given TA’s statistical andsystematic errors, as well as the clear separation in simu-lated (cid:104) X max (cid:105) and σ ( X max ), that the best fit to the data iscompatible to only one element within systematic errors.The situation changes though where statistics are small,as shown in figure 3b. This energy bin shows (cid:104) X max (cid:105) and σ ( X max ) observed for 19 . ≤ log ( E/ eV) < . (cid:104) X max (cid:105) here has increased as pre-dicted due to the relationship between primary particleenergy and X max and observed σ ( X max ) has decreased.But the simulations predict much larger dispersion in the (cid:104) X max (cid:105) and σ ( X max ), causing the data to be indistin-guishable from light single element models such as pro-ton, all the way up to single element nitrogen. Above E ≥ . eV the observed data exhibits this effect tosuch a degree that we impose an additional cut elimi-nating events above this energy for the present analysisto ensure that the tail of the X max distribution used forthe p-air cross section measurement is not significantlycontaminated with heavy elements. We also estimate the
660 680 700 720 740 ) > (g/cm max 720 740 760 780 800 820 840) > (g/cm max 9. The star represents (cid:104) X max (cid:105) and σ ( X max ) observed by TA in the two energy bins,as well as the statistical and systematic uncertainties. EachMonte Carlo chemical element shows the 68.3% (blue ellipse),90% (orange ellipse), and 95% (red ellipse) confidence inter-vals. contamination in the tail due to helium and this estimateis described later in this section.The second additional cut added to this analysis isa zenith angle cut. Because the atmosphere acts as acalorimeter, fluorescence detectors are limited in theirability to reconstruct showers based upon the atmo-spheric overburden available for air showers to develop.Fluorescence detectors far above sea level in general willaccept fewer showers for reconstruction than those closerto sea level, because the higher detector has less atmo-sphere for the shower to reach maximum size before hit-ting the ground. Thus, showers with small zenith angle(close to vertical) traverse less atmosphere than thosewith large zenith angle, reducing the chance that theshower will achieve maximum size before penetrating theground or falling below the field of view of the detector.Figure 4 shows Telescope Array hybrid X max acceptanceof QGSJET II.4 protons. As seen in the figure, eventswith zenith angle less than 30 degrees show a break inacceptance roughly corresponding to the vertical depthof ground level, indicated by the dashed line showing thevertical depth of the Central Laser Facility (CLF) at thecenter of the SD array. Events that have zenith angle 700 750 800 850 900 950 1000) (g/cm max X - - · - · - · - · - · - E ff i c i e n cy (cid:176) < 30 q (cid:176) > 30 q T A C LF D e p t h FIG. 4. Telescope Array hybrid QGSJET II.4 proton X max acceptance in two zenith angle bands for energies 18 . ≤ log ( E/ eV) < . 0. The black shows the reconstructionefficiency of events with zenith angle < 30 degrees and redline is the efficiency with zenith angle > 30 degrees. Smallzenith angle events are more likely to achieve shower maxi-mum below the FD field of view and therefore the acceptancedrops roughly at the vertical depth of ground level at TA.The dashed line shows the vertical depth of the TA CentralLaser Facility (CLF). Inclined events with large zenith angleare more likely to reach shower maximum in the atmosphereand the acceptance does not significantly fall off with slantdepth. greater than 30 degrees, are sufficiently inclined to pro-vide enough slant depth to reach shower maximum inthe atmosphere. Indeed, these events show no significantbreak in the X max acceptance.The proton-air cross section measurement uses infor-mation from the deep tail of the X max distribution underthe assumption that only light primaries such as protons,and possibly some helium contamination, populate thisregion of the distribution. To ensure the highest level ofproton purity in the tail of the X max distribution usedto determine σ inelp − air , we search for the minimum zenithangle cut which results in nearly flat X max acceptance forall X max in the energy range 18 . ≤ log ( E/ ev) < . X max acceptance shows a break in the deep X max re-gion for some range of zenith angles, those events mustbe removed because showers induced by proton primariesmay be lost in the X max distribution tail.Analysis, data, and Monte Carlo used for this work isidentical to that used in [31] except for energy binningand the additional zenith angle cut described above. Theresultant data set contains 1975 events with a resolutionin X max of ∼ 20 g/cm and an average energy of 10 . eV. For further details concerning the hybrid analysisprocedure refer to [31]. IV. ANALYSIS The proton-air inelastic cross section σ is related tointeraction length λ (mean free path) by λ = 1 / ( nσ ) (2)where n is the target particle density. The probability ofan interaction in a slab of target material of thickness dx is P ( x ) = (1 /λ ) dx . Given a “beam” of cosmic rays, beamintensity, I , decreases with increasing number of slabstraversed as dI = − I/λ dx , leading to the expressionof beam intensity, I ( x ) = I exp( − x/λ ), where I is theinitial intensity and x is depth. Therefore for cosmicrays the interaction length can be measured by fitting adistribution of depth of first interaction ( X ) between thecosmic ray primary particle and an air nucleus to find λ .In practice this is not feasible because the starting pointof the upper atmosphere is not well defined due to itsvery low density and there is no appreciable fluorescencesignal generated at first interaction.After the initial inelastic collision an air shower con-tinues to grow in size through production of secondariesmainly by radiative processes of pair production andbremsstrahlung in the electromagnetic portion of theshower. The shower grows until it reaches a maximumsize, X max , dependent primarily on primary particle en-ergy and mass, then decreases as energy loss of secon-daries becomes dominated by collisional processes. X max therefore is a uniquely defined point in the shower pro-file that is observed by fluorescence detectors and canbe used as a proxy for X to determine the interactionlength of a distribution of cosmic ray primaries. Proton-air cross section is therefore measured indirectly for airshowers.In this work, the K-Factor method [23] is used to ob-tain the proton-air cross section. The tail of the X max distribution retains the exponentially falling nature of the X distribution encoded within it and can be parameter-ized as f ( X max ) = exp( − X max / Λ m ), where Λ m is theexponential slope of the tail. The K -Factor method re-lates the slope of the tail of the X max distribution and tothe slope of X distribution, through a constant factor K which is determined by MC simulation and is close tounity. Λ m = Kλ p − air (3)where we now label λ p − air as the proton-air interactionlength and both Λ m and λ p − air are measured in g/cm .Using the relationship between λ and σ expressed inequation 2, cross section can be related to the mean tar-get mass of air as σ inelp − air = (cid:104) m air (cid:105) λ p − air (4)and substituting this into equation 3 we findΛ m = K σ p − air = K . m p σ p − air (5)where (cid:104) m air (cid:105) = 24160 mb g cm − or 14 . m p with theproton mass expressed in g [32]. σ inelp − air is expressed inmb and Λ m is in g/cm . This equation directly linksthe observed X max distribution to the proton-air crosssection. This K -Factor method is the same method used inthe first TA report on the proton-air cross section in2015 [23]. The data analysis here is divided into twoparts. The first part is done by calculating the valueof the attenuation length (Λ m ) of the observed UHECRevents. In the second part, we calculate the inelasticproton-air cross section ( σ inelp − air ) value from the obtainedattenuation length Λ m . A. Measuring the Attenuation Length Λ m The value of attenuation length Λ m , and thereforethe proton-air cross section, can be calculated by fittingthe X max distribution tail to the exponential functionexp( − X max / Λ m ). Here only the tail of the X max dis-tribution is used to obtain Λ m , because it is the mostpenetrating part of the distribution and is assumed to becomposed mostly of protons. UHECR composition cannot be measured on an event by event basis and mustbe inferred from a distribution of events. By restrictingthe determination of Λ m to the tail of the X max distri-bution, potential contamination from heavier elements inthe primary spectrum is reduced.The choice of the starting point of the tail fit ( the loweredge of the fit range) X i for the exponential fit is madeby fitting the X max distribution tail to two exponentialfunctions with separate power indices. The break pointof these two fits ( found to be at 790 g/cm ) describe thebest fit beyond which the distribution can be describedusing a single exponential function. This maximizes thenumber of events in the tail distribution while minimizinginstability in the value of Λ m due to possible detector biasor helium contamination.Figure 5 shows the X max distribution of the data col-lected by the Telescope Array southern most fluorescencedetectors, Black Rock Mesa and Long Ridge, togetherwith the surface detector hybrid events. The distribu-tion includes 1975 events in the energy range between10 . − . eV with an average energy of 10 . eV.The data included here passed the quality cuts describedin the previous section and is fitted to the exponentialfunction exp( − X max / Λ m ) using the unbined maximumlikelihood method.Several systematic checks are applied to test for thestability of the measured attenuation length Λ m . This isdone by dividing the data in two halves based on: thezenith angle, the distance of the shower using the impactparameter R p , and the energy of the event. The dividedsubsets are found to be consistent within statistical fluc-tuations.The final Λ m measured by the Telescope Array detec-tor at an average energy of 10 . eV ( √ s = 73 TeV) in-cluding the statistical checks is found to be Λ m = 55.9 ± . Note that Λ m is directly derived fromthe data and is model independent. Therefore, it can beused at a later time to calculate the inelastic proton-aircross section independent of the method or the UHECR ) Xmax (g/cm600 650 700 750 800 850 900 950 1000 / g ) ( c m m a x X D N u m be r o f E v en t s / -1 FIG. 5. The number of events per X max bin (∆ X max ) vs. X max g/cm for BRM and LR fluorescence detectors and theTelescope Array surface detector in hybrid mode. The line isthe exponential fit to the slope using the unbined likelihoodmethod between 790-1000 g/cm . models used in this paper. B. Proton-Air cross Section Measurement To determine the interaction mean free path of pro-tons in air, λ p − air and therefore the inelastic proton-air cross section σ inelp − air we use the K -Factor technique.Using equation 3, K can be directly computed usingMonte Carlo. K depends on the hadronic model be-ing used in simulations. UHECR simulations rely onthe choice of electromagnetic interaction driver, low en-ergy hadronic generator, and high energy hadronic gen-erator [33]. The most popular high energy hadronicgenerators are SIBYLL2.3 [34, 35], EPOS-LHC [36],QGSJET II.4 [37, 38], and QGSJET01 [39] which, withthe exception of QGSJET01, are tuned to the most re-cent accelerator data at energies accessible to acceleratorsand extrapolated to UHECR energies through theoreti-cal and phenomenological predictions. Hadronic modeldependence is an important and difficult considerationwhen dealing with questions related to the fundamentalproperties of hadronic air showers such as proton-air crosssection or cosmic ray composition. Some of the impor-tant parameters that affect shower development whichare extrapolated from accelerator data to UHECR ener-gies are inelasticity, multiplicity, and cross section. Each E(eV) Log18.2 18.3 18.4 18.5 18.6 18.7 18.8 18.9 19 K - V a l ue / ndf c – c – FIG. 6. The value of K obtained vs. energy in Log (eV) forsimulated data sets using CONEX 6.4 with the high energymodel QGSJET II.4, for the energy range of the data, between10 . and 10 . eV. hadronic generator uses different methods to do this lead-ing to differences in shower development at ultra highenergies. For a summary of these issues refer to [40, 41].For this work we present the results for several differentmodels and report on the systematic uncertainty in theresults of our measurement. K is computed in this work by generating several sim-ulated sets between 10 . − eV for each of the highenergy models. Each generated set contains ten thou-sand events using a one-dimensional air shower MonteCarlo program CONEX 6.4 [42–44]. Figure 6 shows the K -value including the statistical fluctuation calculatedfor each of these simulated sets, using QGSJET II.4 asan example. The value of K is then obtained by fittingthe K -value vs. energy to a horizontal line as shown inFigure 6. It is important to note that the value of Λ m and therefore, the value of K is dependent on the choiceof the lower edge of the tail fit range X i (as shown inFigure 7 ). A consistent procedure needs to be used todetermine X i and therefore the value of K for each en-ergy bin and the high energy model shower simulations.To do so we calculate the difference in slant depth D be-tween the peak of the X max distribution and 790 g/cm ,using a simulated data set at an energy of 10 . eV(equivalent to the mean energy of the data set used inthis work). The same difference in slant depth D is laterused to consistently determine the value of X i from peakof the X max distribution for each of the simulated setsfor each of the high energy models.To confirm the validity of the obtained K values, foreach of the generated data sets, for each of the high en-ergy models, λ p − air is reconstructed and compared to the λ p − air provided by the corresponding high energy model.Figure 8 shows the comparison of the the values of thehigh energy model λ p − air and the obtained λ p − air using Lower Edge of the Fit Range (Xi)700 750 800 850 900 K (E(eV)) 18.2 Log (E(eV)) 18.4 Log (E(eV)) 18.7 Log FIG. 7. The value of K vs. the lower edge in the fit range( X i ) to the tail of the X max distribution for several data sets10 . , . , and 10 . eV simulated using CONEX withthe high energy model QGSJET II.4. Each data set contains10,000 simulated events. the K -Factor technique. Figure 8 shows that the valueof K obtained in this study indeed describes the value of K of the high energy models correctly. E(eV) Log18.2 18.3 18.4 18.5 18.6 18.7 18.8 18.9 19 ) ( g / c m p - a i r l Reconstructed LambdaModel Lambda FIG. 8. The proton-air interaction length λ p − air in g/cm vs. Energy in Log (eV) for the simulated data sets usingCONEX with the high energy model QGSJET II.4, for theenergy range of the data, between 10 . and 10 . eV. Thecircle points are the λ p − air values obtained from the X dis-tribution. Triangle points are the ones determined from re-constructing the λ p − air values using the K -Factor method. The K value is dependent on the high energy modelused. The obtained K value is shown using CONEX6.4 in Table I, together with the corresponding inelasticproton-air cross section σ inelp − air . Note that, the σ inelp − air iscalculated using equation 5 with Λ m obtained from theTA X max distribution and K tabulated for each the highenergy models QGSJET II.4 [37, 38], QGSJET01 [39] ,SIBYLL2.3 [34, 35], and EPOS-LHC [36].Each K listed in Table I is the average value of K over the energy range of 10 . -10 . eV. The value of K ismeasured to be ∼ 20% larger than 1.0 meaning the slopeof the tail of the X max distribution falls more slowly thanthe X tail. This is because the X max distribution resem-bles a convolution of a falling exponential, from the con-tribution of X , and a Gaussian from the growth of theshower and fluctuations of stochastic processes of showerdevelopment [45]. Showers exhibit large intrinsic fluctua-tions in development even for those initiated by particlesof the same mass and energy. If showers did not exhibitthese fluctuations, air shower X max distributions wouldresemble the distribution of X , just shifted to a greaterdepth in the atmosphere by a constant amount. It isimportant to note that the K value model dependenceshown in Table I is on the order of ∼ ± 3% ( K -Valuehistorical improvement is discussed in [23]). This makesthe K -Value method weakly model dependent and thusa reliable method to use in calculating the σ inelp − air . Model K σ inelp − air (mb)QGSJET II.4 1 . ± . 01 505 . ± . . ± . 01 514 . ± . . ± . 01 535 . ± . . ± . 01 527 . ± . K obtained for each of the high en-ergy models and the corresponding inelastic proton-air crosssection for that model. Each K listed is the single averagevalue of K over the energy range of 10 . -10 . eV. Notethat the values of K shows a ∼ ± 3% model uncertainty. In order to quantify the systematic uncertainties in the σ inelp − air measurement, several check were applied. First,systematic uncertainty due to model dependence is re-ported. This is done by quantifying the maximum varia-tion in the σ inelp − air value by each model from the average σ inelp − air obtained from all of the high energy models. Thisuncertainty was found to be equal to ± 15 mb.In addition, the systematic effect of possible energy de-pendent bias in the X max distribution was studied. Thiswas done by shifting the values of X max by their elon-gation rate (the rate of change of shower (cid:104) X max (cid:105) w.r.t.shower energy) prior to fitting. The systematic effectfrom a possible energy bias was found to be negligible.The systematic effects due to detector bias is alsotested. This systematic effect is done by comparing theattenuation length calculated with and without detectoreffects. First, the attenuation length is calculated fromsimulations, where the detector effect is not included.After which, Λ m is calculated from simulation, wherethe events are propagated through the detector, recon-structed, and the quality cuts applied. The value of Λ m was found to be consistent, for all the high energy models,between the thrown events and the reconstructed eventswith quality cuts applied. Therefore, the systematic ef-fect from this test was found to be negligible. For furtherdetails concerning this test procedure refer to [23].0Another systematic check is done by studying the im-pact of contamination from other primaries. The sys-tematic effect of other elements in the tail beside protonincluding photon, CNO, helium and iron is investigated.Only photons and helium introduce a bias in the inelasticproton-air cross section.The upper limit of cosmic-ray photon fraction at theenergy range in this study is found to be ∼ . -10 . eV is measured to notlarger than 43.8% at the 95% c.l. Using this limit, thesystematic uncertainty due to helium contamination isfound to be − 40 mb.Note here that the sign for the systematic uncertaintydue to helium and gamma contamination is negative andpositive respectively. Helium has a larger cross-sectionthan protons. Therefore, helium contamination will re-sult in the observation of a larger cross section than wouldbe the case with pure protons. The opposite occurs dueto gamma contamination. The final systematic uncer-tainty for the σ inelp − air is calculated by adding each of thesystematic uncertainties quadratically.The final proton-air cross section measured by theTelescope Array detector at an average energy of 10 . eV using the K -Factor method and including the statisti-cal and systematic checks is σ inelp − air = 520 . ± . +25 − [Sys.] mb. This result is shown in Figure 9 and iscompared to other experimental measurements [15–23]and current high energy model predictions. Note herethat the current proton-air cross section result includ-ing the error fluctuations is consistent with the high en-ergy models tuned to the LHC (QGSJET II.4 [37, 38],SIBYLL2.3 [34, 35], and EPOS-LHC [36]) shown in Fig-ure 9. C. Proton-Proton Cross Section The analysis to convert from the inelastic proton-aircross section to proton-proton cross section consists oftwo parts.The first part is done by converting the measured in-elastic proton-air cross section to the possible allowedvalues of the proton-proton cross. The conversion isobtained using the Glauber formalism [24] which gives σ inelp − air as a function of σ totpp and B , where B is the for-ward scattering elastic slope. The three curved lines inFigure 10 show the TA measurement of σ inelp − air and its sta-tistical uncertainties allowed region in the ( σ totpp - B ) plane.The second part is done by constraining the relationbetween σ totpp and B using a theoretical model. Themodel used in this work is (Block, Halzen, and Stanev(BHS)) [25], shown as the dashed line in Figure 10. Theintersection of the σ inelp − air allowed region and the theoret-ical constraint (BHS model) gives us σ totpp and B values. FIG. 9. The proton-air cross section result of this work,including the statistical (thin) and systematics (thick) errorbars, in comparison to previous experimental results [15–23]. In addition, the high energy models (QGSJET II.4,QGSJET01, SIBYLL 2.3, EPOS-LHC) cross section predic-tions are shown. Note that the BHS model can be replaced with othermodels or predictions to solve for the σ totpp . Note theBHS model is consistent with the unitarity constraintwhile describing the pp and ¯pp cross section data fromthe Tevatron well [47, 48].The proton-proton cross section in this work is foundto be σ totpp = 139 . +23 . − . [Stat.] +15 . − . [Sys.] mb. This resultis shown in Figure 11 in comparison to previously re-ported values by UHECR experiments [16, 17, 19, 21, 23].The recent result from LHC by TOTEM at √ s = 7and 13 TeV [49, 50] is also shown, in addition to theBHS fit [25]. The best fit of the proton-proton totalcross section data by the COMPETE collaboration is alsoadded [51]. For further details concerning the proton-airto proton-proton procedure refer to [23]. V. SUMMARY AND CONCLUSION Telescope Array has measured the inelastic proton-aircross section of ultra high energy cosmic rays at √ s =73 TeV. This measurement is performed for energies thatare not accessible to accelerator experiments, thereforeprovides an important and unique test of standard modelpredictions about the fundamental nature of matter.The Telescope Array utilizes a large array of surfacedetectors and fluorescence telescopes to record the atmo-spheric depth of maximum size of air showers initiated by1 [mb] p-ptot s 40 60 80 100 120 140 160 180 200 220 ] - B [ ( G e V / c ) FIG. 10. The elastic slope B in ((GeV/c) − ) vs. σ totalp − p inmb. The solid line is the allowed σ totalp − p values from the σ inelp − air and its statistical errors reported in this work. The dashedline is the BHS fit prediction [52]. While the gray shaded areais the unitarity constraint. inelastic collisions of ultra high energy cosmic rays andair molecules in the upper atmosphere. By combining thegeometric and timing information of SDs and the BlackRock Mesa and Long Ridge FDs that observe a hybridevent X max can be determined with a good precision of ∼ 20 g/cm . UHECR X max distributions are related tothe interaction length of cosmic rays in the atmosphere,which in turn depends on the tail of X max distributionsis populated with the deepest penetrating events, pre-dominantly proton initiated events, the slope of which isrelated to the interaction length by a constant, K . UsingMonte Carlo simulations K can be evaluated using MonteCarlo provides access to the depth of first interaction and X max for each event, allowing a direct determination of K . Once K is known the inelastic proton-air cross sectioncan be determined using equation 5. Using nearly nineyears of hybrid data, TA measures σ inelp − air = 520 . ± . +25 − [Sys.] mb for √ s = 73 TeV. Using Glaubertheory and the Block, Halzen, Stanev model The totalproton-proton cross section is determined from σ inelp − air tobe σ totpp = 139 . +23 . − . [Stat.] +15 . − . [Sys.] mb.It is interesting to note that ultra high energy cosmicray model prediction of the proton-air cross section haveconverged closer than was the case prior to tuning toLHC data. This is shown in the K value convergingfrom 7% down to 3%. Most importantly, this is also (GeV)s ( m b ) p - p s Fly's EyeAkenoHiResAugerTelescope Array (2015)Telescope Array (This work)pbar-pppCOMPETEppTOTEM even (QCD-Fit) nn s FIG. 11. A compilation of the proton-proton cross sectionvs. the center of mass energy result of this work, includingthe statistical (thin) and systematics (thick) error bars, inaddition to previous work by cosmic rays detectors [16, 17,19, 21, 23] in addition to, the recent result from LHC byTOTEM at √ s = 7 and 13 TeV [49, 50]. The dashed redcurve is the BHS fit [25] and the dashed black curve is the fitby the COMPETE collaboration [53]. This plot is adaptedand modified from [25]. found to be consistent with results for ultra high energycosmic ray experiments including this work. The datafrom the high energy models and ultra high energy cosmicray experiments continue to show a rising cross sectionwith energy.Future cross section results, using TA × − eV with high statisti-cal power and at several energy intervals. This wouldallow us to make a statement on the functional form ofthe cross section energy dependence. VI. ACKNOWLEDGEMENTS The Telescope Array experiment is supported bythe Japan Society for the Promotion of Science(JSPS)through Grants-in-Aid for Priority Area 431, forSpecially Promoted Research JP21000002, for Sci-entific Research (S) JP19104006, for Specially Pro-moted Research JP15H05693, for Scientific Research2(S) JP15H05741, for Science Research (A) JP18H03705,for Young Scientists (A) JPH26707011, and for Fos-tering Joint International Research (B) JP19KK0074,by the joint research program of the Institute forCosmic Ray Research (ICRR), The University ofTokyo; by the U.S. National Science Foundationawards PHY-0601915, PHY-1404495, PHY-1404502,and PHY-1607727; by the National Research Founda-tion of Korea (2016R1A2B4014967, 2016R1A5A1013277,2017K1A4A3015188, 2017R1A2A1A05071429) ; by theRussian Academy of Sciences, RFBR grant 20-02-00625a(INR), IISN project No. 4.4502.13, and Belgian SciencePolicy under IUAP VII/37 (ULB). The foundations ofDr. Ezekiel R. and Edna Wattis Dumke, Willard L. Ec-cles, and George S. and Dolores Dor´e Eccles all helpedwith generous donations. The State of Utah supportedthe project through its Economic Development Board,and the University of Utah through the Office of the Vice President for Research. The experimental site becameavailable through the cooperation of the Utah School andInstitutional Trust Lands Administration (SITLA), U.S.Bureau of Land Management (BLM), and the U.S. AirForce. We appreciate the assistance of the State of Utahand Fillmore offices of the BLM in crafting the Plan ofDevelopment for the site. Patrick Shea assisted the col-laboration with valuable advice on a variety of topics.The people and the officials of Millard County, Utah havebeen a source of steadfast and warm support for our workwhich we greatly appreciate. We are indebted to the Mil-lard County Road Department for their efforts to main-tain and clear the roads which get us to our sites. Wegratefully acknowledge the contribution from the tech-nical staffs of our home institutions. An allocation ofcomputer time from the Center for High PerformanceComputing at the University of Utah is gratefully ac-knowledged. [1] D. Ivanov, in , International Cosmic Ray Conference,Vol. 36 (2019) p. 298.[2] T. Abu-Zayyad et al. (Telescope Array), Astrophys. J.Lett. , L1 (2013), arXiv:1205.5067 [astro-ph.HE].[3] T. Abu-Zayyad et al. (Telescope Array), Astropart. Phys. , 16 (2013), arXiv:1305.6079 [astro-ph.HE].[4] M. G. Aartsen et al. (IceCube), Phys. Rev. D88 , 042004(2013), arXiv:1307.3795 [astro-ph.HE].[5] R. U. Abbasi et al. (HiRes), Phys. Rev. Lett. , 101101(2008), arXiv:astro-ph/0703099 [astro-ph].[6] T. Abu-Zayyad et al. (HiRes-MIA), Astrophys. J. ,686 (2001), arXiv:astro-ph/0010652 [astro-ph].[7] M. Amenomori, Adv. Space Res. , 467 (2008),arXiv:0803.1005 [astro-ph].[8] D. J. Bird et al. (HiRes), Astrophys. J. , 491 (1994).[9] F. Fenu (Pierre Auger), PoS ICRC2017 , 486 (2018).[10] J. W. Fowler, L. F. Fortson, C. C. H. Jui, D. B. Kieda,R. A. Ong, C. L. Pryke, and P. Sommers, Astropart.Phys. , 49 (2001), arXiv:astro-ph/0003190 [astro-ph].[11] D. Ivanov, Proceedings, 34th International Cosmic RayConference (ICRC 2015): The Hague, The Netherlands,July 30-August 6, 2015 , PoS ICRC2015 , 349 (2016).[12] S. Knurenko, I. Petrov, Z. Petrov, and I. Sleptsov, Pro-ceedings, 18th International Symposium on Very HighEnergy Cosmic Ray Interactions (ISVHECRI 2014):Geneva, Switzerland, August 18-22, 2014 , EPJ WebConf. , 04001 (2015).[13] M. Nagano, M. Teshima, Y. Matsubara, H. Y. Dai,T. Hara, N. Hayashida, M. Honda, H. Ohoka, andS. Yoshida, J. Phys. G18 , 423 (1992).[14] V. V. Prosin et al. , Proceedings, 18th International Sym-posium on Very High Energy Cosmic Ray Interactions(ISVHECRI 2014): Geneva, Switzerland, August 18-22,2014 , EPJ Web Conf. , 04002 (2015).[15] F. Siohan, R. Ellsworth, A. Ito, J. Macfall, R. Streitmat-ter, et al. , J.Phys. G4 , 1169 (1978).[16] R. Baltrusaitis, G. Cassiday, J. Elbert, P. Gerhardy,S. Ko, et al. , Phys.Rev.Lett. , 1380 (1984). [17] M. Honda, M. Nagano, S. Tonwar, K. Kasahara, T. Hara, et al. , Phys.Rev.Lett. , 525 (1993).[18] H. Mielke, M. Foeller, J. Engler, and J. Knapp, J.Phys. G20 , 637 (1994).[19] K. Belov (HiRes Collaboration), Nucl.Phys.Proc.Suppl. , 197 (2006).[20] G. Aielli et al. (ARGO-YBJ Collaboration), Phys.Rev. D80 , 092004 (2009), arXiv:0904.4198 [hep-ex].[21] P. Abreu et al. (Pierre Auger Collaboration),Phys.Rev.Lett. , 062002 (2012), arXiv:1208.1520[hep-ex].[22] Aglietta, M. and others (EAS-TOP Collaboration), Phys.Rev. D , 032004 (2009).[23] R. U. Abbasi et al. (Telescope Array), Phys. Rev. D92 ,032007 (2015), arXiv:1505.01860 [astro-ph.HE].[24] R. Glauber and G. Matthiae, Nucl.Phys. B21 , 135(1970).[25] M. Block and F. Halzen, Phys.Rev. D72 , 036006 (2005),arXiv:hep-ph/0506031 [hep-ph].[26] T. K. Gaisser and A. M. Hillas, International Cosmic RayConference , 353 (1977).[27] T. Z. AbuZayyad, The Energy Spectrum of Ultra HighEnergy Cosmic Rays , Ph.D. thesis, Utah U. (2000).[28] H. Tokuno, Y. Tameda, M. Takeda, K. Kadota,D. Ikeda, et al. , Nucl.Instrum.Meth. A676 , 54 (2012),arXiv:1201.0002 [astro-ph.IM].[29] Y. Tameda et al. , Nucl. Instrum. Meth. A609 , 227(2009).[30] T. Abu-Zayyad et al. (Telescope Array), Nucl. Instrum.Meth. A689 , 87 (2013), arXiv:1201.4964 [astro-ph.IM].[31] R. U. Abbasi et al. (Telescope Array), Astrophys. J. ,76 (2018), arXiv:1801.09784 [astro-ph.HE].[32] R. Ulrich, J. Blumer, R. Engel, F. Schussler,and M. Unger, New J. Phys. , 065018 (2009),arXiv:0903.0404 [astro-ph.HE].[33] D. Heck, G. Schatz, T. Thouw, J. Knapp, and J. Capde-vielle, (1998).[34] R. Fletcher, T. Gaisser, P. Lipari, and T. Stanev,Phys.Rev. D50 , 5710 (1994). [35] E.-J. Ahn, R. Engel, T. K. Gaisser, P. Lipari,and T. Stanev, Phys. Rev. D , 094003 (2009),arXiv:0906.4113 [hep-ph].[36] T. Pierog, I. Karpenko, J. M. Katzy, E. Yatsenko,and K. Werner, Phys. Rev. C92 , 034906 (2015),arXiv:1306.0121 [hep-ph].[37] S. Ostapchenko, Colliders to cosmic rays. Proceedings,2nd International Conference, C2CR07, Granlibakken,Lake Tahoe, USA, February 25-March 1, 2007 , AIP Conf.Proc. , 118 (2007), arXiv:0706.3784 [hep-ph].[38] S. Ostapchenko, Phys. Rev. D83 , 014018 (2011),arXiv:1010.1869 [hep-ph].[39] N. Kalmykov, S. Ostapchenko, and A. Pavlov,Nucl.Phys.Proc.Suppl. , 17 (1997).[40] R. Engel, D. Heck, and T. Pierog, Ann. Rev. Nucl. Part.Sci. , 467 (2011).[41] R. Ulrich, R. Engel, and M. Unger, Phys.Rev. D83 ,054026 (2011), arXiv:1010.4310 [hep-ph].[42] T. Bergmann, R. Engel, D. Heck, N. Kalmykov,S. Ostapchenko, et al. , Astropart.Phys. , 420 (2007),arXiv:astro-ph/0606564 [astro-ph].[43] T. Pierog et al. , Nucl. Phys. Proc. Suppl. , 159 (2006),astro-ph/0411260.[44] G. Bossard et al. , Phys. Rev. D63 , 054030 (2001), hep-ph/0009119.[45] C. J. Todero Peixoto, V. de Souza, and J. A. Bellido,Astropart. Phys. , 18 (2013), arXiv:1301.5555 [astro-ph.HE]. [46] A. Glushkov, I. Makarov, M. Pravdin, I. Sleptsov,D. Gorbunov, G. Rubtsov, and S. Troitsky, Phys. Rev.D , 041101 (2010), arXiv:0907.0374 [astro-ph.HE].[47] J. Dias De Deus, Nucl.Phys. B59 , 231 (1973).[48] A. Buras and J. Dias de Deus, Nucl.Phys. B71 , 481(1974).[49] G. Antchev et al. , EPL , 21002 (2011), arXiv:1110.1395[hep-ex].[50] G. Antchev et al. (TOTEM), Eur. Phys. J. C79 , 103(2019), arXiv:1712.06153 [hep-ex].[51] J. R. Cudell, V. V. Ezhela, P. Gauron, K. Kang,Yu. V. Kuyanov, S. B. Lugovsky, E. Martynov, B. Nico-lescu, E. A. Razuvaev, and N. P. Tkachenko (COM-PETE), Phys. Rev. Lett. , 201801 (2002), arXiv:hep-ph/0206172 [hep-ph].[52] M. Block, F. Halzen, and T. Stanev, Phys.Rev. D62 ,077501 (2000), arXiv:hep-ph/0004232 [hep-ph].[53] J. R. Cudell, V. V. Ezhela, P. Gauron, K. Kang, Y. V.Kuyanov, S. B. Lugovsky, E. Martynov, B. Nicolescu,E. A. Razuvaev, and N. P. Tkachenko (COMPETE Col-laboration), Phys. Rev. Lett. , 201801 (2002).[54] E. Kido, in European Physical Journal Web of Confer-ences , European Physical Journal Web of Conferences,Vol. 210 (2019) p. 06001.[55] R. U. Abbasi et al. (Telescope Array), Astrophys. J.865