Measurement of Ultra-High Energy Cosmic Rays with CHICOS Citation Brobeck, Elina (2009) Measurement of Ultra-High Energy Cosmic Rays with CHICOS. Dissertation (Ph.D.), California Institute of Technology. doi:10.7907/9YH6-N064. https://resolver.caltech.edu/CaltechETD:etd-10192008-143041 Abstract The California HIgh school Cosmic ray ObServatory (CHICOS) is a ground-based scintillator array designed to measure the extended air showers of ultra-high energy cosmic rays. The goal of the project is to gain insight into the origin of ultra-high energy cosmic rays by measuring the energy spectrum and the distribution of arrival directions. The CHICOS array has been in operation since 2003. It consists of 77 pairs of scintillator dectectors deployed at schools in the San Fernando and San Gabriel valleys near Los Angeles, and is designed to observe cosmic ray air showers at energies of 10^18 eV and above. In addition, the Chiquita subarray is designed to observe smaller showers in the energy range of 10 16 - 10 19 eV. We present new descriptions of the air shower lateral distribution function and time distribution function, which have been derived from AIRES-generated simulated air showers. The new functions are specific to the CHICOS altitude and allow for a maximum likelihood shower reconstruction method, which is more appropriate to the CHICOS data than the χ 2 minimization method. We present several analyses of the accuracy of the reconstruction software in the energy ranges available to the Chiquita and CHICOS arrays. The energy spectrum between 10 17 eV and 10 19 eV has been measured by the Chiquita subarray. At the lowest energy range, it is found to agree with previous measurements, while the measured flux falls below previous experiments for energies greater than approximately 10 17.5 eV. The CHICOS energy spectrum above 10 18.4 eV is found to agree with previous results published by AGASA. However, we do not observe the cutoff in the spectrum at 10 20 eV reported more recently by the Auger and HiRes Collaborations. A correlation analysis between CHICOS data and nearby active galactic nuclei (AGN) was performed. No excess of cosmic rays was observed in the vicinity of nearby AGN. The maximum correlation was observed for cosmic ray events with E > 10 20 eV and for AGN with z < 0.009, with P chance = 21%. This is consistent with random correlations from an isotropic distribution, a result also found by HiRes, but in disagreement with Auger. Item Type: Thesis (Dissertation (Ph.D.)) Subject Keywords: cosmic rays; UHECRs Degree Grantor: California Institute of Technology Division: Physics, Mathematics and Astronomy Major Option: Astrophysics Thesis Availability: Public (worldwide access) Research Advisor(s): McKeown, Robert D. Group: Astronomy Department Thesis Committee: Stone, Edward C. (chair) Harrison, Fiona A. Kamionkowski, Marc P. McKeown, Robert D. Readhead, Anthony C. S. Defense Date: 18 September 2008 Non-Caltech Author Email: elina (AT) alum.mit.edu Funders: Funding Agency Grant Number NSF PHY-0244899 NSF PHY-0102502 Weingart Foundation UNSPECIFIED Caltech UNSPECIFIED Record Number: CaltechETD:etd-10192008-143041 Persistent URL: https://resolver.caltech.edu/CaltechETD:etd-10192008-143041 DOI: 10.7907/9YH6-N064 Default Usage Policy: No commercial reproduction, distribution, display or performance rights in this work are provided. ID Code: 4187 Collection: CaltechTHESIS Deposited By: Imported from ETD-db Deposited On: 27 Apr 2009 Last Modified: 05 Nov 2021 22:56 Thesis Files Preview PDF (elina_brobeck_thesis.pdf) - Final Version See Usage Policy. 6MB Repository Staff Only: item control page

Measurement of Ultra-High Energy Cosmic Rays with CHICOS Thesis by Elina Brobeck In Partial Fulfillment of the Requirements for the Degree of Doctor of Philosophy California Institute of Technology Pasadena, California 2009 (Defended September 18, 2008) ii c 2009 ° Elina Brobeck All Rights Reserved iii Acknowledgements The CHICOS project is directed by Robert McKeown (California Institute of Technology). Principal collaborators include Ryoichi Seki (California State University, Northridge), G. B. Yodh (University of California, Irvine), Jim Hill (California State University, Dominguez Hills), and John Sepikas (Pasadena City College). Project Coordinators Michelle Larson and Theresa Lynn oversaw the deployment of sites in the CHICOS array and coordinated among hundreds of teachers and school administrators to keep the array running. The data-acquisition software was designed and implemented by Sandy Horton-Smith and Juncai Gao. Chao Zhang contributed critical parts of the original data-filtering software and a set of diagnostic tools for determining the health of sites in the field. The CHICOS project owes much to Brant Carlson and Chris Jillings for the original development of the shower-reconstruction software. We thank Pat Huber, who has performed the essential task of maintaining the CHICOS server and has also overseen the deployment and maintenance of the CHICOS computers in the field. Many thanks are also given to Bob Gates, who has been crucial in running the CHICOS summer program for high school students. CHICOS has benefitted from the work of summer students Stephen Ho, Keith Chan, Derek Wells, Shawn Ligocki, and Veronica Anderson in developing a CHICOSspecific lateral distribution function and time distribution function. We also wish to acknowledge students Angela Marotta, for her work on the temperature variation of CHICOS data, Clare Kasper, for her analysis of the detector energy spectra, Amanda McAuley, who developed the tools to construct a sky map of the CHICOS cosmic ray data, and Jay Conrod, who designed an interactive graphical user interface to the iv shower reconstruction software. Philipp Bonushkin and Eric Black performed the measurement of the effective area of the CHICOS detectors. We gratefully acknowledge the many high school and middle school teachers and administrators who volunteered their time and energy to help with the deployment and continued operation of the CHICOS array. We thank the many high school students who chose to participate in the CHICOS summer sessions and aided in the commissioning of CHICOS detectors. The work in this paper makes use of unthinned AIRES air shower simulations contributed by Barbara Falkowski. The analysis presented here has benefitted greatly from discussion and feedback offered by Chris Jillings, Theresa Lynn, Robert McKeown, and other CHICOS collaborators. The CHICOS project is grateful to Los Alamos National Laboratory for the donation of scintillator detectors, and IBM for the donation of computer equipment. We acknowledge financial support from the National Science Foundation (grants PHY0244899 and PHY-0102502), the Weingart Foundation, and the California Institute of Technology. v Abstract The California HIgh school Cosmic ray ObServatory (CHICOS) is a ground-based scintillator array designed to measure the extended air showers of ultra-high energy cosmic rays. The goal of the project is to gain insight into the origin of ultra-high energy cosmic rays by measuring the energy spectrum and the distribution of arrival directions. The CHICOS array has been in operation since 2003. It consists of 77 pairs of scintillator dectectors deployed at schools in the San Fernando and San Gabriel valleys near Los Angeles, and is designed to observe cosmic ray air showers at energies of 1018 eV and above. In addition, the Chiquita subarray is designed to observe smaller showers in the energy range of 1016 –1019 eV. We present new descriptions of the air shower lateral distribution function and time distribution function, which have been derived from AIRES-generated simulated air showers. The new functions are specific to the CHICOS altitude and allow for a maximum likelihood shower reconstruction method, which is more appropriate to the CHICOS data than the χ2 minimization method. We present several analyses of the accuracy of the reconstruction software in the energy ranges available to the Chiquita and CHICOS arrays. The energy spectrum between 1017 eV and 1019 eV has been measured by the Chiquita subarray. At the lowest energy range, it is found to agree with previous measurements, while the measured flux falls below previous experiments for energies greater than approximately 1017.5 eV. The CHICOS energy spectrum above 1018.4 eV is found to agree with previous results published by AGASA. However, we do not observe the cutoff in the spectrum at 1020 eV reported more recently by the Auger vi and HiRes Collaborations. A correlation analysis between CHICOS data and nearby active galactic nuclei (AGN) was performed. No excess of cosmic rays was observed in the vicinity of nearby AGN. The maximum correlation was observed for cosmic ray events with E > 1020 eV and for AGN with z < 0.009, with Pchance = 21%. This is consistent with random correlations from an isotropic distribution, a result also found by HiRes, but in disagreement with Auger. vii Contents Acknowledgements iii Abstract v List of Figures x List of Tables xiii 1 Introduction 1 1.1 Historical Background . . . . . . . . . . . . . . . . . . . . . . . . . . 1 1.2 The CHICOS Experiment . . . . . . . . . . . . . . . . . . . . . . . . 5 2 Origins of UHECRs 7 2.1 Diffusive Shock Acceleration . . . . . . . . . . . . . . . . . . . . . . . 7 2.2 The GZK Cutoff . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 9 2.3 Potential UHECR Sources . . . . . . . . . . . . . . . . . . . . . . . . 13 2.3.1 Radio Galaxies and AGN . . . . . . . . . . . . . . . . . . . . 15 2.3.2 Neutron Stars and Magnetars . . . . . . . . . . . . . . . . . . 19 2.3.3 Quasar Remnants . . . . . . . . . . . . . . . . . . . . . . . . . 20 2.3.4 Starburst Galaxies and LIGs . . . . . . . . . . . . . . . . . . . 22 2.3.5 Gamma Ray Bursts . . . . . . . . . . . . . . . . . . . . . . . . 23 2.3.6 Top-Down Models 24 . . . . . . . . . . . . . . . . . . . . . . . . 3 Cosmic Ray Air Showers 3.1 Air Shower Development . . . . . . . . . . . . . . . . . . . . . . . . . 25 25 viii 3.2 3.3 Lateral Distribution Function . . . . . . . . . . . . . . . . . . . . . . 30 3.2.1 AGASA LDF . . . . . . . . . . . . . . . . . . . . . . . . . . . 32 3.2.2 CHICOS LDF . . . . . . . . . . . . . . . . . . . . . . . . . . . 33 Time Distribution Function . . . . . . . . . . . . . . . . . . . . . . . 40 3.3.1 AGASA TDF . . . . . . . . . . . . . . . . . . . . . . . . . . . 40 3.3.2 CHICOS TDF . . . . . . . . . . . . . . . . . . . . . . . . . . . 41 4 The CHICOS Experiment 49 4.1 The Detector Array . . . . . . . . . . . . . . . . . . . . . . . . . . . . 49 4.2 Data-Collection Software . . . . . . . . . . . . . . . . . . . . . . . . . 58 5 Shower Reconstruction Software 5.1 5.2 5.3 62 Overview of libCTShower . . . . . . . . . . . . . . . . . . . . . . . . 62 5.1.1 Data Format . . . . . . . . . . . . . . . . . . . . . . . . . . . 63 5.1.2 Code . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 64 Shower Reconstruction . . . . . . . . . . . . . . . . . . . . . . . . . . 65 5.2.1 Chi-Square Reconstruction Method . . . . . . . . . . . . . . . 68 5.2.2 Maximum Likelihood Reconstruction Method . . . . . . . . . 70 User Interface . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 74 6 Modeling the Array Response 79 6.1 Simulation of Air Showers Using AIRES . . . . . . . . . . . . . . . . 79 6.2 Modeling the Detector Response . . . . . . . . . . . . . . . . . . . . . 81 6.2.1 Energy Deposited in the Scintillator . . . . . . . . . . . . . . . 81 6.2.2 Time-Over-Threshold Measurement . . . . . . . . . . . . . . . 82 6.2.3 Timing Accuracy . . . . . . . . . . . . . . . . . . . . . . . . . 83 Modeling of Photon Interactions . . . . . . . . . . . . . . . . . . . . . 84 6.3.1 Compton Scattering . . . . . . . . . . . . . . . . . . . . . . . 85 6.3.2 Pair Production . . . . . . . . . . . . . . . . . . . . . . . . . . 88 Analysis of Simulated Shower Reconstructions . . . . . . . . . . . . . 95 6.3 6.4 ix 7 Analysis of Low-Energy Data 103 7.1 Low-Energy Data . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 103 7.2 Simulation of Low-Energy Showers . . . . . . . . . . . . . . . . . . . 105 7.3 Estimating the Low-Energy Flux . . . . . . . . . . . . . . . . . . . . 109 8 Analysis of High-Energy Data 116 8.1 Results from Previous Experiments . . . . . . . . . . . . . . . . . . . 116 8.2 High-Energy Data from CHICOS . . . . . . . . . . . . . . . . . . . . 117 8.3 Simulation of High-Energy Showers . . . . . . . . . . . . . . . . . . . 119 8.4 Estimating the High-Energy Flux . . . . . . . . . . . . . . . . . . . . 124 9 Correlation of UHECR Data with AGN 128 9.1 Quantifying the Degree of Correlation . . . . . . . . . . . . . . . . . . 128 9.2 Recent Results from Other Experiments . . . . . . . . . . . . . . . . 130 9.3 Results from CHICOS Data . . . . . . . . . . . . . . . . . . . . . . . 131 10 Conclusions 137 A Site Locations and Parameters 139 B CHICOS Showers Above 1018.4 eV 144 Bibliography 151 x List of Figures 3.1 Diagram of an air shower . . . . . . . . . . . . . . . . . . . . . . . . . 26 3.2 Profile of an air shower . . . . . . . . . . . . . . . . . . . . . . . . . . . 27 3.3 Comparison of CORSIKA and AIRES simulations . . . . . . . . . . . . 30 3.4 Low-energy iron-primary electron LDF . . . . . . . . . . . . . . . . . . 36 3.5 Low-energy iron-primary muon LDF . . . . . . . . . . . . . . . . . . . 37 3.6 High-energy proton-primary electron LDF . . . . . . . . . . . . . . . . 38 3.7 High-energy proton-primary muon LDF . . . . . . . . . . . . . . . . . 39 3.8 High-energy proton-primary electron TDF . . . . . . . . . . . . . . . . 45 3.9 High-energy proton-primary muon TDF . . . . . . . . . . . . . . . . . 46 3.10 Energy invariance of the electron TDF . . . . . . . . . . . . . . . . . . 47 3.11 Energy invariance of the muon TDF . . . . . . . . . . . . . . . . . . . 48 4.1 Sites in the CHICOS array as of July 2005 . . . . . . . . . . . . . . . . 50 4.2 Array size and reporting statistics, 2003–2005 . . . . . . . . . . . . . . 51 4.3 Array size and reporting statistics, 2006–2007 . . . . . . . . . . . . . . 52 4.4 Diagram of a CHICOS array site . . . . . . . . . . . . . . . . . . . . . 53 4.5 Diagram of CHICOS detectors . . . . . . . . . . . . . . . . . . . . . . 53 4.6 CEU schematic . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 54 4.7 Sample decay-constant calibration . . . . . . . . . . . . . . . . . . . . 55 4.8 Sample pulse-energy histogram . . . . . . . . . . . . . . . . . . . . . . 56 4.9 Labview data acquisition program front panel . . . . . . . . . . . . . . 60 4.10 Labview data acquisition program history panel . . . . . . . . . . . . . 60 4.11 Labview data acquisition program satellites panel . . . . . . . . . . . . 61 4.12 Labview data acquisition program energy panel . . . . . . . . . . . . . 61 xi 5.1 Sample shower data . . . . . . . . . . . . . . . . . . . . . . . . . . . . 64 5.2 Transform function f (xang , yang ) used to map (θ, φ) to linear coordinates 67 5.3 Tangent of the vertical angle θ as a function of xang and yang . . . . . . 67 5.4 Sample output of shower reconstructor on raw shower data . . . . . . . 76 5.5 Sample output of shower reconstructor after removal of accidental hits 77 5.6 Sample output of shower reconstructor after removal of accidental hits and removal of inactive sites from the array . . . . . . . . . . . . . . . 78 6.1 Simulated pulse measurement . . . . . . . . . . . . . . . . . . . . . . . 84 6.2 Cross section for Compton scattering . . . . . . . . . . . . . . . . . . . 89 6.3 Distribution of recoil electron energies from Compton scattering . . . . 89 6.4 Cross section for pair production in hydrogen . . . . . . . . . . . . . . 91 6.5 Cross section for pair production in carbon . . . . . . . . . . . . . . . 91 6.6 Pair-produced electron energies . . . . . . . . . . . . . . . . . . . . . . 94 6.7 Distribution of reconstructed energies for simulated showers at 1017 eV 97 6.8 Distribution of reconstructed energies sorted by shower inclination . . . 98 6.9 Distribution of reconstructed angles sorted by shower inclination . . . . 99 6.10 Input core locations of simulated showers . . . . . . . . . . . . . . . . . 100 6.11 Reconstructed core locations of simulated showers . . . . . . . . . . . . 100 6.12 Example of simulated shower data . . . . . . . . . . . . . . . . . . . . 101 6.13 Example reconstruction of simulated shower data . . . . . . . . . . . . 102 7.1 Map of the Chiquita detector sites . . . . . . . . . . . . . . . . . . . . 104 7.2 Energy distribution of showers simulated for Chiquita . . . . . . . . . . 106 7.3 Reconstructed energies of low-energy simulations . . . . . . . . . . . . 107 7.4 Reconstructed core locations of low-energy simulations . . . . . . . . . 107 7.5 Reconstructed vertical angles of low-energy simulations . . . . . . . . . 108 7.6 Reconstructed azimuthal angles of low-energy simulations . . . . . . . 108 7.7 Acceptance of the Chiquita array . . . . . . . . . . . . . . . . . . . . . 109 7.8 Energy distribution of Chiquita showers between 1016 eV and 1019 eV . 110 7.9 Flux times E 3 , as measured by the Chiquita array . . . . . . . . . . . 114 xii 8.1 Core locations of simulated high-energy showers . . . . . . . . . . . . . 120 8.2 Energy distribution of showers simulated for CHICOS . . . . . . . . . 121 8.3 Reconstructed energies of high-energy simulations . . . . . . . . . . . . 121 8.4 Reconstructed core locations of high-energy simulations . . . . . . . . . 122 8.5 Reconstructed core locations of high-energy simulations . . . . . . . . . 122 8.6 Reconstructed vertical angles of high-energy simulations . . . . . . . . 123 8.7 Reconstructed azimuthal angles of high-energy simulations . . . . . . . 123 8.8 Acceptance of the CHICOS array . . . . . . . . . . . . . . . . . . . . . 125 8.9 Relative sky exposure of the CHICOS experiment . . . . . . . . . . . . 125 8.10 Energy distribution of CHICOS showers between 1018.4 eV and 1020.4 eV 126 8.11 Flux times E 3 , as measured by the CHICOS array . . . . . . . . . . . 126 9.1 Exposure-weighted fraction of the sky covered by windows of angular radius θmax centered on nearby AGN . . . . . . . . . . . . . . . . . . . 132 9.2 Cumulative binomial probability P of the observed correlation resulting from a random isotropic distribution . . . . . . . . . . . . . . . . . . . 134 9.3 The set of CHICOS data points used in the correlation search . . . . . 135 9.4 The set of 11 CHICOS data points with E > 100 EeV and the set of 220 AGN with z < 0.009 for which there is a maximum correlation . . 136 xiii List of Tables 4.1 GPS settings for M12 receivers . . . . . . . . . . . . . . . . . . . . . . 57 4.2 GPS settings for UT+ receivers . . . . . . . . . . . . . . . . . . . . . . 57 6.1 Parameters of the AIRES showers used in this work . . . . . . . . . . . 80 6.2 Fraction of showers detected at 1017 eV . . . . . . . . . . . . . . . . . . 95 A.1 Database of site locations, part I . . . . . . . . . . . . . . . . . . . . . 140 A.2 Database of site locations, part II . . . . . . . . . . . . . . . . . . . . . 141 A.3 Database of site parameters, part I . . . . . . . . . . . . . . . . . . . . 142 A.4 Database of site parameters, part II B.1 High-energy showers observed in 2003 . . . . . . . . . . . . . . . . . . 145 B.2 High-energy showers observed in 2004 . . . . . . . . . . . . . . . . . . 146 B.3 High-energy showers observed in 2005, part I . . . . . . . . . . . . . . 147 B.4 High-energy showers observed in 2005, part II . . . . . . . . . . . . . . 148 B.5 High-energy showers observed in 2006 . . . . . . . . . . . . . . . . . . 149 B.6 High-energy showers observed in 2007 . . . . . . . . . . . . . . . . . . 150 . . . . . . . . . . . . . . . . . . . 143 1 Chapter 1 Introduction The energy spectrum of cosmic rays spans over 10 orders of magnitude, from 109 eV to 1020 eV, and perhaps beyond. Over the past 50 years, ground array experiments have increasingly been able to probe the highest energy range of the spectrum, yet much remains uncertain. At the highest energies, the existence of a flux cutoff and the origin of the cosmic ray particles have yet to be determined. The CHICOS experiment has been designed to provide data in this ultra-high energy range of the cosmic ray spectrum. 1.1 Historical Background The phenomenon now known as cosmic radiation was first recognized as having an extraterrestrial origin following the experiments of Victor Hess in 1912 [1]. Hess, intending to show that the pervasive ionizing radiation in the atmosphere emanated from the Earth, measured the intensity change with altitude from a hot air balloon. Counter to expectations, he found that the radiation intensity increased with altitude and must therefore be arriving at the Earth from space. For this discovery, Hess was awarded the Nobel Prize in 1936. Among the explanations put forward to explain cosmic radiation was Millikan’s hypothesis that it was neutral gamma radiation emitted in the process of protons and electrons coming together in space to form atoms [2]. In 1933, however, Arthur Compton carried out a worldwide survey that showed the intensity of cosmic rays 2 varied with latitude, concluding that the radiation must consist of charged particles whose paths were deflected by the Earth’s magnetic field [3]. Assuming the particles were electrons, Compton calculated their energy to be approximately 7 × 109 eV. The energy spectrum was soon to be expanded even further. In 1938, Pierre Auger observed that some particles, separated by as much as 20 m, arrived in time coincidence [4, 5]. This phenomenon was simultaneously discovered by Werner Kohlhörster, working separately in Germany [6]. Additional experiments with more widely spaced counters showed that coincidence events could be observed as far as 200 m apart [7]. Under the assumption that particles arriving in coincidence derived from a single primary source, Auger estimated the energy of the primary cosmic rays to be 1015 eV. This was the beginning of the study of the particle cascades known as extended cosmic ray air showers. The details of our current understanding of air showers will be discussed in chapter 3. The systematic measurement of ultra-high energy cosmic rays (UHECRs) via sampling of extended air showers was first implemented in 1954 at the Agassiz station of the Harvard College Observatory [8]. This experiment was the first to use plastic scintillator detectors to simultaneously measure particle densities and arrival times (from which the arrival direction of the shower can be derived). The array of 15 detectors was operational between 1954 and 1957, and extended the known energy spectrum to above 1018 eV by the observation of a shower with more than 109 particles. The Harvard array was used as a prototype for the larger Volcano Ranch array built in New Mexico. In 1962, Volcano Ranch measured the first particle with an estimated energy greater than 1020 eV [9]. Following the discovery of the cosmic microwave background in 1965, separate theoretical analsyses by Greisen [10] and Zatsepin and Kuz’min [11] predicted a sharp decline, now known as the GZK cutoff, in the cosmic ray spectrum near this energy due to photopion production. It was noted by Greisen that given the total exposure of UHECR experiments at that time, the observation of even one particle above 1020 eV was unexpected. A number of ground array experiments designed to measure cosmic rays in the ultra-high energy range have since been carried out. These experiments include the 3 Haverah Park array in England [12], the Yakutsk array in Russia [13], the Sydney University Giant Airshower Recorder (SUGAR) in Australia [14], the KASCADE experiment in Germany [15, 16, 17, 18] and its follow-up KASCADE-Grande [19], and the Akeno Giant Air-Shower Array (AGASA) in Japan [20, 21]. Of these, AGASA has accumulated the largest data set, and sets a standard against which CHICOS can be compared. The AGASA array was in operation between 1990 and 2003. It covered 100 km2 with an array of 111 scintillation detectors. Over 14 years of observation, AGASA recorded nearly 1000 events above 1019 eV, including 11 events above 1020 eV [22, 23, 24]. Their data showed that the slope of the low-energy spectrum extended up to the highest observed energies. This result was in conflict with the theoretical GZK cutoff, which predicted AGASA would observe only 1.9 events above 1020 eV [24]. In addition, the sky map of the UHECR data collected by AGASA showed evidence of small-scale clustering as early as 1996 [25] and this result has been updated and expanded several times [26, 27, 28, 29, 30]. The final data set finds one triplet and 6 additional doublets in a data set of 67 events above 4 × 1019 eV [24]. Each of the 9 pairs of data points has an angular separation of less than 2.5◦ , corresponding to the angular resolution of the array. The probability of that number of pairs arising from a random isotropic distribution is given as less than 0.1%. The combined data set of 92 events above 4 × 1019 eV from Volcano Ranch, Haverah Park, Yakutsk, and AGASA also showed statistically significant clustering [31]. These two unexpected results from AGASA, if confirmed, would have important implications for both astronomy and physics and have fueled continued research of ultra-high energy cosmic rays. A nondetection of the predicted GZK cutoff would imply either that ultra-high energy cosmic rays are not primarily protons or that there are nearby sources capable of accelerating protons to these energies. The identification of small-scale clustering may be a first step to identifying the astrophysical sources of UHECRs. The potential origins of ultra-high energy particles are still an area of great debate and the possibilities will be discussed in chapter 2. In addition to the ground array experiments, high-energy cosmic rays can be 4 measured via the air fluorescence the showers produce. This method was pioneered by the Fly’s Eye experiment in Utah [32, 33], which has since been replaced by the second-generation experiment, the High Resolution Fly’s Eye (HiRes) [34, 35]. The HiRes collaboration has recently reported that they have observed the GZK cutoff at the expected energy of 6 × 1019 eV [36, 37]. They also note that as of 2006, the AGASA collaboration has revised their energy estimates downward by 10–15%, lowering the observed number of super-GZK events from 11 to 5 or 6. The remaining points no longer have sufficient statistical significance to constitute a nondetection of the cutoff [38, 39]. The air fluorescence technique is being used in conjunction with a ground array at the Pierre Auger Observatory currently under construction in Argentina [40, 41]. The Auger Observatory is the first of two planned sites, the second of which will be located in the northern hemisphere. The ground array at each site is to consist of 1600 water Cherenkov detectors spread out over 3000 km2 . Although still under construction, the exposure of the Auger Observatory is already twice that of HiRes and 4 times that of AGASA [42]. In addition, the combination of a ground array with fluorescence detectors provides a unique advantage in the calibration of their results. Initial data from the southern Pierre Auger Observatory have confirmed the HiRes detection of the GZK cutoff [43, 44, 42]. However, the exact shape of the upper end of the cosmic ray spectrum is still of great interest. Both HiRes and the Pierre Auger Observatory have failed to observe the smallscale clustering reported by AGASA [45, 46]. However, the Pierre Auger Observatory has recently claimed to observe a correlation between ultra-high energy cosmic rays and active galaxies [47, 48]. Using the same methodology, data from HiRes shows only the degree of correlation expected by chance from a random, isotropic distribution of cosmic rays [49]. Given that the nuclei of nearby active galaxies are considered to be likely candidates for UHECR sources, this possibility merits further investigation. A the results of a correlation searach between CHICOS UHECR data and nearby active galaxies is presented in chapter 9. 5 1.2 The CHICOS Experiment The California HIgh School Cosmic Ray Observatory, or CHICOS, is a collaboration between U.C. Irvine, C.S.U. Northridge, and the California Institute of Technology (Caltech). The Project Director is Dr. Robert McKeown of Caltech and the Education Director is Dr. Ryoichi Seki of C.S.U. Northridge. Financially, the project is primarily supported by an NSF grant, with hardware donations from Los Alamos National Laboratory and IBM. The CHICOS project was conceived as a collaboration with Los Angeles-area high schools for the dual purposes of education outreach and UHECR research. The CHICOS array is made up of pairs of solid-scintillator cosmic ray detectors spread throughout the San Gabriel and San Fernando valleys. Each pair of detectors is situated in a high school (or in some cases a middle or elementary school), with the detectors and a GPS antenna typically placed on the roof and a workstation in a nearby science classroom. A major advantage of using secondary schools as detector sites is that the infrastructure needed for power and data transfer is already in place, allowing for a very large array to be built with minimal cost. See chapter 4 for details of the construction and operation of the array. The teachers who are involved with the project are encouraged to integrate it into the science curriculum. All CHICOS data is made available to teachers and students via the project website for this purpose. The project also offers a series of week-long summer programs for students from participating schools. Other cosmic ray detector arrays have used schools as detector sites (for example, ALTA in Alberta and CROP in Nebraska), but the CHICOS array differs from these projects in its greater emphasis of science goals in addition to educational contributions. Much work has gone into the development of user-friendly event reconstruction software. In keeping with the educational mission of CHICOS, this is available in interactive format on the CHICOS webpage.1 . The details of the event reconstruction software are discussed in chapter 5. Chapter 6 describes the methods used to assess 1 www.chicos.caltech.edu 6 the accuracy of the reconstructor software using simulated, unthinned air showers at 1017 eV. The CHICOS array has been designed to observe cosmic ray air showers with energies of about 1018 eV and above. In addition, a more closely spaced subset of the array, nicknamed Chiquita and located on the Caltech campus, is designed to observe showers down to energies of 1016 eV. Data from the smaller array is in the energy range where the spectrum has been more accurately measured, and thus provides a useful calibration of the data reconstruction methods. Chapters 7 and 8 present the data obtained by the Chiquita and CHICOS arrays, respectively. 7 Chapter 2 Origins of UHECRs The flux of cosmic rays appears to fall smoothly over at least 10 orders of magnitude, decreasing approximately as the inverse cube of the energy. There is a slight break at approximately 1015.5 eV, known as the “knee,” where the slope steepens from E −2.7 to E −3 . The spectrum steepens again to E −3.3 at 1017.7 eV, then flattens slightly to E −2.7 at the “ankle,” around 1019 eV [51, 52]. Among the physical processes that may be able to explain the power-law spectrum is diffusive shock acceleration, reviewed briefly in section 2.1. Theoretical considerations, discussed in Section 2.2, predict that the flux of cosmic rays should drop sharply above 6 × 1019 eV, though there remains disagreement over whether this has been observed. Despite the relative uniformity of the spectrum over the measured energy range, cosmic rays are believed to come from a diversity of sources, ranging from solar to galactic to extragalactic. In the ultra-high energy range around the ankle and above, it is believed that extragalactic particles dominate the flux for reasons discussed in section 2.3, although the specific sources are unknown. 2.1 Diffusive Shock Acceleration One process by which cosmic rays may acquire ultra-high energies is diffusive shock acceleration. This is a process in which the the particle repeatedly crosses a shock front, gaining energy at each crossing [53]. This theory is appealing both because shock fronts are a common astrophysical phenomenon and because the output of 8 diffusive shock acceleration is a power-law energy spectrum. Following Malkov [54], consider a shock front in which there is a velocity change across the shock front from u1 to u2 . A particle with velocity v and momentum vector p crossing the shock front at angle θ to the shock normal emerges with momentum p0 . Define the dimensionless velocity change between frames upstream and downstream 2 of the shock to be β = u1 −u . It can then be shown from the transformation between c frames in special relativity that the relationship between p and p0 is µ 0 ¶2 µ ¶ 1 2βc β 2 c2 p 2 2 = 1+ cos θ + 2 − β sin θ . p 1 − β2 v v (2.1) For a nonrelativistic shock, β ¿ 1, and to first order in β we have µ ¶ βc p =p 1+ cos θ , v 0 (2.2) which can equivalently be written in vector notation as ∆p = p0 − p = p · (u1 − u2 ) . v (2.3) The flux of particles that go from a momentum less than p to a momentum greater than p as they cross the shock can be found by integrating p over all possible directions: Z p Φ(p) = Z 0 dp f (p0 ) v · n p02 dΩ Zp−∆p ≈ ∆p f (p) v · n p2 dΩ ¸ · Z p · (u1 − u2 ) 2 p dΩ. ≈ f (p) v · n v (2.4) When ∆p ¿ p and β ¿ 1, the momentum distribution function will be approximately isotropic, f (p) ≈ f (p). Under these assumptions, equation (2.4) simplifies 9 to Z 2 Φ(p) = p f (p) · ¸ p · (u1 − u2 ) 4π 3 v·n dΩ = p f (p)n · (u1 − u2 ) . v 3 (2.5) Particle conservation requires that divergence of the momentum space acceleration flux balance the difference between the upstream and downstream momentum distributions and the source term Q(p) of particles being injected into the shock. This can be written [55] as ∂Φ(p) − n · u1 4πp2 f1 (p) + n · u2 4πp2 f2 (p) = 4πQ(p). ∂p (2.6) Using equation (2.5) with equation (2.6), we obtain the momentum distribution produced by the shock: −q Z p f2 (p) = p (Q(p0 ) + n · u1 f1 (p0 )) p0q−1 dp0 , (2.7) 0 where q= 3r 3n · u1 = . n · (u1 − u2 ) r−1 (2.8) It can be seen from this expression that the output energy spectrum is a power law with slope q determined by r, the compression ratio of the shock. For a strong shock, r = 4, and the output spectrum f (p) ∝ p−4 corresponds to an energy spectrum proportional to E −2 . 2.2 The GZK Cutoff As ultra-high energy cosmic rays travel through space, they interact with the cosmic microwave background (CMB). There are two main types of interactions involving cosmic ray protons: pair production and photo-pion production. Photo-pion production may proceed as p + γ → π 0 + p, p + γ → π + + n, or via the production of multiple pions. Single pion production dominates at energies just above the in- 10 teraction threshold, while the cross section for multiple pion production dominates at higher energies [56]. Photo-pair production, which becomes important at energies below the photo-pion threshold, proceeds as p + γ → p + e+ + e− . To obtain the energy threshold for either interaction (following Schlickeiser [57]), we must work in the Lorentz geometry, where the line segment is defined by ds2 = c2 dt2 − dx2 − dy 2 − dz 2 . We define the four-momentum of a particle to be P = ( c² , p), where P 2 = m2 c2 is an invariant quantity. In general, the energy threshold for particle production occurs when the initial energy of all particles in the center-of-mass frame is equal to the rest mass of all particles following the interaction. The total energy before the interaction in the center-of-mass frame is given by COM 2 ) = (EaCOM + EbCOM )2 = c2 (PaCOM + PbCOM )2 = c2 (Pa + Pb )2 (Etotal (2.9) by virtue of the invariance of P 2 . Using P 2 = m2 c2 and Pa Pb = ²ca ²cb − pa pb , we obtain COM 2 (Etotal ) = m2a c4 + m2b c4 + 2²a ²b − 2pa pb c2 . (2.10) The energy threshold for the interaction is Eth = ma c2 + mb c2 + ∆mc2 , (2.11) where ∆m is the difference in rest mass between the incoming and outgoing particles. COM Setting Etotal = Eth , we have µ 2 4 ²a ²b − pa pb c = ma mb c + ∆mc 4 ∆m ma + mb + 2 ¶ . (2.12) We can simplify equation (2.12) by first rewriting it as ²a ²b − pa pb c2 = 1 + ∆m ma mb c4 µ 1 1 ∆m + + ma mb 2ma mb ¶ . (2.13) The left-hand side of equation (2.13) can now be written in terms of the Lorentz p factor γ = mc² 2 , where we have also used pa pb = pa pb cos θ and p = γ 2 − 1 mc. This 11 produces µ q γa γb − (γa2 − 1)(γb2 − 1) cos θ = 1 + ∆m 1 1 ∆m + + ma mb 2ma mb ¶ . (2.14) In the case of a proton-photon interaction, where particle b is massless, equation (2.12) reduces to ³ ²b γa − p ´ γa2 − 1 cos θ ¶ µ ∆m . = ∆mc 1 + 2ma 2 (2.15) For relativistic cosmic rays, with γa À 1, equation (2.15) becomes ∆mc2 γa = (1 − cos θ)²b µ ∆m 1+ 2ma ¶ . (2.16) The energy required for the interaction is therefore £ ¤ 4 2 2 (∆m + m ) − m c a a E = γa m a c 2 = . 2²b (1 − cos θ) (2.17) The minimum energy for the interaction occurs in a head-on collision, with cos θ = −1. In this case, equation (2.17) becomes £ ¤ (∆m + ma )2 − m2a c4 . Emin = 4²b For photo-pion production, ma = mp and ∆m = P (2.18) mπ , the total mass of pions produced. Using ²b = h²i, the average energy of CMB photons, the minimum energy needed for the proton to initiate pion production is given by ¤ £P ( mπ + mp )2 − m2p c4 Emin = . 4 h²i (2.19) The average energy of CMB photons is approximately h²i = 7 × 10−4 eV. Using 2 2 2 0 m± π = 139.570 MeV/c (mπ = 134.977 MeV/c ), and mp = 938.272 MeV/c , we have Emin = 1.0 × 1020 eV (2.20) 12 for the case of single pion production. The threshold for multiple pion production is correspondingly higher. However, because the blackbody distribution of photons has a tail that extends to higher energies, protons of lower energies can occasionally undergo photo-pion production. Repeated encounters will eventually cause the energy of the proton to fall below the energy threshold for a given interaction with the majority of CMB photons. The attenuation length due to photo-pion production for a proton with energy 1020 eV is approximately 100 Mpc, but drops to 10 Mpc for a proton at 1021 eV [58]. Given these limits, ultra-high energy cosmic rays would only be observable if they originate from a relatively small volume around our location. A volume 10 Mpc in radius would encompass only the Local Group of galaxies. Ultra-high energy cosmic rays that originate farther away would be observed as an accumulation of flux just below the threshold for photo-pion production, beyond which the spectrum would drop quickly. This predicted cutoff in the cosmic ray spectrum is known as the GZK effect after Greisen [10], and Zatsepin and Kuz’min [11], who developed the theory independently in 1966. The cutoff energy for analyses of ultra-high energy cosmic rays is typically taken to be 4 × 1019 eV. This is based on simulations that show UHECRs emitted by relatively nearby sources (z . 0.057) accumulate just above that energy, at approximately 5 × 1019 eV with a steep drop-off around 6 × 1019 eV [56]. Energy loss by pair production begins to dominate below about 3 × 1019 eV [59]. By equation (2.18), the threshold for electron pair production with a photon at the average energy of the CMB is ¤ (2me + mp )2 − m2p c4 . Emin = 4 h²i £ (2.21) Given the electron mass of 0.511 MeV, the threshold energy for this process is Emin = 6.9 × 1017 eV. (2.22) 13 The mean energy loss for this process is only 0.1% per encounter, compared to 20% for photo-pion production, making photo-pair production a less efficient mechanism for energy loss [60]. The attenuation length for pair production reaches a minimum of approximately 1000 Mpc at 2 × 1019 eV [61]. Heavier nuclei are limited in the distance they can travel by photo-disintegration effects [62, 63]. The current theory that cosmic rays at the highest energies are predominantly protons or light nuclei is supported by the data from multiple experiments, including AGASA and HIRES [64]. 2.3 Potential UHECR Sources The observation of cosmic rays above the GZK cutoff raises questions about the origins of these particles. Particles at energies at and below the knee are believed to be galactic in origin, with the primary source being supernova shocks [65]. A secondary source may be OB associations, in which particles are accelerated by turbulent motion and stellar winds [66]. No individual sources have yet been identified, however. Only a few known astrophysical phenomena are plausible sources of UHECRs. These are defined by the “Hillas criterion” [67], which states that a particle accelerated in a magnetic field can only continue gaining energy until its Larmor radius becomes comparable to the size of the acceleration region. Following Longair [68], the Larmor radius of a relativistic particle can be obtained from its equation of motion, d (γm0 v) = Ze (v × B) . dt Using γ = (2.23) p 1 − v·v , this becomes c2 ³v · a´ d dv 3 m0 (γv) = m0 γ + m0 γ v . dt dt c (2.24) For movement in a magnetic field, the acceleration is perpendicular to the parti- 14 cle’s velocity, v · a = 0. Hence γm0 dv = Ze (v × B) . dt (2.25) Considering only the component of v perpendicular to the magnetic field, and equating the acceleration with the centrifugal acceleration, we have v2 ZevB = , γm0 r (2.26) which leads directly to the relativistic Larmor radius γm0 v . ZeB rL = (2.27) For a relativistic particle, E ≃ pc = γm0 vc. Rewriting the Larmor radius in terms of the energy of the particle, we have rL = E . ZeBc (2.28) Expressing the particle’s energy in units of E18 ≡ E/1018 eV and the magnetic field in microgauss, equation (2.28) becomes rL = 1018 eV E18 E18 = 1.08 kpc. −6 ec · 10 G ZBµG ZBµG (2.29) The size L of the region that accelerates the particle must be at least 2rL . Hence Lkpc & 2E18 . ZBµG (2.30) It is necessary to modify this result to take into account the shock speed βc that is causing the acceleration [67], yielding Lkpc & 2E18 . ZBµG βc (2.31) 15 Equivalently, the Hillas criterion for the maximum energy to which a region of size L can accelerate a particle is E18,max ∼ 0.5ZBµG Lkpc βc. (2.32) At higher energies, the particle will move beyond the region permeated by the magnetic field, and will escape from the system. The interstellar magnetic field, for example, is approximately 2–4 µG [69]. Given the disk thickness of the galaxy of approximately 300 pc, protons can be accelerated in the galactic magnetic field to at most ∼ 1018 eV [60]. For this reason, it is speculated that most UHECRs are extragalactic in origin. Speculated extragalactic sources of UHECRs include the following astrophysical phenomena, as well as more exotic possibilities [52]. 2.3.1 Radio Galaxies and AGN The extended lobes of radio galaxies typically contain “hot spots,” which are interpreted to be the shock front of the relativistic jets that emanate from the active galactic nucleus, or AGN. The hot spots contain a magnetic field up to a few hundred µG in an area of a few kpc2 [70]. Under these conditions, the Hillas criterion yields Emax ≈ 1020 eV. This estimate can be refined by taking into account losses due to synchrotron radiation and photon interactions [71]. Balancing the timescale for energy loss against the timescale for accleration yields an upper bound on the energy of the cosmic ray particles that can be produced. To obtain the timescale for acceleration, we first write the momentum-space particle conservation equation [54]. Defining κ1 and κ2 to be the upstream and downstream diffusion coefficients respectively, the number of particles interacting with the shock is µ 4πf (p) κ1 κ2 + u1 u2 ¶ . (2.33) 16 Particle conservation requires the change in particle number to be balanced by the divergence of the momentum-space flux and the “source” term, which in this case represents the downstream flow of particles away from the shock: · µ ¶¸ ∂ κ1 ∂Ω(p) 2 4πp f (p) + κ 2 u2 + = 4πp2 f (p)u2 . ∂t u1 ∂p (2.34) Equation (2.34) can be simplified to µ κ1 κ2 + u1 u2 ¶ ∂f u1 − u2 ∂f + p + u1 f = 0. ∂t 3 ∂p (2.35) As shown by Drury [72], it follows that the mean acceleration time from some momentum p0 to p is 3 htacc (p)i = u1 − u2 Z pµ p0 κ1 κ2 + u1 u2 ¶ dp . p (2.36) The timescale for acceleration of particles of momentum p is therefore 3 τacc = u 1 − u2 µ κ1 κ2 + u1 u2 ¶ . (2.37) For a strong shock, r = u1 /u2 = 4. If the upstream and downstream diffusion lengths are assumed to be equal, the acceleration timescale further simplifies to κ τacc = 20 2 . u1 (2.38) Following Biermann and Strittmater [71], in order to evaluate this timescale in the environment of an active galaxy, we need to evaluate the diffusion coefficient κ. The diffusion coefficient is related to the mean free path λ and to the scattering time τS ∼ λ/c by µ κ∼ 4 3𠶵 2¶ λ . τS (2.39) In the small-angle resonant scattering approximation, where the particle deflection is dominated by Alfvèn waves with wavelength equal to the gyroradius of the particle, 17 the mean free path is given by λ = rg B 2 /8π . I(k)k (2.40) Here I(k) is the magnetic energy density per unit wavenumber k in the magnetic field. The resonant scattering approximation requires k ∼ 1/rg . The mean free path therefore depends on the spectrum of the turbulent magnetic field. If we assume Kolmogorov-type turbulence, I(k) = I0 (k/k0 )β , where β ≃ 5/2, we have (B 2 /8π) λ = rg = rg I(k)k µ B 2 /8π k0 I0 ¶µ k k0 ¶β−1 . (2.41) The factor k0−1 corresponds to the outer scale of turbulence, or equivalently, to rg,max , the gyration radius of the most energetic particles. This can be simplified by introducing b, the ratio of turbulent to ambient magnetic energy density: Z ∞ I0 k0 I0 k0 = . 2 (B /8π) (β − 1)(B 2 /8π) b= k0 (2.42) Inserting equation (2.42) into the expression for λ in equation (2.41), we have · rg λ= b(β − 1) ¸µ rg,max rg ¶β−1 . (2.43) From equation (2.38), we can now write the acceleration timescale as 80 τacc ∼ 3π µ c u21 ¶· rg b(β − 1) ¸µ rg,max rg ¶β−1 . (2.44) The timescale for proton energy loss to synchrotron radiation is τsyn = 6πm3p c , σT m2e γp B 2 (2.45) where mp is the proton mass, σT is the Thompson cross section, and γp is the Lorentz factor of the accelerated proton. 18 The general expression for energy loss due to proton-photon interactions is 1 = τpγ Z ∞ c d² n(²) 2 2 2γp ² ²rmth /2γp Z 2γp ² kp (²0 )σ(²0 )²0 d²0 , (2.46) ²th where n(²) is the number density of photons per unit energy interval, ²th is the energy threshold for inelastic collisions, kp (²) is the inelasticity, and σ(²) is the cross section for interaction in the relativistic proton frame. The number density of photons is assumed to have the form n(²) =   (N0 /²0 )(²/²0 )−2 , ²0 ≤ ² ≤ ²∗ ,  0, (2.47) otherwise. where ²0 and ²∗ correspond to radio and γ-ray energies respectively. The integral in equation (2.46) can then be evaluated to be ¸µ 2 ¶ · B 1 a σγp = γp , ∗ τpγ 6π ln (² /²0 ) mp c (2.48) where a is the ratio of photon to magnetic energy density, given by a= N0 ²0 ln (² ∗ /²0 ) . (B 2 /8π) (2.49) The total energy loss timescale for protons is therefore 1 1 1 1 = + = (1 + Aa), τp τp,sy τpγ τp,syn (2.50) σγp σγp (mp /me )2 ≈ 1.6 × 105 ≈ 200. A= ∗ σT ln (² /²0 ) σT (2.51) where This leads to a maximum Lorentz factor for accelerated protons of · 27πb e γp,max = (β − 1)1/2 2 320 r0 B ¸1/2 ³ ´ µ ¶ µ ¶1/2 u mp 1 , c me 1 + Aa (2.52) 19 where r0 is the classical electron radius. Given typical hotspot parameters (β ≃ 5/3, a ∼ 0.1, b ∼ 0.5, u ∼ 0.3c, and B ∼ 3 × 10−4 G) [52], the corresponding maximum energy to which a proton can be accelerated is Ep,max = γp,max mp c2 ∼ 2 × 1020 eV. (2.53) Particles can also be accelerated to ultra-high energies within the jets or within the AGN itself. For example Knot A in the M87 jet has linear dimension LM87 ∼ 2 × 1020 cm and magnetic field B ∼ 300 µG [73]. A typical active galactic nucleus can have L ∼ 1015 cm and B ∼ 1 G [74]. It should be noted, however, that there is limited number of AGN within 100 Mpc of our location, and none are clear candidate sources for the 1020 eV AGASA events. Associations between UHECR data and BL Lac objects have been investigated [75, 76, 77, 78] but the claims of a correlation are contested [79]. More recently, the Auger Collaboration has claimed to observe a correlation between their UHECR data and nearby AGN [47, 48]. The HiRes Experiment has failed to reproduce this result [49]. The details of these correlation searches are presented in section 9.1. 2.3.2 Neutron Stars and Magnetars Given the constraints of the GZK cutoff, it is attractive to consider nearby phenomena that might produce the observed cosmic ray events above 1020 eV. Unfortunately there are very few plausible possibilities within our own galaxy. One suggestion is that neutron stars may transfer their rotational kinetic energy to the kinetic energy of heavy nuclei via relativistic magnetohydrodynamic wind [80]. A young neutron star may have a rotation rate of Ω ∼ 3000 rad s−1 and a surface magnetic field of up to BS & 1013 G at RS = 106 cm. The field strength decreases as B(R) = BS (RS /r)3 . The light cylinder of the star (the maximum radius at which the dipole field can be sustained), is located at RLC = c/Ω. The magnetic field at the light cylinder is 20 therefore µ BLC = BS RS c/Ω ¶3 = 1010 B13 Ω33k G, (2.54) where B13 ≡ B/1013 G and Ω3k ≡ Ω/3000 rad s−1 . The maximum energy of particles that can be contained in the system out to the radius of the light cylinder is Emax = ZeBLC RLC c ≃ 8 × 1020 Z26 B13 Ω23k eV, (2.55) where Z26 ≡ Z/26. Magnetars are neutron stars with unusually high magnetic fields, in the range of 1015 G. A “fast magnetar” may have a rotational frequency of 104 rad s−1 . Using these values in equation (2.55), we find the maximum energy is Emax = ZeBLC RLC c ≃ 3 × 1022 ZB15 Ω24 eV, (2.56) where B15 ≡ B/1015 G and Ω4 ≡ Ω/104 rad s−1 [81]. 2.3.3 Quasar Remnants A quasar remnant is the end-stage evolution of a luminous quasar: a spinning supermassive black hole, threaded by magnetic fields generated by currents flowing in a disc around it. We appear to live in an epoch where luminous quasars are rare. However, extrapolating from the number of luminous quasars at high redshift, the number of quasar remnants nearby may be large and these have been postulated to be a source of UHECRs [82]. The relatively dormant supermassive black holes found in many giant elliptical galaxies are likely examples of such “dead” quasars. A Kerr black hole whose event horizon is threaded by an external magnetic field can act as a battery [83], and the EMF generated would potentially be sufficient to accelerate a proton to ultrahigh energies. If B is the strength of the ordered poloidal magnetic field near the hole, then V ∼ aB, where a is the hole’s specific angular momentum [84]. (For a black hole of mass M , a ≤ M .) In appropriate astrophysical 21 units, the EMF generated is 20 ∆V ∼ 9 × 10 ³a´ M M9 B4 V, (2.57) where M9 ≡ M/109 M¯ and B4 ≡ B/104 G. In the case of an advection-dominated accretion flow (ADAF) onto the black hole, the strength of the magnetic field near the event horizon is related to the accretion rate Ṁ (in units of M¯ yr−1 ) by B4 = 1.33M9−1 Ṁ 1/2 , (2.58) under the assumption that the energy density of the magnetic field is in equipartition with the rest mass of the accreting matter [85]. The combination of equation (2.57) and equation (2.58) yields a maximum possible EMF of ∆V = 1.2 × 1021 Ṁ 1/2 V, (2.59) where we have taken a ≃ M for a maximally rotating black hole. The maximum obtainable energy, however, is less than this quantity because energy is lost to curvature radiation [86]. For an average curvature radius ρ, the rate of energy loss by a particle of energy E = γmc2 is P = 2 Z 2 e2 cγ 4 . 3 ρ2 (2.60) The energy change per unit distance for a particle with mass µmp is eZ∆V P dE = − , ds h c (2.61) where h is the gap height of the black hole. Integrating over s from 0 to h yields the 22 maximum energy to which the particle can be accelerated: 19 Emax = 3 × 10 µZ −1/4 1/2 1/4 M9 B4 µ 2 ¶1/4 ρh eV. Rg3 (2.62) This can be simplified by assuming h ≈ Rg and r ≈ Rg . For a proton (µ = 1 and Z = 1), we can use then equation (2.58) to obtain 1/4 Emax = 1.0 × 1020 Ṁ10 M9 eV, (2.63) where Ṁ10 ≡ Ṁ /10M¯ yr−1 . 2.3.4 Starburst Galaxies and LIGs Starbursts are galaxies undergoing a period of intense star formation. Due to numerous supernovae, a cavity of hot gas can be created in the center of an active region. Given that the cooling time of the gas is longer than the expansion timescale, the hot gas will expand and form a shock front as it contacts the cooler interstellar medium. Ions such as iron nuclei can be accelerated to super-GZK energies in these conditions by Fermi’s mechanism [87]. The acceleration of nuclei in this scenario is a two-stage process beginning with diffusive acceleration to energies of 1014 −−1015 eV at supernova shock fronts [88]. The ions are then injected into the galactic-scale wind created by the starburst region [89, 90, 91]. The maximum particle energy that can be obtained from this process is 1 2 τ, Emax = ZeBvsh 4 (2.64) where vsh is the shock velocity and τ is the age of the starburst [87]. The shock velocity is related to the kinetic energy flux of the superwind, Ėsw , and the mass flux, Ṁ , generated by the starburst as 1 2 . Ėsw = Ṁ vsh 2 (2.65) 23 Substituting this into equation (2.64), we have 1 Ėsw Emax = ZeB τ. 2 Ṁ (2.66) Two nearby starburst galaxies that are candiates for UHECR production are M82 and NGC253. NGC253, for example, has a kinetic energy flux of 2 × 1042 erg s−1 and a mass flux of 1.2 M· yr−1 [92], and a magnetic field strength of B ∼ 50 µG [93]. This leads to an estimated maximum energy for iron nuclei of Fe Emax = 3.4 × 1020 eV. (2.67) In an axisymmetric (ASS) galactic field model, the arrival directions of the 4 highest-energy cosmic rays observed as of 2003 were found to be associated with starburst galaxies [94]. However, in a bisymmetric (BSS) galactic field model, smaller cosmic ray deflections result in an absence of correlation. Luminous infrared galaxies (LIGs), which may form after a collision between galaxies, are similar to starburst galaxies on a larger scale [95]. LIGs have luminosities above 1011 L¯ , and are the dominant extragalactic objects in the local universe in that luminosity range. The triplet event observed by AGASA [25, 26, 27] is potentially associated with the LIG Arp299 [96]. 2.3.5 Gamma Ray Bursts Gamma ray bursts (GRBs) are short bursts of high-energy radiation [97]. They are among the most energetic phenomena in the universe; a single gamma ray burst may be brighter than all other gamma ray sources combined. The most popular theory of the origin of GRBs is the “fireball” model: GRBs are believed to arise from the dissipation of the kinetic energy of a relativistically expanding wind, the cause of which remains unknown [98]. Gamma ray bursts feature a rapid rise time and short duration (∼ 1 ms), which implies a compact source. The 24 detection of afterglows has allowed the measurement of the redshifts of some GRB host galaxies, and confirmed that GRBs originate at cosmological distances [99, 100]. The compactness and high gamma ray luminosity result in a high optical depth to pair creation. This creates a thermal plasma, the radiation pressure of which drives relativistic expansion. Conditions within the fireball may accelerate protons to energies greater than 1020 eV, provided the magnetic field is close to equipartition with electrons [52]. The principal difficulty with the GRB theory of cosmic ray origins is the cosmological distances involved. If the GRB redshift distribution follows that of the star formation rate in the universe, which increases with redshift, the flux of ultra-high energy cosmic rays is predicted to be attenuated by the GZK cutoff at energies above 3 × 1019 eV [101, 102]. 2.3.6 Top-Down Models Due to the difficulty in finding physical phenomena that can accelerate particles to ultra-high energies, many alternative models have been proposed in which ultra-high energy cosmic rays originate in the decay of massive unstable particles. This idea originated with Georges Lemaı̂tre [103], who in 1931 proposed that all material in the universe originated in the decay of a “primeval atom.” In top-down models, massive particles (generically known as “X” particles) with mass mX > 1011 GeV are generated from high energy processes in the early universe, and their decay continues in the present time. UHECRs emitted by such decays avoid the GZK attenuation experienced by particles with a cosmological origin. A wide variety of specific mechanisms involving theories such as string/M theory, supersymmetry (SUSY), grand unified theories (GUTs), and TeV-scale gravity have been invoked as possible origins of ultra-high energy cosmic rays [52, 51, 61]. 25 Chapter 3 Cosmic Ray Air Showers When an ultra-high energy cosmic ray enters the atmosphere and precipitates an air shower, much of the information describing the incident particle is lost. Properties of interest include the species, energy and incident angle of the primary particle. In order to extract this information from the ground data, we require a reliable model of air shower development. The CHICOS project has used extensive simulations of air showers to construct analytical descriptions of the shape of the air shower front. The components of the particle cascade are examined in section 3.1. The measured intensity of the air shower is characterized by the lateral distribution function (LDF) and the time distribution function (TDF). The CHICOS-specific LDF is presented in section 3.2 and the CHICOS-specific TDF is presented in section 3.3. 3.1 Air Shower Development An ultra-high energy cosmic ray incident on the Earth will eventually collide with an atom in the atmosphere. The output of such a collision will include protons, neutrons, smaller atomic nuclei, and mesons [64]. Some of these particles will go on to interact with other atoms in the atmosphere, forming a hadronic cascade that makes up the core of an air shower (figure 3.1). Large numbers of pions are produced in the hadronic interactions. The main decay 26 Primary Particle Molecule in the Atmosphere K+ K− π+ π − K0 K0 π0 K+ K− µ+ µ− γ γ e+ νµ _ νµ νµ µ+ µ− π+ π − e− e+ e− _ νµ µ+ µ− Muonic Component p, n, K , π + Nuclear Fragments Hadronic Cascade γ e+ γ e− γ e+ γ e− Electromagnetic Component Figure 3.1. Diagram of an air shower. An air shower comprises a hadronic core, a muonic component, and an electromagnetic cascade. Decay paths leading to the these three main components are shown. mode of π 0 particles is π0 → γ + γ (τ = 0.83 × 10−16 s). (3.1) The high energy photons produced by this decay initiate an electromagnetic cascade via alternating electron-positron pair production and bremsstrahlung. This process is interrupted when the electrons fall below the critical energy for air of ∼ 81 MeV, at which point more energy is lost to ionization than to bremsstrahlung, and the intensity of the electromagnetic cascade begins to attenuate [55]. In addition to the hadronic and electromagnetic components of the air shower, there is also a muonic component. Muons are created by the decays π ± → µ± + νµ (ν µ ) (τ = 2.063 × 10−8 s) (3.2) 27 Shower Axis Curvature Delay Spread Shower Plane Detectors θ Ground r perp Figure 3.2. Profile of an air shower. The shower front is a curved surface with a finite thickness. Moving away from the center of the shower, the particle intensity decreases, while the spread in the depth of the shower front increases. Distance from the core of the shower is measured along r⊥ , perpendicular to the shower axis. and K ± → µ± + νµ (ν µ ). (τ = 1.237 × 10−8 s) (3.3) Muons are the most penetrating component of the air shower, and reach the ground with little attenuation and only slight energy loss to ionization. They do contribute somewhat to the electromagnetic cascade via the decay µ± → e± + νe (ν e ). (τ = 2.197 × 10−6 s) (3.4) The electromagnetic component dominates the air shower, comprising about 90% of shower particles. The muonic component accounts for most of the remaining 10%, with the hadronic core making up less than 1% of the total shower. The resulting particle front of the air shower is a thin curved surface, traveling close to the speed of light, which spreads out from the axis of the primary particle’s trajectory. The width of this shower front increases with distance from the shower axis (figure 3.2). 28 The precise evolution of a cosmic ray air shower can be modeled with codes such as AIRES (AIRshower Extended Simulations) [104, 105, 106] and CORSIKA (COsmic Ray SImulations for KAscade) [107]. The CHICOS project is currently using AIRES version 2.6.0, which is freely available from the Universidad Nacional de La Plata, Argentina [108]. The simulation code in turn depends on specific models of hadronic interactions; AIRES uses the SIBYLL and QGSJET models. We have made use of a series of AIRES simulations in order to accurately model the air showers observed by CHICOS. The simulations were divided into two groups: low-energy (for the Chiquita subarray) and high-energy (for the CHICOS array). Protons and iron nuclei were used as the primary particles, and the resulting showers were measured at the CHICOS average altitude of 250 meters above sea level. The low-energy simulations cover the energy range between 1016 eV and 1017.5 eV. Ten showers were simulated at each primary energy (log (E/eV) = 16.0, 16.5, 17.0, 17.5) and each zenith angle (cos θ = 0.75, 0.85, 0.95), for each type of primary particle (proton or iron nucleus). The iron showers were used as the basis for the low-energy LDF, based on evidence that heavy nuclei predominate at those energies [109, 110, 111]. The high-energy simulations cover the energy range between 1018 eV and 1020.5 eV. Ten showers were simulated at each primary energy (log (E/eV) = 18.0, 18.5, 19.0, 19.5, 20.0, 20.5) and each zenith angle (cos θ = 0.55, 0.65, 0.75, 0.85, 0.95), for each type of primary particle (proton or iron nucleus). The proton showers were used as the basis for the high-energy LDF, based on evidence that protons predominate at those energies [112]. Tracking all particles generated in a simulated ultra-high energy air shower is beyond the computational resources available. (An air shower with primary energy 1020 eV will generate approximately 1011 secondary particles.) All simulations have therefore employed statistical thinning, beginning at an energy threshold Eth = 10−7 Eprimary . When an interaction within the shower generates particles with energy below this threshold, only a subset of the secondary particles with E < Eth will continue to be tracked by the simulation. The accepted particle is assigned a 29 statistical weight equal to the number of particles it represents in the simulation. AIRES employs the Hillas thinning algorithm [105]. When a particle with energy E ≥ Eth generates a set of secondary particles with energies Ei , each secondary particle is individually tested against the thinning energy and accepted with probability Pi =   1, Ei ≥ Eth   Ei , Ei < Eth . Eth (3.5) If the primary particle has E < Eth , only one secondary particle will be conserved. It is selected from the set of secondary particles with probability Ei P i = Pn j=1 Ei . (3.6) The weight of the accepted secondary particle is equal to the weight of the primary multiplied by the inverse of Pi . AIRES provides an optional statistical weight factor, Wf , which limits the particle weights that may be assigned. Given a value for Wf , AIRES sets two internal parameters Wr = 14 GeV−1 Eth Wf (3.7) Wy = Wr /8. (3.8) and In an interaction that generates 3 or fewer secondary particles, if the weight of the primary is w > Wy or if wE/min(E1 , ..., Ei ) > Wr , then all of the secondary particles will be kept; otherwise the standard Hillas algorithm is used. If more than 3 particles are generated, then the Hillas algorithm is always used, but if the weight w0 of the single selected secondary is larger than Wr , then m copies of the secondary are kept (each with weight equal to the weight of the secondary particle divided by m. The integer m is adjusted to ensure that Wy < w0 /m < Wr . For the CHICOS simulations, the AIRES statistical weight factor was set to Wf = 1. 30 (a) (b) Figure 3.3. Comparison of CORSIKA and AIRES simulations. A set of 10 showers was generated with each program, using the input parameters E = 1017 eV and cos θ = 0.95. The lateral distributions are shown for (a) electrons and (b) muons. The discrepancy in the first bin is due to a difference in the inner radial cutoff of the simulations. Shower particles were tracked down to Ee± ,γ = 1 MeV and Eµ± = 20 MeV. A 5 MeV cutoff corresponding to the detector sensitivity threshold was applied to particles reaching the ground. The 5 MeV energy threshold applied to ground particles is the same as that used by the KASCADE experiment [113, 114]. The number of electrons reaching the ground with energy between 1 and 5 MeV is approximately 20% of the total. The total number above threshold, however, is not a sensitive function of the cutoff energy in the 5 MeV range. The accuracy of our AIRES simulations has been confirmed by performing a small series of simulations at 1017 eV using the CORSIKA code (figure 3.3). The lateral distributions of electrons and muons generated by the two codes were found to agree well. 3.2 Lateral Distribution Function An air shower front is characterized by its lateral distribution function (LDF), which describes the intensity of particles ρ(r⊥ ; E, θ) as a function of perpendicular distance r⊥ from the shower core and is an implicit function of the energy and angle of incidence 31 of the shower. The CHICOS reconstruction software originally used the LDF obtained empirically by the AGASA experiment as a first approximation to the LDF at our altitude. After completing a representative set of simulated showers, a new CHICOSspecific LDF was developed. For a pure electromagnetic cascade, the lateral distribution function is given by the Nishimura-Kamata-Greisen (NKG) function, Ne ρ(r⊥ ) = C 2 RM µ r⊥ RM ¶s−2 µ r⊥ 1+ RM ¶s−4.5 , (3.9) where Ne is the number of particles in the shower, and s is the “age parameter” of the shower [115, 116]. The Molière unit, RM , characterizing the scattering length,1 is equal to 91.6 m at the altitude of AGASA, and 85 m at the altitude of CHICOS. In a cosmic ray air shower, the electromagnetic component is a combination of electromagnetic cascades initiated by the π 0 particles produced in successive interactions of the central hadronic cascade. Thus the electromagnetic component near the center of the shower consists of “younger” (less developed) showers than the electromagnetic component far from the shower axis. In this case the lateral distribution of charged particles becomes flatter than for a single electromagnetic cascade. This distribution can described by the generalized NKG function [118] as µ ρ(r⊥ ) ∝ r⊥ RM ¶−α µ ¶−(η−α) r⊥ 1+ . RM (3.10) This formula is the basis for both the AGASA and CHICOS lateral distribution functions. 1 The Molière unit is defined by RM = XR ES /EC , where the radiation length XR is the scale length for energy losses from electron bremsstrahlung, the critical energy EC is the energy at which bremsstrahlung and ionization losses are equal, and the scattering energy ES relates the meansquare scattering angle to the distance x traversed by an electron in the multiple-scattering formula hθ2 i = (ES /EC )2 x/XR [117]. 32 3.2.1 AGASA LDF The AGASA LDF is given by the modified NKG function µ ρ(r⊥ ) = C r⊥ RM ¶−α µ r⊥ 1+ RM ¶−(η−α) · ³ 1+ r⊥ ´2 1000 m ¸δ , (3.11) where r is the distance in meters from the core of the shower, and C is a proportionality constant related to the energy of the primary particle. The parameters α and δ are found to be 1.2 and 0.6, respectively [119]. The parameter η depends on the incident angle θ, measured from the vertical: η = (3.97 ± 0.13) − (1.79 ± 0.62)(sec θ − 1), (3.12) for incident angles θ ≤ 45◦ . No energy dependence of η has been observed, so it is assumed that this formula for the LDF can be used to describe even the highest-energy showers [23]. The measured intensity S(r) is a function of the LDF and the detector response. For scintillating detectors, the signal is determined by the average energy loss in the scintillator of electrons, photons, and muons. This function can be expressed in units of the energy loss of vertically penetrating muons, Ce , a convenient measure because they determine the peak of the spectrum of single-particle events. Thus, the measured intensity of a vertical shower is given by µ S0 (r) = Ne Ce r⊥ RM ¶−α µ r⊥ 1+ RM ¶−(η−α) · ³ r⊥ ´2 1.0 + 1000 m ¸δ . (3.13) This function has been shown to be valid between 500 m and 3 km from the core of the shower, at energies up to 1020 eV [120]. Using Monte Carlo simulations [121], AGASA finds that for vertical showers, the energy of the incident cosmic ray is related to S0 (600), the measured intensity at a 33 distance of 600 meters from the core, by the formula E0 = (2.03 ± 0.10) × 1017 eV · S0 (600)1.02±0.02 . (3.14) A shower that enters the atmosphere with an inclined trajectory passes through a greater air depth, and the shower development is correspondingly affected. To determine the energy of an air shower at incident angle θ, the measured intensity Sθ (600) must first be converted to an equivalent value of S0 (600) by the formula ¸ X0 X0 2 Sθ (600) = S0 (600) exp − (sec θ − 1) − (sec θ − 1) , Λ1 Λ2 · (3.15) 2 Here X0 = 920 g/cm2 , Λ1 = 500 g/cm2 , and Λ2 = 594+268 −120 g/cm . This conversion formula is valid for θ ≤ 45◦ [119]. 3.2.2 CHICOS LDF Each CHICOS LDF (low-energy and high-energy) was fit separately to the distributions of muons and electrons. For each species (muons and electrons), the AIRES simulations were used to fill histograms of particle intensity as a function of r⊥ ; lowenergy showers were fit between 25 m and 1000 m using 10 m bins, while high-energy showers were fit between 25 m and 4000 m using 50 m bins. The histograms were averaged over the 10 runs at each energy and zenith angle and the standard deviation of the runs was used as the uncertainty in the histogram. The scintillator detectors used by CHICOS do not distinguish between electrons and muons, therefore the measured intensity must be compared with the sum of the electron and muon LDFs: ρtot (r⊥ , E, θ) = ρe (r⊥ ; E, θ) + ρµ (r⊥ ; E, θ). (3.16) Each particle LDF is given by a modified NKG formula similar to that used by 34 AGASA: ρe,µ (r⊥ ; E, θ) = µ ¶−αe,µ µ ¶−(ηe,µ −αe,µ ) · ³ r ´2 ¸δe,µ r⊥ r⊥ ⊥ Ce,µ (E) 1+ 1+ . (3.17) (RM )e,µ (RM )e,µ 1000 m In this function, the parameter C is explicitly a function of energy. Thus no conversion to S0 (600) is necessary in order obtain the energy of a shower after it has been fit to the CHICOS LDF. The Molière radius has in this expression been replaced by an effective Molière radius, which was fit simultaneously with the other parameters. In addition, the constant α has been replaced by a parameterized function. For the lowenergy (iron-primary) LDF, αe,µ = αe,µ (E). For the high-energy (proton-primary) LDF, αe = αe (θ), while αµ remains a constant. The parameters of the low-energy electron LDF are as follows: RM e = 82.0 m δe = 0.4 αe = 1.429 + 0.6220(log(E/eV) − 17.0) (3.18) ηe = 0.307 + 3.656 cos θ log10 (Ce ) = 1.88 + 1.0(log(E/eV) − 17.0) + 5.0(cos θ − 0.85) Similarly, the parameters of the low-energy muon LDF are as follows: RM µ = 102.5 m δµ = −0.9 αµ = 0.5647 + 0.06972(log(E/eV) − 17.0) ηµ = 1.247 + 0.8214 cos θ log10 (Cµ ) = 0.78 + 0.9(log(E/eV) − 17.0) + 1.2(cos θ − 0.85) (3.19) 35 The parameters of the high-energy electron LDF, expressed in a slightly different format are: RM e = 2477 m δe = 0.03107 αe = 2.774 + 1.326(sec θ − 1) (3.20) ηe = 7.794 − 2.404(sec θ − 1) log10 (Ce ) = −0.015 + 0.95(log(E/eV) − 19.0) − 0.56(sec θ − 1) Similarly, the parameters of the high-energy muon LDF are: RM µ = 2560 m δµ = 0.01939 αµ = 0.7701 (3.21) ηµ = 9.020 + 2.552(sec θ − 1) log10 (Cµ ) = 1.2 + 0.97(log(E/eV) − 19.0) − 0.72(sec θ − 1) The CHICOS low-energy LDF is considered valid for energies approximately between 1016 eV and 1019 eV, and for zenith angles out to 45◦ . The high-energy LDF is considered valid for energies of 1018 eV and above, and for zenith angles out to approximately 60◦ . Figure 3.4 shows the behavior of the low-energy (iron-primary) electron LDF over a range of energies and zenith angles, compared with AIRES simulations of particle density. Figure 3.5 shows the same series of plots for the muon component of the showers. Figure 3.6 shows the behavior of the high-energy (proton-primary) electron LDF over a range of energies and zenith angles, compared with AIRES simulations of particle density. Figure 3.7 shows the same series of plots for the muon component of the showers. ρ [m-2] 10 ρ [m-2] ρ [m-2] 36 10 10 1 1 1 10-1 10-1 10-1 10-2 10-2 10-2 -3 10 -3 10 0 100 200 300 400 500 600 700 800 900 1000 r [m] -3 0 100 200 300 400 500 600 700 800 900 1000 r [m] 10 0 100 200 300 400 500 600 700 800 900 1000 r [m] 100 200 300 400 500 600 700 800 900 1000 r [m] 102 ρ [m-2] 103 ρ [m-2] ρ [m-2] (a) 102 102 10 10 10 1 1 10-1 10-2 0 1 10-1 100 200 300 400 500 600 700 800 900 1000 r [m] 10-2 0 10-1 100 200 300 400 500 600 700 800 900 1000 r [m] 0 (b) Figure 3.4. Low-energy iron-primary electron LDF. The behavior of the simulated electron/positron density as a function of r⊥ for iron primaries of energy (a) E = 1016 eV, and (b) E = 1017 eV. Within each set at a given energy, results are shown (from left to right) for zenith angles cos θ = (0.75, 0.85, 0.95). Points with error bars are AIRES output (mean and standard deviation of 10 runs). The solid curve overlay shows the electron LDF parameterization defined in equation (3.18). 1 10-1 ρ [m-2] ρ [m-2] ρ [m-2] 37 1 10-1 10-2 10-1 10-2 0 100 200 300 400 500 600 700 800 900 1000 r [m] 1 10-2 0 100 200 300 400 500 600 700 800 900 1000 r [m] 0 100 200 300 400 500 600 700 800 900 1000 r [m] 100 200 300 400 500 600 700 800 900 1000 r [m] 10 1 ρ [m-2] ρ [m-2] ρ [m-2] (a) 10 1 10-1 0 1 10-1 100 200 300 400 500 600 700 800 900 1000 r [m] 10 0 10-1 100 200 300 400 500 600 700 800 900 1000 r [m] 0 (b) Figure 3.5. Low-energy iron-primary muon LDF. The behavior of the simulated muon density as a function of r⊥ for iron primaries of energy (a) E = 1016 eV, and (b) E = 1017 eV. Within each set at a given energy, results are shown (from left to right) for zenith angles cos θ = (0.75, 0.85, 0.95). Points with error bars are AIRES output (mean and standard deviation of 10 runs). The solid curve overlay shows the muon LDF parameterization defined in equation (3.19). 102 ρ [m-2] ρ [m-2] ρ [m-2] 38 104 103 103 10 102 1 10 10-1 1 -2 -1 102 10 1 10 10 10 10-3 10-4 -5 10 10-1 -2 0 500 1000 1500 2000 2500 3000 3500 10-2 10-3 10 10-4 10-4 4000 r [m] -3 0 500 1000 1500 2000 2500 3000 3500 4000 r [m] 0 500 1000 1500 2000 2500 3000 3500 4000 r [m] 500 1000 1500 2000 2500 3000 3500 4000 r [m] 500 1000 1500 2000 2500 3000 3500 4000 r [m] ρ [m-2] ρ [m-2] ρ [m-2] (a) 5 10 104 102 104 3 10 103 10 102 1 10 102 10 1 1 10-1 10-1 -2 -2 10-1 10-2 10 10-3 10 10-3 10-3 0 500 1000 1500 2000 2500 3000 3500 4000 r [m] 0 500 1000 1500 2000 2500 3000 3500 4000 r [m] 0 104 ρ [m-2] ρ [m-2] ρ [m-2] (b) 6 10 105 5 10 10 3 104 104 102 103 103 10 102 1 10 10-1 1 102 10 1 -1 10 10-2 10-2 10-2 -3 10 10-1 0 500 1000 1500 2000 2500 3000 3500 4000 r [m] 10-3 0 500 1000 1500 2000 2500 3000 3500 4000 r [m] 0 (c) Figure 3.6. High-energy proton-primary electron LDF. The behavior of the simulated electron/positron density as a function of r⊥ for proton primaries of energy (a) E = 1018 eV, (b) E = 1019 eV, and (c) E = 1020 eV. Within each set at a given energy, results are shown (from left to right) for zenith angles cos θ = (0.55, 0.75, 0.95). Points with error bars are AIRES output (mean and standard deviation of 10 runs). The solid curve overlay shows the electron LDF parameterization defined in equation (3.20). 10 ρ [m-2] ρ [m-2] ρ [m-2] 39 102 10 10 1 1 1 10-1 10-1 10-2 10-2 10-3 10-3 0 500 1000 1500 2000 2500 3000 3500 4000 r [m] 10-1 10-2 10-3 0 500 1000 1500 2000 2500 3000 3500 4000 r [m] 0 500 1000 1500 2000 2500 3000 3500 4000 r [m] 500 1000 1500 2000 2500 3000 3500 4000 r [m] 500 1000 1500 2000 2500 3000 3500 4000 r [m] 102 ρ [m-2] ρ [m-2] ρ [m-2] (a) 103 102 102 10 10 1 10-1 10 1 1 10-1 10-1 10-2 10-2 10-2 0 500 1000 1500 2000 2500 3000 3500 4000 r [m] 0 500 1000 1500 2000 2500 3000 3500 4000 r [m] 0 3 10 104 ρ [m-2] ρ [m-2] ρ [m-2] (b) 3 10 102 10 103 102 102 10 10 1 1 1 10-1 10-1 10-1 0 500 1000 1500 2000 2500 3000 3500 4000 r [m] 0 500 1000 1500 2000 2500 3000 3500 4000 r [m] 0 (c) Figure 3.7. High-energy proton-primary muon LDF. The behavior of the simulated muon density as a function of r⊥ for proton primaries of energy (a) E = 1018 eV, (b) E = 1019 eV, and (c) E = 1020 eV. Within each set at a given energy, results are shown (from left to right) for zenith angles cos θ = (0.55, 0.75, 0.95). Points with error bars are AIRES output (mean and standard deviation of 10 runs). The solid curve overlay shows the muon LDF parameterization defined in equation (3.21). 40 3.3 Time Distribution Function In a ground array, the angle of incidence of a cosmic ray shower is determined by fitting the relative particle arrival times at sites in the array to the shape of the shower front. This is complicated by the fact that the particle front of a cosmic ray air shower is not planar, but rather curves back from the center of the shower. In addition, the width of the particle front varies; it is narrow close to the shower axis and wider toward the edges. In general, the front edge of the shower is a steep rise to a maximum particle intensity followed by a longer tail of particles trailing behind the shower front. The time distribution function (TDF) describes the time (relative to a plane perpendicular to the shower axis) at which particles in the shower front will reach a detector at a given distance from the core. There is no well-motivated model for the TDF similar to the NKG formula for the LDF. The TDF developed by AGASA was obtained experimentally and was originally used as a first approximation to the CHICOS TDF. A complete description of the CHICOS-specific TDF has since been developed by fitting a parameterized function to AIRES-generated simulated showers. 3.3.1 AGASA TDF AGASA divided the time distribution function into two separate parts: the average time delay, Td , due to curvature of the shower front, and the average time spread, TS , which characterized the width of the shower front [55]. The average time delay, Td , from a plane perpendicular to the shower axis, at given distance from the core, is given by ³ r ´1.5 ρ(r)−0.5 ns, Td (ρ, r) = 2.6 1 + 30 where r is in meters. (3.22) 41 The AGASA formula for the width of the shower front, Ts , is given by ³ r ´1.5 Ts (ρ, r) = 2.6 1 + ρ(r)−0.3 ns. 30 (3.23) The time-delay and time-spread formulae were modified for CHICOS by removing the ρ(r) term in Td , and replacing ρ(r)−0.3 with ρ(r)−0.5 (i.e., pure counting statistics) in Ts . This was done because the CHICOS detectors have a much shorter time constant than AGASA detectors; hence the CHICOS detectors can generally resolve individual particles (sufficiently far from the core of the shower), while AGASA measurements integrated all particles in the shower front in a single pulse. The equations for Td and Ts in their original form describe the time delay and spread of the first particle to hit the detector, whereas it is more appropriate for CHICOS to use the average time delay and overall spread of all incident particles. 3.3.2 CHICOS TDF The AGASA TDF was designed to be used with chi-square fit methods. Such parameterizations have traditionally taken the form of a time delay function combined with a Gaussian uncertainty in the arrival time of particles within the shower front. Detailed shower simulations show that this is not an accurate model on timescales measureable by CHICOS; the shape of the particle distribution within the shower front is decidedly non-Gaussian, with a steep initial rise and a broad tail. The greater resolution of CHICOS hardware makes it more appropriate and desirable to use a maximum likelihood method in conjuction with a more complete description of the time distribution at all distances from the shower core. The CHICOS TDF, P (t; r⊥ , E, θ) describes the distribution of particles hitting the ground as a function of time at a given distance, r⊥ , from the core of the shower. As with the lateral distribution function, we have derived separate models for the electron and muon TDFs. The AIRES simulations used in this process is the set of high-energy, proton-primary showers used to construct the high-energy LDF. In the case of the TDF, however, it was observed that the shape of the arrival time 42 distribution has very little dependence on energy; thus the set of simulations was averaged over energy before proceeding. For each species (muons and electrons), the AIRES simulations were used to fill histograms of particle intensity as a function of r⊥ and t, using 50-m and 50-ns bins, respectively. The histograms were averaged over the 10 runs at each energy and zenith angle and the standard deviation of the runs was used as the uncertainty in the histogram. Because the scintillator detectors do not distinguish between electrons and muons, the two TDFs must be combined as a weighted average before being compared to a given measurement. Thus P (t) = ρe (r⊥ )Pe (t) + ρµ (r⊥ )Pµ (t) , ρe (r⊥ ) + ρµ (r⊥ ) (3.24) where ρe,µ (r⊥ ) is the appropriate electron or muon LDF. It has been found that the electron arrival time has a structure that is best fit as the sum of two similar distributions: a narrow curve with a fast rise time superimposed on a slower, broader one. The two curves have the same parameterized form and are related by the weight function w: Pe (t) = wPe1 (t) + (1 − w)Pe2 (t), (3.25) w = 1 − exp [(0.0021 − 0.0043 cos θ) r⊥ ]. (3.26) where The form of the curves Pe1 and Pe2 , including the curvature delay, are given by Pe1,e2 (t; r⊥ , E, θ) =   Ne1,e2 (t − ae1,e2 )be1,e2 exp[−ce1,e2 (t − ae1,e2 )], t ≥ ae1,e2 ,  0, t < ae1,e2 . (3.27) It has been found that, over the energy range observable by CHICOS, the parameters a, b, and c can be described as functions of r⊥ and θ only. The result is an 43 energy-independent time distribution formula that is presumed to be valid over the entire energy range of CHICOS data. The parameters of the electron TDF are as follows: ¡ ¢ 2 ae1 = 1.3 × 107 cos θ − 1 × 10−8 r⊥ ae2 = ae1 be1 = 1.5 cos θ − 0.825 + r⊥ exp (−57.38 cos3 θ) be2 = 10 be1 ³ r ´−1.5(cos θ−0.13) ⊥ ce1 = 0.8 2200 (3.28) ce2 = 10 ce1 µ be1 +1 ¶ ce1 Ne1 = Γ(be1 + 1) µ be2 +1 ¶ ce2 Ne2 = Γ(be2 + 1) The muon TDF has been found to have a simpler description. Only one term is necessary, with Pµ (t; r⊥ , E, θ) =   Nµ (t − aµ )bµ exp[−cµ (t − aµ )], t ≥ aµ ,  0, (3.29) t < aµ . The parameters of the muon TDF are as follows: ¡ ¢ 2 aµ = 1.3 × 107 cos θ − 2.7 × 10−8 r⊥ bµ = 1.8 cos θ − 0.000175r⊥ cµ = 0.68 + exp (6.1 − 3.58 cos θ − 0.0011r⊥ ) à ! b +1 cµµ Nµ = Γ(bµ + 1) (3.30) All TDFs are normalized to integrate to 1 (with time in units of microseconds). The TDF is considered valid out to zenith angles of approximately 60◦ and distances 44 r⊥ of approximately 5000 km. Figure 3.8 shows the behavior of the electron TDF over a range of zenith angles and radial distances from the shower core, compared with AIRES simulations of particle arrival time distributions. Figure 3.9 shows the same series of plots for the muon component of the showers. The insensitivity of the electron and muon TDFs to the energy of the shower primary is demonstrated in figure 3.10 and figure 3.11. These figures show simulated arrival time histograms at a given energy, zenith angle, and distance from the core, with the energy-invariant muon or electron TDF superimposed on the AIRES output. Results for other values of the zenith angle and distance from the core show similar agreement. While the particular parameterization of the TDF presented here is specific to the CHICOS altitude and experimental design, the form of the TDF should be applicable to the analysis of cosmic ray data from other ground array experiments. P [µ s-1] P [µ s-1] 45 6 2.2 2 1.8 5 1.6 4 1.4 1.2 3 1 0.8 2 0.6 1 0.4 0.2 0 0 1 2 3 4 0 0 5 t [µ s] 1 2 3 4 5 t [µ s] 1 2 3 4 5 t [µ s] 1 2 3 4 5 t [µ s] P [µ s-1] P [µ s-1] (a) 1.8 1.6 1 0.8 1.4 1.2 0.6 1 0.8 0.4 0.6 0.4 0.2 0.2 0 0 1 2 3 4 0 0 5 t [µ s] P [µ s-1] P [µ s-1] (b) 1.2 0.6 0.5 1 0.4 0.8 0.3 0.6 0.2 0.4 0.1 0.2 0 0 1 2 3 4 0 0 5 t [µ s] (c) Figure 3.8. High-energy proton-primary electron TDF. Sample electron/positron arrival-time distributions, averaged over primary energy, are shown for zenith angle (a) cos θ = 0.55, (b) cos θ = 0.75, and (c) cos θ = 0.95. Arrival-time distributions are shown (from left to right) at r⊥ = (1.0 km, 2.0 km). Points with error bars are AIRES output (mean and standard deviation of 10 runs at each of 6 energies). The solid curve overlay shows the electron TDF defined in equation (3.28). P [µ s-1] P [µ s-1] 46 6 2.4 2.2 2 1.8 5 1.6 4 1.4 1.2 3 1 0.8 2 0.6 1 0.4 0.2 0 0 1 2 3 4 0 0 5 t [µ s] 1 2 3 4 5 t [µ s] 1 2 3 4 5 t [µ s] 1 2 3 4 5 t [µ s] P [µ s-1] P [µ s-1] (a) 3.5 3 1.2 1 2.5 0.8 2 0.6 1.5 0.4 1 0.2 0.5 0 0 1 2 3 4 0 0 5 t [µ s] P [µ s-1] P [µ s-1] (b) 2 0.8 0.7 0.6 1.5 0.5 1 0.4 0.3 0.5 0.2 0.1 0 0 1 2 3 4 0 0 5 t [µ s] (c) Figure 3.9. High-energy proton-primary muon TDF. Sample muon arrival-time distributions, averaged over primary energy, are shown for zenith angle (a) cos θ = 0.55, (b) cos θ = 0.75, and (c) cos θ = 0.95. Arrival-time distributions are shown (from left to right) at r⊥ = (1.0 km, 2.0 km). Points with error bars are AIRES output (mean and standard deviation of 10 runs at each of 6 energies). The solid curve overlay shows the muon TDF defined in equation (3.30). P [µ s-1] P [µ s-1] 47 1.8 1.6 2 1.4 1.2 1.5 1 0.8 1 0.6 0.4 0.5 0.2 0 0 0 1 2 3 4 5 t [µ s] 0 1 2 3 4 5 t [µ s] 1 2 3 4 5 t [µ s] 1 2 3 4 5 t [µ s] P [µ s-1] P [µ s-1] (a) 1.8 1.4 1.6 1.2 1.4 1.2 1 1 0.8 0.8 0.6 0.6 0.4 0.4 0.2 0.2 0 0 1 2 3 4 0 0 5 t [µ s] P [µ s-1] P [µ s-1] (b) 1.6 1.4 1.2 1 1.2 0.8 1 0.8 0.6 0.6 0.4 0.4 0.2 0.2 0 0 0 1 2 3 4 5 t [µ s] 0 (c) Figure 3.10. Energy invariance of the electron TDF. Sample electron/positron arrivaltime distributions are shown for zenith angle cos θ = 0.75 and energy (a) E = 1018 eV, (b) E = 1019 eV, and (c) E = 1020 eV. Within each set at a given energy, results are shown (from left to right) for r⊥ = (1.0 km, 2.0 km). Points with error bars are AIRES output (mean and standard deviation of 10 runs at each of 6 energies). The solid curve overlay shows the energy-invariant electron TDF defined in equation (3.28). P [µ s-1] P [µ s-1] 48 4 3.5 1.4 1.2 3 1 2.5 0.8 2 0.6 1.5 0.4 1 0.2 0.5 0 0 0 1 2 3 4 5 t [µ s] 0 1 2 3 4 5 t [µ s] 1 2 3 4 5 t [µ s] 1 2 3 4 5 t [µ s] P [µ s-1] P [µ s-1] (a) 3.5 3 1.2 1 2.5 0.8 2 0.6 1.5 0.4 1 0.2 0.5 0 0 0 1 2 3 4 5 t [µ s] 0 P [µ s-1] P [µ s-1] (b) 3.5 3 1.2 1 2.5 0.8 2 0.6 1.5 1 0.4 0.5 0.2 0 0 1 2 3 4 0 0 5 t [µ s] (c) Figure 3.11. Energy invariance of the muon TDF. Sample muon arrival-time distributions are shown for zenith angle cos θ = 0.75 and energy (a) E = 1018 eV, (b) E = 1019 eV, and (c) E = 1020 eV. Within each set at a given energy, results are shown (from left to right) for r⊥ = (1.0 km, 2.0 km). Points with error bars are AIRES output (mean and standard deviation of 10 runs at each of 6 energies). The solid curve overlay shows the energy-invariant muon TDF defined in equation (3.28). 49 Chapter 4 The CHICOS Experiment The CHICOS experiment consists of a ground array of scintillator detectors networked to data-processing computers at Caltech. The first CHICOS sites were deployed in 2001, and data collection of cosmic ray air showers has continued since 2003. The design of the array and the detector hardware are discussed in section 4.1. Section 4.2 discusses the data handling and shower search methods. 4.1 The Detector Array The CHICOS array was originally envisioned to consist of 90 sites located in the San Gabriel and San Fernando valleys. The maximum nuber of sites simultaneously installed since the project’s inception has been 77, in 2005 (figure 4.1). The array covers a total area of approximately 400 km2 . The San Gabriel valley contained as many as 44 sites spread over an area of about 150 km2 . Included in the San Gabriel array is the set of 12 pairs of detectors on or near the Caltech campus that make up the Chiquita array. The sparser San Fernando array comprised a maximum of 33 sites over an area of approximately 250 km2 . Array size and reporting statistics are shown in figures 4.2 and 4.3. See appendix A for site locations and parameters. The Chiquita array consists of the detectors at Pasadena City College and the Polytechnic School, as well as 5 sites on the Caltech campus. (One of the Caltech sites is a closely spaced set of 6 individual pairs.) The Chiquita array is sensitive to lower-energy events than the larger CHICOS array; the range of sensitivity lies 50 y [m] CHICOS Array (July 2005) 20000 15000 10000 5000 0 -5000 -10000 -50000 -40000 -30000 -20000 -10000 0 10000 x [m] Figure 4.1. Sites in the CHICOS array as of July 2005. The origin is placed on one of the Caltech sites, Caltech 03, at longitude −118.12522◦ and latitude 34.13657◦ . The western group of sites is located in the San Fernando valley; the eastern group of sites is located in the San Gabriel valley. Los Angeles is located to the south. approximately between 1016 eV and 1019 eV. Each site in the array is equipped with two scintillation detectors, a GPS antenna and receiver (Motorola UT+ or M12), and a workstation with two National Instruments PCI-6602 data acquisition cards (figure 4.4). The scintillator detectors were donated to the CHICOS project by the Los Alamos National Laboratory, where they were previously used in the CYGNUS cosmic ray project. Each detector consists of a sheet of plastic scintillator approximately 1 m2 in area and 5–10 cm thick housed in a lightproof fiberglass enclosure (figure 4.5). The top of the unit contains either a 3” or 5” photomultiplier tube. Some of our photomultiplier tubes were recycled from the Palo Verde Neutrino Oscillation Experiment, while others were purchased new. The effective area of the scintillator has been measured and found to be at least 90% of the real area. A time-over-threshold discriminator circuit (figure 4.6), built at Caltech, measures the length of the exponential PMT pulse from each detector. The decay constant, τ , of the photocurrent is approximately 80 ns. This decay constant is significantly greater than the time spread of the scintillator pulse, so the output can be accurately modeled as an exponentially decaying pulse, V (t) = V0 e−t/τ . (4.1) 51 Fraction of Sites Reporting (2003) 90 Sites Reporting/Installed Sites Number of Installed Sites Installed Sites (2003) 80 70 60 50 40 30 20 50 100 150 200 Sites Reporting/Installed Sites Number of Installed Sites 0.5 0.4 0.3 50 100 150 200 250 300 350 Day (Jan 1, 2003 = 1) 200 250 300 350 Day (Jan 1, 2004 = 1) 200 250 300 350 Day (Jan 1, 2005 = 1) Fraction of Sites Reporting (2004) 80 70 60 50 40 30 20 1 0.9 0.8 0.7 0.6 0.5 0.4 0.3 0.2 10 0.1 50 100 150 200 0 0 250 300 350 Day (Jan 1, 2004 = 1) Installed Sites (2005) 50 100 150 Fraction of Sites Reporting (2005) 90 Sites Reporting/Installed Sites Number of Installed Sites 0.6 0 0 250 300 350 Day (Jan 1, 2003 = 1) 90 80 70 60 50 40 30 20 1 0.9 0.8 0.7 0.6 0.5 0.4 0.3 0.2 10 0 0 0.7 0.1 Installed Sites (2004) 0 0 0.8 0.2 10 0 0 1 0.9 0.1 50 100 150 200 250 300 350 Day (Jan 1, 2005 = 1) 0 0 50 100 150 Figure 4.2. Array size and reporting statistics, 2003–2005. Information regarding the health and configuration of the array is updated daily based on data quality and diagnostic information received from the array sites. 52 Fraction of Sites Reporting (2006) 90 Sites Reporting/Installed Sites Number of Installed Sites Installed Sites (2006) 80 70 60 50 40 30 20 0.8 0.7 0.6 0.5 0.4 0.3 0.2 10 0 0 1 0.9 0.1 50 100 150 200 250 300 350 Day (Jan 1, 2006 = 1) 0 0 50 100 150 200 250 300 350 Day (Jan 1, 2006 = 1) Figure 4.3. Array size and reporting statistics, 2006–2007. Information regarding the health and configuration of the array is updated daily based on data quality and diagnostic information received from the array sites. 53 Detector A GPS Antenna Detector B Roof Classroom with data station CEU Figure 4.4. Diagram of a CHICOS array site. Each site is equipped with two detectors and a GPS unit on the roof, and a workstation for data acquisition in a nearby science classroom. 21.0" 20.75" 52.5" 51.5" 24.0" 24.25" 3.0" 4.5" 2.75" 3.75" Square Detector Round Detector 4.5" 13" 41.5" 44.5" 13" 42.25" 46" Figure 4.5. Diagram of CHICOS detectors. Detectors may have either a round or a square base. The total scintillator area in each type is approximately 1 m2 . 54 Figure 4.6. CEU schematic. The CEU is the CHICOS Electronics Unit, which controls the voltages applied to the detector PMTs, and contains a discriminator circuit through which the detector signal is passed. As described in section 6.2, the number of particles passing through the detector is measured in units of the energy, EMIP , deposited by a minimum ionizing particle, defined as a muon propagating through the scintillator at the average zenith angle. The total energy deposited in the scintillator is proportional to the integral of the PMT pulse. A discriminator threshold is applied to the pulse, and the time over threshold tth is measured. The discriminator voltage is typically set to −10 or −15 mV. The energy in MIP is related to the time over threshold by 100 EMIP = Cetth /τ , (4.2) where C is the pulse-height calibration constant. The precise values of C and τ vary 55 Figure 4.7. Sample decay-constant calibration. The decay constant τ is measured by collecting pulse traces and performing a linear fit to the logarithm of the pulse integral as a function of the time over threshold. A sample calibration is shown for site 1001, detectors A (left) and B (right). from site to site, depending on the hardware (see appendix A). The decay constant τ is measured by collecting pulse traces and performing a linear fit to the logarithm of the pulse integral as a function of the time over threshold (figure 4.7). The pulseheight constant C is then set such that the mean of the distribution of pulse energies C tth /τ e is approximately 1 MIP (figure 4.8). 100 There is evidence that the detector response varies with temperature. This has been quantified by fitting the pulse-energy histogram with a Landau distribution curve, and measuring the correlation between the most probable value of the distribution with the daily temperature reported by Burbank airport. This analysis allowed a correction factor to be retroactively applied to the detector calibration. The time over threshold is measured by two NI PCI-6602 Timer/Counter cards. These data acquisition cards compare the discriminator pulse with an 80 MHz oscillator, for an internal timing precision of approximately 12.5 ns. Variation in the oscillator frequency is adjusted for on a second-by-second basis. 56 Mip: Site 1006B, 10/15/2005 hour 0 Mip: Site 1006B, 10/15/2005 hour 0 Entries Mean RMS ×10 2 1200 40 97.91 67.67 1000 800 600 400 200 0 2 10 10 3 4 10 Measured Charge Figure 4.8. Sample pulse-energy histogram. Charge is measured in 1/100 MIP. The calibration constant C is defined such that the mean of the pulse-energy histogram is 1 MIP. This histogram represents 1 hour of accumulated data. Clocks at widely separated locations are kept synchronized by the GPS signal. When the GPS receiver is used set in “stationary” mode (i.e., assuming fixed location), the accuracy of the timing pulse used to resynchronize the clock is accurate to better than ±50 ns. The receiver uses the Time Receiver Autonomous Integrity Monitoring (T-RAIM) algorithm, which allows it to isolate and remove faulty satellite data from its time calculation. A minimum of 3 satellite signals is required to maintain timing accuracy. No data are recorded if there are fewer satellites in range. Each time the data acquisition software is restarted, the GPS receiver is reinitialized with the settings listed in Tables 4.1 and 4.2. At the beginning of 2004, all existing sites were upgraded. A revised version of the data collection software was installed, some hardware improvements were made, and all detectors were recalibrated. 57 Table 4.1. GPS settings for M12 receivers Command Name Setting @@Cj @@Ge @@Gf @@Gc @@Aq @@Aw @@Gb @@Ay @@As @@Gd @@Az @@Ha Poll Receiver T-RAIM Select T-RAIM Alarm Limit PPS Control Atmospheric Mode Time Mode GMT Offset 1-PPS Time Offset Position-Hold Position Position Control 1-PPS Cable Delay Position/Status/Data <none> T-RAIM on Alarm limit 1400 ns PPS on continuously Ionospheric model only UTC 0 hours 0 ns <site coordinates> Position-hold enabled <site cable delay time> Update once per second Table 4.2. GPS settings for UT+ receivers Command Name Setting @@Cj Poll Receiver @@En T-RAIM Setup/Status Message @@Aq @@Aw @@Ab @@Ay @@As @@At @@Az @@Ea Atmospheric Mode Time Mode GMT Offset 1-PPS Time Offset Position-Hold Position Position-Hold Option 1-PPS Cable Delay Position/Status/Data <none> T-RAIM on Alarm limit 1400 ns PPS on continuously Ionospheric model only UTC 0 hours 0 ns <site coordinates> Position-hold enabled <site cable delay time> Update once per second 58 4.2 Data-Collection Software The pair of detectors at each site in the array are separated by approximately 3– 5 meters. Each detector records approximately 200 single-particle events per second. Hits in both detectors are considered “coincident” if they are separated by 50 ns or less. A coincidence event is considered a “trigger” if both hits have an intensity equal to two minimum ionizing particles, or 2 MIP. All singles events at each detector are written to a daily file. In addition, all triggers recorded at a site are written to a daily trigger file. Each site also records a number of history files containing diagnostic information, including: • minute-by-minute average singles rate per second for each detector • minute-by-minute coincidence rate • minute-by-minute trigger rate • minute-by-minute oscillator frequency on each data acquisition card • minute-by-minute average number of satellites tracked by the GPS unit • hour-by-hour accumulation of pulse-energy histogram for each detector • an alarm log file noting singles, coincidence, or trigger rates outside of expected ranges • a log file noting software restarts and changes in GPS status • a summary file of site parameters and reporting statistics The software that coordinates the data collection at each site is written in Labview. It is designed to have a user-friendly interface (in order to be accessible to teachers and students) as well as to operate with minimal oversight. The workstations are designed to take data continuously, without any human intervention; data transfer and recovery from hardware glitches are done automatically. Several views of the 59 Labview data acquisition software are shown in figure 4.9, figure 4.10, figure 4.11, and figure 4.12. Each data station maintains a 7-day archive of relevant data on site. This archive includes the raw data files as well as the history and log files. Since these are useful for diagnosing problems remotely, the file handler sends the daily history files with the first data transfer the following day at 12:45 am. These files are available through a web interface to assess the status of any site on any given date. Given the large quantity of data collected at every site (approximately 100 MB per detector per day), it would be impractical to transfer all the data files to Caltech every day. The following system of data handling is designed to minimize the amount of data that must be sent over the internet, while keeping the reliablity of data transfer as high as possible. Each day at 12:45 am, every site sends a file of timestamped trigger events via FTP or SCP to a server at Caltech. If not immediately successful, the local site will continue trying to send the file for up to 3 hours. At 4:00 am the trigger events from all sites in the array are combined into a “master trigger file” by an automatic shell script. At 5:00 am, the data stations download the master trigger file. Once the local site has downloaded the master trigger file, it compares the arraywide list of triggers with the local lists of single events for each detector. A single hit occurring within ±50 µs of a trigger is considered a “match.” The matching process typically takes about 45 minutes, after which the site uploads the match files to Caltech. At 8:20 am the match files from all CHICOS sites are automatically assembled into candidate shower files. At 9:00 am, the daily shower file is passed to event filtering software. This program selects candidate events that meet certain criteria. In order to be considered a “candidate shower,” a set of matching hits must include 3 of the 5 sites closest to the trigger site. Further, the trigger and at least one set of 3 matching hits must be fit by a plane with an RMS time residual of less than 10 µs. The set of candidate showers is designed to accept events liberally; the data set is pruned again during the analysis phase based on the quality of the reconstructor fits. 60 Figure 4.9. Labview data acquisition program front panel. The rate of individual hits per second is shown in the upper left. The rates of coincidences and triggers are shown in the upper right. The time difference between A and B detector coincident hits is is histogrammed in the lower right. Figure 4.10. Labview data acquisition program history panel. Information characterizing the site’s performance can be viewed for any day in the preceeding week. 61 Figure 4.11. Labview data acquisition program satellites panel. The number and relative strenth of GPS satellite signals is displayed. A minimum of 3 signals is required to maintain timing accuracy. Figure 4.12. Labview data acquisition program energy panel. The distributions of time over threshold (top) and the corresponding energy (bottom) of the hits recorded by each detector are histogrammed. 62 Chapter 5 Shower Reconstruction Software The reconstruction software used by the CHICOS project was largely written by Brant Carlson and Chris Jillings. It is based on the methods used by AGASA, and is implemented in C++. The C++ library, libCTShower, uses the ROOT class library developed by CERN [122]. Section 5.1 contains an overview of the data format and handling. The first version of the CHICOS reconstructor used a χ2 minimization method. More recently, this has been replaced with a maximum likelihood method, which is more appropriate for the type of data collected by CHICOS. These and other reconstruction options are described in section 5.2. An online user interface has been developed through which any shower data can be reconstructed with any of several versions of the reconstructor code. An interactive form allows the user to select or deselect data and to guide the reconstructor software. Numerical output is accompanied by a graphical representation of the shower data and reconstruction results. Details of the user interface are discussed in section 5.3. 5.1 Overview of libCTShower Raw data files are filtered daily before being transferred to Caltech, where trigger and match events from all sites are organized into showers as described in chapter 4. The shower is the basic unit of data that is given to the reconstruction software. The purpose of the reconstructor software is to extract an estimate of the orignal energy and trajectory of the primary cosmic ray associated with a given shower. 63 5.1.1 Data Format Shower information is stored in the form of daily shower log files. Each shower in the file has two header lines. The first contains the date and time (to the nearest nanosecond). The second line lists the number of trigger events, the number of hits (triggers plus matches), a shower ID (SID) number (which is unique within a given day but not in the data set as a whole), and a set of 8 flags, 7 of which are used. The meanings of the flags are as follows, beginning with the least significant (rightmost) bit: 1. At least 10 hits in the shower, no conditions. 2. At least 2 triggers in the shower. 3. At least 6 hits (including triggers) within 20 km of any trigger in the shower. 4. At least 4 sites (including the trigger site) within 20 km of any trigger in the shower. 5. Caltech 03 - Polytechnic - Pasadena High triplet. 6. LaSalle - Alverno - Pasadena High triplet. 7. Blair - Mayfield - Westridge triplet. The “triplet” flags identify showers that hit all three of a group of nearby sites. These flags been useful in the statistical analysis of the CHICOS data [123]. The data in each shower event is organized into 6 columns: site ID number, detector type (A or B), time, A detector intensity, B detector intensity, and site name. Since each line is associated with a hit in either an A or a B detector, only one of the intensity columns will be filled. The data of a sample shower is shown in figure 5.1. The first two lines of shower data are the hits that form the shower trigger. The hit in the A detector is listed first, identified by an uppercase “T” in the detector type column. The B detector hit is listed second, identified by a lowercase “t.” The 64 shower 01-12-07 08:10:13 PST UTC: 3251463013 ns: 629881825 Ntrig: 1 Nhit: 9 SID: 81013 Trig: 00001100 0063 T 0.000 2.11 WashingtonMS 0063 t -0.061 6.70 WashingtonMS 0010 A -11.683 0.83 SylmarHS 0065 A -3.031 1.09 ChandlerSchool 0131 B -2.249 0.68 NewHorizon 0103 A 2.923 0.80 LongFellowElem 0103 B 3.133 0.53 LongFellowElem 0044 B 3.969 0.72 EliotMS 0092 A 36.555 0.44 SanFernandoHS Figure 5.1. Sample shower data. The columns represent the site ID number, the detector label, the hit time (in milliseconds, relative to the A-detector trigger), the hit intensity (in MIP), and the site name. time of all events in the shower are normalized to the time of the A detector of the trigger, which is set to time t = 0. Times are listed in microseconds, and intensities in units of minimum ionizing particles (MIP), defined in section 4.2. 5.1.2 Code The libCTShower code contains a “shower” class that stores the information about the event being reconstructed. Member data of the shower class include an array of sites involved in the shower, an array of hits, and an array of triggers. Each site object in the shower’s site array is also associated with the hits and triggers at that site. Hence one can access the set of hits in a shower, triggers in a shower, hits at a site, or triggers at a site [124]. Each hit and trigger event is treated as object whose member data is the detector label (A or B), the hit intensity (in MIP), the time of the hit (in microseconds), and a pointer to the site object with which the hit is associated. Each site object stores information such as its (x, y, z) coordinates, the size and thickness of the scintillator in each detector, its intallation date and calibration results, as well as the set of hits and triggers associated with it for a given shower. The shower reconstruction is done by the “reconstructor” class. This class contains 65 the bulk of the code that calls the ROOT MINUIT minimization routine to fit the shower. The reconstructor class contains a “fitshower” object that stores the shower data. The fitshower is modified from the original shower object as needed by the reconstructor. For example, when using the χ2 minimization method, hits and triggers are averaged by site before being stored. The CIT 1000 array, which includes 6 individual pairs of detectors, is treated as one site and averaged together. When using the log likelihood minimization, no averaging is done and all hits are treated independently. Shower visualization is done by the “ui” class. This class contains tools for plotting the shower data and the reconstruction results. The user can interact with the reconstructor class through the user interface. Other classes include the “latdist” class, which holds all options for the form of the lateral distribution, and a “timedist” class, which holds the options for the time distribution function. An “array” class acts as a container of “array site” objects, and is used to model the state of the array on a given day. There is also a shower simulator class, which contains code for generating simulated showers. Simulated showers are based on the lateral distribution function and time distribution function, with statistical noise added. A shower analyzer class is provided for working with the data in ROOT. 5.2 Shower Reconstruction The libCTShower code makes use of the ROOT MINUIT minimization package. The parameters fit are core location (x, y, z), hit time of the shower core (t0 ), energy, and direction (θ, φ). Core location and hit time are not independent parameters; hence z is set to the average vertical height of sites involved in the shower and held fixed in the minimization. The angles θ and φ do not work well as parameters in the fit routine, for several reasons. First, the periodicity can cause the reconstruction to become confused. Second, when θ is near zero (a vertical shower), φ becomes poorly constrained. To avoid these problems, the angles θ and φ are mapped to Cartesian coordinates (xang , yang ) 66 by the introduction of a transform function f (xang , yang ). This is accomplished by defining a vector v with direction angles (θ, φ) and extending this vector until it intersects the surface of the transform function at the point (xang , yang , f (xang , yang )). Following each minimization step, it is necessary to reverse the transformation to obtain θ and φ from the coordinates xang and yang . The transformations between these coordinates and the shower angles are Ãp θ = arctan 2 x2ang + yang f (xang , yang ) and µ φ = arctan yang xang ! (5.1) ¶ . (5.2) The function f (xang , yang ) is chosen to be smooth and continuous over the region where it is defined. The function currently in use is q ¡ ¢ Nq 2 2 2 2 f (xang , yang ) = 1 + xang 1 + yang 2 − x2ang − yang + x2ang yang , 2 (5.3) with −1 < xang < 1 and −1 < yang < 1. The proportionality constant N determines the maximum angle of θ. A function of this form is constant along the boundaries of xang and yang , where it evaluates to p 2 √ 2 xang + yang 2 = . tan θ = f (xang , yang ) N (5.4) In other words, the maximum allowed vertical angle θ does not depend on the azimuthal angle φ. For N = 1, the maximum shower inclination is θ = 54.7◦ , and for p N = 2/3, the maximum shower inclination is θ = 60◦ . The value of N is chosen to be consistent with the LDF model being used, as different versions of the LDF were based on sets of simulations with different maximum inclinations. See figures 5.2 and 5.3 for illustrations of the transform function and its tangent. In the case where showers are being reconstructed using the AGASA LDF, the energy is not directly determined by the fit. Instead, the natural log of the proportion- 67 Figure 5.2. Transform function f (xang , yang ) used to map (θ, φ) to linear coordinates. Figure 5.3. Tangent of the vertical angle θ as a function of xang and yang . The maximum allowed θ is constant along the boundary of the allowed region. 68 ality constant C in equation (3.11) for ρ(r) is used as a fit parameter. The value of the lateral distribution function 600 m from the core can then be used in equation (3.14) to determine the energy of the incident particle. For all other CHICOS-specific LDFs, the proportionality constant in ρ(r) is explicitly a function of log(E), so the log of the energy can be fit directly. 5.2.1 Chi-Square Reconstruction Method The original version of the reconstruction code works by performing a χ2 minimization of the shower data to a model LDF parametrized by incident angle and energy. As mentioned above, this fit is done on a shower object in which the data has been averaged by site. First, multiple events at a single detector are summed, and the results are then averaged together (weighted by the area of the scintillator in each detector). Events with measured intensity greater than 1000 MIP are removed from the shower data before the average intensity is calculated, since it is known that the detector hardware is unreliable beyond that range. An error bar is assigned to the intensity measurement based on the variance of the data. Given a mean intensity µI = 1 NX events Nevents i=1 Ii , Ai (5.5) where A is the detector area, the spread is calculated as either v u u err1 = t 1 NX events µ Nevents i=1 Ii − µI Ai ¶2 (5.6) or q err2 = P events max(Ii , 1) Nevents × N i=1 PNevents Ai i=1 (5.7) 69 and the error is assigned to be max(err1, err2). (5.8) Since there is very often only one hit in each detector, a time error bar is calculated based on the theoretical spread of the shower front at the location of the detector, as given in equation (3.23). The procedure used to fit a set of shower data is as follows [124]: 1. The core location is given a first estimate of 200 m or 1000 m from the trigger site trigger site in one of the cardinal directions.a (The minimization procedure is repeated 8 times, with different starting points; the fit that produces the lowest χ2 is used.) 2. Each set of 3 hits in the shower can be fit with a plane perpendicular to the trajectory of the shower front. A first estimate of of the incident angle is found by averaging the results from each trio of events in the data. Trios that form a straight line are rejected, as are sets that contain a time difference greater than the light travel time between the sites. 3. No attempt is made at a first estimate of the energy. Instead, the logarithm of the proportionality constant C in ρ(r) is set to 10, which corresponds roughly to a starting energy estimate of 1018 eV. 4. The direction of propagation is fit by χ2 minimization of the hit times with the time distribution function, which may be either a CHICOS-specific TDF or the AGASA time delay and time spread functions. 5. The energy is fit by χ2 minimization of the hit intensities with the LDF, ρ(r). 6. The core location and the energy are fit simultaneously, again by χ2 minimization of the intensities given by LDF. a site. For Chiquita showers, the initial core location estimates are 40 m and 200 m from the trigger 70 7. Steps 4 through 6 are repeated until the fit converges, up to a maximum of 10 iterations. 8. The minimization procedure is repeated with a different initial core location and the best result of 8 is chosen. The reconstruction returns two separate χ2 values, one for the angle/time fit, and one for the energy/core location fit. The output of the reconstructor is given as the (x, y, z) coordinates of the core location in meters (with the origin centered on Caltech 03, at latitude 34.13657◦ and longitude −118.12522◦ ), the zenith angle θ and the azimuthal angle φ (measured counterclockwise from east), the time of impact t of the core relative to the first measured event, and the energy E. It is also possible to fit shower data by assigning initial values to the parameters, or forcing a reconstruction with a given set of parameters. 5.2.2 Maximum Likelihood Reconstruction Method The χ2 fit code has been superseded by a maximum likelihood method, which is better suited to the type of data collected by CHICOS. A χ2 fit is subject to loss of information when data is binned or averaged. It is also only strictly valid if the data follows a Gaussian distribution. This is not the case with data collected from a shower front, where the distribution of particles has a sharp rise and a long tail, more closely approximating a te−t distribution. The maximum likelihood reconstructor uses the MINUIT package to minimize the negative of the log of the likelihood function. The two detectors making up each array site are treated independently. The reconstruction procedure proceeds analogously to the χ2 method, with the exception that the energy of the shower enters directly as a parameter, and is given the initial value of log(E) = 1018 eV. The log likelihood function is given by log Ptotal (x, y, z, t, xang , yang , log(E)). The likelihood function Ptotal is defined to be the product of the likelihood functions for each site in the array, which in turn are considered to be the product of the likelihood functions for each of the two detectors within the site. Hence the overall likelihood 71 function is the product of the likelihood functions for all individual detectors in the array: Ptotal = Y Pdetector . (5.9) detectors The total log likelihood function is the sum of the individual log likelihood functions: log Ptotal = X log Pdetector . (5.10) detectors In the current version of the likelihood minimization method, the likelihood is calculated as the product of two probability functions: the likelihood of measuring the total intensity at a given detector, PI (Idetector ), and the likelihood of measuring each individual event at the observed time, PT (tevent ). Thus Pdetector = PI (Idetector ) × Y PT (tevent ), (5.11) events where Idetector = X Ievent (5.12) events and the set of events under consideration is the set of individual hits at the given detector. We can now substitute equation (5.11) into equation (5.9) to obtain the total likelihood function Ptotal = Y PI (Idetector ) × detectors Y PT (tevent ) (5.13) events and the total log likelihood function X log Ptotal = detectors log PI (Idetector ) + X log PT (tevent ). (5.14) events The set of events considered in equations (5.13) and (5.14) is the entire set of individual hits across the array. The likelihood function PT (tevent ) is given by the integral of the TDF over an inter- 72 val surrounding tevent , plus a term that accounts for the background noise (accidental hits across the array, unrelated to the shower): Z tevent −t∆ PT (tevent ) = (PTDF (t0 ) + Background(t0 )) dt0 . (5.15) tevent +t∆ The background noise is assumed to be a constant 200 MIP/s in each detector, hence Z tevent −t∆ PT (tevent ) = PTDF (t0 )dt0 + 2(2 × 10−4 µs−1 )t∆ . (5.16) tevent +t∆ The actual noise background varies slightly between detectors because of differing area and depth of scintillator, and pulse threshold sensitivity. The interval of integration 2t∆ is currently set to 25 ns. The likelihood function PI (Idetector ) is given by the probability of measuring the total intensity from a Poisson distribution about the total expected intensity ILDF : PI (Idetector ) = PPoisson (Idetector ; ILDF ). (5.17) When the total measured intensity at a site is greater than 10 MIP, this function is modified to reflect the greater uncertainty in the measurement: PI (Idetector > 10 MIP) = 1 − exp (−ILDF ). (5.18) This allows large hits to fit to large expected intensities, even if the exact values differ. An alternative method for calculating Pdetector uses the K-S test, which calculates the Kolmogorov-Smirnov probability of obtaining the observed time distribution of hits at a given detector from the time distribution function at the detector’s location [125]. The detector likelihood function is Pdetector = PI (Idetector ) × PKS (t1 , ..., tn ), (5.19) where t1 , ..., tn is the set of time measurements for the hits at that detector. This 73 gives the following expression for the total likelihood function: Ptotal = Y PI (Idetector ) × detectors X log Ptotal = Y PKS (t1 , ..., tn ), (5.20) detectors log PI (Idetectors ) + detectors X log PKS (t1 , ..., tn ) (5.21) detectors The probability PKS is a function of the quantity D, which is the maximum difference between two cumulative distributions, in this case the total charge measured at the detector before time t, and the integrated TDF: D = Max |SN (t) − P (t)| over − ∞ < t < ∞. (5.22) The quantity SN (t) is the sum of hits at the detector before time t, normalized by the total charge measured at the detector: SN (t) = X Ievent (tn )/Idetector . (5.23) tn <t The quantity P (t) is the integral of the TDF from 0 to t: Z t P (t) = PTDF (t0 , r, θ, log E)dt0 . (5.24) 0 The probability of D being greater than the observed value is given by PKS (D > observed) = QKS ³h√ √ i ´ N + 0.12 + 0.11/ N D (5.25) where N is the number of data points, and ∞ X 2 2 QKS (λ) = 2 (−1)j−1 e−2j λ . (5.26) j=1 The function that calculates QKS in ROOT calculates the number of terms to sum based on the value of λ. 74 Because the integral of measured charge at a site over time is a step function, the maximum value of D will occur either immediately before or immediately after a hit. These are the only times checked by the code, to avoid having to step through all values of t. The K-S test was not found to produce a significant improvement in the quality of shower reconstructions, so it is not the default method for calculating the likelihood. 5.3 User Interface The CHICOS reconstruction software is integrated with an online user interface. An online search function is available to filter showers based on date, time, size, or sites involved. Once a shower is loaded, the online reconstructor allows the user to select which data points to include in the fit. This is useful for removing data points that are temporally associated with the shower window, but are likely to be noise hits. In addition a control string may be fed to the reconstructor to specify starting values of the parameters, special reconstruction methods, or sites to be ignored in the array.b Each time a shower event is loaded into the reconstruction software, the site data files and history files listed in section 4.2 are used to assess the “health” status of each site that was active in the array at the time of the shower (to the nearest second). Criteria for identifying a site as “healthy” during a shower event include: • the existence of a trigger file and both A and B match files, each with a minimum size of 500 bytes • A and B match files that each contain at least one entry during a 30-minute time window centered on the shower Sites that did not transmit their data that day, or that were not recording data at the time of the shower, are designated “unhealthy.” Such sites are shown on the shower b Removing a data point by de-selecting it from the list changes the measured intensity at the site, although the site is still included in the array. By contrast, the “ignore” function tells the software to delete a site – and all its data – from the array completely. 75 display as yellow triangles labeled with the site number. If such sites are not manually ignored, the reconstructor will treat them as sites that measured zero intensity.c The output of the reconstructor is a list of the estimated parameters, their estimated errors, and a graphical representation of the fit. (See figures 5.4, 5.5, and 5.6.) A map of the array is plotted, showing all sites in the array (black dots) near the shower core. Sites involved in the shower are shown as circles, color coded according to hit time, from red (earliest) to blue (latest). Sites not involved in the shower, but that have a best-fit intensity of 0.1 MIP or greater are circled in gray. The estimated core location of the shower is shown with an X surrounded by an error box. The estimated direction and uncertainty is indicated by a wedge. The second plot in the reconstructor output is an overlay of the lateral distribution function on the data. The third plot shows measured intensity vs. time. Individual data points are shown as solid circles, the site average is shown as an X with an associated error box, and the estimated fit value is shown as an open circle. It is important to note that all of our shower reconstruction software is made available to participating schools via the CHICOS website. It is part of the philosophy of the project that both data and analysis tools be made available for educational purposes. The reconstruction tools are accompanied by detailed documentation as well as a tutorial aimed at high-school students written by Theresa Lynn. c Rarely it happens that only one detector is malfunctioning while the other continues recording data. In this instance, a site may record a shower hit yet still be marked as unhealthy. In this case it is at the discretion of the user to determine whether the hit at that site should be removed from the data. All such events have been removed from the data sets being considered in this paper in order to simplify the task of modeling the data with simulated air showers. 76 Figure 5.4. Sample output of shower reconstructor on raw shower data. 77 Figure 5.5. Output of shower reconstructor after removal of accidental hits. 78 Figure 5.6. Sample output of shower reconstructor after removal of accidental hits and removal of inactive sites from the array. 79 Chapter 6 Modeling the Array Response This chapter describes the method used to characterize the response of the Chiquita subarray to cosmic-ray air showers at energies close to 1017 eV. This work was undertaken in order to assess the accuracy of the CHICOS shower-reconstruction software. Unthinned shower simulations are combined with a detailed model of the detector response to produce simulated data. The simulated data, when fed back into the reconstruction software, allows us to gauge any bias in the distribution of reconstructed energies. The generation of unthinned simulated air showers is discussed in section 6.1. Section 6.2 describes the details of the detector response model, and section 6.3 describes the analytical and numerical methods used to model particle interactions with the detector scintillator. The results of this exercise are given in section 6.4. 6.1 Simulation of Air Showers Using AIRES The simulated air showers used in this analysis were created with AIRES, a simulation program made available by the Universidad Nacional de La Plata, Argentina [108]. Air showers were generated with energies of 1017 eV and 1017.5 eV, and zenith-angle cosines of 0.60, 0.70, 0.80, and 0.90. The depth of first interaction and the polar angle were randomly generated for each shower. No particle thinning was used. Because of the computation time involved in generating unthinned high-energy air showers, only 2 or 3 showers were generated at each energy and zenith angle (ZA). 80 Table 6.1. Parameters of the AIRES showers used in this work Energy Cosine ZA Polar Angle Depth of First Interaction Reference Number 1e17 eV 1e17 eV 1e17 eV 1e17 eV 1e17 eV 1e17 eV 1e17 eV 1e17 eV 1e17 eV 1e17 eV 1e17 eV 0.60 0.60 0.70 0.70 0.70 0.80 0.80 0.80 0.90 0.90 0.90 12.38◦ 97.04◦ 121.88◦ -20.52◦ -119.68◦ 155.76◦ 73.38◦ 2.74◦ -131.64◦ -90.42◦ 62.86◦ 4.5451 g/cm2 27.1158 g/cm2 10.9821 g/cm2 4.9164 g/cm2 14.9705 g/cm2 0.2366 g/cm2 0.1710 g/cm2 22.9553 g/cm2 7.5142 g/cm2 6.0402 g/cm2 16.1197 g/cm2 E17-ZA60 01 E17-ZA60 04 E17-ZA70 01 E17-ZA70 05 E17-ZA70 06 E17-ZA80 01 E17-ZA80 05 E17-ZA80 06 E17-ZA90 01 E17-ZA90 06 E17-ZA90 07 Only showers generated at 1017 eV were used in this work. The shower parameters for the data set are shown in Table 6.1. To simulate CHICOS data, each AIRES shower is dropped on the array, and the number of particles landing on each detector is recorded. Each AIRES shower can be used multiple times by shifting the core location relative to the array and rotating the polar angle of the shower. The core locations have been chosen randomly over an area of 1 km2 , centered on the Chiquita subarray. Only Chiquita sites are included in the array used to reconstruct the showers. The exact area of the scintillator in each detector is drawn from the site database (see appendix A). The relative locations of the sites are known to within several meters. Each site consists of two detectors, which are modeled as being separated by 2.8 m in an arbitrary direction. The difference in measured intensity between the two detectors at a site should be negligible; the separation is necessary to obtain independent statistics. 81 6.2 Modeling the Detector Response CHICOS uses plastic-scintillator detectors, with an area of approximately 1 m2 . The output pulse of the photomultiplier tube (PMT) is characterized by a time constant τ . The pulse is routed through a discriminator circuit, which measures the time over threshold. The relation between the time over threshold and the energy deposited in the scintillator is calibrated for each detector, as described in section 4.1. 6.2.1 Energy Deposited in the Scintillator The energy of leptons (electrons, antielectrons, and muons) created in the air shower can be measured by the energy released in the scintillator of the CHICOS detectors. The energy deposited by a particle moving through the scintillator is assumed to be 2.2 MeV/cm. For electrons, antielectrons, and muons that pass completely through the scintillator of the detector, the energy in MeV released is given by Emax (d, θ) = 2.2 MeV/cm × d cos θ (6.1) where d is the depth of the scintillator in cm, and θ is the zenith angle of the particle. If the energy of the particle is less than Emax , it is assumed that all of the kinetic energy of the particle is released. The number of particles is measured by converting the energy to units of minimum ionizing particles (MIP), with 1 MIP being equivalent to the energy deposited by a muon traveling through the scintillator at the average cosine zenith angle. Assuming that the distribution of incident particles goes as I(θ) ∝ cos2 θ, the average cos θ is R1 R1 I(θ) cos θ d cos θ cos3 θ d cos θ = R0 1 = 0.75. hcos θi = 0 R 1 I(θ) d cos θ cos θ d cos θ 0 0 (6.2) It is important to note that the energy equivalent of 1 MIP in eV varies from one detector to another depending on the thickness of the scintillator. 82 The energy contained in 1 MIP is given by E1MIP (d) = 2.2 MeV/cm × d . hcos θi (6.3) Combining equation (6.1) with equation (6.3), the energy in MIP of a particle passing completely through the scintillator is Emax (θ) = hcos θi . cos θ (6.4) The thickness of the scintillator can be ignored except in determining whether a particle is stopped in the detector. 6.2.2 Time-Over-Threshold Measurement The output pulse of the PMT is assumed to be a decaying exponential with time constant τ and height V0 : V (t) = V0 e−t/τ . (6.5) The energy of the pulse in MIP is proportional to the integral of the pulse, Z ∞ EMIP = V0 e−t/τ dt = τ V0 , (6.6) 0 where V0 includes the proportionality constant that relates the voltage to the energy, and is expressed in units of MIP/ns. The discriminator threshold is typically set to −10 or −15 mV. The time-overthreshold is measured by a counter card in units of 12.5 ns clock ticks. For a threshold of Vth , the pulse is over threshold for a time tth defined by Vth = V0 e−tth /τ . (6.7) EMIP = τ V0 = τ Vth etth /τ . (6.8) Thus 83 The conversion of time over threshold in ns to pulse energy in MIP is defined by 100 EMIP = Cetth /τ , (6.9) where C and τ are calibrated and stored in the database for each detector (see appendix A). The proportionality constant C from the database is therefore related to the threshold Vth by C = 100 τ Vth . (6.10) Hence the threshold can be expressed in units of MIP/ns as Vth = C . 100 τ (6.11) Expressing the threshold and pulse height in MIP/ns allows the simulation to ignore the variation of PMT sensitivity and threshold values between sites, as these effects are calibrated out by the choice of C and τ . For a typical detector, τ ≃ 80 ns and C ≃ 20 MIP. This gives Vth ≃ 0.0025 MIP/ns. A 1-MIP pulse would have V0 ≃ 0.0125 MIP/ns and a time over threshold of tth ≃ 135 ns. The minimum pulse with non-zero time over threshold would, in this case, be 0.2 MIP. In reality, the minimum length of the discriminator output pulse is ∼ 50 ns; the minimum measured pulse is therefore typically ∼ 0.37 MIP. 6.2.3 Timing Accuracy The time resolution of a CHICOS detector is assumed to be 12.5 ns. Random time offsets are introduced in order to model uncertainty in the detector timing. An offset of up to 6.25 ns is allowed between the A and B detectors at a given site, to reflect hardware differences between the detectors. Up to 25 ns is allowed between the timing at any two sites, to reflect the precision of the GPS timing. The clock rates on the counter cards are nominally 80 MHz, but may deviate 84 pulse intensity measurement point 12.5 ns time bin time Figure 6.1. Simulated pulse measurement. The intensity of the decaying exponential pulse is measured at 12.5 ns intervals. The measurement point is offset with respect to the start of the time bin by a random amount 0 < toffset < (1/4)× bin width. from this rate slightly, causing the time bins of the A and B data to be offset. In the simulation, the length of the time bins is taken to be exactly 12.5 ns, and the timescale for the bins to become misaligned is assumed to be long compared to the timescale of the shower (which is typically < 1 µs). The misalignment of the A and B time bins is therefore taken to be constant over any given shower. Figure 6.1 shows a diagram of sample pulse timing. 6.3 Modeling of Photon Interactions AIRES output identifies the species of each particle in the air shower. The signal contribution from charged particles (electrons and muons) is easily calculated as described above. However, the number of photons in the air shower at ground level is approximately 10 times the number of electrons, antielectrons, and muons. The photons are not detected directly, but can create secondary electrons in the scintillator via Compton scattering or pair production. These secondary electrons make a significant contribution to the signal. In order to accurately model the detector 85 response to air showers, care has been taken to estimate the signal contribution from photon interactions. Both processes (Compton scattering and pair production) have been modeled using Monte Carlo methods. This includes the interaction cross section (as a function of the photon energy) and the distribution of secondary electron energies. The distribution of polar angles of the secondary particles is ignored; the electrons and antielectrons are assumed to follow the same path as the incident photon. 6.3.1 Compton Scattering Compton scattering is the scattering of a photon off an electron in the scintillator. This is the dominant interaction process for photons below about 20 MeV. Compton scattering produces a recoil electron, which deposits energy in the scintillator. The differential cross section for Compton scattering is given by the Klein-Nishina formula [126] r2 1 dσ = e dω 2 [1 + γ0 (1 − cos θ)2 ] µ γ 2 (1 − cos θ)2 1 + cos θ + 0 1 + γ0 (1 − cos θ) 2 ¶ , (6.12) where γ0 is the incident photon energy, re is the classical electron radius, and θ is the angle through which the photon is scattered. The integral of this function is the total scattering probability per electron: ½ σc = 2πre2 · ¸ ¾ 1 + γ0 2(1 + γ0 ) 1 1 1 + 3γ0 − ln (1 + 2γ0 ) + ln (1 + 2γ0 ) − . γ02 1 + 2γ0 γ0 2γ0 (1 + 2γ0 )2 (6.13) The photon cross section for Compton scattering is shown in figure 6.2. The mean free path of a photon of energy γ0 in the scintillator can be calculated from the cross section as λ= 1 , nσ (6.14) where n is the electron density of the scintillator. For the plastic scintillator used by the CHICOS detectors, n ≃ 3.37 × 1023 electrons per cm3 . 86 Given λ, the probability of a photon penetrating x cm into the scintillator before scattering is P (x) = e−x/λ . (6.15) The probability that the photon will scatter before traveling x cm, Pint , is therefore defined by P (x) = 1 − Pint (x) = e−x/λ . (6.16) The cumulative interaction probability, Pint , is used as the Monte Carlo rejection function. Given a random number 0 < r < 1, the scattering depth xint of the photon is xint = −λ ln (1 − r). (6.17) If xint is less than the path length of the photon in the scintillator (taking into account the incident angle), the particle is assumed to undergo compton scattering after traveling distance xint . For a photon with incident zenith angle θ, the vertical depth dint is given by dint = xint cos θ. (6.18) The recoil electron is assumed to travel along the path of the incident photon until it exits the scintillator, or until it loses all of its kinetic energy. The maximum depth of scintillator traversed by the recoil electron, in terms of the total scintillator depth d, is therefore dmax = d − dint = d − xint cos θ. (6.19) The recoil electron can therefore be treated as a particle incident on a scintillator of depth dmax with zenith angle θ. The distribution of the recoil energy of the electron, T , is given by · ¶¸ µ πre2 s2 s 2 dσ = 2+ 2 + , s− dT me c2 γ02 γ0 (1 − s)2 1 − s γ0 (6.20) where s = T /γ0 [126]. The recoil photon energy γ1 , is more conveniently expressed as a fraction of the 87 intial energy: ² = γ1 /γ0 . The minimum photon recoil energy, ²min , is determined by the kinematics to be ²min = me c2 . me c2 + 2γ0 (6.21) Following the Monte Carlo procedure outlined in the Geant4 Physics Reference Manual [127], the probability function for the recoil energy of the photon can be approximated by · ¸· ¸ 1 ² sin2 θ P (²) ≃ +² 1− . ² 1 + ²2 (6.22) This formula can be expressed as P (²) = f (²) · g(²) = [α1 f1 (²) + α2 f2 (²)] · g(²), where α1 = ln ( 1 1 , (α1 ²) ² f2 = . α2 ), (6.23) f1 = ²min (1 − ²min ) α2 = , 2 (6.24) The functions f1 (²) and f2 (²) are probability density functions defined for ²min < ² < 1. The function · g(²) = 1 − ² sin2 θ 1 + ²2 ¸ (6.25) is the rejection function. For ²min < ² < 1, we have 0 < g(²) < 1. The Monte Carlo sampling of the recoil energy requires 3 random numbers 0 < r1 , r2 , r3 < 1. The first step is to select one of the two probability density functions: if r1 < α1 /(α1 /α2 ), select f1 (²); otherwise select f2 (²). The recoil energy ² is then sampled from the distribution corresponding to the probability density function: for f1 , for f2 , 2 ² = ²rmin , ²2 = ²2min + (1 − ²2min )r2 . (6.26) 88 The term sin2 θ in g(²) is given by sin2 θ = t(2 − t), (6.27) where t ≡ (1 − cos θ) = me c2 (1 − ²) . γ0 ² (6.28) The rejection function is tested against the third random number: if g(²) ≥ r3 , then the recoil energy of the photon is ²; otherwise the calculation should be repeated with a new set of r1 , r2 , r3 . The kinetic energy of the recoil electron is simply given by T = γ0 (1 − ²). (6.29) The distribution of electron energies given by this formula is shown in figure 6.3. 6.3.2 Pair Production Pair production is the transformation of a photon into an electron and an antielectron. The process must involve a second body (generally a nucleus in the scintillator) in order to conserve momentum. In general, numerical integration is required to calculate the pair production crosssection. This simulation uses a combination of two approximate formulae [126]. At low photon energies, me c2 ¿ hν ¿ 137me c2 Z −1/3 , σpp = 4Z 2 αre2 ¶ ¸ · µ 7 2hν 109 ln − f (Z) − , 9 me c2 54 (6.30) where f (Z) ≃ a2 [(1 + a2 )−1 + 0.20206 − 0.0369a2 + 0.0083a4 − 0.002a6 ] for a = αZ. For hν À 137me c2 Z −1/3 , · σpp = 4Z 2 αre2 ¸ ¢ ¢ 1 7¡ ¡ −1/3 ln 183Z − f (Z) − . 9 54 (6.31) The low-energy approximation assumes no electron screening, while the highenergy approximation assumes complete screening. The simulation chooses the former 89 2 σc × Re-2 [cm ] Cross Section for Compton Scattering in Scintillator 10 8 6 4 2 0 0 20 40 60 80 100 Incident Photon Energy [MeV] Figure 6.2. Cross section for Compton scattering. The cross section is a function of incident photon energy. Distribution of Recoil Electron Energies 600 ctest1 hν = 10 MeV Entries 10000 h ν = 25 MeV Mean 6.473 RMS 3.181 hν = 50 MeV hν = 100 MeV 500 400 300 200 100 0 0 10 20 30 40 50 60 70 80 90 100 Electron Recoil Energy [MeV] Figure 6.3. Distribution of recoil electron energies from Compton scattering. The distribution of recoil electron energies from 10000 Compton-scattered photons is shown for each of 4 different photon energies. 90 for photons of energy 6 MeV ≤ γ0 < 600 MeV and the latter for γ0 ≥ 600 MeV. The cross section is forced to go smoothly to 0 at 2 me c2 by the linear function σpp (γ0 < 6 MeV) = γ0 − 2me c2 σpp (γ0 = 6 MeV). 6 − 2me c2 (6.32) The pair-production cross sections for Z = 1 and Z = 6 are shown in figures 6.4 and 6.5. The composition of the plastic scintillator is 5.23 × 1022 H/cm3 and 4.74 × 1022 C/cm3 . The element participating in the interaction is chosen according to the probability defined by nk · σpp (Zk , γ0 ) , P (Zk , γ0 ) = P i [ni · σpp (Zi , γ0 )] (6.33) where ni is the number of atoms per volume of the ith element of the scintillator. The Monte Carlo calculation of the electron kinetic energy described in the Geant4 Physics Reference Manual [127] is based on the Bethe-Heitler cross section: dσpp (Z, ²) = d² ½ αre2 Z[Z + ξ(Z)] · ¸ · ¸¾ F (Z) 2 F (Z) [² + (1 − ²) ] Φ1 (δ(²)) − + ²(1 − ²) Φ2 (δ(²)) . 2 3 2 2 2 (6.34) In this expression, ² = E/γ0 , where E is the total energy carried by one particle in the the electron/antielectron pair. The kinematics constrain ² to the range me c2 = ²0 ≤ ² ≤ 1 − ²0 . γ0 (6.35) The screening variable δ(²) is given by δ(²) = 136 ²0 . Z −1/3 ²(1 − ²) (6.36) 91 Cross Section for Pair Production in Hydrogen 2 σpp × Re-2 [cm ] 0.08 0.07 0.06 0.05 0.04 0.03 0.02 0.01 0 0 20 40 60 80 100 Incident Photon Energy [MeV] Figure 6.4. Cross section for pair production in hydrogen. The cross section is a function of incident photon energy. σpp × Re-2 [cm^2] Cross Section for Pair Production in Carbon 3 2.5 2 1.5 1 0.5 0 0 20 40 60 80 100 Incident Photon Energy [MeV] Figure 6.5. Cross section for pair production in carbon. The cross section is a function of incident photon energy. 92 The two screening functions, Φ1 (δ) and Φ2 (δ), are defined by ) Φ1 (δ) = 20.867 − 3.242δ + 0.625δ 2 δ ≤ 1, Φ2 (δ) = 20.209 − 1.930δ − 0.086δ 2 (6.37) Φ1 (δ) = Φ2 (δ) = 21.12 − 4.184 ln (δ + 0.952) δ > 1. The function F (Z) is a Coulomb correction function: F (Z) = 8/3 ln Z for γ0 < 50 MeV, F (Z) = 8/3 ln Z + 8f (Z) for γ0 ≥ 50 MeV, (6.38) where f (Z) is the correction to the Born approximation given above. The contribution to pair production from the electron cloud of the atom is given by ξ(Z) = ln (1440/Z 2/3 ) . ln (183/Z 1/3 ) − f (Z) (6.39) The Bethe-Heitler formula is symmetric under ² ↔ (1−²), so it suffices to consider the range ² ∈ [²0 , 1/2]. (6.40) However, the requirement that the cross section be positive produces an addtional constraint: ¸ 42.24 − F (Z) − 0.952. δ ≤ δmax (²1 ) = exp 8.368 · Therefore, 1 1 ² ≥ ²1 = − 2 2 where µ 1 δmin = δ ² = 2 r δmin , δmax (6.42) 136 4²0 . Z 1/3 (6.43) 1− ¶ = (6.41) 93 The range of ² is therefore ² ∈ [²min = max(²0 , ²1 ), 1/2]. (6.44) The screening functions Φ1 (δ) and Φ2 (δ) can be combined with F (Z) to form two additional screening functions that are used in the factorization of the cross section: F1 (δ) = 3Φ1 (δ) − Φ2 (δ) − F (Z), 3 1 F2 (δ) = Φ1 (δ) − − F (Z). 2 2 (6.45) The functions F1 (δ) and F2 (δ) are decreasing functions of δ. They reach a maximum at δmin : F10 = max F1 (δ) = F1 (δmin ), (6.46) F20 = max F2 (δ) = F2 (δmin ). Using these definitions, the Bethe-Heitler formula can be written · ¸ dσpp (Z, ²) 2 1 2 = αre Z[Z + ξ(Z)] − ²min × [N1 f1 (²)g1 (²) + N2 f2 (²)g2 (²)] , (6.47) d² 9 2 where · 1 − ²min N1 = 2 3 N2 = F20 . 2 ¸2 F10 , (6.48) The probability density functions f1 (²) and f2 (²) are f1 (²) = £ 3 · ¤ 3 1 − ²min 2 1 −² 2 ¸2 1 ¤. f2 (²) = const = £ 1 − ²min 2 (6.49) The rejection functions are given by g1 (²) and g2 (²), defined on the interval 0 < 94 Distribution of Pair-Produced Electron Energies ptest1 hν = 10 MeV Entries 10000 hν = 25 MeV Mean 1.971 RMS 1.293 hν = 50 MeV 300 hν = 100 MeV 250 200 150 100 50 0 0 10 20 30 40 50 Electron Energy [MeV] Figure 6.6. Pair-produced electron energies. The distribution of pair-produced electron energies from 10000 photons is shown for each of 4 different photon energies. gi (²) < 1: F1 (²) , F10 F2 (²) g2 (²) = . F20 g1 (²) = (6.50) As for Compton scattering, the Monte Carlo sampling of the electron energy requires 3 random numbers 0 < r1 , r2 , r3 < 1. For r1 < N1 /(N1 + N2 ), f1 (²) and g1 (²) are used, f2 (²) and g2 (²) otherwise. The electron energy is then sampled according to the appropriate probability density: for f1 , for f2 , µ ¶ 1 1/3 − ²min r2 , 2 ¶ µ 1 2 − ²min r2 . ² = ²min + 2 1 ²= − 2 (6.51) The appropriate rejection function is then tested against the third random number: if gi (²) ≥ r3 , the energy of the electron (and of the antielectron) is ². For γ0 < 2 MeV, ² is sampled uniformly on ² ∈ [²min , 1/2], without rejection. The distribution of electron energies is shown in figure 6.6. 95 Table 6.2. Fraction of showers detected at 1017 eV 6.4 Energy Cosine ZA Number detected Fraction of total 1017 eV 1017 eV 1017 eV 1017 eV 0.90 0.80 0.70 0.60 962 644 266 86 48.1% 32.3% 13.3% 4.3% Analysis of Simulated Shower Reconstructions At each shower inclination (cos θ = 0.60, 0.70, 0.80, 0.90), two AIRES-generated showers were each dropped on the simulated Chiquita array 1000 times. For each drop, the AIRES shower was rotated through a randomly chosen azimuthal angle, and displaced to a random core location (x, y), with −500 < x < 500 and −500 < y < 500 (with the origin centered on the site Caltech 03.) The results were passed through the CHICOS filtering software; only a fraction of the showers at each angle produced a trigger and enough matching hits to be considered an event (Table 6.2). The number of showers detected falls steeply with increasing inclination. The resulting shower data was fed back to the CHICOS reconstruction software (using the χ2 method). Since all simulated showers were given an input energy of 1017 eV, the expected distribution of reconstructed energies should be centered near that value. As shown in figure 6.7, the average reconstructed energy was 1017.11 eV, or 1.3 × 1017 eV. The peak of the distribution is higher than the input energy; however the distribution extends farther below the peak than above it. There is a selection effect at work in that showers that fluctuate up in the vicinity of a site are more likely to be found than showers that fluctuate down. This may partially account for the shift in the peak energy. The average reconstructed energy did not vary except for showers with the greatest inclination (figure 6.8). Selection effects are amplified for highly inclined showers, which are more strongly attenuated by the atmosphere. The distributions of reconstructed angles are plotted in figure 6.9. The average 96 reconstructed inclination is very good in general, with the most inclined showers again being the least well reconstructed. The software does not allow showers at this energy to reconstruct with an angle greater than 54.7◦ (cos θ = 0.58). (The software has more recently been upgraded to allow high-energy reconstructions out to cos θ = 0.5). There is a slight correlation between the reconstructed angle and the reconstructed energy. If a shower is reconstructed at an inclination greater than the input angle, the reconstructed energy tends to be greater than the input energy. Finally, it is interesting to look at the distribution of shower core locations. Showers were dropped randomly over a 1 km2 area centered on Caltech 03. The core locations of showers that were picked up as CHICOS events are shown in figure 6.10. There is an overdensity near the center because showers that land near the array sites are more likely to generate enough hits to pass the filtering software. The reconstructed core locations are shown in figure 6.11. There is a strong tendency for the reconstructor to place the core near, but not directly over, sites in the array. Error in the core location is slightly correlated with error in the energy reconstruction. If the reconstruction moves the shower closer to the array, it will also lower the energy estimate. Figure 6.12 shows an example of the simulated shower data. The input parameters of the shower are shown. Figure 6.13 shows the reconstructor ouput for the same shower. 97 htemp Entries 1958 Mean 17.11 RMS 0.3506 Reconstructed Log(E [eV]) 180 160 140 120 100 80 60 40 20 0 13 14 15 16 17 18 Log(E [eV]) Figure 6.7. Distribution of reconstructed energies for simulated showers at 1017 eV. Showers of all zenith angles are included in the histogram. The peak of the reconstructed distribution is shifted upward from the input energy to 1.3 × 1017 eV. 98 Reconstructed Log(E) for input Cosθ = 0.6 htemp Entries 86 Mean 17.18 RMS 0.2768 5 Reconstructed Log(E) for input Cosθ = 0.7 htemp Entries 266 Mean 17.11 RMS 0.3175 18 16 4 14 12 3 10 8 2 6 4 1 2 0 16.4 16.6 16.8 17 17.2 17.4 17.6 Reconstructed Log(E) for input Cosθ = 0.8 17.8 18 Log(E [eV]) htemp Entries 644 Mean 17.11 RMS 0.3395 45 0 15 80 70 30 60 25 50 20 40 15 30 10 20 5 10 15.5 16 16.5 17 17.5 18 Log(E [eV]) 16.5 17 17.5 0 13 14 15 16 18 Log(E [eV]) htemp Entries 962 Mean 17.1 RMS 0.3711 90 35 15 16 Reconstructed Log(E) for input Cosθ = 0.9 40 0 15.5 17 18 Log(E [eV]) Figure 6.8. Distribution of reconstructed energies sorted by shower inclination. All showers had a simulated energy of 1017 eV. An equal number of showers were simulated at each zenith angle, but fewer are detected at large inclinations. 99 Reconstructed Cosθ for input Cosθ = 0.6 htemp Entries 86 Mean 0.6223 RMS 0.04948 24 22 Reconstructed Cosθ for input Cosθ = 0.7 htemp Entries 266 Mean 0.705 RMS 0.05549 18 20 16 18 14 16 14 12 12 10 10 8 8 6 6 4 4 2 2 0 0.6 0.65 0.7 0.75 Reconstructed Cosθ for input Cosθ = 0.8 0.8 Cosθ htemp Entries 644 Mean 0.8038 RMS 0.05539 40 0 0.55 0.65 0.7 0.75 0.8 0.85 0.9 Reconstructed Cosθ for input Cosθ = 0.9 60 30 50 0.95 1 Cosθ htemp Entries 962 Mean 0.9005 RMS 0.04227 70 35 25 0.6 40 20 30 15 20 10 10 5 0 0.55 0.6 0.65 0.7 0.75 0.8 0.85 0.9 0.95 1 Cosθ 0 0.55 0.6 0.65 0.7 0.75 0.8 0.85 0.9 0.95 1 Cosθ Figure 6.9. Distribution of reconstructed angles sorted by shower inclination. The range of cos θ is limited by the reconstruction software to cos θ ≥ 0.58. All showers had a simulated energy of 1017 eV. An equal number of showers were simulated at each zenith angle, but fewer are detected at large inclinations. 100 Input core locations y [m] 1000 500 0 -500 -1000 -1000 -500 0 500 1000 x [m] Figure 6.10. Input core locations of simulated showers. Only showers that are detected as candidate events are shown. Locations of sites in the Chiquita array are shown as blue triangles. Simulated core locations were evenly distributed over 1 km2 ; showers that landed near the array sites were more likely to pass the filtering software. Reconstructed core locations y [m] 1000 500 0 -500 -1000 -1000 -500 0 500 1000 x [m] Figure 6.11. Reconstructed core locations of simulated showers. Only showers that are detected as candidate events are shown. Locations of sites in the Chiquita array are shown as blue triangles. 101 Figure 6.12. Example of simulated shower data. Input parameters are shown. 102 Figure 6.13. Example reconstruction of simulated shower data. 103 Chapter 7 Analysis of Low-Energy Data Cosmic ray airshowers in the energy range 1016 –1019 eV have been measured by a subarray of CHICOS, as described in section 7.1. Comparison of this analysis with other experiments was used to refine the CHICOS data analysis before examining the ultra-high energy data. The simulations used to characterize the Chiquita subarray are detailed in section 7.2. Section 7.3 describes the method used to estimate the low-energy cosmic ray flux. 7.1 Low-Energy Data As described in section 4.1, a subsection of the CHICOS array called Chiquita has been in operation on the Caltech campus to capture the lower-energy range of the UHECR spectrum. The Chiquita array consists of 12 neighboring pairs, 6 of which are grouped in a closely spaced set. The remaining 6 sites cover an area of approximately 1 km2 (figure 7.1). The denser spacing between sites results in the Chiquita setup being sensitive to lower-energy events than the larger CHICOS array; data from Chiquita spans the energy range 1016 to 1019 eV, with peak sensitivity approximately between 1017 and 1018 eV. The flux of cosmic rays in this range has been measured by previous experiments such as Haverah Park [12, 128, 129], Yakutsk [130, 131, 132], Akeno [133, 134], Fly’s Eye [135, 136], and HiRes [137, 138]. The data collected in this energy range can therefore serve as a useful calibration of the larger CHICOS array. 104 y [m] Map of Chiquita sites 800 600 400 200 0 -200 -400 -1000 -500 0 500 1000 x [m] Figure 7.1. Map of the Chiquita detector sites. The origin is placed on the site Caltech 03, at longitude −118.12522◦ , latitude 34.13657◦ . The larger marker at approximately (−100 m, −50 m) represents the six closely spaced CIT 1001–1006 sites. Incoming shower events are classified as either “CHICOS” or “Chiquita” based on the sites involved in the shower. All sites that are not part of the Chiquita array are part of the CHICOS array. Within the Chiquita array, 9 sites are members of the Chiquita array only: CIT 1001–1006, CIT 2001, CIT 3001, and CIT 4001. The remaining three sites in the Chiquita array are also considered to be members of the CHICOS array: Caltech 03, Polytechnic 096, and PCC 105. A CHICOS shower satifies the following criteria: it has a trigger at a site that is a CHICOS site and it has hits at 3 out of the 5 nearest-neighbor CHICOS sites. Similarly, a Chiquita shower must have a trigger at a Chiquita site and hits at 3 out of the 5 nearest-neighbor sites from the Chiquita array. (CIT 1001–1006 are treated as a single site for this purpose.) Note that because some sites belong to both arrays, some showers satisfy the requirements for both categories. All showers must additionally pass an initial quality test on the time distribution of the hits. The trigger plus any 3 of the additional sites are fit to a planar shower 105 front passing through the array. The quality requirement is given as sX TRMS = (texpected − tmeasured )2 < 10, (7.1) sites where texpected is the fit time in µs at that site, and tmeasured is either the average hit time in µs (for CIT 1001–1006) or the measured time in µs of the hit closest to the trigger time (for all other sites). The data set of 18,403 Chiquita showers used for this analysis was drawn from a selection of time periods totalling ∼ 3743 hours between August 1, 2003 and May 31, 2004 during which all sites were operational. Selecting only data recorded while all 12 sites were running made it possible to avoid the problem of time dependence when calculating the effective aperture of the array. There are 15 showers in the data set that could also be classified as CHICOS showers. In order to automate the reconstruction of this very large data set, an attempt is made to remove accidental hits from the shower data (as such spurious hits have a tendancy to result in poor fits to the data). All hits at non-Chiquita sites are removed. All remaining hits that fall more than 10 µs away from the trigger are assumed to be accidentals. Hits more than 1 µs after the first hit at a given detector are assumed to be caused by electronic afterpulsing and are also removed. Reconstruction of the data was done using a χ2 method and the low-energy ironprimary LDF described in section 3.2. Following reconstruction, a final cut on the data is applied such that only showers with θ ≤ 45◦ are included in the analysis. The final data set consists of 16,592 showers. 7.2 Simulation of Low-Energy Showers The effective aperture of the array is a function of the physical size of the array, the spacing between sites, the sensitivity of the detectors, and the methods used to identify shower events. The response of the array is modeled by evaluating the appropriate LDF and time-delay formula at the location of each site in the array and 106 Number of Showers Input Spectrum 5 10 104 103 16 16.5 17 17.5 18 18.5 19 19.5 Log(E) Figure 7.2. Energy distribution of showers simulated for Chiquita. This set of simulated showers was used in characterizing the performance of the Chiquita array. The distribution is not intended to reflect the actual power-law flux of showers, but rather to generate an approximately equal number of detected showers in each bin. applying Poisson statistics to produce simulated data. The resulting shower data is then filtered through the same cuts as the real data to determine whether a real shower with those parameters would be detected. To evaluate the aperture of the Chiquita array as a function of energy, a series of showers with energies ranging between 1016.0 eV and 1019.6 eV was simulated, covering a (2 km)2 area centered on the origin. A Monte Carlo routine was used to select core location and shower angle. The vertical angle cos θ was chosen to lie in the range of 0.5 and 1. Shower data is divided into bins of width 0.2 in log-energy space. To obtain good simulated statistics in each energy bin, showers were generated within that bin until 1000 showers were found that passed the event filtering software. (This number was further reduced after the cut on reconstructed angle θ ≤ 45◦ was applied.) The energy distribution of showers within each bin follows a power-law spectrum with spectral index −3. The number of showers generated in each energy bin is shown in figure 7.2. The energy distribution of showers that are found by the filtering software is shown in figure 7.3. Also shown are the core locations (figure 7.4) and angles (figure 7.5 and figure 7.6). 107 Simulated Log(E) Reconstructed Log(E) 900 1200 800 1000 700 600 800 500 600 400 300 400 200 200 100 0 14 15 16 17 18 19 20 0 14 21 22 Log(E [eV]) 15 16 17 (a) 18 19 20 21 22 Log(E [eV]) (b) Simulated core locations Reconstructed core locations 1000 1000 y [m] y [m] Figure 7.3. Reconstructed energies of low-energy simulations. (a) Distribution of simulated energies for showers that are detected as events. Showers were generated with an E −3 distribution within each bin, although the distribution of detected showers within each bin tends to be flatter because the acceptance is smaller at the lowerenergy end. In each bin, 1000 showers passed the data filtering software; some fraction of these were subsequently removed from the simulated data set if the reconstructed vertical angle was ≥ 45◦ . Highly inclined showers were more likely to be detected at higher energies. (b) Distribution of reconstructed energies for the same set of showers. 500 500 0 0 -500 -500 -1000 -1000 -1000 -500 0 (a) 500 1000 x [m] -1000 -500 0 500 1000 x [m] (b) Figure 7.4. Reconstructed core locations of low-energy simulations. (a) Distribution of simulated core locations for showers that are detected as events. Simulated core locations were evenly distributed over the area shown; showers that landed near the array sites were more likely to pass the filtering software. (b) Distribution of reconstructed core locations for the same set of showers. 108 Simulated Cos(θ) Reconstructed Cos(θ) 1400 1400 1200 1200 1000 1000 800 800 600 600 400 400 200 200 0 0.5 0.55 0.6 0.65 0.7 0.75 0.8 0.85 0.9 0.95 1 Cos(θ) 0 0.5 0.55 0.6 0.65 (a) 0.7 0.75 0.8 0.85 0.9 0.95 1 Cos(θ) (b) Figure 7.5. Reconstructed vertical angles of low-energy simulations. (a) Distribution of the cosine of the simulated vertical angle θ for showers that are detected as events. Simulated shower angles were evenly generated over the surface of a sphere down to cos(θ) = 0.5. (b) Distribution of the cosine of the reconstructed angle θ for the same set of showers. The cut on reconstructed angle requiring cos θ ≤ 45◦ has been applied. Simulated φ Reconstructed φ 600 600 500 500 400 400 300 300 200 200 100 100 0 0 1 2 3 (a) 4 5 6 φ [rad] 0 -3 -2 -1 0 1 2 3 φ [rad] (b) Figure 7.6. Reconstructed azimuthal angles of low-energy simulations. (a) Distribution of simulated azimuthal angle φ for showers that are detected as events. (b) Distribution of reconstructed azimuthal angle φ for the same set of showers. Spread in reconstructed energies 22 Reconstructed Log(E [eV]) Acceptance N(passed)/N(in) 109 1 10-1 250 21 20 200 19 150 18 17 10-2 100 16 50 15 -3 10 16 16.5 17 17.5 18 18.5 19 19.5 Log(E [eV]) 14 16 16.5 (a) 17 17.5 18 18.5 19 19.5 Simulated Log(E [eV]) 0 (b) Figure 7.7. Acceptance of the Chiquita array. (a) Acceptance of the array as a function of energy. (b) Reconstructed energy vs. simulated shower energy. The acceptance of the array at a given energy is defined as the ratio of the number of showers that trigger the array and pass the filters to the number of showers thrown on the array at that energy. The effective aperture is the acceptance ratio multiplied by the area and solid angle covered by the simulated data set. To fit the data, it is also necessary to characterize the accuracy of the reconstruction software, which introduces a spread in the reconstructed energies, and in some cases an offset in the average reconstructed energy. Figure 7.7 shows the acceptance of the Chiquita array and a plot of reconstructed energy vs. simulated energy. 7.3 Estimating the Low-Energy Flux In order to correctly plot the cosmic ray flux, it is necessary to adjust the observed flux for the fact that not all showers are detected, and for the fact that some showers will be erroneously reconstructed in a higher or lower energy bin. Given an assumed cosmic ray flux J0 (E), this function can be multiplied by the acceptance, and convolved with the spread in reconstructed energy to model the observed shower flux. We work backward from this point to derive an expression for the flux (evaluated at the central point of each energy bin) as a function of the observed number of showers per bin. The Chiquita data set of 16592 showers, shown in figure 7.8, was used. Only data 110 Chiquita Data 3000 2500 2000 1500 1000 500 0 16 16.5 17 17.5 18 18.5 19 Log(E [eV]) Figure 7.8. Energy distribution of Chiquita showers between 1016 eV and 1019 eV. points above 1017 eV (where the aperture drops to ∼ 10% of its maximum value) were included in the fit. We begin by assuming that the number flux (per eV) of cosmic rays follows the power law J0 (E) = dN0 = CE −γ . dE (7.2) The measured number flux per interval in energy, dN , is the integral of the actual dE 0 number flux, dN , multiplied by the function R(E, E 0 ), which describes the distribudE 0 tion of reconstructed energies for showers of energy E 0 . Thus the observed flux J(E) is related to the real flux by dN = J(E) = dE Z dE 0 A(E 0 )R(E, E 0 ) dN0 . dE 0 (7.3) The function A(E 0 ) is the acceptance of the array (the fraction of showers that pass all cuts on the data at a given energy). 111 The histogram of Chiquita data is binned in log space, with bin width log (Eupper ) − log (Elower ) = 0.2. (7.4) In terms of the central energy EC of the bin, the bin width is Eupper − Elower = 100.1 EC − 10−0.1 EC = 0.46EC . (7.5) Each bin in the histogram therefore contains a number of showers given by Z 100.1Ei N (Ei ) = dN dE, 10−0.1Ei dE (7.6) where Ei is the central energy (in log space) of the ith bin. To relate the number of Chiquita showers in a given bin to the incoming flux, we in equation (7.3) into equation (7.6) above: substitute the expression for dN dE Z Z N (Ei ) = dE 0 A(E 0 )R(E, E 0 ) dE Ei bin all E0 dN0 . dE 0 (7.7) We have the function R(E, E 0 ) in the form of a data table for a fixed set of values of E 0 , so instead of a continuous integral, a sum would be more appropriate. We make the approximation that the entire flux in the bin centered around Ej is being transferred to the bin centered around the observed energy Ei with the efficiency R(Ei , Ej ).a This gives N (Ei ) = X Ej0 Z A(Ej0 )R(Ei , Ej0 ) dE 0 . E0j bin dN0 dE 0 (7.8) a The data table of values for R(Ei , Ej ) was calculated using simulated showers with an E −3 distribution over the width of the bin. It is not the efficiency for moving showers from the central energy of one bin to the central energy of another bin, but rather an approximation to the fraction of showers in one bin that get moved to another bin due to error in the reconstructed energies. This is why substituting R(Ei , Ej ) for R(E, E 0 ) takes care of the integral over the Ei bin. Note, however, that the farther γ is from 3, the more error is introduced through R. (This also applies to A(E), which is tabulated from the same set of simulated showers.) 112 The integral of the flux over the Ej bin is given by Z CE Ej bin −γ Z 100.1Ej dE = 10−0.1Ej CE −γ dE ¯100.1Ej CE −γ+1 ¯¯ = −γ + 1 ¯10−0.1Ej h¡ ¢−γ+1 ¡ −0.1 ¢−γ+1 i 1 = 100.1 − 10 CEj−γ+1 . −γ + 1 (7.9) (7.10) (7.11) If we assume that γ is close to 3, we can reasonably approximate the coefficient in the above equation by the value at γ = 3. For γ = 3, Z Ej bin CE −γ dE = 0.477CEj−γ+1 . (7.12) The function that should be used to fit the data is therefore N (Ei ) ≃ X 0.477A(Ej )R(Ei , Ej )CEj−γ+1 . (7.13) Ej We want to plot flux times E 3 : J0 (E)E 3 = CE 3−γ . (7.14) When this expression is evaluated at the central energy of each bin, it yields the set of data points F1 (Ei ) ≡ J0 (Ei )Ei3 = CEi3−γ . (7.15) Equation (7.13) relates the parameters C and γ to the observed number per bin, so we begin by isolating these in the sum: N (Ei ) ≃ X 0.477A(Ej )R(Ei , Ej )CEj3−γ Ej−2 . (7.16) Ej We make the assumption that γ is close to 3, so Ej3−γ is approximately constant over the range of energies being considered. We also assume that R(Ei , Ei ) dominates 113 the sum, and approximate Ej3−γ by its value at Ei , allowing us to extract it from the sum. This gives F1 (Ei ) = CEi3−γ = P N (Ei ) −2 . Ej 0.477 A(Ej )R(Ei , Ej ) Ej (7.17) We now have a set of points that is an approximation to the flux times E 3 at the central point of each bin. We would like to improve this, however, by estimating the error and dividing it out. To do this, we fit the points to a power law, then use this power law to generate a set of artificial data. By then plotting the results we would obtain from this artifical data and comparing it to the input power law, we can estimate the offset in each bin that results from the assumptions in this analysis and divide it out. The first step is to fit a line to the data, and get values for C and γ. Let the fit values be C1 and γ1 . To perform one iteration on this result, we first calculate the number of showers per bin that would correspond to the fitted values of C1 and γ1 : N2 (Ei ) = X 0.477 A(Ej )R(Ei , Ej ) C1 Ej−γ1 +1 . (7.18) Ej This generated data set is used to calculate a new histogram of flux times E 3 : N2 (Ei ) −2 . Ej 0.477 A(Ej )R(Ei , Ej ) Ej F2 (Ei ) = P (7.19) In the absence of error introduced by A(Ej ) and R(Ei , Ej ), and the assumption that γ close to 3, the set of points generated in this way would lie on the fit line C1 E 3−γ1 . The flux is now estimated as C1 Ei3−γ1 3 . J(Ei )Ei = F1 (Ei ) F2 (Ei ) (7.20) This is plotted in figure 7.9. Error in each bin is estimated by allowing a shift in 114 Flux*E^3 [eV^2/m^2/sr/s] 1025 1024 Chiquita Spectrum HiRes-II 1023 Fly’s Eye Stereo Akeno 1022 17 17.2 17.4 17.6 17.8 18 18.2 18.4 18.6 18.8 19 Log(E [eV]) Figure 7.9. Flux times E 3 , as measured by the Chiquita array. An error band based on an assumed energy reconstruction error of 20% is shown. Data from the following experiments is shown for comparison: HiRes [112], Fly’s Eye [135], and Akeno [133, 134]. the reconstructed shower energies of up to ±20%. This is assumed to dominate over other sources of error. The Chiquita data agrees with other experiments at 1017 eV, but results in a lower flux estimate at higher energies. Results at the lower end of the energy range are likely to be more reliable. Above 1018 eV, the flux estimates are based on only a few showers per bin. Also, showers at 1017 eV are significantly smaller than showers at 1018 or 1019 eV, and are thus unlikely to be detected by the Chiquita array unless they fall within the central cluster of sites. However, showers that do trigger the array by falling near the center can be well measured by surrounding sites. Higher-energy showers are capable of triggering the array even when the shower core lies outside the perimeter of the array. This reduces the constraints on the reconstruction, and increases the uncertainty associated with the core location estimate and the closely related energy estimate. In addition, the relatively linear configura- 115 tion of the Chiquita array reduces the accuracy with which the incident angle can be determined. It is likely that the assumed 20% systematic error depicted in figure 7.9 is an underestimate for this energy range. 116 Chapter 8 Analysis of High-Energy Data The goal of the CHICOS experiment has been to measure the energy spectrum and sky distribution of ultra-high energy cosmic rays (above approximately 1019 eV). As described in section 8.1, results from previous experiments have been conflicting, although with better quality data emerging from Auger and HiRes, consensus is gradually being reached. Section 8.2 describes the data collection, filtering, and quality selection used in producing the CHICOS data set. The simulations used to characterize the CHICOS array are detailed in section 8.3. Section 8.4 describes the method used to estimate the high-energy cosmic ray flux. 8.1 Results from Previous Experiments Above 1019 eV, the flux of cosmic rays is less than 1 per km2 per year. Only the largest detector arrays have been in a position to make accurate measurements of the flux in this energy regime. Early experiments that measured the ultra-high energy spectrum include Yakutsk [139] and AGASA [23]. The results obtained from these two experiments were inconsistent with a GZK cutoff at the expected energy of approximately 1020 eV. However, the Yakutsk data had very few super-GZK events, with only 4 data points above 1019.9 eV. The AGASA data initially appeared more definitive, with 11 events above 1020 eV, but a reanalysis of the data has since reduced this number to 6 [38, 39], and AGASA no longer claims a statistically significant nondetection of the 117 GZK cutoff. More recently, HiRes [36, 37] and the Pierre Auger Observatory [43, 44, 42] have measured the ultra-high energy spectrum with a sky exposure exceeding that of AGASA. The results from both experiments are consistent with a GZK cutoff at approximately 1020 eV. HiRes claims to have observed the cutoff with a statistical significance of 5 standard deviations [37]. Auger similarly rejects a uniform continuation of the spectrum with a statistical significance of 6 standard deviations [42]. A major goal of the CHICOS experiment has been to obtain an ultra-high energy cosmic ray spectrum that can be compared to these results. 8.2 High-Energy Data from CHICOS The CHICOS array has been in operation since 2003. The design of the array is described in section 4.1. A shower event must trigger at least one detector in the CHICOS array and register hits at 3 out of the 5 sites nearest to the trigger (within a ±50 µs window). Some sites on the Caltech campus are excluded from the set of nearest neighbors (CIT 1001–1006, CIT 2001, CIT 3001, and CIT 4001), but hits at these sites will appear on showers that meet the criteria. All CHICOS showers must additionally pass the quality test given in equation (7.1). With data generated by the Chiquita array, it was possible to select only showers that occured when all sites were operational, thus simplifying the analysis. This was not practical with the larger CHICOS array, firstly because the configuration of the array (including the total number of sites) changed with time, and secondly because some subset of sites in the array are in need of maintenance at any given time. All shower events between January 1, 2003 and December 31, 2007 are included in the initial data set, and the changing acceptance is accounted for in the simulations of the array response. There were 1164 candidate showers detected during this time period. Prior to reconstructing the showers in this data set, showers with two triggers separated by more than 10 µs were removed from the data set, as the possibility exists for two nearly simultaneous low-energy showers to mimic the spread of a single 118 high-energy event, contaminating the data. In the remaining showers, hits that fell more than 25 µs before or after the trigger, or were recorded more than 8 km from the trigger site, were assumed to be accidentals and were removed. Shower reconstruction was done using the maximum-likelihood analysis described in section 5.2, with the high-energy proton-primary LDF described in section 3.2. Following reconstruction, several additional quality cuts were made on the data set. In order to be considered a valid reconstruction, a shower must have hits at a minimum of 4 sites, including the trigger. All hits at shower match sites in a reconstruction must fall within 5 µs of the fit time at the hit site.a (This criterion is not applied to the trigger site, but this has little impact on the data set, since the trigger site usually drives the reconstruction.) Hits that fail this criterion are removed from the shower, and the reconstruction is iterated. Hits at sites that were reported as unhealthy (as defined in section 5.3) during the shower are also removed. If no fit to the remaining data can be found before the number of sites in the shower drops below 4, the shower is discarded. Next, a cut on the quality of the lateral distribution function fit to the intensity data is made by requiring that at least 2 of the shower sites be fit with an intensity of 0.1 MIP or greater, and within 2 standard deviations of the average intensity at the site. This is a relatively weak criterion, but it serves the purpose of eliminating a specific kind of reconstructor problem where the the core is placed directly over the trigger site and the LDF is fit to that site only. Finally, the reconstruction is required to have a core location not more than 2500 m from at least one of the sites in the array. This eliminates the possibility that the reconstructor will attempt to fit the time information by moving the core of the shower outside the array (so that the time distribution function associated with the shower sites will be flatter) while simultaneously increasing the estimated energy. After these cuts are made, 275 showers remain in the CHICOS data set. A final a A variable allowed time discrepancy based on the reported time error from the reconstructor was considered, but was rejected on the basis of being unphysically small for sites near the trigger (because the intrinsic width of the shower front is not taken into account) and excessively large for sites at the edge of the shower (i.e., larger than the predetermined cut to eliminate accidental hits). 119 cut requiring a vertical angle of θ < 1 rad is applied, bringing the final number of showers down to 208. The maximum vertical angle allowed by the reconstructor for CHICOS showers is 60◦ ; this cut corresponds to requiring θ < 57.3◦ and eliminates all showers that may have a real inclination of greater than 60◦ but were reconstructed with the maximum allowed angle. Of the 208 showers in the final data set, 145 reconstruct with energy above 1018.4 eV, which is the lower limit of the energy range used in the analysis. 8.3 Simulation of High-Energy Showers As with the low-energy analysis, the effective aperture of the CHICOS array is a function of the physical structure of the array as well as the methods used to identify and select showers. However, since any given site in the CHICOS array may be offline at a given time, the aperture of the CHICOS array changes in a complex way on small timescales. The response of the array in a given configuration is modeled by evaluating the appropriate LDF and time-delay formula at the location of each site in the array and applying Poisson statistics to produce simulated data. The resulting shower data is then filtered through the same cuts as the real data to determine whether a real shower with thoese parameters would have been detected at that particular point in time, given the state of the array.. To evaluate the aperture of the CHICOS array as a function of energy, a series of showers with energies ranging between 1018.4 eV and 1021.0 eV was simulated. Each shower was randomly assigned a UTC time (to the nearest second) between January 1, 2003 and December 31, 2007, and the configuration of the array at that time was used for the simulation and the reconstruction. Between the energies of 1018.4 eV and 1018.8 eV, showers were thrown over an area of 360 km2 covering the San Gabriel array. Showers in this energy range are not detected in the sparser San Fernando array. Above 1018.8 eV, showers were thrown over this area plus an area of 810 km2 covering the San Fernando array (figure 8.1). Showers were generated by energy bin and year until 1000 simulated showers 120 Simulated core locations, Log(E [eV]) < 18.8 Simulated core locations, Log(E [eV]) > 18.8 400 1200 15000 y [m] 20000 y [m] 20000 350 15000 300 1000 10000 10000 250 5000 200 800 5000 600 0 400 150 0 100 -5000 -10000 -50000 200 -40000 -30000 -20000 -10000 0 10000 x [m] 0 -5000 50 -10000 -50000 -40000 (a) -30000 -20000 -10000 0 10000 x [m] 0 (b) Figure 8.1. Core locations of simulated high-energy showers. (a) Distribution of core locations of simulated showers, for energies below 1018.8 eV. (b) Distribution of core locations for simulated showers, for energies above 1018.8 eV. generated detectable events in each energy bin (of width 0.2 in log-energy space) and each year. Over 5 years and 13 energy bins, this totals 65,000 simulated shower events. Following cuts on reconstruction quality, as described in section 8.2, the number of simulated shower events was reduced to 38,973. As with the low energy simulations, the simulated shower angle was randomly drawn from a uniform distribution on a sphere, with cos θ ranging from 0.4 to 1, and the energy of each shower was drawn from a power-law distribution with spectral index of −3. Following a cut on the reconstructed vertical angle, the number of accepted events was further reduced to 30,180. The number of showers generated in each energy bin is shown in figure 8.2. The energy distribution of showers that are found by the filtering software is shown in figure 8.3. Also shown are the core locations of detected shower events (figure 8.4 and figure 8.5) and their angles (figure 8.6 and figure 8.7). The acceptance of the array is calculated as a time average over the 5 years during which data was taken. Figure 8.8 shows the acceptance of the CHICOS array as a function of energy, and a plot of reconstructed energy vs. simulated energy. The method used to generate simulated showers (by assigning a random arrival time and angle) results in a set of simulated showers (those that are detected as 121 Number of Showers Input Spectrum 106 105 104 18.5 19 19.5 20 20.5 21 Log(E [eV]) Figure 8.2. Energy distribution of showers simulated for CHICOS. This set of simulated showers was used in characterizing the performance of the CHICOS array. This distribution is not intended to reflect the actual power-law flux of showers, but rather to generate an approximately equal number of detected showers in each bin. There is a jump in the number of showers generated at 1018.8 eV, where the area covered by the simulations increases. Simulated Log(E) Reconstructed Log(E) 3000 2500 2500 2000 2000 1500 1500 1000 1000 500 0 17 500 18 19 20 (a) 21 22 23 Log(E [eV]) 0 17 18 19 20 21 22 23 Log(E [eV]) (b) Figure 8.3. Reconstructed energies of high-energy simulations. (a) Distribution of simulated energies for showers that are accepted as events. Showers were generated with an E −3 distribution within each bin, although the distribution of detected showers within each bin tends to be flatter because the acceptance is smaller at the lower-energy end. In each bin, 5000 showers passed the data filtering software; some fraction of these were subsequently removed from the simulated data set if the reconstructed vertical angle was ≥ 1 rad. (b) Distribution of reconstructed energies for the same set of showers. 122 Simulated core locations, Log(E [eV]) < 18.8 Reconstructed core locations, Log(E [eV]) < 18.8 y [m] 120 15000 120 20000 y [m] 140 20000 100 15000 100 10000 10000 80 5000 60 0 40 -5000 20 80 5000 60 0 40 -5000 -10000 -50000 20 -40000 -30000 -20000 -10000 0 10000 x [m] 0 -10000 -50000 -40000 -30000 (a) -20000 -10000 0 10000 x [m] 0 (b) Figure 8.4. Reconstructed core locations of high-energy simulations. (a) Distribution of simulated core locations for showers with energy less than 1018.8 eV that are accepted as events. (b) Distribution of reconstructed core locations for the same set of showers. Simulated core locations, Log(E [eV]) > 18.8 Reconstructed core locations, Log(E [eV]) > 18.8 y [m] 140 15000 180 20000 y [m] 160 20000 160 15000 140 120 10000 100 10000 120 100 5000 80 5000 80 60 0 0 60 40 -5000 40 -5000 20 -10000 -50000 -40000 -30000 -20000 (a) -10000 0 10000 x [m] 0 20 -10000 -50000 -40000 -30000 -20000 -10000 0 10000 x [m] 0 (b) Figure 8.5. Reconstructed core locations of high-energy simulations. (a) Distribution of simulated core locations for showers with energy greater than 1018.8 eV that are accepted as events. (b) Distribution of reconstructed core locations for the same set of showers. 123 Simulated Cos(θ) Reconstructed Cos(θ) 2000 2000 1800 1800 1600 1600 1400 1400 1200 1200 1000 1000 800 800 600 600 400 400 200 0 0.4 200 0.5 0.6 0.7 0.8 0.9 1 Cos(θ) 0 0.4 0.5 0.6 (a) 0.7 0.8 0.9 1 Cos(θ) (b) Figure 8.6. Reconstructed vertical angles of high-energy simulations. (a) Distribution of the cosine of the simulated vertical angle θ for showers that are accepted as events. Simulated shower angles were evenly generated over the surface of a sphere down to cos(θ) = 0.4. (b) Distribution of the cosine of the reconstructed angle θ for the same set of showers. The cut on reconstructed angle requiring cos θ ≤ 1 rad has been applied. Simulated φ Reconstructed φ 1400 1400 1200 1200 1000 1000 800 800 600 600 400 400 200 200 0 0 1 2 3 (a) 4 5 6 φ [rad] 0 0 1 2 3 4 5 6 φ [rad] (b) Figure 8.7. Reconstructed azimuthal angles of high-energy simulations. (a) Distribution of simulated azimuthal angle φ for showers that are accepted as events. (b) Distribution of reconstructed azimuthal angle φ for the same set of showers. 124 events and pass all quality cuts) whose distribution on the sky is proportional to the exposure of the array. A map of the relative sky exposure for CHICOS is shown in figure 8.9. 8.4 Estimating the High-Energy Flux The method for finding the high-energy cosmic ray flux as function of energy proceeds as described for the low-energy data in section 7.3. The CHICOS data set consists of 145 showers above 1018.4 eV, the lowest energy for which simulated data was generated. The distribution of measured energies is shown in figure 8.10. The estimated spectrum is plotted in figure 8.11. Error in each bin is estimated by allowing a shift in the reconstructed shower energies of up to ±20%. This is assumed to dominate over other sources of error. The lowest-energy bin is unreliable, as the set of simulated showers did not extend below this energy. This makes it impossible to estimate what fraction of showers may actually be lower-energy showers erroneously reconstructed at this energy, and leaving out this effect will cause the aperture to be underestimated (and the flux to be overestimated). There is no observable change in spectral index out to the limit of the measurements. The flux above 1019 eV is in close agreement with the AGASA results [23], but lies above the more recent HiRes [37] and Auger [43] measurements. CHICOS and AGASA are both ground arrays, which must use simulated shower fronts as the basis for their calibration. HiRes and Auger base their calibration on the total air fluorescence yield. The difference in results between the two types of experiments indicates a possible systematic discrepancy in the two calibration methods. A 40% systematic reduction in the CHICOS energy estimates would yield a spectrum which reproduces the shape measured by HiRes and Auger. The angular reconstruction of showers depends in large part on the relative time of hits at sites involved in the shower, and this part of the shower reconstruction is believed to have a smaller systematic uncertainty than the energy estimate. The estimated angular resolution of CHICOS is approximately 5◦ in the vertical direction 125 Acceptance N(passed)/N(in) Reconstructed Log(E [eV]) Spread in reconstructed energies 10-1 -2 10 23 900 800 22 700 21 600 20 500 19 400 300 18 200 17 -3 10 18.5 19 19.5 20 16 20.5 21 Log(E [eV]) 100 18.5 (a) 19 19.5 20 20.5 21 Simulated Log(E [eV]) 0 (b) Figure 8.8. Acceptance of the CHICOS array. (a) Acceptance of the array as a function of energy. The drop in acceptance at 1018.8 eV is an artifact of the change in area covered by the simulations at that energy. A plot of acceptance multiplied by the area under consideration would be smoothly increasing with energy. (b) Reconstructed energy vs. simulated shower energy. CHICOS Sky Exposure Galactic Latitude 100 180 80 160 60 140 40 120 20 100 0 80 -20 60 -40 -60 40 -80 20 -100 -200 -150 -100 -50 0 50 100 150 200 Galactic Longitude 0 Figure 8.9. Relative sky exposure of the CHICOS experiment. This is the cumulative exposure beginning January 1, 2003 and ending December 31, 2007. The map is based on the reconstructed coordinates of simulated showers that were generated over all visible regions of the sky during that time period. 126 CHICOS Data 30 25 20 15 10 5 0 18.5 19 19.5 20 20.5 21 Log(E [eV]) Figure 8.10. Energy distribution of CHICOS showers between 1018.4 eV and 1020.4 eV. Flux*E^3 [eV^2/m^2/sr/s] CHICOS Energy Spectrum 25 10 1024 18.5 19 19.5 20 20.5 21 Log(E [eV]) Figure 8.11. Flux times E 3 , as measured by the CHICOS array. An error band based on an assumed energy reconstruction error of 20% is shown. The first bin may be artificially high due to the fact that it is at the lower energy limit of the simulations, which will will cause the aperture to be underestimated. 127 and 8◦ in the azimuthal direction. This is sufficient to allow a search for associated point sources, as will be explored in chapter 9. 128 Chapter 9 Correlation of UHECR Data with AGN There has recently been substantial interest recently in the possibility of discovering a correlation between ultra-high energy cosmic rays and nearby active galactic nuclei (AGN). The mechanism by which AGN might accelerate particles to the observed energies is discussed in section 2.3. AGN are among the most powerful phenomena in the universe, and the number in relative proximity to our location makes them an attractive target for correlation searches. In section 9.1, we present the method of correlation analysis used by both the Auger and HiRes Collaborations to search their data for excess alignments with AGN. The detailed results from these experiments are covered in section 9.2. In section 9.3, an analogous correlation search is performed on the CHICOS ultra-high energy data. 9.1 Quantifying the Degree of Correlation A correlation is said to exist between events in a set of cosmic ray data and a collection of AGN if an excess of cosmic ray events falling close to one or more AGN is observed, compared to the expected result from a random, isotropic distribution of data. In order to make a quantitative comparison, 3 parameters must be chosen: Emin , the energy threshold for the data being considered, zmax , the maximum redshift of AGN included in the search, and θmax , the angular window around each AGN that will be searched for correlated cosmic rays. The set of parameters that yields the greatest 129 correlation for a given data set is found by scanning all possible combinations. To evaluate the degree of correlation for a given set of parameters, it is first necessary to calculate p, the exposure-weighted fraction of the sky covered by all windows of angular radius θmax around the set of AGN with z < zmax . As zmax and θmax become large, p approaches unity. Given a random, isotropic data set of size N , the probability P that k or more data points will lie within angle θmax of an AGN is given by the cumulative binomial distribution N µ ¶ X N j P = p (1 − p)N −j . j j=k (9.1) Minimizing P results in a set of paramters zmax , Emin , and θmax for which there is a maximum alignment between the cosmic ray data and the set of AGN. However, Pmin cannot be directly interpreted as the probability that the data set is isotropically distributed, because the values that minimize P depend on the characteristics of the data set. To determine the true probability Pchance that a data set has only random correlations, it is necessary to generate a large number of simulated data sets (with the same distribution of energy as the real data set, and drawn from an isotropic distribution proportional to the array’s exposure pattern on the sky). The fraction of simulated data sets with a smaller Pmin than the real data set gives the estimated probability Pchance that the observed correlation in the data set is due to chance alone [30]. The AGN catalogue used by CHICOS and other experiments for this purpose is the Véron-Cetty and Véron (VCV) Catalogue [140]. This catalogue lists 85221 quasars, 1122 BL Lac objects, and 21737 active galaxies. Of these objects, 694 have redshift z ≤ 0.024, the maximum redshift considered. This redshift corresponds to a distance of approximately 100 Mpc for a Hubble constant of H0 = 71 km s−1 Mpc−1 . The VCV catalogue is a compilation of individual survey results and is therefore not a homogeneous sample. It is notably incomplete near the plane of the galaxy as well as at distances greater than approximately 100 Mpc. 130 9.2 Recent Results from Other Experiments The Auger Collaboration performed several parameter searches on their data and concluded that a statistically significant correlation was observed [47, 48]. An exploratory scan yielded a maximum correlation for the parameter values Emin = 56 eV, zmax = 0.018, and θmax = 3.1◦ . Of the 15 events in this set, 12 were correlated with at least one AGN, where only 3.2 correlated events would be expected from an isotropic flux. The parameters found in the exploratory scan were then applied to an independent data set. Given these parameters, 8 out of 13 cosmic ray events were correlated with at least one AGN in the independent data set, where 2.7 alignments would be expected by chance. The probability for this degree of correlation to be observed in an isotropic flux is given as Pchance = 1.7 × 10−3 . Having determined, according to these results, that their data set shows correlated events, the Auger Collaboration performed a parameter search over their entire data set, consisting of 81 events with energy above 40 EeV and zenith angle smaller than 60◦ . Their scan covered the parameter ranges 1◦ ≤ θmax ≤ 8◦ , 0 ≤ zmax ≤ 0.024, and Emin ≥ 40 EeV. Maximum correlation was found for the parameter set Emin = 57 EeV, zmax = 0.017, and θmax = 3.2◦ . For this paramter set, 20 of 27 cosmic ray events above Emin are aligned with at least one AGN, where only 5.6 would be expected by chance. The cumulative binomial probability of observing that degree of correlation from an isotropic distribution is Pmin = 4.6 × 10−9 . After performing a similar analysis on a large number of simulated data sets, they report that the true probability of the observed correlation arising from an isotropic distribution is approximately Pchance = 10−5 . The HiRes experiment has followed a similar procedure [49] to search for correlations between their data and AGN. They begin by evaluating the correlation of their data using the parameters found in the Auger exploratory scan. For the parameter set Emin = 56 eV, zmax = 0.018, and θmax = 3.1◦ , they observed 2 of 13 events that were aligned with at least one AGN. They evaluated the chance probability of this result by observing the fraction of simulated data sets (out of 5000) that yielded a greater 131 number of alignments for the same set of parameters. They found that this result has a chance probability of Pchance = 0.82 and is consistent with random correlations from an isotropic flux. Next they divided their data into two independent sets. The first set yielded a maximum correlation for the parameter set Emin = 15.8 EeV, zmax = 0.020, and θmax = 1.7◦ , with 20 correlated events from a total of 97. Applying these parameters to the second data set, they found 14 correlated events from a total of 101. The estimated chance proability of that level of correlation was Pchance = 0.15, again consistent with random correlations. Finally, the HiRes Collaboration performed a parameter search over their entire data set. The maximum correlation was found at Emin = 15.8 EeV, zmax = 0.016, and θmax = 2.0◦ , with 36 correlated events out of 198. The probability (based on the probability for a simulated data set to show greater correlation for this particular set of paramters) is Pmin = 1.8 × 10−3 . The true probability of the observed correlation is determined by conducting a full parameter search on each simulated data set, and observing the fraction that yield a greater degree of correlation for any parameter combination. This process results in an estimated probability of the observed correlation arising by chance of Pchance = 0.24. It is concluded that the HiRes observations are consistent with an isotropic cosmic ray flux, showing only random correlations with AGN. 9.3 Results from CHICOS Data Given that the Auger and HiRes conclusions are incompatible, more data are clearly needed in order to confirm or reject the hypothesis that nearby AGN are cosmic ray sources. We use the CHICOS data to perform an analysis similar to that done by these two experiments. In order to evaluate equation (9.1), we must tabulate the value of p, the exposureweighted fraction of the sky that falls within angle θ of any AGN with z < zmax . As shown in figure 9.1, p increases linearly with θ for small zmax , but for large zmax 132 p Fraction of the sky within θ deg of an AGN 0.9 0.8 0.7 0.6 0.5 0.4 0.3 0.2 0.1 0 1 2 3 4 5 6 7 8 9 10 θ [deg] Figure 9.1. Exposure-weighted fraction of the sky covered by windows of angular radius θmax centered on nearby AGN. The curves correspond to different maximum redshifts, in increments of 0.001, from zmax = 0.000 to 0.024. (when more AGN are included), p approaches unity as the angular windows begin to overlap. The value of p(θ, z) has been calculated by observing the fraction of simulated showers (distributed proportionally to the sky exposure) that fall within the specified windows on the sky. We first examine the CHICOS data using the results of the Auger exploratory scan. For the parameter set Emin = 56 eV, zmax = 0.018, and θmax = 3.1◦ , the CHICOS data contains 3 correlated events out of 23. The expected number of correlations is 4.2, and the cumulative binomial probability is P = 0.81. Next we examine the CHICOS data using the results of the Auger scan of their complete data set. For the parameter set Emin = 57 eV, zmax = 0.017, and θmax = 3.2◦ , the CHICOS data contains 3 correlated events out of 22. The expected number of correlations is 3.8, and the cumulative binomial probability is P = 0.77. Both of these results are consistent with an isotropic flux. Finally, we scan the CHICOS data set over energy cutoff, angular window, and redshift, to determine the parameter set that minimizes P . The energy cutoff Emin 133 was scanned from 10 EeV to 100 EeV in increments of 1 EeV. For the smallest value of Emin , the CHICOS data set contains 84 events. There are 11 events above the highest value of Emin . The angular window θmax was scanned from 1◦ to 10◦ in increments of 0.1◦ . The redshift cutoff zmax was scanned from 0 to 0.024 in increments of 0.001. The maximum correlation was found for Emin = 100 eV, zmax = 0.009, and θmax = 1.9◦ . For this configuration, 4 out of 11 events were aligned with an AGN, where the expected number of correlations was 0.5. The cumulative binomial probability is Pmin = 8.6 × 10−4 . The parameter scans for this result are shown in figure 9.2. In order to calculate the true probability of observing this result from a random isotropic distribution, we performed the same analysis on 1000 simulated data sets. The simulated data sets contained a set of events that were matched in energy to the CHICOS data set. Galactic latitude and longitude were drawn from the set of simulated showers used to estimate the aperture of the array, described in section 8.3. It was found that the chance probability of the observed degree of correlation was approximately Pchance = 0.21. This is consistent with random correlations from an isotropic flux. A sky map the CHICOS data above 1019 eV is shown in figure 9.3. The sky map of the set of AGN and CHICOS data for which there is a maximum correlation is shown in figure 9.4. 134 Log10P(θ) P(θ), for E = 100 EeV, z = 0.009 0 -0.5 -1 -1.5 -2 -2.5 -3 1 2 3 4 5 6 7 8 9 10 θ [deg] 50 60 70 80 90 100 E [eV] Log10P(E) P(E), for θ = 1.9 deg, z = 0.009 -1 -1.5 -2 -2.5 -3 ×10 18 10 20 30 40 Log10P(z) P(z), for E = 100 EeV, θ = 1.9 deg 0 -0.5 -1 -1.5 -2 -2.5 -3 0 0.005 0.01 0.015 0.02 0.025 z Figure 9.2. Cumulative binomial probability P of the observed correlation resulting from a random isotropic distribution. The probability P is a function of the angular window θmax , the energy threshold Emin , and the redshift cutoff zmax . P is plotted vs. each of these variables, with the other two being held constant at the values that result in the minimum probability Pmin (Emin = 100 EeV, θmax = 1.9◦ , zmax = 0.009). 135 Sky Map of CHICOS Data Galactic Latitude 100 180 80 160 60 140 40 120 20 100 0 80 -20 60 -40 -60 40 -80 20 -100 -200 -150 -100 -50 0 50 100 150 200 Galactic Longitude 0 Figure 9.3. The set of CHICOS data points used in the correlation search. Data points are shown as open circles superimposed on a map of the relative exposure on the sky. 136 Sky Map of AGN and CHICOS data Galactic Latitude 100 80 60 40 20 0 -20 -40 -60 -80 -100 -200 -150 -100 -50 0 50 100 150 200 Galactic Longitude Figure 9.4. The set of 11 CHICOS data points with E > 100 EeV and the set of 220 AGN with z < 0.009 for which there is a maximum correlation. Data points are shown as open circles; AGN are shown as blue dots. There are 4 events that fall within the angular window θmax = 1.9◦ of at least one AGN. (One event falls within 1.9◦ of 2 AGN.) The galactic coordinates (l, b) of the correlated cosmic ray events are (-67.37◦ , 67.77◦ ), (168.24◦ , 43.89◦ ), (80.13◦ , 19.48◦ ), and (132.87◦ , 44.20◦ ). 137 Chapter 10 Conclusions The CHICOS project has been in active operation since 2003 and continues to record cosmic ray air showers. During that time, CHICOS has fulfilled many scientific and educational goals, summarized below. New lateral distribution and time distribution functions have been developed for use with the CHICOS reconstruction software. The time distribution function in particular represents a significant improvement over previous descriptions of particle arrival times. These particular parameterizations of the distribution functions describe the shape of an air shower front at the altitude of the CHICOS experiment, but the general form is expected to be applicable to other ground-array experiments. Event reconstruction software was written for use with CHICOS. The accuracy of this software has been assessed with unthinned simulated air showers at 1017 eV. At higher energies, the accuracy of event reconstructions and the acceptance of the array were estimated using simulated showers based on the lateral distribution function with added statistical noise. The Chiquita subarray was used to measure the energy spectrum between 1017 and 1019 eV. At the lower end of this range, the results were found to be comparable to previous experiments. At the upper end of this range, the measured flux fell below previously measured values. This is considered to be due to limitations imposed by the configuration of the Chiquita array, which made it difficult to estimate the shower core location and incidence angle with sufficient precision. The high-energy spectrum above 1018.4 eV was measured using data from the entire 138 CHICOS array between January 1, 2003 and December 31, 2007. A total of 145 showers were included in the final data set (appendix B). The resulting spectrum agreed very closely with results originally published by AGASA, but fell significantly above more recent results from HiRes and the Pierre Auger Observatory. This suggests that the calibration of ground-array data varies systematically from the calibration of air-fluorescence data. Finally, CHICOS data above 1019 eV was used in a correlation analysis with AGN within 100 Mpc. The Pierre Auger Observatory has previously reported a significant correlation between their ultra-high energy data and AGN in the southern hemisphere, while HiRes has found their data to be constistent with no correlation. In the CHICOS data, no significant excess of correlated events were observed for any values of energy threshold and redshift cutoff. The maximum correlation was observed for events above 1020 eV, AGN with z ≤ 0.009, and an angular window of 1.9◦ . In addition to scientific goals, CHICOS has encorporated a strong educational component. CHICOS data is made available to participating schools so that students and teachers can access it for individual projects. Through methods including teacher workshops, week-long summer sessions for high school students, and online tutorials, CHICOS has been active in bringing scientific research to the community. 139 Appendix A Site Locations and Parameters The following series of tables lists site locations and parameters for all sites in the CHICOS array. Tables A.1 and A.2 list each site location in Cartesian coordinates (relative to the location of Site 03, at longitude -118.12522◦ , latitude 34.13657◦ , and an elevation of 245 m above sea level), the installation and activation dates, the date of the most recent calibration and (where applicable) the deactivation date. Note that two sites appear twice in the list. Site 74 was renamed from MaranathaHS to Judson on 09.01.05, and site 86 was moved to a new location on 03.27.07. Tables A.3 and A.4 list the following parameters characterizing each site: scintillator thickness (in cm), scintillator area (in m2 ), the pulse-height calibration constant C, and the decay constant τ . The parameters C and τ are used in the conversion between time-over-threshold (in ns) and pulse energy (in MIP) as shown in equation (6.9). Site name Caltech03 SylmarHS MonroeHS HarvardWestlake AlhambraHS AlvernoHS ArcadiaHS BirminghamHS BlairHS LAMissionCollege CanogaParkHS Hillcrest ChaminadeHS SGChristian ChatsworthHS CrespiHS CSUN FultonMS FrancisPolyHS EliotMS AlemanyHS SanFernandoMS GranadaHillsHS GrantHS NorthHollywoodHS HolmesMS HuntingtonMS BuckleySchool EmeryPark WashingtonMS JohnMuirHS ChandlerSchool LaSalle Louisville MadisonMS MaranathaHS Judson MarkKeppelHS MayfieldSchool WilsonMS TowneandCountry Site number 0003 0010 0011 0012 0016 0017 0018 0021 0022 0023 0029 0031 0032 0034 0035 0040 0041 0042 0043 0044 0045 0050 0051 0052 0055 0058 0059 0060 0061 0063 0064 0065 0068 0071 0073 0074 0074 0075 0076 0083 0085 0.00 -29086.00 -32518.28 -26482.67 -300.86 5380.42 8018.06 -35123.67 -1935.43 -27100.87 -44246.97 -35917.81 -47043.76 3015.51 -42429.31 -36459.06 -37144.68 -30633.11 -25591.59 -476.63 -30958.86 -28615.25 -36747.99 -26740.62 -24248.20 -34968.52 1066.52 -29391.13 -2727.18 -2009.10 -3256.36 -3513.45 5235.92 -44660.61 -26968.61 737.61 737.61 927.72 -3307.05 4125.64 2256.93 x [m] 0.00 18971.27 11246.86 429.92 -4929.60 3236.27 -796.37 5778.39 -1113.53 19808.57 6611.15 15944.09 8040.20 -3610.19 12989.22 2859.58 11389.41 7870.40 9669.08 5219.36 15386.12 16300.61 13896.07 4750.24 3321.95 11599.53 -2120.88 384.20 -5267.53 3962.67 4637.47 2455.53 2688.17 1526.93 6602.33 4043.46 4043.46 -7334.41 -667.21 550.00 735.73 y [m] 0.00 42.44 -70.76 -83.93 -90.92 63.91 -114.08 -118.23 -4.39 91.62 -151.75 -23.02 -177.39 -111.74 -104.24 -98.73 -63.25 -86.36 -46.64 133.84 -31.68 2.98 -51.94 -92.77 -61.41 -89.34 -72.44 -69.66 -97.77 94.45 70.48 27.56 19.28 -76.37 -81.38 80.67 80.67 -125.30 10.11 -45.36 -32.44 z [m] 08.21.01 09.10.01 09.10.01 09.20.01 12.09.02 01.18.02 03.04.02 10.11.01 05.15.02 06.06.03 03.27.03 11.26.02 01.10.03 03.30.04 06.25.03 11.21.02 02.15.02 12.02.03 06.01.05 06.11.02 05.26.04 11.19.03 03.10.03 10.11.01 05.11.05 03.26.03 05.23.02 07.02.03 03.03.05 04.06.04 06.19.02 11.06.03 11.22.02 09.20.01 03.18.04 03.13.03 09.01.05 09.03.03 04.18.03 03.05.04 12.18.03 Installation Date [mm.dd.yy] 08.21.01 09.10.01 09.10.01 09.20.01 12.09.02 01.18.02 03.04.02 10.11.01 05.15.02 06.06.03 03.27.03 11.26.02 01.10.03 03.30.04 06.25.03 11.21.02 02.15.02 12.02.03 06.01.05 05.17.04 05.26.04 11.19.03 03.10.03 10.11.01 05.11.05 03.26.03 05.23.02 07.02.03 03.03.05 04.06.04 06.19.02 11.06.03 11.22.02 09.20.01 03.18.04 03.13.03 09.01.05 09.03.03 04.18.03 03.05.04 01.13.04 Activation Date [mm.dd.yy] Table A.1. Database of site locations, part I 12.01.03 01.14.04 02.23.04 01.07.04 02.17.04 01.13.04 02.25.04 02.23.04 02.04.04 01.15.04 01.20.04 01.22.04 02.12.04 03.30.04 01.20.04 02.12.04 01.26.04 01.28.04 06.01.05 00.00.00 05.26.04 11.09.03 03.08.04 02.10.04 05.11.05 02.10.04 01.08.04 02.19.04 03.03.05 04.06.04 12.03.02 01.08.04 12.03.02 02.05.04 03.18.04 01.12.04 01.12.04 02.02.04 01.06.04 03.05.04 02.09.04 Calibration Date [mm.dd.yy] – – – – – – – – – – – – – – – – – – – – – – – – – – – – 12.12.06 – 02.03.05 – – – – 09.01.05 07.28.06 – – – 08.04.05 Deactivation Date [mm.dd.yy] 140 Site name NotreDameHS NotreDameHS MarshallFundamental RamonaElementary PacoimaMS SanFernandoHS PasadenaHS PolytechnicSchool RamonaConvent HaleMS SGAcademy LongfellowElem PCC SanGabrielMissionHS SanMarinoHS ShermanOaksCES SouthPasadenaHS Southwestern MadisonElementary NorthridgeMS ChaminadeMS PierceCollege VanNuysHS VanNuysMS TempleCityHS WestridgeSchool HugoReid NewHorizon MilkenCommunityHS CenturyHS BarnhartSchool WilsonElementary CIT1001 CIT1002 CIT1003 CIT1004 CIT1005 CIT1006 CIT2001 CIT3001 CIT4001 Site number 0086 0086 0087 0088 0090 0092 0093 0096 0099 0100 0101 0103 0105 0108 0109 0111 0112 0114 0116 0121 0122 0123 0124 0125 0128 0129 0130 0131 0132 0133 0152 0179 1001 1002 1003 1004 1005 1006 2001 3001 4001 -28135.08 -28177.20 1289.12 -325.83 -27653.77 -29251.57 3175.77 -392.96 -1378.16 -48089.79 4465.35 -331.94 454.71 1703.49 2432.67 -37996.85 -2594.60 -545.27 -1276.27 -36702.08 -40532.05 -41628.90 -30274.51 -30143.39 5679.41 -2913.30 5674.48 -2927.05 -32566.27 -1505.24 7866.41 3392.37 -96.48 -100.91 -92.00 -89.77 -96.48 -103.12 -325.61 -198.32 -17.52 x [m] 2414.80 2514.76 3016.28 -6653.37 12533.68 15080.04 2916.87 -457.02 -7073.70 4654.33 -4066.62 3650.74 731.03 -4181.80 -1195.49 5380.23 -3002.47 -2353.87 2504.88 9750.21 13273.90 5485.78 5883.57 3621.27 -1972.79 -1168.74 115.97 2125.87 -1146.19 -5133.60 777.56 -2071.51 -51.03 -48.03 -48.03 -52.58 -54.02 -52.58 39.94 -337.22 394.91 y [m] -104.45 -104.67 50.15 -118.49 -33.24 -15.89 8.54 -17.03 -102.09 -156.80 -121.86 72.94 -1.06 -94.60 -45.58 -123.33 -43.24 -64.46 35.38 -117.96 -0.51 -149.77 -89.50 -109.18 -97.83 -9.77 -75.52 23.97 34.85 -82.25 -94.89 -88.24 -18.50 -18.50 -18.50 -18.50 -18.50 -18.50 -4.01 -15.01 -1.01 z [m] 03.18.03 03.27.07 11.13.03 02.04.04 03.31.03 09.11.03 06.20.02 11.26.01 05.24.02 06.25.04 04.27.04 10.31.05 02.07.03 12.04.02 04.17.02 09.20.01 04.11.02 07.19.02 01.05.05 01.06.05 04.16.03 05.06.05 10.11.01 11.20.02 04.24.03 01.01.01 09.27.04 04.07.05 03.04.03 12.05.02 11.24.04 10.08.04 05.15.03 05.15.03 05.15.03 05.17.03 05.17.03 05.17.03 05.24.03 05.29.03 05.28.03 Installation Date [mm.dd.yy] 03.18.03 03.27.07 11.13.03 02.05.04 03.31.03 09.11.03 06.20.02 11.26.01 05.24.02 06.25.04 04.27.04 10.31.05 02.07.03 12.04.02 04.17.02 09.20.01 04.11.02 07.19.02 01.05.05 07.07.05 04.16.03 05.11.05 10.11.01 11.20.02 04.24.03 01.01.01 09.27.04 04.07.05 03.04.03 12.05.02 11.24.04 10.08.04 05.15.03 05.15.03 05.15.03 05.17.03 05.17.03 05.17.03 05.24.03 05.29.03 05.28.03 Activation Date [mm.dd.yy] Table A.2. Database of site locations, part II 02.19.04 02.19.04 03.08.04 02.04.04 02.03.04 12.04.03 01.09.04 02.09.04 02.18.04 06.25.04 04.27.04 07.18.05 02.07.03 01.28.04 01.29.04 01.07.04 02.17.04 01.29.04 01.05.05 07.07.04 02.12.04 05.11.05 01.27.04 01.27.04 02.02.04 01.06.04 08.03.04 04.07.05 02.03.04 01.21.04 07.20.04 08.10.04 10.28.03 10.28.03 05.19.03 05.19.03 10.29.03 10.29.03 12.05.03 12.08.03 12.05.03 Calibration Date [mm.dd.yy] 03.27.07 – – – – – – – – – – – – – – – – – – – – – – – – – – – – – – – – – 12.12.06 – – 12.12.06 12.12.06 – 06.03.04 Deactivation Date [mm.dd.yy] 141 Site name Caltech03 SylmarHS MonroeHS HarvardWestlake AlhambraHS AlvernoHS ArcadiaHS BirminghamHS BlairHS LAMissionCollege CanogaParkHS Hillcrest ChaminadeHS SGChristian ChatsworthHS CrespiHS CSUN FultonMS FrancisPolyHS EliotMS AlemanyHS SanFernandoMS GranadaHillsHS GrantHS NorthHollywoodHS HolmesMS HuntingtonMS BuckleySchool EmeryPark WashingtonMS JohnMuirHS ChandlerSchool LaSalle Louisville MadisonMS MaranathaHS MarkKeppelHS MayfieldSchool WilsonMS TowneandCountry Site number 0003 0010 0011 0012 0016 0017 0018 0021 0022 0023 0029 0031 0032 0034 0035 0040 0041 0042 0043 0044 0045 0050 0051 0052 0055 0058 0059 0060 0061 0063 0064 0065 0068 0071 0073 0074 0075 0076 0083 0085 10.16 5.00 10.20 10.20 5.08 7.60 7.60 5.00 8.90 8.89 7.60 6.35 5.08 5.08 5.08 10.20 10.20 8.90 5.20 5.10 7.50 10.00 10.16 10.00 7.00 10.16 10.20 10.16 7.60 4.00 5.10 5.08 5.08 5.08 5.08 6.35 7.62 7.62 5.00 9.00 Scintillator Thickness A [cm] 0.86 1.00 0.85 1.00 1.03 0.60 1.05 1.03 0.69 0.60 0.56 0.53 1.04 1.03 1.03 0.70 0.54 0.82 1.05 1.03 0.83 1.05 1.03 1.00 0.81 1.03 1.03 1.03 0.62 1.03 1.03 1.03 1.03 1.03 1.03 0.59 0.62 0.62 1.08 1.07 Scintillator Area A [m2 ] 15.00 24.00 22.50 20.20 17.30 23.00 24.00 24.00 22.00 15.00 18.50 20.00 20.00 20.00 22.00 18.00 22.00 13.00 28.00 11.00 19.00 18.00 23.00 24.00 17.50 18.00 27.00 3.00 33.00 20.00 24.00 16.00 21.50 40.50 21.50 20.00 20.00 18.00 26.00 15.00 Pulse-Height Calibration A 83.30 92.50 100.60 81.00 85.60 98.30 105.50 88.90 94.80 79.10 84.90 68.60 85.80 81.80 81.00 81.60 81.90 75.70 91.10 77.00 77.70 89.70 93.40 73.80 84.20 86.70 97.20 31.80 82.30 83.30 100.90 91.80 94.30 103.50 88.20 82.50 76.10 87.70 85.20 80.90 Decay Constant A [ns] 7.62 5.00 7.60 7.60 7.62 7.60 5.10 7.00 8.30 7.62 5.10 7.62 10.16 5.08 5.08 8.90 10.20 8.90 5.20 8.30 10.00 5.00 10.16 8.00 8.00 10.16 7.60 5.08 7.60 5.00 7.60 5.08 5.08 10.20 5.08 7.62 7.62 7.62 5.00 4.30 Scintillator Thickness B [cm] Table A.3. Database of site parameters, part I 0.82 0.62 0.58 0.80 0.82 0.82 0.84 0.82 0.62 0.80 1.03 0.61 0.64 1.03 1.03 0.60 1.03 0.82 1.03 0.81 0.84 1.03 1.03 0.85 0.81 1.03 0.81 1.04 0.82 1.05 0.81 1.03 0.84 0.62 1.03 0.61 0.82 0.62 1.08 1.02 Scintillator Area B [m2 ] 20.00 17.30 18.40 22.00 19.00 22.00 22.40 27.00 15.00 20.00 19.00 16.00 20.00 23.00 25.00 27.00 24.00 19.00 28.00 6.00 16.00 25.00 25.00 20.00 21.00 19.00 25.00 8.00 15.00 27.00 15.00 16.00 27.00 25.00 18.00 19.00 20.00 22.00 22.00 23.00 Pulse-Height Calibration B 79.30 71.20 83.10 77.10 66.90 91.50 96.50 98.90 93.80 88.00 79.80 75.80 78.80 85.80 82.50 100.40 84.20 79.10 92.60 77.00 71.20 83.70 84.60 78.70 83.80 90.50 98.40 45.10 74.80 89.50 92.50 86.60 97.20 95.80 90.00 79.80 91.40 85.80 88.70 85.10 Decay Constant B [ns] 142 Site name NotreDameHS MarshallFundamental RamonaElementary PacoimaMS SanFernandoHS PasadenaHS PolytechnicSchool RamonaConvent HaleMS SGAcademy LongfellowElem PCC SanGabrielMissionHS SanMarinoHS ShermanOaksCES SouthPasadenaHS Southwestern MadisonElementary NorthridgeMS ChaminadeMS PierceCollege VanNuysHS VanNuysMS TempleCityHS WestridgeSchool HugoReid NewHorizon MilkenCommunityHS CenturyHS BarnhartSchool WilsonElementary CIT1001 CIT1002 CIT1003 CIT1004 CIT1005 CIT1006 CIT2001 CIT3001 CIT4001 Site number 0086 0087 0088 0090 0092 0093 0096 0099 0100 0101 0103 0105 0108 0109 0111 0112 0114 0116 0121 0122 0123 0124 0125 0128 0129 0130 0131 0132 0133 0152 0179 1001 1002 1003 1004 1005 1006 2001 3001 4001 9.90 9.00 5.00 3.81 7.62 10.20 8.90 5.10 8.00 7.00 8.00 7.62 5.08 5.10 3.00 10.16 10.20 10.00 5.00 7.62 10.50 5.08 5.10 5.08 10.16 5.00 7.80 8.89 5.08 5.00 10.50 5.08 5.08 6.99 7.91 5.08 9.31 10.20 7.62 7.62 Scintillator Thickness A [cm] 0.73 0.68 1.05 0.95 0.82 0.62 0.71 1.01 1.09 0.68 0.59 0.78 1.03 1.02 0.99 1.03 1.03 1.04 1.06 0.82 1.00 1.03 1.03 1.03 1.04 1.01 0.61 0.61 1.03 1.12 0.53 1.03 1.03 1.03 0.73 1.03 0.77 0.62 0.82 0.82 Scintillator Area A [m2 ] 12.00 18.00 17.00 22.00 20.00 18.00 24.00 24.30 17.00 36.00 24.00 9.00 20.00 25.00 26.00 21.50 18.20 22.00 21.00 15.00 30.00 50.00 35.00 20.00 20.00 18.00 29.00 20.00 20.00 18.00 25.00 23.00 25.00 24.00 17.00 29.00 25.00 14.00 10.00 10.00 Pulse-Height Calibration A 75.40 79.10 79.60 82.60 96.30 102.90 104.80 100.10 81.30 83.50 83.30 78.50 81.10 101.50 95.30 102.60 96.70 85.30 77.40 81.10 87.40 86.20 99.50 90.00 105.20 87.00 84.80 79.60 81.50 92.30 83.60 87.00 86.60 83.90 86.76 85.60 85.50 83.80 79.00 78.60 Decay Constant A [ns] 7.62 4.80 6.00 10.20 7.62 7.60 7.00 10.20 6.50 5.50 7.25 7.62 8.89 7.60 10.00 7.62 5.10 5.00 9.00 7.62 8.35 7.62 7.62 8.26 7.62 4.00 7.80 7.62 5.08 5.00 10.50 7.79 5.08 7.62 5.08 5.08 5.08 6.70 10.20 10.16 Scintillator Thickness B [cm] Table A.4. Database of site parameters, part II 0.82 1.02 0.80 0.64 0.82 0.81 1.02 1.04 1.01 1.00 0.82 0.60 0.85 0.83 0.65 0.91 0.81 1.04 0.55 0.82 0.54 0.82 0.80 0.58 0.82 0.99 0.61 0.85 0.89 1.04 0.53 0.75 1.03 0.77 1.03 1.03 1.03 0.89 0.63 0.54 Scintillator Area B [m2 ] 14.00 19.00 18.00 17.00 18.00 18.00 24.00 25.60 13.00 24.00 24.00 8.00 12.00 32.00 26.00 31.50 27.00 22.00 25.00 20.00 30.00 38.00 35.00 22.00 22.00 25.00 24.00 20.00 22.40 22.50 32.00 13.00 30.00 18.00 22.00 26.00 24.50 12.00 12.00 10.00 Pulse-Height Calibration B 82.00 77.80 85.70 78.50 105.70 89.20 99.40 100.40 88.30 85.20 78.80 77.40 79.10 82.80 89.00 94.40 85.40 80.00 79.20 82.00 75.50 93.80 101.30 81.10 99.20 88.10 84.00 83.20 82.00 87.50 89.90 91.40 85.00 80.10 89.10 85.90 89.70 77.00 82.30 75.20 Decay Constant B [ns] 143 144 Appendix B CHICOS Showers Above 1018.4 eV The following series of tables lists all CHICOS showers used in the analysis presented in Chapters 8 and 9. Each shower is identified by a date, time, UTC second, and nanosecond, corresponding to the time of the trigger hit in the A detector. Reconstructed shower parameters include energy (in eV); x, y, and z (in meters, with the origin located on site 03, at longitude −118.12522◦ , latitude 34.13657◦ , and 250 m above sea level; the zenith angle θ and the azimuthal angle φ (in radians). Three measures of the quality of the fit are given: the angle/time fit parameter χ2AT , the energy/core-location fit parameter χ2EL , and the value of the negative logarithm of the likelihood function, − log(l). Time 14:42:14 PST 16:47:48 PST 05:36:01 PDT 18:22:08 PDT 14:55:35 PDT 11:01:02 PDT 23:54:04 PDT 11:56:35 PST 08:27:21 PST 15:04:24 PST 15:24:11 PST Date 03-03-03 03-24-03 05-11-03 05-17-03 06-27-03 08-31-03 09-12-03 11-07-03 11-14-03 12-18-03 12-23-03 3129576134 3131398068 3135501361 3136065728 3139595735 3145197662 3146280844 3151079795 3151672041 3154633464 3155066651 UTC Second 521862986 559630966 525599534 567265938 354206868 468585258 640992158 964607040 769742722 952186254 663406676 ns 8.57e+18 5.36e+19 2.59e+18 7.33e+18 2.59e+19 7.55e+19 9.22e+19 2.04e+20 3.98e+19 2.21e+19 6.10e+18 Energy [eV] -2147 -1933 -1064 -1240 -1516 2580 1240 -31341 -30656 -36972 1339 x [m] -560 627 -387 -5448 1150 3998 1222 4319 1337 10289 -1355 y [m] -10 -16 -11 -91 -10 28 -5 -108 -61 -83 -76 z [m] 0.80 0.61 0.17 1.00 0.64 0.63 0.91 0.76 0.51 0.66 0.36 θ [rad] Table B.1. High-energy showers observed in 2003 1.62 -3.00 -0.45 1.24 2.07 -2.49 2.58 -2.98 -0.43 -2.28 0.41 φ [rad] 10.27 0.68 0.04 23.62 3.66 214.26 11.00 0.06 1.58 0.96 16.82 χ2AT 0.52 1.68 0.20 0.75 12.66 0.25 4.29 0.25 1.13 0.35 4.22 χ2EL 40.86 28.11 42.08 34.74 66.35 45.02 119.14 47.07 30.49 57.71 57.60 − log(l) 145 Time 12:44:06 PST 05:24:13 PST 01:06:25 PST 14:30:55 PST 12:44:55 PST 13:00:59 PST 22:47:31 PST 13:34:17 PST 16:29:23 PST 21:02:44 PST 15:28:11 PST 02:05:39 PST 23:43:18 PST 14:25:45 PST 02:31:59 PDT 16:49:51 PDT 04:13:00 PDT 15:32:03 PDT 09:14:21 PDT 19:11:46 PDT 00:14:54 PDT 07:45:03 PDT 03:10:45 PDT 03:46:17 PDT 07:40:32 PDT 03:13:26 PDT 23:49:46 PDT 10:32:37 PDT 02:58:31 PDT 04:43:53 PDT 23:10:56 PST 21:27:58 PST 00:50:02 PST 19:00:28 PST Date 01-01-04 01-05-04 01-08-04 01-26-04 02-05-04 02-05-04 03-01-04 03-02-04 03-03-04 03-07-04 03-21-04 03-23-04 03-23-04 03-26-04 04-10-04 04-28-04 05-25-04 05-25-04 06-04-04 06-04-04 06-12-04 07-16-04 07-19-04 08-24-04 09-02-04 09-24-04 10-06-04 10-12-04 10-21-04 10-30-04 11-27-04 11-28-04 12-05-04 12-13-04 3155834646 3156153853 3156397585 3158001055 3158858695 3158859659 3161054851 3161108057 3161204963 3161566964 3162756491 3162881139 3162958998 3163184745 3164434319 3166040991 3168328380 3168369123 3169210461 3169246306 3169869294 3172833903 3173076645 3176189177 3176980832 3178865606 3179976586 3180447157 3181197511 3181981433 3184470656 3184550878 3185081402 3185838028 UTC Second 882497223 304615070 34546149 849112349 782397209 316862151 829980376 779467455 977150250 755540156 889019314 159671525 274684900 209067216 315704765 966772223 600733938 260252132 670923718 98644908 564464646 123707640 344216719 153983197 703114214 791256207 320257994 172214226 212748564 379820630 679434822 409207174 620271417 781316889 ns 9.77e+18 5.68e+19 1.34e+19 8.37e+18 5.17e+18 5.35e+19 1.52e+19 8.67e+18 3.84e+18 9.06e+18 1.06e+19 1.34e+20 1.11e+19 2.23e+19 2.68e+19 3.36e+18 3.40e+19 6.39e+18 4.41e+18 8.27e+18 4.57e+19 1.76e+19 5.69e+18 3.66e+18 8.07e+18 1.19e+19 1.66e+20 1.61e+20 4.01e+19 4.12e+19 3.29e+18 2.91e+18 3.55e+18 1.53e+19 Energy [eV] -3071 1663 -458 -2190 320 -675 -606 -208 -297 -1973 1042 1971 -1228 1922 3887 797 5430 -2339 -145 5184 -1395 -2658 -1602 -444 -3302 87 -34809 -32209 2692 2092 -312 265 1258 2791 x [m] 333 5097 -5582 -1641 -630 1783 -7288 -962 -6881 -648 848 -6085 333 818 -2625 894 3097 -510 -929 2026 473 2684 -458 -985 -1948 -6634 9748 14118 -585 3957 -619 455 2521 -3880 y [m] 16 61 -95 -8 -12 -10 -118 -17 -118 -10 -8 -118 -3 -32 -98 -2 61 -5 -18 13 -17 24 -5 -17 -10 -118 -86 -32 -36 50 -17 -1 46 -95 z [m] 0.61 0.99 0.54 0.89 0.80 0.17 0.83 0.77 0.34 0.90 0.92 0.52 0.52 0.75 0.48 0.92 0.85 0.48 0.75 0.73 0.74 1.00 0.37 0.59 0.93 0.13 0.29 0.12 0.75 0.87 0.81 0.82 0.44 0.55 θ [rad] Table B.2. High-energy showers observed in 2004 1.80 0.68 1.43 1.55 1.23 1.81 -1.44 -0.45 -1.86 0.25 1.28 2.84 -0.03 0.28 1.07 2.18 -3.14 -0.78 0.54 0.29 -2.01 1.54 2.99 2.83 0.67 1.81 2.13 1.88 -2.99 -2.95 0.11 -0.58 -2.19 -0.83 φ [rad] 25.99 1.91 0.63 3.84 6.00 9.93 18.61 10.06 1.72 4.43 57.69 0.79 1.93 3.26 235.47 11.19 15.88 0.30 2.56 0.89 52.22 62.03 1.53 7.99 1.79 0.49 107.69 0.64 4.97 2.42 4.28 665.04 7.20 0.70 χ2AT 1.05 0.81 15.03 0.36 4.22 55.13 0.65 10.21 1.81 3.12 6.08 0.44 16.66 113.96 0.11 2.95 0.87 2.65 0.99 2.23 97.49 0.56 74.01 3.21 0.24 1.37 0.06 0.11 3.13 2.12 14.26 392.59 0.39 0.84 χ2EL 41.99 28.86 24.30 64.62 80.23 354.52 89.86 157.44 21.37 80.46 65.02 44.69 110.15 79.45 54.47 32.89 63.91 57.36 77.49 20.82 362.40 29.48 49.19 130.52 20.79 24.56 41.63 17.64 101.42 48.19 62.90 397.79 35.61 21.73 − log(l) 146 Time 19:01:09 PST 07:07:28 PST 01:29:25 PST 15:43:25 PST 07:50:45 PST 13:05:08 PST 20:13:07 PST 09:58:21 PST 10:53:42 PST 15:31:10 PST 23:51:43 PST 16:24:20 PST 18:46:54 PDT 02:31:59 PDT 21:57:15 PDT 01:36:47 PDT 03:28:32 PDT 07:41:20 PDT 16:32:53 PDT 05:56:15 PDT Date 01-03-05 01-04-05 01-06-05 01-13-05 01-21-05 01-29-05 02-03-05 02-06-05 02-22-05 03-03-05 03-07-05 04-02-05 04-03-05 04-05-05 04-15-05 04-27-05 05-18-05 05-24-05 06-01-05 06-21-05 3187652469 3187696048 3187848565 3188504605 3189167445 3189877508 3190335187 3190557501 3191943222 3192737470 3193113103 3195332660 3195424014 3195538319 3196472235 3197435807 3199256912 3199790480 3200513573 3202203375 UTC Second 936948536 337624576 771148700 201000887 441977521 593985267 376768045 213946237 276357220 578094031 219651862 847091825 402424237 603402274 34878905 176466607 343114478 629585008 254318701 798893723 ns 1.44e+19 2.68e+18 4.44e+19 1.68e+19 1.54e+19 3.68e+19 9.72e+18 2.13e+19 1.80e+19 1.19e+19 9.89e+18 3.40e+18 3.30e+19 1.06e+19 4.72e+18 7.77e+18 2.75e+18 1.95e+19 4.01e+19 1.11e+20 Energy [eV] -1259 593 -30742 6277 -26737 5837 6241 -599 1392 7115 815 -429 -36086 -1132 -1358 2170 -2708 -1924 -2161 -39547 x [m] -3395 4436 2603 2085 1641 934 927 5100 4687 -96 -2683 644 13348 -429 2095 798 2063 3037 -6986 12181 y [m] -41 74 -106 10 -89 -45 -70 126 50 -76 -56 -5 -52 -4 35 -31 28 35 -100 -14 z [m] 0.78 0.42 0.92 0.86 0.42 0.83 0.52 0.83 0.55 0.49 0.95 0.43 0.84 0.66 0.97 0.94 0.90 0.38 0.94 0.90 θ [rad] 0.10 2.51 -1.64 0.02 -2.34 2.86 -2.86 0.17 2.72 -0.23 0.03 0.14 0.63 1.88 -0.72 2.38 2.22 3.03 -1.65 1.40 φ [rad] Table B.3. High-energy showers observed in 2005, part I 4.71 1.67 890.96 1.86 0.18 3.29 1.27 1.67 8.30 0.03 7.24 3.36 1.45 6.81 11.42 8.08 0.07 1.44 3.39 233.89 χ2AT 2.80 0.55 0.25 0.32 0.38 9.59 0.98 5784.22 1.14 0.34 2.16 118.58 1.30 553.27 1.65 299.37 0.65 1.48 3.03 0.21 χ2EL 181.55 181.55 47.31 22.57 20.11 139.32 41.71 27.42 47.88 7.89 38.80 131.08 18.56 141.23 30.67 22.27 27.07 69.63 29.54 46.43 − log(l) 147 Time 12:12:40 PDT 08:02:18 PDT 22:22:07 PDT 03:56:41 PDT 08:45:18 PDT 09:42:15 PDT 20:55:44 PDT 16:14:22 PDT 03:55:55 PDT 14:45:18 PDT 02:22:45 PDT 06:02:06 PDT 01:20:25 PDT 15:53:05 PDT 14:06:43 PDT 01:14:01 PDT 00:10:08 PST 15:14:30 PST 05:24:18 PST 09:10:10 PST 04:20:51 PST 04:25:36 PST 02:37:22 PST 03:56:16 PST 18:14:10 PST 12:59:00 PST 03:46:07 PST 17:45:31 PST 07:38:44 PST Date 07-26-05 07-29-05 08-17-05 09-20-05 09-20-05 09-20-05 09-22-05 09-23-05 09-27-05 09-30-05 10-03-05 10-07-05 10-11-05 10-14-05 10-15-05 10-17-05 11-03-05 11-08-05 11-11-05 11-22-05 11-30-05 12-03-05 12-15-05 12-17-05 12-18-05 12-21-05 12-22-05 12-24-05 12-26-05 3205249960 3205494138 3207187327 3210058601 3210075918 3210079335 3210292544 3210362062 3210663355 3210961518 3211176165 3211534926 3211863625 3212175185 3212255203 3212381641 3213850208 3214336470 3214560258 3215524210 3216198051 3216457536 3217487842 3217665376 3217803250 3218043540 3218096767 3218319931 3218456324 UTC Second 50163864 715240186 664168834 772244089 728609320 875578744 923107009 354083715 90813075 641149903 397295727 912899013 626503801 556059825 428772600 263761630 372906239 760940724 445117735 741145311 725973382 893569734 318694028 358833368 10904031 816517131 974834114 982756520 989180794 ns 5.91e+18 3.34e+19 7.79e+18 3.69e+19 1.45e+19 2.89e+19 5.91e+19 5.40e+18 4.70e+18 9.52e+19 4.34e+19 2.64e+18 1.51e+19 3.11e+18 6.11e+18 9.60e+18 3.37e+18 2.89e+19 3.65e+19 2.00e+20 5.41e+18 8.46e+19 9.27e+18 1.03e+19 2.53e+18 1.58e+19 1.55e+20 8.51e+19 7.77e+18 Energy [eV] 989 -654 -35400 -885 -735 253 -1487 -3499 -3687 -29831 -608 -1416 -3923 -1425 -1957 -676 119 -4010 -31382 -43457 11 -28207 -777 -587 915 2625 55 -2175 -1007 x [m] 32 4959 10714 -3571 -4636 -3323 -4016 -1545 2702 7069 -5778 -1309 3634 -6421 -356 5567 -6768 -719 7089 14476 -6944 4700 -1070 3000 3583 -3458 -8650 -5400 -750 y [m] -1 133 -77 -83 -65 -65 -98 -5 28 -75 -91 -11 22 -118 -15 0 -117 -2 -88 -103 -118 -93 -15 74 81 -112 -119 -98 -18 z [m] 0.21 0.61 0.46 0.73 0.96 0.76 0.88 0.93 0.61 0.98 1.00 0.59 0.51 0.45 0.17 0.64 0.78 0.68 0.97 0.16 0.36 0.64 0.30 0.80 0.39 0.43 0.80 0.65 0.92 θ [rad] 3.05 -2.54 2.50 -2.94 -2.15 -1.75 0.31 -1.85 2.89 -2.31 2.28 1.99 -1.58 -2.53 -0.15 2.53 -2.92 0.68 -0.18 2.97 -1.26 -1.39 -0.72 -2.06 -0.09 -1.51 -0.53 0.07 0.50 φ [rad] Table B.4. High-energy showers observed in 2005, part II 3.80 0.30 0.39 254.18 82.77 66.92 5.29 1.10 0.94 3.72 3.93 14.35 0.13 1.69 2.67 0.00 1.20 1.30 1.72 0.37 0.49 368.18 6.36 1.02 1.19 2.40 3.86 0.99 7.52 χ2AT 2.51 1.31 0.22 9.03 2.41 6.70 4.07 1.18 0.76 1.06 2.09 0.33 0.18 0.24 1.96 2.33 1.29 2.64 1.68 0.60 0.20 1.40 6.64 1.13 7.78 0.13 2.86 1.06 1.91 χ2EL 111.61 28.93 32.20 85.19 70.13 92.03 146.99 56.53 38.38 34.88 64.16 41.45 0.00 84.95 56.80 42.52 24.50 39.39 38.84 12.11 18.66 44.56 169.62 17.73 48.03 71.30 95.49 62.49 59.39 − log(l) 148 Time 22:48:33 PST 16:15:42 PST 17:52:39 PST 23:41:41 PST 19:16:58 PST 06:26:11 PST 12:38:08 PST 14:27:38 PST 06:04:28 PST 00:36:52 PST 01:43:19 PST 20:24:16 PST 01:59:23 PST 11:22:50 PDT 08:44:29 PDT 20:29:52 PDT 15:15:50 PDT 23:26:25 PDT 12:20:45 PDT 17:31:41 PDT 12:55:26 PDT 07:51:47 PDT 14:26:23 PDT 22:49:35 PDT 21:30:15 PDT 19:59:08 PDT 02:26:05 PDT 09:11:12 PDT 12:50:39 PDT 09:14:39 PDT 13:08:24 PDT 19:52:35 PDT 09:33:33 PST 17:08:03 PST 15:18:42 PST 22:03:31 PST 03:31:03 PST Date 01-02-06 01-13-06 02-05-06 02-05-06 02-13-06 02-14-06 02-23-06 02-26-06 03-06-06 03-08-06 03-10-06 03-11-06 03-23-06 04-09-06 04-16-06 04-16-06 04-19-06 04-20-06 05-13-06 05-17-06 05-18-06 05-19-06 05-25-06 06-04-06 06-14-06 06-17-06 07-02-06 07-10-06 07-10-06 08-10-06 08-18-06 08-30-06 11-03-06 11-03-06 11-29-06 12-11-06 12-16-06 3219115713 3220042542 3222035559 3222056501 3222731818 3222771971 3223571888 3223837658 3224498668 3224651812 3224828599 3224982256 3225952763 3227451770 3228047069 3228089392 3228329750 3228445585 3230392845 3230757101 3230826926 3230895107 3231437183 3232331375 3233190615 3233444348 3234677165 3235392672 3235405839 3238071279 3238776504 3239837555 3245420013 3245447283 3247687122 3248748211 3249113463 UTC Second 654684605 471682219 881340359 565110461 897526412 523174082 275950367 804468076 954325652 970034770 669536912 890329074 741541067 701018423 479369026 393290312 463305936 715724547 491206832 57764195 268460827 265661284 836991963 628865504 56319450 849449265 974970639 443594233 430240110 585201792 483758782 535481479 692405330 439622076 838847374 489856685 579478785 ns 1.43e+19 8.26e+19 4.51e+18 1.40e+20 1.39e+19 1.28e+20 7.50e+18 5.51e+18 6.40e+18 1.01e+20 2.08e+19 1.25e+20 5.24e+18 2.62e+19 7.75e+18 4.75e+18 1.04e+19 6.59e+18 2.87e+18 3.39e+18 2.04e+19 1.81e+19 9.97e+19 4.50e+18 6.36e+19 3.66e+18 2.69e+19 7.85e+18 5.43e+18 1.31e+19 3.44e+19 5.95e+19 5.51e+19 4.16e+19 9.22e+18 1.41e+19 1.35e+19 Energy [eV] 3605 -718 1470 -30534 -3483 8470 3286 739 -1016 -35361 -657 1186 -922 1737 -737 -793 942 648 -380 4771 6284 -1182 -43388 -867 -31312 1072 3755 2077 -880 -1028 1070 -1182 8606 4591 -1218 -1786 -2440 x [m] 3314 2322 2502 13432 2559 -188 -2146 23 -5594 11996 1263 156 -6158 181 -754 4331 2262 4057 3580 2412 2596 2755 4420 -519 3886 3804 -2247 -1999 4224 825 4109 -669 -1846 1907 2677 3370 4995 y [m] 8 53 52 -39 28 -111 -88 -0 -82 -89 -18 -1 -117 -25 -10 93 69 81 73 18 41 36 -142 -17 -105 81 -88 -53 73 -18 48 -10 -111 42 35 94 89 z [m] 0.99 1.00 0.67 0.29 0.69 0.58 0.85 0.91 0.19 0.65 0.23 0.99 0.86 0.92 0.84 0.34 0.63 0.74 0.50 0.66 0.96 0.86 0.86 0.70 0.89 0.81 0.42 0.25 0.69 0.54 0.16 0.78 0.55 0.36 0.58 0.70 0.58 θ [rad] Table B.5. High-energy showers observed in 2006 2.04 -2.46 -2.71 1.35 -2.81 0.72 -2.88 -0.36 -3.02 1.64 0.70 -1.59 -1.08 1.94 2.82 2.69 0.24 2.29 -0.05 -0.54 3.06 -2.01 -1.38 1.40 1.85 -2.97 1.30 -1.64 -0.34 -2.01 1.72 1.83 -2.57 -2.15 2.49 1.13 -1.76 φ [rad] 7.98 21.78 0.79 0.42 0.33 0.01 1.09 9.47 0.95 2.04 136.77 14.46 2.13 22.86 8.97 0.33 1.95 1.14 0.96 89.48 2.80 1.44 2.21 33.81 11.13 1.60 1.94 1.59 0.59 35.14 0.10 21.26 0.77 5.89 1.08 1.15 2.16 χ2AT 4.84 12.13 0.72 0.85 0.09 0.12 8.05 1.24 5.45 5.13 23.51 480.02 1.01 3.60 2.02 0.36 0.99 1047.41 1965.55 0.12 0.63 1.06 0.11 3.73 1.20 1.33 0.33 0.50 0.14 55.51 0.16 20.79 0.27 1.86 495.47 0.49 0.11 χ2EL 16.14 157.98 24.76 28.23 42.46 42.46 14.99 85.67 68.53 27.24 390.85 134.53 34.44 61.15 77.57 77.57 90.34 19.47 10.86 25.60 17.35 27.91 26.85 108.49 32.66 18.11 71.38 58.61 58.61 258.98 55.62 294.04 26.06 52.43 20.46 36.57 36.57 − log(l) 149 Time 08:10:13 PST 10:07:12 PST 08:23:49 PST 02:48:37 PST 02:05:03 PST 21:14:47 PST 13:36:38 PDT 03:08:27 PDT 03:33:50 PDT 12:10:44 PDT 14:42:56 PDT 13:17:01 PST 18:10:31 PST 17:47:54 PST Date 01-12-07 01-13-07 02-08-07 03-05-07 03-10-07 03-18-07 05-02-07 06-29-07 07-04-07 09-09-07 09-29-07 10-31-07 10-31-07 11-22-07 3251463013 3251556432 3253796629 3255936517 3256365903 3257126087 3260982998 3265956507 3266390030 3272209844 3273946976 3276710221 3276727831 3278627274 UTC Second 629881825 279798382 461152642 403219213 125097721 100150769 333400124 437648942 665612095 529921593 771410499 518559147 131262702 508377239 ns 1.46e+19 3.92e+18 3.03e+19 1.20e+19 4.06e+18 9.27e+18 2.79e+18 1.48e+19 1.42e+19 1.75e+19 6.57e+18 6.42e+18 1.71e+19 9.00e+19 Energy [eV] -2754 -510 -28288 -327 -3348 749 202 1229 -155 -2838 5197 -2928 6053 -37638 x [m] 4517 901 2705 3519 -1458 -275 -675 1201 -961 186 2558 -1130 1096 13647 y [m] 85 -1 -104 72 -8 -10 -1 23 -18 10 20 -8 -78 -50 z [m] 0.60 0.21 0.95 0.93 0.69 0.59 0.57 0.73 0.87 0.41 0.63 0.72 0.50 0.80 θ [rad] Table B.6. High-energy showers observed in 2007 -2.57 -0.72 -1.09 1.78 0.12 -2.41 3.05 2.16 -2.60 -1.42 -1.76 1.36 2.73 0.47 φ [rad] 1.32 6.65 16.03 2.23 0.66 80.54 2.24 10.70 5.42 2.24 0.69 6.86 3.53 1.29 χ2AT 0.69 1.48 0.31 185.16 0.55 110.97 1.48 0.66 16.82 2.45 0.24 352833.35 17.17 37.68 χ2EL 29.88 206.88 20.72 25.01 30.91 231.22 51.80 36.20 71.52 39.79 22.63 29.97 82.67 18.98 − log(l) 150 151 Bibliography [1] V. F. Hess. Observation of penetrating radiation in seven balloon flights. Z. Phys., 13:1084, 1912. [2] R. A. Millikan and G. H. Cameron. The origin of the cosmic rays. Phys. Rev., 32(4):533–557, October 1928. [3] A. H. Compton. A geographic study of cosmic rays. Phys. Rev., 43(6):387–403, 15 March 1933. [4] P. Auger, R. Maze, and T. Grivet-Mayer. Heavy electrons in cosmic ray showers originating in the atmosphere. C. R. Acad. Sci. Ser. B, 206:1721, 1938. [5] P. Auger and R. Maze. C. R. Acad. Sci. Ser. B, 207:288, 1938. [6] W. Kohlhörster, I. Matthes, and E. Weber. Naturwissenschaften, 26:576, 1938. [7] P. Auger, P. Ehrenfest, R. Maze, J. Daudin, Robley, and A. Fréon. Extensive cosmic-ray showers. Rev. Mod. Phys., 11:288–291, July–October 1939. [8] G. W. Clark, J. Earl, W. L. Kraushaar, J. Linsley, B. B. Rossi, F. Scherb, and D. W. Scott. Cosmic-ray air showers at sea level. Phys. Rev., 122(2):637–654, 15 April 1961. [9] J. Linsley. Evidence for a primary cosmic-ray particle with energy 1020 eV. Phys. Rev. Lett., 10(4):146–148, 15 February 1963. [10] K. Greisen. End to the cosmic-ray spectrum? Phys. Rev. Lett., 16(17):748–750, 25 April 1966. 152 [11] G. T. Zatsepin and V. A. Kuz’min. Upper limit of the spectrum of cosmic rays. Zh. Eskp. Teor. Fiz., 4:114–117, 1 August 1966. JETP Lett. 4(3):78–80, 10 August 1966. [12] M. A. Lawrence, R. J. Reid, and A. A. Watson. The cosmic ray energy spectrum above 4 × 1017 eV as measured by the Haverah Park array. J. Phys. G: Nucl. Part. Phys., 17(5):733–758, May 1991. [13] A. A. Ivanov, S. P. Knurenko, and I. Y. Sleptsov. The energy spectrum of cosmic rays above 1015 eV derived from air Cherenkov light measurements in Yakutsk. Nucl. Phys. B Proc. Suppl., 122:226–230, July 2003. arXiv:astro-ph/0305053. [14] M. M. Winn, J. Ulrichs, L. S. Peak, C. B. A. McCusker, and L. Horton. The cosmic-ray energy spectrum above 1017 eV. J. Phys. G: Nucl. Phys., 12(7):653– 674, July 1986. [15] T. Antoni et al. [KASCADE-Collaboration]. The cosmic-ray experiment KASCADE. Nucl. Instr. Meth. A, 513(3):490–510, November 2003. [16] T. Antoni et al. [The KASCADE Collaboration]. anisotropy with KASCADE. Large-scale cosmic-ray Astrophys. J., 604(2):687–692, 1 April 2004. arXiv:astro-ph/0312375. [17] T. Antoni et al. [The KASCADE Collaboration]. Search for cosmic-ray point sources with KASCADE. Astrophys. J., 608(2):865–871, 20 June 2004. arXiv:astro-ph/0402656. [18] T. Antoni et al. KASCADE measurements of energy spectra for elemental groups of cosmic rays: Results and open problems. Astropart. Phys., 24(1– 2):1–25, September 2005. arXiv:astro-ph/0505413. [19] K.-H. Kapert et al. [KASCADE-Grande Collaboration]. Recent results from KASCADE-Grande. In F. Giovannelli and G. Mannocchi, editors, Proc. of Vulcano 2006. Italian Physical Society, 2006. arXiv:astro-ph/0608340. 153 [20] N. Chiba et al. Akeno Giant Air Shower Array (AGASA) covering 100 km2 area. Nucl. Instr. Meth. A, 311(1–2):338–349, 1 January 1992. [21] H. Ohaka et al. Further development of data acquisition system of the Akeno Giant Air Shower Array. Nucl. Instr. Meth. A, 385(2):268–276, 21 January 1997. [22] M. Takeda et al. Extension of the cosmic-ray energy spectrum beyond the predicted Greisen-Zatsepin-Kuz’min cutoff. Phys. Rev. Lett., 81(6):1163–1166, 10 August 1998. arXiv:astro-ph/9807193. [23] M. Takeda et al. Energy determination in the Akeno Giant Air Shower Array experiment. Astropart. Phys., 19(4):447–462, July 2003. arXiv:astro-ph/0209422. [24] K. Shinozaki and M. Teshima for the AGASA Collaboration. AGASA results. Nucl. Phys. B Proc. Suppl., 136:18–27, November 2004. [25] N. Hayashida et al. Possible clustering of the most energetic cosmic rays within a limited space angle observed by the Akeno Giant Air Shower Array. Phys. Rev. Lett., 77(6):1000–1004, 5 August 1996. [26] M. Takeda et al. Small-scale anisotropy of cosmic rays above 1019 eV observed with the Akeno Giant Air Shower Array. Ap. J., 522(1):225–237, 1 September 1999. arXiv:astro-ph/9902239. [27] N. Hayashida et al. Updated AGASA event list above 4 × 1019 eV. 2000. arXiv:astro-ph/0008102. [28] M. Takeda et al. In Proc. 27th Int. Cosmic Ray Conf., Hamburg, page 345, 2001. [29] M. Teshima et al. In Proc. 28th Int. Cosmic Ray Conf., Tsukuba, page 437, 2003. 154 [30] C. B. Finley and S. Westerhoff. On the evidence for clustering in the arrival directions of AGASA’s ultrahigh energy cosmic rays. Astropart. Phys., 21(4):359– 367, July 2004. arXiv:astro-ph/0309159. [31] Y. Uchihori, M. Nagano, M. Takeda, M. Teshima, J. Lloyd-Evans, and A. A. Watson. Cluster analysis of extremely high energy cosmic rays in the northern sky. Astropart. Phys., 13(2–3):151–160, May 2000. arXiv:astro-ph/9908193. [32] R. M. Baltrusaitis et al. The Utah Fly’s Eye detector. Nucl. Instr. Meth. A, 240(2):410–428, 15 October 1985. [33] R. M. Baltrusaitis, G. L. Cassiday, R. Cooper, B. R. Dawson, J. W. Elbert, B. E. Fick, D. F. Liebing, E. C. Loh, P. Sokolsky, and D. Steck. The Fly’s Eye detector: Present and future. Nucl. Instr. Meth. A, 264(1):87–92, 1 February 1988. [34] S. C. Corbató, H. Y. Dai, J. W. Elbert, D. B. Kieda, E. C. Loh, P. V. Sokolsky, P. Sommers, and J. K. K. Tang. HiRes, a High Resolution Fly’s Eye detector. Nucl. Phys. B Proc. Suppl., 28(2):36–39, November 1992. [35] T. Abu-Zayyad et al. The prototype high-resolution Fly’s Eye cosmic ray detector. Nucl. Instr. Meth. A, 450(2–3):253–269, 11 August 2000. [36] P. Sokolsky and G. B. Thompson. Highest energy cosmic rays and results from the HiRes experiment. J. Phys. G: Nucl. Part. Phys., 34(11):R401–R429, November 2007. arXiv:astro-ph/0706.1248. [37] R. U. Abbasi et al. [High Resolution Fly’s Eye Collaboration]. First observation of the Greisen-Zatsepin-Kuzmin suppression. Phys. Rev. Lett., 100:101101, 10 March 2008. [38] M. Teshima et al. Overview of AGASA results. In Proc. of CRIS 2006, Catania, Sicily, 2006. 155 [39] K. Shinozaki. In Proceedings of Quarks, 14th Int. Seminar on High Energy Physics, St. Petersburg, 2006. [40] The Auger Collaboration. The Pierre Auger Observatory Design Report, 14 March 1997. http://www.auger.org/admin/DesignReport. [41] J. Abraham et al. [Auger Collaboration]. Properties and performance of the prototype instrument for the Pierre Auger Observatory. Nucl. Instr. Meth. A, 523(1–2):50–95, 1 May 2004. [42] J. Abraham et al. [The Pierre Auger Collaboration]. Observation of the suppression of the flux of cosmic rays above 4 × 1019 eV. Accepted to Phys. Rev. Lett., 2008. arXiv:astro-ph/0806.4302. [43] T. Yamamoto for the Pierre Auger Collaboration. The UHECR spectrum measured at the Pierre Auger Observatory and its astrophysical implications. In Proc. 30th Int. Cosmic Ray Conf., Méridia, 2007. arXiv:astro-ph/0707.2638. [44] J. Abraham et al. [The Pierre Auger Collaboration]. Upper limit on the cosmic-ray photon flux above 1019 eV using the surface detector of the Pierre Auger Observatory. Astropart. Phys., 29(4):243–256, May 2008. arXiv:astroph/0712.1147. [45] R. U. Abbasi et al. [The HiRes Collaboration]. Study of small-scale anisotropy of ultra-high energy cosmic rays observed in stereo by the High Resolution Fly’s Eye detector. Astrophys. J., 610(2):L73 – L76, 1 August 2004. [46] S. Mollerach for the Pierre Auger Collaboration. Studies of clustering in the arrival directions of cosmic rays detected at the Pierre Auger Observatory above 10 EeV. In Proc. 30th Int. Cosmic Ray Conf., Méridia, 2007. arXiv:astroph/0706.1749. [47] J. Abraham et al. [The Pierre Auger Collaboration]. Correlation of the highestenergy cosmic rays with nearby extragalactic objects. Science, 318(5852):938– 943, 9 November 2007. arXiv:astro-ph/0711.2256. 156 [48] J. Abraham et al. [The Pierre Auger Collaboration]. Correlation of the highestenergy cosmic rays with the positions of nearby active galactic nuclei. Astropart. Phys., 29(3):188–204, April 2008. arXiv:astro-ph/0712.2843. [49] R. U. Abbasi et al. [The High Resolution Fly’s Eye Collaboration]. Search for correlations between HiRes stereo events and active galactic nuclei. 2008. arXiv:astro-ph/0804.0382. [50] J. W. Cronin. Cosmic rays: the most energetic particles in the universe. Rev. Mod. Phys., 71(2):S165–S172, March 1999. [51] L. Anchordoqui, T. Paul, S. Reucroft, and J. Swain. Ultrahigh energy cosmic rays: The state of the art before the Auger Observatory. Int. J. Mod. Phys., 18(13):2229–2366, 20 May 2003. [52] D. F. TOrres and L. A. Anchordoqui. Astrophysical origins of ultrahigh energy cosmic rays. Rep. Prog. Phys., 67:1663–1730, July 2004. arXiv:astroph/0402371. [53] E. Fermi. On the origin of the cosmic radiation. Phys. Rev., 75(8):1169–1174, 15 April 1949. [54] M. A. Malkov and L. O’C. Drury. Nonlinear theory of diffusive acceleration of particles by shock waves. Rep. Prog. Phys., 64(4):429–481, April 2001. [55] M. Takeda. Extremely high energy cosmic rays observed with AGASA: Energy spectrum and arrival direction distribution. PhD thesis, Tokyo Institute of Technology, March 1998. [56] S. Yoshida and M. Teshima. Energy spectrum of ultra-high energy cosmic rays with extra-galactic origin. Prog. Theor. Phys., 89(4):833–845, April 1993. [57] R. Schlickeiser. Cosmic Ray Astrophysics. Springer-Verlag, 2002. 157 [58] F. W. Stecker. Effect of photomeson production by the universal radiation field on high-energy cosmic rays. Phys. Rev. Lett., 21(14):1016–1018, 30 September 1968. [59] V. S. Berezinksy and S. I. Grigor’eva. A bump in the ultra-high energy cosmic ray spectrum. Astron. Astrophys., 199(1–2):1–12, June 1988. [60] M. Nagano and A. A. Watson. Observations and implications of the ultrahighenergy cosmic rays. Rev. Mod. Phys., 72(3):689–732, July 2000. [61] P. Bhattacharjee and G. Sigl. Origin and propagation of extremely high energy cosmic rays. Phys. Rept., 327:109–247, 2000. arXiv:astro-ph/9811011. [62] J. L. Puget, F. W. Stecker, and J. J. Bredekamp. Photonuclear interactions of ultrahigh energy cosmic rays and their astrophysical consequences. Astrophys. J., 205(1):638–654, 15 April 1976. [63] L. N. Epele and E. Roulet. Comment on “On the origin of the highest energy cosmic rays”. Phys. Rev. Lett., 81(15):3295–3295, 12 October 1998. arXiv:astroph/9806251. [64] A. Haungs, H. Rebel, and M. Roth. Energy spectrum and mass composition of high-energy cosmic rays. Rep. Prog. Phys., 66(7):1145–1206, July 2003. [65] A. R. Bell. The acceleration of cosmic rays in shock fronts. I. Mon. Not. R. Astron. Soc., 182:147–156, January 1978. [66] A. M. Bykov and G. D. Fleishman. On non-thermal particle generation in superbubbles. Mon. Not. R. Astron. Soc., 255:269–275, 15 March 1992. [67] A. M. Hillas. The origin of ultra-high energy cosmic rays. Ann. Rev. Astron. Astrophys., 22:425–444, September 1984. [68] M. S. Longair. High Energy Astrophysics. Cambridge University Press, 1992. 158 [69] J. Alvarez-Muñiz and T. Stanev. The large scale structure of the galactic magnetic field and high energy cosmic ray anisotropy. J. Phys.: Conf. Ser., 47:126– 131, 2005. arXiv:astro-ph/0507273. [70] M. C. Begelman, R. D. Blandford, and M. J. Rees. Theory of extragalactic radio sources. Rev. Mod. Phys., 56(2):255–351, April 1984. [71] P. L. Biermann and P. A. Strittmatter. Synchrotron emission from shock waves in active galactic nuclei. Astrophys. J., 322(2):643–649, 15 November 1987. [72] L.Ø’C. Drury. An introduction to the theory of diffusive shock acceleration of energetic particles in tenuous plasmas. Rep. Prog. Phys., 46(8):973–1027, August 1983. [73] J. T. Stocke, G. H. Rieke, and M. J. Lebofsky. New observational constraints on the M87 jet. Nature, 294:319–322, 26 November 1981. [74] J. R. P. Angel and H. S. Stockman. Optical and infrared polarization of extragalactic objects. Ann. Rev. Astron. Astrophys., 18:321–361, September 1980. [75] P. G. Tinyakov and I. I. Tkachev. BL Lacertae are probable sources of the observed ultra-high energy cosmic rays. JETP Lett., 74(9):445–448, 10 November 2001. arXiv:astro-ph/0102476. [76] P. G. Tinyakov and I. I. Tkachev. Tracing protons through the Galactic magnetic field: a clue for charge composition of ultra-high-energy cosmic rays. Astropart. Phys., 18(2):165–172, October 2002. arXiv:astro-ph/0111305. [77] P. G. Tinyakov and I. I. Tkachev. Cuts and penalties: Comment on “Clustering of ultrahigh energy cosmic rays and their sources”. Phys. Rev. D, 69(12):128301, 30 June 2004. arXiv:astro-ph/0301336. [78] D. S. Gorbunov et al. Testing the correlations between ultra-high-energy cosmic rays and BL Lac type objects with HiRes stereoscopic data. JETP Lett., 80(3):145–148, 10 August 2004. arXiv:astro-ph/0406654. 159 [79] D. F. Torres, S. Reucroft, O. Reimer, and L. A. Anchordoqui. On the crosscorrelation between the arrival direction of ultra-high-energy cosmic rays, BL Lacertae objects, and EGRET detections: a new way to identify EGRET sources? Astrophys. J., 595:L13–L16, 20 September 2003. arXiv:astro- ph/0307079. [80] P. Blasi, R. I. Epstein, and A. V. Olinto. Ultra-high energy cosmic rays from young neutron star winds. Astrophys. J., 533(2):L123–L126, 20 April 2000. arXiv:astro-ph/9912240. [81] J. Arons. Magnetars in the metagalaxy: An origin for ultra high energy cosmic rays in the nearby universe. Astrophys. J., 589(2):871–892, 1 June 2003. arXiv:astro-ph/0208444. [82] E. Boldt and P. Ghosh. Cosmic rays from remnants of quasars? Mon. Not. R. Astron. Soc., 307(3):491–494, August 1999. arXiv:astro-ph/9902342. [83] R. D. Blandford and R. L. Znajek. Electromagnetic extraction of energy from Kerr black holes. Mon. Not. R. Astron. Soc., 179:433–456, May 1977. [84] R. L. Znajek. The electric and magnetic conductivity of a Kerr hole. Mon. Not. R. Astron. Soc., 185:833–840, 1978. [85] E. Boldt and M. Loewenstein. Cosmic ray generation by quasar remnants: constraints and implications. Mon. Not. R. Astron. Soc., 316(3):L29–L33, August 2000. arXiv:astro-ph/0006221. [86] A. Levinson. Particle acceleration and curvature TeV emission by rotating, supermassive black holes. Phys. Rev. Lett., 85(5):912–915, 31 July 2000. arXiv:astro-ph/0002020. [87] L. A. Anchordoqui, G. E. Romero, and J. A. Combi. Heavy nuclei at the end of the cosmic-ray spectrum? Phys. Rev. D, 60(10):103001, 25 October 1999. arXiv:astro-ph/9903145. 160 [88] P. O. Lagage and C. J. Cesarsky. The maximum energy of cosmic rays accelerated by supernova shocks. Astron. Astrophys., 125(2):249–257, September 1983. [89] J. R. Jokipii and G. R. Morfill. On the origin of high-energy cosmic rays. Astrophys. J., 290:L1–L4, 1 March 1985. [90] J. R. Jokipii. Rate of energy gain and maximum energy in diffusive shock acceleration. Astrophys. J., 313:842–846, 15 February 1987. [91] P. L. Biermann. Cosmic rays I. The cosmic ray spectrum between 104 GeV and 3109 GeV. Astron. Astrophys., 271:649–661, 1993. arXiv:astro-ph/9301008. [92] T. M. Heckman, L. Armus, and G. K. Miley. On the nature and implications of starburst-driven galactic superwinds. Astrophys. J. Suppl., 74:833–868, 1 December 1990. [93] T. A. D. Paglione, A. P. Marscher, and J. M. Jackson. Diffuse gamma-ray emission from the starburst galaxy NGC 253. Astrophys. J., 460:295–302, 20 March 1996. [94] L. A. Anchordoqui, H. Goldberg, and D. F. Torres. Anisotropy at the end of the cosmic ray spectrum? Phys. Rev. D, 67(12):123006, 26 June 2003. arXiv:astroph/0209546. [95] D. B. Sanders and I. F. Mirabel. Luminous infrared galaxies. Ann. Rev. Astron. Astrophys., 34:749–792, September 1996. [96] A. Śmialkowski, M. Giller, and W. Michalak. Luminous infrared galaxies as possible sources of UHE cosmic rays. J. Phys. G: Nucl. Part. Phys., 28(6):1359– 1374, June 2002. arXiv:astro-ph/0203337. [97] G. J. Fishman and C. A. Meegan. Gamma-ray bursts. Ann. Rev. Astron. Astrophys., 33:415–58, September 1995. 161 [98] E. Waxman. Gamma-ray bursts: The underlying model. In K. W. Weiler, editor, Supernovae and Gamma Ray Bursters, volume 598 of Lecture Notes in Physics, page 393. Springer-Verlag, 2003. arXiv:astro-ph/0303517. [99] C. A. Meegan, G. J. Fishman, R. B. Wilson, W. S. Paciesas, G. N. Pendleton, J. M. Horack, M. N. Brock, and C. K. Kouveliotou. Spatial distribution of γ-ray bursts observed by BATSE. Nature, 355:143–145, 9 January 1992. [100] M. R. Metzger, S. G. Djorgovski, S. R. Kulkarni, C. C. Steidel, K. L. Adelberger, D. A. Frail, E. Costa, and F. Frontera. Spectral constraints on the redshift of the optical counterpart to the γ-ray burst of 8 May 1997. Nature, 387(6636):878– 880, 26 June 1997. [101] F. W. Stecker. Can gamma ray bursts produce the observed cosmic rays above 1020 eV? Astropart. Phys., 14(3):207–210, November 2000. arXiv:astroph/9911269. [102] S. T. Scully and F. W. Stecker. On the spectrum of ultrahigh energy cosmic rays and the γ-ray burst origin hypothesis. Astropart. Phys., 16(3):271–276, January 2002. arXiv:astro-ph/0006112. [103] A. G. Lemaı̂tre. Contributions to a British Association discussion on the evolution of the universe. Nature, 128(3234):704–706, 1931. [104] S. J. Sciutto. Air shower simulations with the AIRES system. In Proc. 26th Int. Cosmic Ray Conf., Salt Lake City, 1999. arXiv:astro-ph/9905185. [105] S. J. Sciutto. AIRES: A system for air shower simulations (Version 2.2.0), 1999. arXiv:astro-ph/9911331. [106] S. J. Sciutto. The AIRES system for air shower simulations. an update. In Proc. 27th Int. Cosmic Ray Conf., Hamburg, 2001. arXiv:astro-ph/0106044. [107] Corsika webpage: http://www-ik.fzk.de/corsika/. [108] Aires webpage: http://www.fisica.unlp.edu.ar/auger/aires/. 162 [109] M. T. Dova, A. G. Mariazzi, and A. A. Watson. The mass of cosmic rays above 1017 eV. In Proc. 29th Int. Cosmic Ray Conf., Pune, volume 7, pages 275–279, 2005. arXiv:astro-ph/0512408. [110] A. A. Watson. The mass composition of cosmic rays above 1017 eV. Nucl. Phys. B Proc. Suppl., 136:290–300, November 2004. arXiv:astro-ph/0312475. [111] M. A. K. Glasmacher et al. The cosmic ray composition between 1014 and 1016 eV. Astropart. Phys., 12(1–2):1–17, October 1999. [112] R. U. Abbasi et al. [The High Resolution Fly’s Eye Collaboration]. A study of the composition of ultra-high-energy cosmic rays using the High-Resolution Fly’s Eye. Astrophys. J., 622(2):910–926, 1 April 2005. arXiv:astro-ph/0407622. [113] T. Antoni et al. [KASCADE Collaboration]. Electron, muon, and hadron lateral distributions measured in air showers by the KASCADE experiment. Astropart. Phys., 14(4):245–263, January 2001. arXiv:astro-ph/0004233. [114] T. Antoni et al. [KASCADE Collaboration]. A non-parametric approach to infer the energy spectrum and the mass composition of cosmic rays. Astropart. Phys., 16(3):245–263, January 2002. arXiv:astro-ph/0102443. [115] K. Greisen. In J. G. Wilson, editor, Progress in Cosmic Ray Physics III. Interscience, 1956. [116] K. Kamata and J. Nishimura. The lateral and the angular structure functions of electron showers. Prog. Theor. Phys. Suppl., 6:93–155, 1 January 1958. [117] C. J. Crannell, H. R. Allan, and M. P. Sun. Widths of electron-photon cascades. Phys. Rev. D, 18(1):216–219, July 1978. [118] J. Linsley, L. Scarsi, and B. Rossi. Energy spectrum and structure of large showers. J. Phys. Soc. Japan, 17(Suppl. A-III):91–102, 1 January 1962. 163 [119] S. Yoshida et al. Lateral distribution of charged particles in giant air showers above 1 EeV observed by AGASA. J. Phys. G: Nucl. Part. Phys., 20(4):651– 664, April 1994. [120] N. Sakaki et al. In Proc. 26th Int. Cosmic Ray Conf., Salt Lake City, volume 1, page 353, 1999. [121] H. Y. Dai et al. On the energy estimation of ultra-high-energy cosmic rays observed with the surface detector array. J. Phys. G: Nucl. Phys., 14(6):793– 806, June 1988. [122] Root webpage: http://root.cern.ch/. [123] B. E. Carlson et al. Search for correlated high energy cosmic ray events with CHICOS. J. Phys. G: Nucl. Part. Phys., 31(5):409–416, May 2005. arXiv:astroph/0411212. [124] B. Carlson. libCTShower/Root, and CGI documentation. SURF Project Report, August 2003. [125] W. H. Press, S. A. Teukolsky, W. T. Vetterling, and B. P. Flannery. Numerical Recipes in C. Cambridge University Press, 1999. [126] W. R. Leo. Techniques for Nuclear and Particle Physics Experiments. SpringerVerlag, 1994. [127] Geant4. Physics Reference Manual. http://geant4.cern.ch/G4UsersDocu- ments/UsersGuides/PhysicsReferenceManual. [128] A. A. Watson. The highest energy cosmic rays. Nucl. Phys. B Proc. Suppl., 22(2):116–137, July 1991. [129] M. Ave, J. Knapp, J. Lloyd-Evans, M. Marchesini, and A. A. Watson. The energy spectrum of cosmic rays in the range 3 × 1017 –4 × 1018 eV as measured with the Haverah Park array. Astropart. Phys., 19(1):47–60, April 2003. arXiv:astro-ph/0112253. 164 [130] N. N. Efimov et al. In M. Nagano and F. Takahara, editors, Proc. Int. Symp. on Astrophysical Aspects of the Most Energetic Cosmic Rays and Particle Physics, page 20, Singapore, 1991. World Scientific. [131] M. I. Pravdin et al. In Proc. 26th Int. Cosmic Ray Conf., Salt Lake City, volume 3, page 292, 1999. [132] M. I. Pravdin et al. In Proc. 28th Int. Cosmic Ray Conf., Tuskuba, page 389, 2003. [133] M. Nagano et al. Energy spectrum of primary cosmic rays between 1014.5 and 1018 eV. J. Phys. G: Nucl. Phys., 10(9):1295–1310, September 1984. [134] M. Nagano et al. Energy spectrum of primary cosmic rays above 1017 eV determined from extensive air shower experiments at Akeno. J. Phys. G: Nucl. Part. Phys, 18(2):423–442, February 1992. [135] D. J. Bird et al. Evidence for correlated changes in the spectrum and composition of cosmic rays at extremely high energies. Phys. Rev. Lett., 71(21):3401– 3404, 22 November 1993. [136] D. J. Bird et al. The cosmic ray energy spectrum observed by the Fly’s Eye. Astrophys. J., 424(1):491–502, 20 March 1994. [137] T. Abu-Zayyad et al. Measurement of the cosmic-ray energy spectrum and composition from 1017 to 1018 eV using a hybrid technique. Astrophys. J., 557(2):686–699, 20 August 2001. arXiv:astro-ph/0010652. [138] R. U. Abbasi et al. [The High Resolution Fly’s Eye Collaboration]. Monocular measurement of the spectrum of UHE cosmic rays by the FADC detector of the HiRes experiment. Astropart. Phys., 23(2):157–286, March 2005. arXiv:astroph/0208301. [139] V. P. Egorova et al. Spectrum of cosmic rays with energy above 1017 eV. Int. J. Mod. Phys., 20(29):6878–6880, 20 November 2005. arXiv:astro-ph/0411484. 165 [140] M.-P. Véron-Cetty and P. Véron. A catalogue of quasars and active nuclei: 12th edition. Astron. Astrophys., 455(2):773–777, August 2006. Catalogue is available online via the VizieR catalogue access tool, CDS, Strasbourg, France, at http://vizier.u-strasbg.fr/viz-bin/VizieR.