Detection of horizontal air showers and neutrino induced showers with the Pierre Auger observatory

UNIVERSIDAD DE SANTIAGO DE COMPOSTELA Departamento de F´ısica de Part´ıculas DETECTION OF HORIZONTAL AIR SHOWERS AND NEUTRINO INDUCED SHOWERS WITH THE PIERRE AUGER OBSERVATORY In´es Vali˜no Rielo Santiago de Compostela, Diciembre 2007. UNIVERSIDAD DE SANTIAGO DE COMPOSTELA Departamento de F´ısica de Part´ıculas DETECTION OF HORIZONTAL AIR SHOWERS AND NEUTRINO INDUCED SHOWERS WITH THE PIERRE AUGER OBSERVATORY Memoria presentada para optar al Grado de Doctor en F´ısica por: In´es Vali˜no Rielo Santiago de Compostela, Diciembre 2007 Fdo: In´es Vali˜no Rielo Jaime ´ Alvarez Mu˜niz, Investigador Ram´on y Cajal de la Universidad de Santiago de Compostela CERTIFICA: que la memoria titulada “Detection of Horizontal Air Showers and Neutrino induced Showers with the Pierre Auger Observatory” ha sido realizada, bajo su direcci´on, por In´es Vali˜no Rielo en el Departamento de F´ısica de Part´ıculas de la Universidad de Santiago de Compostela, y constituye el trabajo de Tesis que presenta para optar al grado de Doctor en F´ısica. Santiago de Compostela, Octubre de 2007 Fdo: Jaime ´ Alvarez Mu˜niz Agradecimientos Mi especial agradecimiento a mi director de tesis Jaime, por su inestimable ayuda, apoyo y amistad. A todos mis compa˜neros del Grupo de Astropart´ıculas: ´ Angeles, Dave, Enrique, Gonzalo P., Gonzalo, Iago, Javier, Lorenzo, Patricia, Pedro, Quique, Ricardo, Ver´onica y V´ıctor, por su ayuda durante estos a˜nos de trabajo. Quiero tambi´en expresar mi agradecimiento al profesor Alan Watson, por su hospitalidad durante mis estancias en Leeds, por su ayuda y sus inestimables ense˜nanzas durante todos estos a˜nos. Quiero agradecerles a todos los miembros del proyecto “Intersecci´on de Aspectos No-perturbativos y Perturbativos en F´ısica de Solitones y Materia de Alta Densidad y Fenomenolog´ıa de Astropart´ıculas” dirigido por el catedr´atico Carlos Pajares por la ayuda econ´omica prestada en la elaboraci´on de este trabajo. Finalmente quiero dar las gracias a mis compa˜neros de la facultad de F´ısica, por todos los buenos momentos que hemos pasado juntos: Diego, Elis, Javier, Jose, N´estor y Teresa. A mis amigas In´es, Inma, Isa, Luz, Patricia, Sonia y Sonia V´azquez por todo el ´animo y cari˜no que siempre me han mostrado. A mi familia y a Juan, por todo su cari˜no y apoyo. A mis padres y a Juan. Contents 1 Introduction 1 2 Ultra High Energy Cosmic Rays 3 2.1 An overview of the study of Ultra High Energy Cosmic Rays . 3 2.1.1 A brief history of Cosmic Rays . . . . . . . . . . . . . 3 2.1.2 The Cosmic Ray Spectrum and Composition . . . . . . 4 2.1.3 Propagation and interactions of UHECRs . . . . . . . . 6 2.1.4 Origin of the bulk of cosmic rays . . . . . . . . . . . . 8 2.2 Extensive Air Showers and measurement techniques . . . . . . 10 2.2.1 General features of air showers . . . . . . . . . . . . . . 11 2.2.2 Detection Techniques . . . . . . . . . . . . . . . . . . . 18 2.3 Extremely High Energy Neutrinos and their detection . . . . . 20 2.3.1 Candidate production mechanisms of EeV neutrinos . . 20 2.3.2 Neutrino detection . . . . . . . . . . . . . . . . . . . . 21 3 The Pierre Auger Observatory 27 3.1 The concept of a Hybrid Detector . . . . . . . . . . . . . . . . 27 3.2 The Surface Detector . . . . . . . . . . . . . . . . . . . . . . . 29 3.2.1 Calibration of the Surface Detector . . . . . . . . . . . 30 3.2.2 The Surface Detector Trigger System . . . . . . . . . . 33 3.3 Reconstruction of vertical showers with the Surface Detector . 36 3.4 Reconstruction of inclined showers with the Surface Detector . 37 3.5 Energy spectrum with vertical showers . . . . . . . . . . . . . 39 3.6 Energy spectrum with inclined showers . . . . . . . . . . . . . 39 4 Study of the signals in the Surface Detector stations of the Pierre Auger Observatory 43 4.1 S(1000) USC code: an alternative method to simulate the TankResponse .......................... 43 4.1.1 Description of the method . . . . . . . . . . . . . . . . 43 4.1.2 Corrections to the muonic signal . . . . . . . . . . . . . 53 ix 5.19 Illustration of the two types of deflection in the muon trajectories due to the effect of the geomagnetic field B⊥. . . . . . . 105 5.20 The components Bg ⊥and Bs ⊥of B⊥plotted as a function of the azimuth direction of the shower for different zenith angles. . . 106 5.21 Muon and electromagnetic signal maps in the shower plane for 10 EeV proton showers with a zenith angle of 86◦and different azimuth angles under the effect of the geomagnetic field . . . . 109 5.22 The relative differences between the lateral distributions of the ratio of the electromagnetic to muon signals without and with geomagnetic field effect for 10 EeV proton showers with θ= 70◦and 80◦...........................110 5.23 The relative differences between the lateral distributions of the ratio of the electromagnetic to muon signals without and with geomagnetic field effect for 10 EeV proton showers with θ= 82◦, 86◦and 88◦........................111 6.1 Schematic picture illustrating the dependence of the apparent transmission speed of the signal on the zenith angle of the shower................................121 6.2 Footprint of an aligned event where the candidate stations after the station selection are too separated. . . . . . . . . . . 123 6.3 Comparison between the signal map in the transverse plane of proton-induced showers and νµ−induced showers in CC interactions. ..............................127 6.4 Muon and electromagnetic contributions to the tank signal in VEM as a function of the distance from the shower axis in the shower plane for proton-induced showers and νµ−induced showers. ..............................128 6.5 Sketch of an inclined shower reaching the ground. . . . . . . . 131 6.6 The depths crossed by the early and late planes that hit a tank at a distance of r= 4.5 km from the core as a function of the depth crossed by the shower axis for different zenith angles . . 132 6.7 The difference in depth crossed by the early and late planes as a function of the distance from the core on the ground when the depth crossed by the shower core is ∆X∼1500 g cm−2. 133 6.8 Schematic representation of the development of a neutrino inducedshower ...........................134 6.9 Maps of ratios SEM /Stotal and SM/Stotal for a 1 EeV proton induced shower at 85◦and ∆X= 1500 g cm−2........135 6.10 Maps of ratios SEM /Stotal and SM/Stotal for 1 EeV proton induced shower at 89◦and ∆X= 1500 g cm−2.........136 xvi 6.11 Ratio of the electromagnetic signal to the total signal as a function of the distance from the core of early and late tanks for a 1 EeV proton shower with ∆X= 1500 g cm−2. . . . . . 137 6.12 Ratio of the electromagnetic signal to the total signal as a function of ∆Xin tanks located at two distances from the core for a shower at 1 EeV proton induced shower at different zenithangles. ...........................139 6.13 Ratio of the electromagnetic signal to the total signal as a function of ∆Xfor the earliest tank in 1 EeV proton induced showers at θ= 85◦and 89◦. ...................140 6.14 Simulated FADC traces of an electromagnetic and muonic showerfront. ...........................142 6.15 Average RT +FT as a function of the normalized start-time averaged over several simulated events of neutrino showers and real inclined events . . . . . . . . . . . . . . . . . . . . . . . . 145 6.16 Example of simulated events produced by neutrino showers and of the time structure of their signals. . . . . . . . . . . . . 146 6.17 Example of an inclined event produced by an ordinary nucleonic shower and of the time structure of their signals . . . . . 147 6.18 Example of an ordinary inclined event whose earliest station has a FADC trace with a double peak because of an accidental muon entering the tank in coincidence with a muon from the shower................................147 6.19 Distributions of risetime in the earliest tank in real inclined events and simulated deep events at 10 EeV, 80◦and different ∆Xintervals. ...........................150 6.20 Distributions of falltime in the earliest tank in real inclined events and in simulated events at 10 EeV, 80◦and different ∆Xintervals. ...........................151 6.21 Distributions of risetime and falltime in the earliest tank in simulated events at 10 EeV, and ∆X∈(1500,2000) gcm−2 and different zenith angles. . . . . . . . . . . . . . . . . . . . . 152 6.22 Distributions of risetime and falltime in the earliest tank in simulated events at 80◦, ∆X∈(1500,2000) gcm−2and different shower energies. . . . . . . . . . . . . . . . . . . . . . . . . 152 6.23 Distributions of risetime and falltime for simulated neutrino showers with E∈(0.1,1.) EeV, θ∈(75◦,89◦) and ∆X∈ (0,2500) g cm−2compared with reconstructed real events. . . 153 xvii 6.24 Trigger efficiency of 3-fold events or larger for down-going neutrinos as a function of the slant crossed by the shower measured from the ground for different zenith angles and different showerenergies...........................156 6.25 Trigger efficiency of 5-fold events or larger for down-going neutrinos as a function of the slant crossed by the shower measured from the ground for different zenith angles and different showerenergies...........................157 6.26 Average number of triggered tanks per event as a function of the slant injection depth measured from the ground, for neutrino induced showers at different shower energies and zenith angles................................158 6.27 Fraction of simulated triggered events (5-fold or more) selected as having θrec ≥75◦as a function of the slant crossed by the shower measured from the ground, for different zenith angles and shower energies from 1 to 10 EeV. . . . . . . . . . . . . . 161 6.28 Example of the reasons of losing efficiency of selection after the event reconstruction. . . . . . . . . . . . . . . . . . . . . . 162 6.29 Difference between the simulated zenith angle (θsim) and the reconstructed zenith angle (θrec) as a function of the slant injection depth measured from the ground, for different primary shower energies. Each panel corresponds to a different value of the primary zenith angle. . . . . . . . . . . . . . . . . . . . 163 6.30 Examples of two neutrino induced events of shower energy Esh = 1 EeV at θsim = 75◦. One even reconstruted with θrec > θsim and other with θrec < θsim. .................164 6.31 The difference between the reconstructed zenith angle (θrec) and the primary zenith angle (θsim) performed by the standard and the aligned reconstruction methods for θsim = 89◦as a function of the slant injection depth measured from the ground and for 1 EeV shower energy.. . . . . . . . . . . . . . . . . . . 165 6.32 Selection efficiencies for down-going neutrinos as a function of the slant crossed by the shower for different zenith angles and different shower energies. . . . . . . . . . . . . . . . . . . . . . 167 6.33 Fraction of simulated events after passing the reconstruction and the deep shower cuts as a a function of the slant injection depth measured from the ground, for different zenith angles and shower energies. . . . . . . . . . . . . . . . . . . . . . . . 168 6.34 Sketch of the effect of the zenith angle on the identification efficiencies for down-going neutrinos for a fixed energy and slant injection depth from the ground. . . . . . . . . . . . . . 169 xviii 6.35 Identification efficiencies for down-going neutrinos as a function of the slant injection depth measured from the ground, for different zenith angles and different shower energies. . . . . 170 6.36 Effective depth for down-going neutrino identification as a function of the shower zenith angle for different shower energies. ...............................171 7.1 The SD event with ID=1452015 that passes all the cuts for selection of neutrino candidates . . . . . . . . . . . . . . . . . 176 7.2 FADC traces of the two earliest stations of the event 1452015 . 177 7.3 The SD event labeled as event 1956182 that passes all the cuts for selection of neutrino candidates . . . . . . . . . . . . . . . 179 7.4 FADC traces of the two earliest stations of the event 1956182 . 179 7.5 FADC traces of two stations with double peak of the event 1956182 ..............................180 7.6 The exposure of the Surface Detector of the Pierre Auger Observatory for down-going neutrino showers for one year and assuming a constant geometrical area A= 3000 km2as a function of shower energy. . . . . . . . . . . . . . . . . . . . . 183 7.7 The neutrino-nucleon (and antineutrino-nucleon) cross-section in CC and NC interactions obtained using the CTEQ6 set of parton distribution functions. . . . . . . . . . . . . . . . . . . 184 7.8 Sensitivity of the Surface Detector of the Pierre Auger Observatory to an E−2diffuse neutrino flux at 90% C.L. . . . . . . 186 7.9 Upper limits at 90% C.L. for an E−2diffuse neutrino flux assuming all νflavours ......................188 7.10 The upper limit at 90% C.L. for an E−2diffuse neutrino flux integrating the event rate from Emin ν= 5 ×1017 eV up to Eν as function of Eν..........................189 xix List of Tables 4.1 Results of the fitted parameters for the fraction of total electromagnetic tracklength contained inside the tank. . . . . . . . 67 4.2 Results of the fitted parameters for the proportional constant between electromagnetic tracklength and energy . . . . . . . . 72 4.3 Muon signal in VEM in an Auger tank as obtained in Geant4 and the S(1000) USC code for different kinetic energies and angles of incidence. . . . . . . . . . . . . . . . . . . . . . . . . 80 4.4 Signal produced by an electron in an Auger tank in Geant4 and the S(1000) USC code for different kinetic energies and angles of incidence. . . . . . . . . . . . . . . . . . . . . . . . . 81 4.5 Signal produced by an positron in an Auger tank in Geant4 and the S(1000) USC code for different kinetic energies and angles of incidence. . . . . . . . . . . . . . . . . . . . . . . . . 82 4.6 Signal produced by a gamma in an Auger tank as obtained in Geant4 and the S(1000) USC code for different kinetic energies and angles of incidence. . . . . . . . . . . . . . . . . . . . . . . 82 7.1 Number of events surviving the cuts for identifying neutrino candidates.............................174 xxi Chapter 1 Introduction The Pierre Auger Observatory is currently the largest cosmic ray observatory. Its goal is to characterize the properties of Ultra High Energy Cosmic Rays (UHECR) with energies above 1018 eV in order to understand their origin, mass composition and energy spectrum. The observatory is a hybrid detector combining an array of surface particle detectors and fluorescence telescopes to measure extensive air showers initiated by UHECR at energies greater than 1018 eV. A high energy cosmic ray typically initiates an air shower soon after entering the upper part of the atmosphere, achieving shower maximum at ∼800 g cm−2. The atmosphere has just the adequate matter depth (∼1000 g cm−2) so that a vertical shower results in a shower front containing a large number of electrons, positrons and photons (the electromagnetic component) at the ground. As the arrival direction of the cosmic ray increases with the zenith angle, the atmospheric slant depth crossed by the shower rises approximately in proportion with sec θ. Beyond θ= 60◦the atmospheric slant depth at Auger level increases from 1760 g cm−2to ∼31000 g cm−2at θ= 90◦. As a result, most of the electromagnetic component of showers with θ > 60◦, namely horizontal showers, is rapidly absorbed in the atmosphere, and only muons arrive at the ground accompanied by an electromagnetic halo that is mainly due to muon decay. On the contrary, high energy neutrinos might induce an horizontal shower deep into the atmosphere easily identifiable by a significant electromagnetic component at ground. This results in the idea of identifying neutrino showers in the background of horizontal showers initiated by nucleonic cosmic rays. This thesis is devoted to the study of horizontal (inclined) showers and the capability of the Surface Detector of the Pierre Auger Observatory to detect ultra high energy neutrinos using horizontal down-going showers. The present work is organized as follows: In Chapter 2 we give a brief introduction 1 to cosmic ray physics and we review some features of air showers, including showers induced by neutrinos. In Chapter 3 we give an introduction to the Pierre Auger Observatory. We describe the reconstruction techniques of the surface detector of the Pierre Auger Observatory and the methods to obtain the cosmic energy spectrum using vertical and horizontal showers. In Chapter 4 we describe a new alternative and fast method to simulate the response of the surface detector. In Chapter 5 we apply the method to the study of the signals in horizontal showers. We study the ratio of the electromagnetic to muonic contributions to the signal in an Auger tank (SEM /Sµ). This ratio is used for the energy reconstruction of inclined events. We also study the asymmetries in the ratio SEM /Sµin absence and in presence of the geomagnetic field. In Chapter 6 we develop a criterion to identify neutrino candidates in the data recorded by the surface detector. We describe the algorithms used to select and reconstruct inclined events. We compute the identification efficiencies of the surface detector for neutrino induced down-going showers assuming an ideal infinite array. In Chapter 7 we search for neutrino candidates in the data recorded by the surface detector. We study the potential of the surface detector of the Pierre Auger Observatory to constrain the diffuse flux of UHE neutrinos and we present a prospective upper limit to the diffuse flux of UHE neutrinos assuming a constant with time geometrical area of the surface detector A= 3000 km2and one year of operation. In Chapter 8 we summarize this thesis and present the main conclusions of this work. 2 Chapter 2 Ultra High Energy Cosmic Rays 2.1 An overview of the study of Ultra High Energy Cosmic Rays Cosmic Rays are relativistic particles that are continuously bombarding the Earth’s atmosphere from all directions and spanning over a wide range of energies from 109eV to beyond 1020 eV. 2.1.1 A brief history of Cosmic Rays The cosmic ray radiation was discovered almost 100 years ago. At that time, at the beginning of the 20th century, several scientists were very interested in the ionization phenomena and in understanding why a heavily shielded ion chamber still recorded ionization. It was assumed that this was some ionizing radiation associated with the earth’s radioactivity, so the detected radiation should be reduced at increasing heights above the ground. However, when Victor Hess and collaborators, in 1912, took ionization chambers in a balloon flight to an altitude about 5 kilometers, it was observed that the amount of radiation increased as the balloon climbed, discovering evidence of a very penetrating radiation coming from outside our atmosphere. This radiation was named “Cosmic Rays” by R.A. Millikan in 1925, and in those days the cosmic rays were supposed to be gamma rays. However, during the 1930s it was found that cosmic rays must be electrically charged particles because of the East-West asymmetry observed in their arrival directions, which is due to the effect of the earth’s magnetic field. During the years before man-made particle accelerators, cosmic rays served 3 Hillas-plot (100 EeV) (1 ZeV) Γ Neutron star max E ~ ZBL max E ~ ZBL White dwarf Protons (candidate sites for E=100 EeV and E=1 ZeV) GRB Galactic disk halo galaxies Colliding jets nuclei lobes hot-spots SNR Active galaxies Clusters (Fermi) (Ultra-relativistic shocks-GRB) 1 au 1 pc 1 kpc 1 Mpc -9 -3 3 9 15 3 6 9 12 15 18 21 log(Magnetic field, gauss) log(size, km) Fe (100 EeV) Protons Figure 2.4: The Hillas plot shows the size and magnetic field strength of astrophysical objects that are candidate sites of cosmic ray acceleration at 1020 eV. Objects below the diagonal line can not be sources of ultra high energy cosmic rays. these particles must be sufficiently massive mX>> 1020 eV; (c) the number density and rate of decay of Xparticles must be large enough to produce a detectable flux of UHECRs. There are basically two ways of generating X particles that decay at the present time: (1) producing them in the decays of topological defects; (2) making them quasi-stable in the early universe. 2.2 Extensive Air Showers and measurement techniques Direct observation of cosmic rays is only possible from space by flying detectors with balloons or spacecrafts. Such detectors are very limited in size and because of the steeply falling energy spectrum, direct observations run out 10 of statistics typically around 1014 eV. Above 1014 eV the flux of cosmic rays decreases so much that cosmic rays must be detected indirectly, by observing the shower of secondary particles created in the inelastic collision of the primary cosmic ray with the atmosphere and subsequent interactions. In the collision of a single high energy particle with an atmospheric nucleus of an air molecule (such as nitrogen and oxygen), its energy is distributed among the secondary particles. Then, these products and the remnant cosmic ray continue to propagate and produce after several generations an extensive air shower (EAS). Extensive air showers can be electromagnetic or hadronic depending on the nature of the primary particle. 2.2.1 General features of air showers A helpful tool to visualize the main features of an extensive air shower development was given by Heitler [22] through a ’Toy Model’. He introduced it in the context of a discussion of purely electromagnetic showers, but its basic structure also applies to air showers initiated by hadrons [23]. In Heitler’s approach the particle cascade is seen as a sequence of generations via branching processes. At each generation, each particle undergoes a splitting process into two other particles after traveling a path length (λ), each of them carrying half of the progenitor energy. The splitting continues until the average particle energy is reduced to the critical energy Ec, where the number of particles is maximum (Nmax) and no more interactions take place. After this, the particles only lose energy or get absorbed. The model displays the two most important features of air showers: the depth of the shower maximum, Xmax, depends on the primary energy in a logarithmic way: Xmax =λ ln(E0/Ec)/(ln2) (2.4) and the number of particles at shower maximum, Nmax, is proportional to the primary energy: Nmax =E0/Ec(2.5) Showers induced by proton or nucleus In an hadronic shower induced by a baryon, typically more than 80% of the particles produced in the first interaction are pions (the rest of particles are kaons, other mesons, hyperons and nucleon-antinucleon pairs). If the secondary hadrons are sufficiently energetic they will themselves initiate new hadronic interactions, produce secondaries and build up a hadron cascade that forms the core of the extensive air shower. Unstable particles such as 11 pions, kaons and another particles will sometimes decay depending on their energy. The neutral pions πos (roughly a third of all the pions produced) have a mean lifetime of 10−16 s and so will nearly always decay, except at the most extreme energies (above ∼1018 eV). The most common decay mode is into 2 photons (and electron-positron pairs for a small fraction of the decays), these photons produce an electromagnetic subshower through two processes: photons undergo pair production and electrons/positrons radiate bremsstrahlung photons. The size of the shower grows until the mean energy of the electrons reaches the critical energy (∼84 MeV in air) at which the energy losses by ionization and bremsstrahlung are equal. At this shower stage, approximately 90% of the total energy is carried in the electromagnetic cascade. Below the critical energy, the ionization losses overcome bremsstrahlung, and the electromagnetic cascade will begin to die out. Electrons and positrons in electromagnetic showers suffer multiple scattering, which is mostly going to determine the main features of the transverse structure of these cascades. Charged mesons, because of a larger mean lifetime (10−8s), not only decay but also interact strongly with atmospheric nuclei. The competition between the two processes depends essentially on the balance between the interaction mean free path (dependent on the cross-section and the density of the medium traversed) and the mean decay length. Both vary substantially with energy and become equal at an energy of ∼115 GeV for charged pions and ∼850 GeV for kaons [24]. Thus, at lower energies than these the decay probability is larger than the interaction probability. Charged pions and kaons give rise to muons and muon-neutrinos in the shower mostly through the following decay modes: π±→µ±+νµ(99.9%) K±→µ±+νµ(63.5%) →π±+π0(21.2%) (2.6) Neutrinos are weakly interacting particles that escape carrying roughly ∼2% of the primary energy. Muons are nearly relatively and have a small cross-section for interactions, so they are very penetrating. This component increases its size as the shower develops to reach a plateau that slowly attenuates, because muons mainly lose energy gradually by ionization (∼2 MeV/g cm−2in air), bremsstrahlung, electromagnetic and hadronic interactions with nuclei and pair production at very small rate compared to electrons. The radiative processes are only dominant at high energy (>500 GeV). 12 Muons are affected by decay in flight when their energies have become quite low (typically below tens of GeV) through the following modes: µ−→e−+νe+νµ µ+→e++νe+νµ(2.7) Muon decay is another source of secondary neutrinos. So, an air shower induced by a baryon can be understood as a core of high energy hadrons that is continuously feeding an electromagnetic component (electrons, positrons, and photons) mainly through π0decay, and both a muonic and a neutrino component through charged pion decay. This is schematically presented in figure 2.5 Figure 2.5: Schematic representation of an hadronic extensive air shower. A simplified view of the interaction of a cosmic ray nucleus with the atmosphere is given by the superposition model [25]. A shower induced by a 13 nucleus with atomic number Adevelops like a superposition of Aindependent nucleon air showers all starting at the same point, each carrying 1/A of the primary energy. These showers from lower energy primaries do not penetrate as deeply. So, the nucleus shower reaches its maximum size higher in the atmosphere than a proton shower of the same total energy. The resulting shower has more muons than a proton shower at the same total primary energy because pions are produced higher and they are more likely to decay before interacting. For example, an iron shower has ∼1.8 times as many muons as a proton shower of the same energy and its Xmax is higher than proton showers by ∼150 g cm−2at all energies [23]. Showers induced by gamma rays A shower induced by gamma rays shows slightly different features than a shower induced by a baryon. As figure 2.6 illustrates, it is a purely electromagnetic shower where the dominant processes are pair production and bremsstrahlung in the manner described previously and its behaviour can be accurately predicted from quantum electrodynamics. At energies above 1019 eV, there are other important processes that need to be taken into account. The LPM effect (Landau-Pomeranchuk-Migdal) becomes important reducing the cross-sections for pair production and bremsstrahlung. Additionally, photon interactions with the geomagnetic field induce pair production before entering the atmosphere what effectively reduces the energy of the particles that interact in the atmosphere, which to a large degree compensates for the LPM effect in the final of shower observables. Shower longitudinal profiles The longitudinal profile of a shower is the number of charged particles as a function of the atmospheric depth. The longitudinal profile of an electromagnetic shower is quite accurately given by the Greisen formula [26]: Ne(E0, t) = 0.31 √tmax exp [t(1 −1.5 ln s)] (2.8) where tis the atmospheric slant depth measured in radiation lengths (t= X/X0), tmax = ln(E0/Ec), and sis the shower age: s≈3t t+2tmax . Many shower properties are well parametrized by the shower age. For a given initial energy, the number of shower particles increases with depth when s < 1, reaches a maximum when s= 1 and declines when s > 1. In the case of hadronic air showers, it is very hard to describe the shower development using an analytical approach. Monte Carlo simulations can be 14 Figure 2.6: Schematic representation of an electromagnetic extensive air shower. performed to model it, but the lack of empirical knowledge of the physical processes which occur at high energies, leads to discrepancies between different models. Laboratory experiments have studied particles collisions (cross sections, inelasticity and multiplicity) only at centre-of-mass energies equivalent to fixed target energies of 1015 eV (in the rest frame of one particle), so the results must be extrapolated to the energies of interest 1020 eV and assumptions must still be made. Another problem, even at lower energies, is that the interactions are primarily ’soft’ interactions, with a low transfer of transverse momentum (Pt) and studies made at accelerators deal primarily with high Ptparticles, where the collision fragments are deflected at large angles into the detectors. An experiment is currently under construction at LHC with capability to measuring very forward particles (TOTEM [27, 28]) and may provide important data to help to refine the current models. The Gaisser-Hillas functional form [29], based on Monte Carlo simulations using the scaling model for nuclear interactions, has proved to be effective in fitting the longitudinal profile of simulated air shower developments resulting from various hadronic models with variable primary masses. The Gaisser- Hillas functional form is: N(X) = Nmax X−X0 Xmax −X0Xmax−X0 λ exp Xmax −X λ(2.9) The four parameters (Nmax, Xmax, X0, λ) provide ample size and shape 15 102 103 104 105 106 107 108 109 1010 1011 0 500 1000 1500 2000 2500 3000 Ne X (g cm-2) Eo = 1019 eV Greisen Gaisser-Hillas Figure 2.7: Longitudinal profile of a purely electromagnetic air shower using the Greisen function 2.8 (continuous line) and profile of an hadronic air shower using the model of Gaisser-Hillas 2.9 (dashed line) . freedom for fitting longitudinal profiles. As illustrated Fig. 2.7, the development of a hadronic shower at high energy tends to be faster than that of an electromagnetic shower due to the high inelasticity and multiplicity of hadronic interactions that distribute the primary energy among many particles. Moreover, after the maximum the hadronic shower has a slower attenuation because the electromagnetic component is being fed continuously into the shower by the hadronic core. Lateral distribution of shower particles The extensive air shower also develops transversally mainly due to electromagnetic and muonic particles spreading away from the shower axis. Energetic secondary hadrons have transverse momenta that are typically very small compared to their longitudinal momentum. They travel close to the shower axis and essentially are confined in a cylinder around the axis because of decay. In the case of pions, the cylinder radius is less than ∼22 m 16 [30]. Electromagnetic particles spread away from the axis primarily by multiple Coulomb scattering of electrons and positrons, and also because of the spreading angles in pair production and bremsstrahlung, which become negligible at high energies. The spread due to Coulomb scattering is given in terms of the Moli`ere radius (rM), which varies inversely with the density in the medium and it is of the order 100 m at the Auger altitude. For pure electromagnetic vertical showers, Nishimura and Kamata, and later Greisen, obtained the well-known NKG formula [31, 26] which gives the charged particle density as a function of the distance from the shower axis depending on shower age, sNKG, ρe=Ne π r2 M Γ(4.5−sNKG) Γ(sNKG)Γ(4.5−2sNKG)r rMsNKG−21 + r rMsNKG−4.5 (2.10) where Neis the total number of electrons. The NKG formula may also be extended to describe the electromagnetic component of hadronic induced showers by modifying the exponents in Eq. (2.10). Fits to the lateral distribution functions (LDF) of electrons and positrons obtained from simulations as a function of depth (t) yield an age parameter given by s=3t t+2β, where the floating parameter βtakes into account the deviations from the electromagnetic shower theory in where it is simply the age sNKG. The modified NKG formula provides a good description of the electromagnetic lateral distribution at all stages of shower development for distances sufficiently far from the hadronic core. Muons are relatively unaffected by multiple Coulomb scattering, and so their lateral distribution function retains information on the primary interactions in the shower. Muons are distributed in a broader lateral region than electromagnetic particles, and their number does not decrease as rapidly as the shower grows old. The lateral spread of muons is determined by the properties of the hadronic interactions, decays, distances to production point and geomagnetic effects. There is no standard functional form for the lateral distribution of the the muonic component. One of the earliest parameterizations of the muon LDF in vertical showers was empirically derived by Greisen [26], ρµ(r) = Nµ(t)fµ(r)≈Nµ(t)r rG−0.75 1 + r rG−2.5 (2.11) 17 where fµ(r) is a structure function describing the lateral shape of the shower, and rG= 320 m is analogous to the Moli´ere radius. Later, Vernov et al proposed an analytical form of the structure function: fµ(r)≈r r0−Γ exp r r0(2.12) with Γ = 0.4 and r0= 80 m. The LDFs are used to fit experimental data. However, neither function reproduces the whole radial range (distances from the core) of an extensive air shower. For the case of inclined (horizontal) showers, the muonic lateral distribution is not azimuthally symmetric about the shower axis because of geomagnetic deviations and geometrical and attenuation effects. For very inclined showers the geomagnetic field effect in the muon LDF becomes dominant. A quantitative description of this effect can be found in [32]. 2.2.2 Detection Techniques The classical method of detection of extensive air showers is to use a number of particle detectors distributed over the ground surface to sample the flux of secondary particles at different points of the shower front. This procedure is based on developments of the technique used by P. Auger and his collaborators in their pioneering work [1] leading to the discovery of air showers. Surface arrays include arrays of muon detectors (e.g. SUGAR), scintillators (e.g. Volcano Ranch, Yakutsk and AGASA), and water Cherenkov tanks (e.g. Haverah Park and Auger). Surface arrays determine the arrival direction of the incoming cosmic ray by recording the relative time at which each detector triggers. The directional precision is limited by the accuracy of the timing measurement, by the sampling area of the detector and by intrinsic fluctuations. The signals collected in the detectors (lateral distribution) can be used to estimate the energy from comparisons with detailed Monte Carlo simulations. Simulations predict the relation between energy and particle density. A surface array has sensitivity to the primary mass through direct or indirect measurement of the muon and electromagnetic content of the shower and/or indirect measurement of Xmax (depth of the electromagnetic shower maximum) [33]. Muon counters placed underground can be used to measure directly the muon component. Scintillator arrays and arrays of water Cherenkov tanks differ in their methods for studying the primary mass distribution. The first are essentially 18 sensitive to the electrons and positrons of the electromagnetic shower which dominate the charged particles. It can be used to estimate Xmax by measuring the shape of the LDF. In general a deeper Xmax will result in a steeper LDF. This gives the opportunity of measuring the fluctuations in Xmax by measuring the fluctuations in the observed LDF. Its sensitivity is limited by statistical fluctuations of the signal, due to the finite number of incident particles on a limited detector area. An array of water Cherenkov tanks is roughly equally sensitive to both muons and electromagnetic particles. The most promising mass indicator is the time structure of the signal. Since the muons suffer less Coulomb scattering, they tend to arrive earlier than the electromagnetic component at large core distance. Heavy nucleus showers have a larger muon component, and a fast shower development that leads to less electromagnetic tails [34]. The predicted muon component is very dependent on the multiplicity of hadronic interactions, so this leads to strongly discrepancies between different models ( ∼40% in number of muons for the same primary). The second class of air shower detectors are those that record radiation from the shower front as it traverses the atmosphere. Those include fluorescence detectors (e.g. Fly’s Eye, HiRes, Auger), air Cherenkov detectors (e.g. HEGRA [35]) and advanced radio frequency antenna arrays (e.g. the LOPES array). The fluorescence detectors record the fluorescence light (λ∼300-400 nm) emitted by deexcitation of nitrogen molecules previously excited by the electromagnetic particles traversing the atmosphere. The shower development appears as a rapidly moving spot of light describing a great circle path across the night sky. The fluorescence light is emitted isotropically with an intensity that is proportional to the number of charged particles in the shower. The efficiency of production is very low (about 4 photons per meter of track of ionizing particle), hence only high energy cosmic rays (>1017 eV) can be observed from large distances. Furthermore, observations can only be done in clear moonless nights, resulting in an average 10% duty cycle. A fluorescence detector consists of a light collector system (mirrors) used to concentrate the fluorescence and several photomultipliers (PMT) that detect the light focused by the mirrors. The timing information and amplitude from the signals together with the pointing direction of the PMTs are used to reconstruct the arrival direction and the longitudinal shower profile. The integral of the longitudinal profile is a direct measurement of the energy deposited by the electromagnetic component of the shower in the atmosphere. The attenuation of the light beam in the atmosphere must be taken into account to determine the energy resolution. The beam attenuation may be due to a combination of absorption and scattering, such as Rayleigh and 19 Chapter 3 The Pierre Auger Observatory Understanding the origin, mass composition and spectrum of the most energetic cosmic rays is one of the foremost issues in Astroparticle physics today. The cosmic ray spectrum with energies exceeding 4 ×1019 eV (above the so-called Greisen-Zatsepin-Kuzmin cutoff) is not very well known, due to the poor statistics and the large systematic errors of the few events detected with those energies. The main experimental difficulty to measure the properties of the cosmic rays with these energies is the extremely low flux of cosmic rays at these energies (of the order of 1 particle km−2sr−1yr−1for energies around 1019 eV). Only detectors that cover vast areas (thousands of kilometers) could collect a significant number of events. The Pierre Auger Observatory was conceived to detect thousands of events in the energy region from 1019 eV to 1021 eV, reconstruct their energy spectrum with unprecedented precision, measure their arrival direction distribution and study the mass composition of the incident cosmic rays over the whole sky. To achieve this coverage, it was decided to build a hybrid twosite observatory, one in the Northern and one in the Southern Hemispheres. The chosen locations are Malarg¨ue in Argentina and Colorado in the USA. The Southern Observatory is currently under construction and is located at the “Pampa Amarilla” at a mean altitude of 879 g cm−2(∼1400 m), near Malarg¨ue in Mendoza Province, Argentina. The site is relatively flat and near the base of the Andes mountains. The weather is classified as “arid” with clear skies and soft temperatures. 3.1 The concept of a Hybrid Detector The Auger Observatory is a hybrid detector, designed to be fully efficient for showers with energies above 3 EeV combining the strengths of two detection 27 techniques: an array of surface detectors and 4 fluorescence telescopes. The hybrid detector has important advantages over either surface detectors or fluorescence detectors operating alone. Observing showers simultaneously with the two different detectors allows to identify the sources of systematic uncertainty in each technique, and to measure independently the properties of the showers. The main benefits of the Surface Detector are a 100% duty cycle, a well defined aperture independent of the energy above 3 ×1018 eV and a high sensitivity to showers arriving at large zenith angles. On the other hand, the Fluorescence Detector provides a direct measurement of the longitudinal profile of the extensive air showers and a calorimetric energy measurement (the small unseen fraction of the total energy carried by muons and neutrinos introduces a small systematic uncertainty (<4%) due to lack of knowledge of the composition of the primary particle and the hadronic interaction model). In this chapter we will only describe the surface detector since it is the most relevant for the work done in this thesis. More information on the FD detector can be found in [60]. The design of the surface array of the Southern Observatory consists of 1600 water Cherenkov detector stations on a hexagonal grid of 1.5 km spacing spreading over an array of 3000 km2, overlooked by four fluorescence detector eyes (figure 3.1). Each eye contains 6 fluorescence telescopes allocated inside a building on the edge of the array. The different types of events that can be detected at the Pierre Auger: •SD events: events only detected by the surface array. •FD events: –Mono events: 1 FD eye –Stereo events: 2 or more FD eyes. •Hybrid events: –Simple hybrid events: 1 FD eye + 1 SD tank or a few SD tanks, but not enough to perform an independent SD reconstruction. –Golden events: 1 FD eye + nSD tanks, with nlarge enough to allow an independent SD reconstruction. –Platinum events or Stereo-hybrid events: 2 or more FD eyes + information from SD. In the following sections, we describe the Surface Detector, the reconstruction of cosmic air showers from the SD, and the most relevant results obtained so far. 28 Figure 3.1: The southern site of the Pierre Auger Observatory near Malarg¨ue, Argentina. The dots represent the current and planned positions of the 1600 Surface Detector tanks. The yellow labels correspond to the localization of the four Fluorescence Detector buildings. The lines mark the 30◦azimuthal field of view of each of the fluorescence telescopes (6 in each eye) of the Fluorescence Detector. 3.2 The Surface Detector The surface detectors used in the Pierre Auger Observatory of the Southern Hemisphere are deep water Cherenkov detectors [61] such as the one shown in Fig 3.2. Each detector unit consists of a cylindrical polyethylene tank, 3.6 m in diameter and 1.55 m in height, enclosing a liner filled with 12000 l of exceptionally pure water. The liner is a plastic cylindrical bag with a height of 1.2 m, which is black in the outside to seal out the external light while it is coated with Tyvek on the inside to diffuse and reflect Cherenkov light. Above the tank, there are three 9” photomultiplier tubes (PMTs) located within the space between the top of the liner and the top of the tank, which are in optical contact with the volume of water through three plastic windows. Each PMT provides two signals: from the last dynode and from the anode. The last dynode signal is amplified 32 times to match the dynamic range. The anode is used for high signals such as seen when the station is near the core of the shower. This six signals are digitized in time slots of 25 ns by a Flash Analog to Digital Converter (FADC) running at 40 MHz. The signals are sent to a Programmable Logic Device, which is used to implement the 29 local trigger conditions as described in the following section. Figure 3.2: Picture of a SD tank installed in the site. The electronics include a commercial GPS unit that provides the event time with ∼8 ns resolution. This was checked by studying the trigger times of two pairs of stations located a few meters from each other. The wireless LAN communication between the tanks and the Central Data Acquisition System (CDAS) is made by conventional radio systems and each tank has its own antenna. Each tank is a standalone system. There are two solar panels and two 12 V batteries that supply power to the electronic read-out system and to the high voltage PMTs. The total power consumption is less than 10 W. A schematic view of the main components of the SD tank is shown in Fig. 3.3. 3.2.1 Calibration of the Surface Detector The tank FADCs measure the light generated by shower particles crossing the water volume of the tanks by sampling the current generated at the PMT. However, the fact that particles crossing different detectors generate equal light does not result in an equal count in the FADCs. This is due to several factors such as differences in the PMT gains, in the water quality, the 30 Antena GPS Antena de comunicaciones Panel solar Bolsa Tanque 1.20 m 3.60 m Muon que atraviesa el tanque Caja de baterias fotomultiplicador Tubo electronica Caja de Communications antenna GPS antenna Solar panels Battery box Photomultiplier tube Electronics enclonsure Plastic tank 3.60 m 1.20 m the tankCharged particle crossing Figure 3.3: Schematic view of a Surface Detector tank, with the main components labeled. Tyvek reflectivity, etc. Therefore, the signal measured by each tank must be normalized to a common calibration unit to cancel out the detector parameter dependence. This normalization factor, called the Vertical Equivalent Muon (VEM or QV EM ), is the signal produced by a vertical muon traveling along the axis of the tank and crossing the entire depth of water. The goal of the calibration procedure is to measure and monitor with good accuracy the VEM unit for each PMT in electronics units. The calibration is carried out in three steps. Firstly the absolute calibration is determined from a sequence of measurements. Secondly, the PMTs are matched in gain. Finally, the evolution with time of the gains is monitored and inserted into the data flow (see [62] for more details). To achieve the absolute calibration of the VEM unit, in Auger we use the flux of atmospheric muons which has roughly a constant value (a rate in a tank of ∼2.5 KHz) producing a peak in a charge histogram. This histogram is understandable as the convolution of distributions of four different classes of incoming particles: (a) muons entering through the top and exiting through the bottom, (b) muons entering through the top and exiting through the side, 31 (c) muons entering and exiting through the side, and (d) small showers that produce the first peak. The peak produced by the first class of incident muons is the VEM. Unfortunately, vertical and central muons can not be selected. The way to relate the peak value to the VEM unit requires measurements of vertical muons using external triggers by means of a muon telescope. This telescope consists of two pairs of scintillator paddles centered, one on the top and the other underneath the tank. Coincidence between scintillators indicates that a vertical muon crosses the tank. In Fig. 3.4 we show an example of the charge histogram produced in a SD tank under the flux of atmospheric muons, and the histogram corresponding to the external calibration by a muon telescope. The second peak of the histogram due to atmospheric muons is found to be very stable with a peak charge equivalent approximately to 1.03 VEM for each PMT (1.09 for the sum of the 3 PMTs), allowing to convert the charge measured in any FADC channel to VEM units. All the surface detectors are calibrated remotely with an overall 5% precision with respect to their absolute VEM value. Figure 3.4: Charge histogram of signals (3 PMTs summed) in a surface detector under the flux of atmospheric muons (black). The first hump is an artifact due to the triggering (3 fold). The second hump corresponds to the signal of single muons going through the tank. The dashed histogram corresponds to events triggered by a muon telescope (see text). The muon peak occurs at 1.09 VEM. Taken from [62]. Besides the charge distribution, other two histograms are stored in the 32 tank calibration for each PMT: one with the value of the first bin before the signal to obtain the baseline and another with the maximum values of the measured FADC traces, called the peak distribution. The mean value of this last one is called VEM peak (Ipeak V EM ), and it is used as the common reference unit for trigger issues. Each station is calibrated online matching the photomultipliers gain by adjusting the voltage on each PMT to get the expected trigger rate for a given VEM threshold. The calibration is operated online every minute, and sent to CDAS every 6 minutes for monitoring, and in addition every 4 hours a charge histogram of the atmospheric muons is made to compute the position of the muon peak. 3.2.2 The Surface Detector Trigger System The SD trigger system is used to select high quality extensive air showers from the background of atmospheric muons. This is a hierarchical system with low level triggers (T1 and T2) implemented by the local tank software, the following level trigger (T3) is formed at the central system (at the observatory campus) based on the spatial and temporal correlation of the level T2 triggers. Additional high levels of trigger are implemented offline to select physical events (T4) and finally quality events which can be well reconstructed (T5). Low level triggers Currently, there are two different triggers implemented at the T1 level. The first is a simple threshold trigger that requires the 3-fold coincidence of signals exceeding 1.75 Ipeak V EM threshold. This trigger with a rate of 100 Hz is used to detect fast signals (<200 ns) corresponding to muons. This trigger is noisier and its rate is used to calibrate the gains of the PMTs (see previous section). The second is a Time over Threshold (ToT) trigger that requires that 13 bins of the FADC trace in a 120 bin window are above a threshold of 0.2 Ipeak V EM in coincidence of 2 PMTs. This trigger with a rate of 1.6 Hz is very efficient to select small spread-out signals, like those produced by distant showers of high energy or close low energy showers. All the ToT triggers are directly promoted to the second level trigger T2, whereas the T1 threshold triggers are requested to pass a higher threshold of 3.2 Ipeak V EM in coincidence of 3 PMTs to be promoted to T2 triggers. The total rate of T2 is close to 20 Hz. Whenever a station fulfills one of the two T2 trigger conditions, the trigger timestamp (start-time) and the type of the trigger are sent to CDAS. The 33 central trigger receives the T2s, which are used to check if the following level trigger (T3) is fulfilled. High level triggers The higher level triggers are intended to select real events and distinguish them from random coincidences. The third level trigger at the CDAS has been designed to have a trigger efficiency close to 1 above energies ∼1018.5 eV (T3). An offline hierarchy of two additional trigger levels are implemented to reject random coincidences (T4) and ensure a good reconstruction (T5). At this level, the trigger nomenclature is based on crowns of stations around any given tank among the triggered stations (see Fig. 3.5). We will refer to it as the “central station”. The six first neighbours around the central station form the first “crown” with hexagonal shape, named C1. The next crown is named C2. Therefore, the mth crown around the central station is named Cm. As the T3 triggers are requirements on the number of triggered stations in each crown, the number of required triggered stations (n) contained within certain number of crowns (m) is denoted as nCm. Once we have introduced the nomenclature, we can present the different trigger levels. Figure 3.5: Topology of the concentric crowns-hexagons of tanks around the central station (red) used for the T3 trigger decision. C1 in blue, C2 in green, C3 in magenta and C4 in cyan. The T3 trigger is implemented at the CDAS where a search is made for 34 the coincidence of at least 3 tanks with T2 trigger together with compactness requirements. This trigger requires at least one of the following two conditions: •TOT −2C1 and 3C2: this trigger requires a 3-fold coincidence of tanks passing the T2 ToT condition. One of the tanks must have a neighbouring tank in the first crown and another one within the 2 first crowns. This trigger is extremely relevant since 90% of the events selected with it are showers and is most efficient for vertical showers. •2C1 and 3C2 and 4C4: this trigger is more permissive, and requires a 4-fold coincidence of any T2 condition with a moderate compactness: one neighbouring tank, 2 tanks inside 2 crowns from the central one and a further tank within 4 crowns. Such a trigger is needed for the detection of horizontal showers but has a lot of noise. From the events selected by this trigger, only ∼2% are real showers. The search of all 3 stations that make up a event is completed as follows. Whenever a station gets a T2 trigger, the trigger time and the type of trigger are sent to CDAS. Placing a 50 µs window around a given T2 ( 25 µs earlier and 25 µs later), all the stations that have a T2 trigger within this time window are examined the patterns required for T3 trigger are searched for. If a pattern is found, the search stops and the T3 trigger flag is assigned to the event. For every T3, all the stations in the array that had a trigger of any level including T2 in coincidence with the central station (of the crown patterns) are recorded. A final timing criteria is imposed, the trigger times must be within (6 + 5n)µs of the central one, where nindicates the crown number. All the FADC traces of the stations that fulfill that later condition are stored in the event file. The two offline higher level triggers described were developed for vertical showers θ < 60◦in accordance with the two main characteristics expected in vertical showers: compactness of the pattern of triggered tanks and FADC traces sufficiently spread in time to satisfy the ToT condition. In principle these conditions are not suitable for the inclined showers that are regularly being selected, because their compactness requirements are too restrictive for the wide-spread topological patterns of inclined showers to some extent because the signals of horizontal showers are typically short in time. The T4 trigger, also known as the “physics” trigger, has been developed to select actual showers from the set of stored T3 data. This trigger requires that the event has at least 3 stations forming a triangle of first neighbours (a 3C1TOT event) or a compact configuration of any local trigger called 4C1. 35 Chapter 4 Study of the signals in the Surface Detector stations of the Pierre Auger Observatory 4.1 S(1000) USC code: an alternative method to simulate the Tank Response The study of the response of the Auger tank to the passage of shower particles (mainly muons, electrons and photons) is a fairly complex task that requires the use of quite sophisticated simulation techniques to model the behaviour of the detector. The response of the tank can be simulated with a number of packages, the most sophisticated one being the well-known Geant4 [72]. This package consists on tools to accurately simulate the passage of particles through matter. It provides routines to describe the behaviour of the Auger tank, and the relevant physical processes suffered by e±,µ±and γinside it. As these detailed simulations typically require a large CPU time, we have developed an alternative fast method (S(1000) USC) to calculate the response of the Auger tanks. An early version of this method is described in Ref. [73]. In this chapter, we describe the physical basis of the method and we compare it to the output of Geant4. 4.1.1 Description of the method This approach is based on parameterizations of the response of the tank to the passage of shower particles. The method stems from two basic ideas: Muons produce signals approximately proportional to their track inside the tank. Electrons, positrons and photons typically induce small electromag- 43 netic showers which give a signal approximately proportional to the energy deposited by the secondary electrons and positrons, which in turn is proportional to their tracklength inside the tank. The signal computed from the tracklength constitutes a first approximation to the average signal produced by a particle entering the tank. In the following we will account for a number of physical effects in the muonic and electromagnetic components and we will determine the signal produced by a shower reaching ground with larger accuracy. Besides the primary observables of the shower such as zenith and azimuth angles (θ,φ), the method uses as input from the shower simulator the following information about the particles reaching ground: •Type: muons, electrons, positrons and photons. •Statistical weight. •Kinetic energy (in GeV). •Distance from the shower core. •Arrival direction of the particle : θp,φp. This information can be provided by Monte Carlo codes that perform the simulation of extensive showers in the atmosphere such as AIRES [74] and CORSIKA [75]. Although the default reference plane to record the particle information is the ground plane, it is useful to work in the plane transverse to the shower axis. In this case, the particle positions from the ground are projected onto the shower plane by means of a simple rectangular projection. The results in this study are based on a library of proton showers simulated with AIRES 2.6.0 with a thinning level of 10−6. Showers were generated with an energy E = 10 EeV, and with θranging from 0◦to 88◦for the hadronic model QGSJET01. A total of 100 showers were simulated for each zenith angle. The simulations were performed in the conditions of the southern site of the Pierre Auger Observatory. In this study, the geomagnetic field effect is neglected but will be accounted for later in the following chapter. We describe step by step the procedure to compute the signal in the tank from the number densities of particles (ρ) and energy densities of particles (ǫ) given by the simulations. We will work in the shower plane unless otherwise indicated. 44 Unthinning procedure The number of particles that are produced in an air shower at the energies relevant for Auger can be very large and the computing time needed to follow all of them becomes excessively large. A way out is to use a statistical sampling algorithm (thinning algorithm) which allows to propagate only a small representative fraction of the total number of particles. Statistical weights (Wi) are assigned to the sampled particles in order to compensate for the rejected ones [76]. As the output of the simulations is a ground particle file with weighted entries, we need to perform a unthinning procedure that allows us to extract a set of unweighted particles entering a given tank. The standard procedure [77] consists of selecting all the particles in the simulation that fall inside a sampling region. Then, their weight needed to calculate the signal produced inside the tank, can be computed as a first approximation as follows: wi=Wi Atank Asampling (4.1) Here Atank is the total area of the tank projected onto the shower plane (Eq. 4.11) and Asampling is the area of the sampling region projected onto the shower plane. Later in this chapter, we will take into account the zenith angle of the particle entering the tank (θp) and project the areas onto the plane transverse to particle direction instead of onto the plane transverse to the shower axis. We can use different sampling regions depending on the results that we want to calculate. For instance in this chapter, we study the Lateral Distribution Functions (LDF), and we consider particles falling within concentric rings in the shower plane. Using polar coordinates (r, ξ), a ring limited by r−δr and r+δr has a sampling area: Aring = 2πr 2δr (4.2) If we wanted to calculate signals maps on the transverse plane (xvs y), we could consider square cells of area Acell =l×las sampling regions. In any case the sampling region has to be large enough so that a significant amount of particles falls inside it, but at the same time it should be small enough so that the properties (energy, etc...) of the particles are representative of their expected properties in the particular region in the ground in which the sampling area is located. The unthinning procedure may induce biases and artificial fluctuations. Also note that the sampling ratio may be abnormally large if the zenith angle used in the projection onto the shower plane is close to 90◦. To avoid 45 this problem, particles arriving at θ= 90◦are not accounted for in our calculations. Energy lost by a particle inside the wall of the Auger tank Since what we get from the shower simulation is the kinetic energy of the shower particles reaching the ground (Kpart 0), we must calculate the energy lost by them inside the walls of the tank in order to obtain their energies just before they enter inside the instrumented volume of water. To implement this effect, the detector has been modeled as a cylinder of 1.2 m height, 1.8 m radius and with mean wall thickness d∼1.27 cm. The energy lost by muons, electrons and positrons when crossing the walls of the tank is given by: ∆Kwall =αion ρwall < twall >(4.3) where ρwall = 0.94 g cm−3[61] is the density of the wall material (polyethylene C2H4), αion is the average energy loss in that medium, assumed to be ∼2.079 MeV g−1cm2for muons [83] and ∼1.655 MeV g−1cm2for electrons, and < twall >is the tracklength of the particle inside the wall of the tank averaged over the impact parameter. The mean tracklength is obtained taking into account that particles can enter through the top or the side of the tank (Fig. 4.1): < twall >=Ptop < ttop >+Pside < tside >(4.4) where < ttop >and < tside >are the mean tracklengths inside the top and side of the tank wall, respectively and Ptop (Pside) is the probability that a particle crosses through the top (side) wall of the tank, The mean tracklengths are calculated as: < ttop >=Veff top /Aeff top (θ)< tside >=Veff side /Aeff side(θ) where Veff and Aeff are the effective volume and area of the side and top walls of the tank. By effective we mean that if a particle hits the effective area it will enter the water volume and produce a signal. In fact, a particle could cross the wall without entering inside the tank and therefore, it would not contribute to the signal inside the tank. Therefore, to calculate the average track of the particles through the wall that enter inside the tank, we must apply a correction to account for this effect. This involves subtracting from the area and volume of the side and top walls, the area and volume of the sectors painted in ’cyan’ colour in Fig.4.1. 46 side h d R− top d d −d R Figure 4.1: Schematic view of the wall tank: top and side regions. The areas of the sectors in the top and side walls as viewed by the incoming particle are: Atop sector =πR2−2 (R−d)pR2−(R−d)2−2R2arcsin R−d R Aside sector = 2 (h−d)pR2−(R−d)2(4.5) The corresponding volumes are: Vtop sector =Atop sector d Vside sector =Atop sector (h−d) (4.6) Therefore, the effective area and the effective volume of the side wall projected onto the shower plane are: Aeff side(θ) = [2R(h−d)−Aside sector] sin θ(4.7) Veff side =π[R2−(R−d)2] (h−d)−Vside sector (4.8) The effective area and the effective volume of the top wall projected onto the shower plane are: 47 θsec 1 1.5 2 2.5 3 ] 2 [m θ A 7.5 8 8.5 9 9.5 10 10.5 11 Graph Figure 4.2: Area of an Auger tank projected onto the traverse plane to the θdirection. Aeff top (θ) = (πR2−Atop sector) cos θ+ 2Rd sin θ(4.9) Veff top =πR2d−Vtop sector (4.10) On the other hand, the probabilities that the particle crosses through the top and side walls are calculated using the total areas projected onto the direction of the incoming particle: Ptop =Atop/Atank Pside =Aside/Atank where the total projected area of the wall (see Fig.4.2) is the sum: Atank =Aθ=Atop +Aside =πR2cos θ+ 2Rh sin θ(4.11) Once we have calculated the mean tracklength of the particle in the tank wall, the average energy of a charged particle after entering inside the water volume is obtained as: Kpart =Kpart 0−∆Kwall (4.12) 48 Energy cuts on the input particles The Auger tank is a water Cherenkov detector, hence we only account for particles with energies inside the instrumented volume above the Cherenkov threshold in water. The minimum kinetic energy that a charged particle must have to produce Cherenkov light in water is: Kpart th =mpart 1 p1−(1/nw)2−1!∼264 keV for e± 54.6 MeV for µ±(4.13) where mpart is the mass of the particle and nw= 1.33 is the refractive index of water for optical wavelengths. For photons the minimum energy is chosen that so that the photon produces an electron-positron pair with at least one of them having energy above Ke th. The threshold energy for photons is: Eγ th =Ke+e− min +Ke th = 1.286 MeV (4.14) First estimate of S(r) As a first approximation for the average value of the signal Sat a distance rfrom the shower core, S(r), we calculate the so-called uncorrected S(r), separating the contributions from the muonic (µ) and the electromagnetic (EM) components of the shower: S(r) = Sµ(r) + SEM (r) (4.15) In this preliminary estimate we have assumed the following approximations: •Muons give signals proportional to their tracks. •The EM component induces typically small electromagnetic showers that give a signal approximately proportional to the energy deposited inside the tank. •All particles travel parallel to the shower axis at the speed of light β=v/c = 1. •As a result of β= 1, all particles have maximum and equal Cherenkov emission efficiency. 49 [deg]θ 0 10 20 30 40 50 60 70 80 90 > [m] θ µ <L 1.2 1.4 1.6 1.8 2 2.2 2.4 2.6 2.8 ltrack Figure 4.3: Average tracklength of muons as a function of the zenith angle. Under these assumptions, the average signal in VEM units induced by a muon number density ρµin an Auger tank is given by: Sµ(r) = ρµAθ < Lθ µ> LV EM µ =ρµAθ < Lθ µ> h=ρµ V h=ρµA0(4.16) where LV EM µis the average tracklength of a vertical muon (a VEM unit) which under the above approximations is equal to the height hof water, and < Lθ µ>=V/Aθis the average tracklength of a muon entering at zenith angle θaveraged over impact parameter. Here, Aθis the area of the tank projected onto the direction θ(Eq. 4.11), and V=πR2his the tank volume. In particular A0=πR2≃10 m2is the area of the tank as seen by a particle entering at θ= 0◦. In Fig. 4.3 we show < Lθ µ>as a function of the zenith angle. It has a minimum at around θ∼23◦. Therefore, the signal induced by the muonic component of the shower is simply: Sµ(r) = ρµ(r)A0=A0 Aring(r) Nµ X i=1 wi(4.17) with Aring(r) the sampling area in the shower plane. The sum runs over all muons inside the ring area Aring. 50 The electromagnetic particles induce subshowers in the tank which produce signals assumed to be proportional to the total tracklength of electrons and positrons (LEM ), which at the same time and as a first approximation is assumed to scale linearly with energy [78]. Under the assumption that all the energy of the electromagnetic shower is deposited inside the tank, the EM tracklength is proportional to the total energy of the EM particles entering the tank: LEM =kEEM =kǫEM Aθ(4.18) where ǫEM is the electromagnetic energy density and kis a proportionality constant. Therefore, the average signal induced by the EM component in VEM units in an Auger tank is given in this approximation by: SEM (r) = LEM h=Aθ k hǫEM (r) = Aθ Aring(r) k h NEM X i=1 Eiwi(4.19) The sum runs over all electromagnetic particles inside the ring area Aring. Photons are included in the sum since they induce EM subshowers inside the tank. The value of the proportionality constant was obtained by performing shower simulations in water using the ZHS Monte Carlo code [78], giving k= 5.25 m GeV−1for the proportionality constant 1. We have applied our method to the calculation of the signal at r= 1000 m from the core S(1000), and in particular we compute its dependence on sec θ. This dependence was chosen as an example to illustrate the relative importance on the muonic and electromagnetic signals of the different corrections with respect to the first estimate that we will introduce later in this chapter. As we mentioned in Section 3.5, S(1000) has been chosen as the energy estimator for the air showers detected at the surface detector of the Pierre Auger Observatory. In Fig. 4.4 we show the dependence of S(1000) on sec θfor 10 EeV proton showers. The relative contributions to the signal from the electromagnetic and muonic components of the shower are shown in the same figure. of S(1000) increases by ∼8% from 50.9 VEM for a vertical shower up to ∼55.0 VEM at sec θ≃1.06 (θ≃20◦) and then decreases monotonically. The increase from 0◦to 20◦is mainly due to the fact that the transverse area of the tank also rises in this range of θ(see Fig. 4.2). However, the increase of the area does not affect the average muon signal in Eqs. 4.16 and 4.17 although it 1It is worth remarking that the ZHS code shows a very good agreement with Geant4 [79] at the 9% level in the total tracklength of electromagnetic showers in water. 51 θsec 1 1.2 1.4 1.6 1.8 2 2.2 2.4 2.6 2.8 3 [VEM] µ S 6 10 20 (uncorrected) µ S decay)µ ( µ S 1)≠β ( µ S -rays)δ1 + ≠β ( µ S 1 + pp)≠β ( µ S 1 + DL)≠β ( µ S Figure 4.8: Corrections specific to the muonic signal Sµ(1000). (1) The upward triangles correspond to the uncorrected signal. (2) The stars correspond to Sµafter including the corrections due to muon decay inside the tank. (3) The full circles show how Sµis reduced after accounting for Cherenkov efficiencies and muon energy loss (β6= 1). All the remaining curves include this first correction. (4) The squares correspond to Sµafter including the corrections due to δ−rays. (5) The downward triangles correspond to Sµincluding Cherenkov efficiency and pair production. (6) The empty circles correspond to Sµincluding direct light hitting the PMTs in the tank. In this plot, we are assuming that particles travel parallel to shower axis. The simulation was performed for 10 EeV proton showers, with AIRES and the QGSJET hadronic model. 58 Energetic knock-on electrons (δ-rays) Muons also produce secondary electrons, called δ-ray electrons, along their paths inside the tank. They are mainly produced by the incident muon interacting primarily with a single atomic electron which is ejected from the atom with a considerable kinetic energy K≫I(Iis the mean excitation energy). The kinetic energy distribution of secondary δ-rays per unit thickness xis given by [7, 82]: d2N dKdx =C1 β2K21−β2K Kmax +K2 2E2(4.26) for I≪K≤Kmax. Here C= 2πr2 emec2NAz2Z/A is a constant that depends on the material (C= 0.08445 MeV cm−1for water), βis the particle velocity, Eis the total electron energy, and Kmax is the maximum kinetic energy carried by δ-ray and given by: Ke max =2mec2β2γ2 1 + 2γ(me/M) + (me/M)2(4.27) where Mis the muon mass. The secondary δ-rays with kinetic energies above Ke th = 0.264 MeV will produce Cherenkov light, and will contribute to the total signal. The number of δ-rays produced by a muon track that will contribute to the signal is given by the integral of Eq. 4.26 from Ke th to Ke max: dN dx =C β2 1 Ke th −1 Ke max −β2 Ke max ln Ke max Ke th +Ke max −Ke th 2E2(4.28) The effect of the δ-rays contribution is equivalent to an increase in the effective muon tracklength, Leff ∗ µ=Leff µ(1 + rδ) (4.29) Using Geant4 simulations, a parameterization of this contribution as a function of muon kinetic energy (Kµ) and the real muon tracklength (Eq. 4.24) was obtained in [68] with the result: rδ(Kµ, Lreal µ) = 0.135 (0.8qLreal µ+ 1) Kµ 0.8qLreal µ+Kµ (4.30) which is valid in the kinetic energy range [0.5,1000] GeV and in the tracklength range [0.3,3.8] m. For muons below 0.5 GeV, which stop inside the tank, we use an approximate constant correction: rδ(Kµ= 0.5 GeV) ≃0.09 59 Figure 4.9: Fraction of Cherenkov light produced by δ-rays taken from [68] which produces an increase in the effective muon tracklength of Leff ∗ µ= Leff µ+0.09 (Kµ/0.5). In Fig. 4.9 (taken from [68]), rδis plotted as a function of the muon kinetic energy for different muon tracklengths. As one can see δ-ray production increases with the muon energy. The behaviour with the muon tracklength depends on the muon energy, so for muons below 1 GeV, rδdecreases as the muon path rises due to muon energy loss. However, above this energy the larger the muon path, the bigger the increase because the muon can produce more energetic δ−ray electrons before leaving the tank. The correction to the muonic signal due to δ-rays must be taken into account in both the tracklength of the incident muon and the tracklength of a vertical muon, LV EM µ, which is needed to express both Sµand SEM in VEM units. The effective tracklength of vertical muons heff is obtained using the parameterization in Eq. 4.30 assuming a single spectrum of monochromatic muons of 1.05 GeV passing through the tank. We obtain a contribution of rδ= 0.13 for vertical muons, which corresponds to heff ∼1.35 m. It is important to note that in this correction we have included the previous correction accounting for reduced Cherenkov efficiency and stopping muons. The effect of the δ-rays correction on the dependence of Sµwith θis shown in Fig. 4.8. One can see the increase of Sµin all the sec θrange with respect to the correction accounting only for the Cherenkov efficiency. For vertical 60 showers, the increase is smaller due to the large fraction of muons below 1 GeV that have Leff ∗ µ< heff . As zenith angle increases, the growth is larger because the average muon energy in the shower increases (see bottom panel of Fig. 4.7) and also the muon tracklengths inside the tank are larger. Hard muon interactions Muons can suffer hard interactions in the tank namely, bremsstrahlung (bs), pair production (pp) and nuclear interactions via photo-nuclear processes (ni). Eq. 4.21 is the corresponding expression for the energy loss per unit thickness, where a(E) represents the ionization losses (ion) and b(E) accounts for the hard processes, b(E) = bbs(E) + bpp(E) + bni(E) (4.31) In Fig. 4.10 we show the rate of muon energy loss in water as a function of its kinetic energy [80]. The critical energy of muons in water is ∼1 TeV, at this energy the ionization losses are equal to the losses due to hard processes. Pair production becomes the most relevant mechanism of energy loss followed by bremsstrahlung and finally photo-nuclear interactions. A detailed study of the contribution of the hard processes to the total signal was performed in [68]. On one hand, it was found that the pair production contribution (rpp) increases with the muon energy and with its tracklength, a parameterization of this contribution is: rpp(Kµ, Lreal µ) = 2.1×10−4L0.88 µKµ 1 + 3.7×10−4L−0.16 µKµ (4.32) In Fig. 4.11 taken from [68], we show the fraction rpp as a function of the muon kinetic energy for different muon tracklengths. It was also found in [68] that the contribution due to bremsstrahlung and photonuclear interactions can be neglected. Therefore, the effective muon tracklength corrected by hard muon interactions is: Leff ∗ µ=Leff µ(1 + rpp) (4.33) Note that this correction also includes through Leff µthe correction due to Cherenkov efficiency and stopping muons. This correction must be also taken into account in the estimate of the VEM tracklength (heff ) although it is practically negligible (about 0.02%). The effect of the correction on the dependence of Sµwith zenith angle is shown in Fig. 4.8. We observe that the effect is completely negligible for the range of zenith angles considered because this correction starts to be 61 important at Kµ≃250 GeV and as one can see in the bottom panel of Fig. 4.7, the mean muon energy is much smaller than this value in this angular range. This correction will be only relevant if the muon energy is sufficiently high, i.e., close to the shower axis and for very large zenith angles. (GeV) µ K -2 10 -1 10 1 10 2 10 3 10 4 10 (GeV) µ K -2 10 -1 10 1 10 2 10 3 10 4 10 ) -1 g 2 ( MeV cm dX dE - 1 10 2 10 total bs+pp+ni pp bs ni ion Figure 4.10: Rate of muon energy loss in (liquid) water. The blue line shows the energy loss by ionization. The green line is the energy loss by pair production. The red line indicates the energy loss by bremsstrahlung. The pink line is the energy loss by nuclear interactions. The cyan dashed line is the energy loss by all the hard processes. The black line is the total energy loss. Muon decay inside the tank A muon decays into an electron and two neutrinos (see Eq. 2.7) with a branching ratio ≃100%. The probability that the muon decays inside the tank is: Pdecay = 1 −exp −< Lθ µ> λ!= 1 −exp −< Lθ µ> γ cτ !(4.34) where γ=Ei/mµand cτ = 658.654 m is the mean decay length of the muon. This probability is taken into account in this work as a correction to the muonic signal. If a muon decays inside the tank its tracklength is 62 Figure 4.11: Fraction of Cherenkov light produced by pair production processes taken from [68]. reduced and hence the muon signal diminishes. We follow here a very simple treatment of this reduction and write that on average the reduction of the muon signal is: Sµ(r) = Sµ(r)(1 −Pdecay) (4.35) This correction is very small. Assuming that the muon decays at rest (Eµ=mµ) and assuming the maximum possible tracklength ∼2.7 m (see Fig. 4.3), the correction (1 −Pdecay)≃(1 −4×10−3) at most. When a muon decays, the energy distribution of the resulting electron is known as the Michel spectrum, which has an endpoint at 53 MeV and an average electron energy of 37 MeV. The Michel electron may produce a signal in the tank if its energy is larger than the Cherenkov threshold, which we take into account as a correction to the muonic signal. Following this simplistic treatment we assume that a muon decays only when it stops inside the tank. From the moment the muon becomes subthreshold Eth = 160.3 MeV until it stops (Eth = 105.7 MeV), the muon crosses a distance of about l= 0.274 m inside the tank. Therefore, the muon 63 needs to travel a distance Ldecay =Leff µ+linside the tank to decay. We assume that the stopping muon decays into a Michel electron of energy < EMichel >= 37 MeV and two neutrinos after crossing a depth of Ldecay. The subshower initiated by the Michel electron might not always be fully contained inside the tank. In fact there is a maximum distance available for the shower to develop given by: dMichel =< Lθ µ>−Ldecay In this case, we correct the signal of each muon that stops inside the tank adding the signal produced by the Michel electron: Sµ=Sµ(1 −Pdecay) + SMichel EM (4.36) where the Michel electron signal is: SMichel EM =k(EMichel) heff (1 + rhδ +rhpp)EMichel fcont (4.37) with k= 5.15 m GeV−1. In the treatment of the signal produced by the Michel electron, we have accounted for all the corrections that affect the VEM unit (heff ) and the corrections that affect the electromagnetic component of the shower that will be described in the following sections, namely: Fraction of the shower contained inside the tank (fcont) and departure from linearity of the relation between energy and tracklength k(EMichel). For a fully contained subshower initiated by a Michel electron of < E >= 37 MeV the signal is SMichel EM ≃0.14 VEM Two important remarks: The effect of this correction is expected to be more relevant the lower the energy of the muons (they have a larger probability of stopping inside the tank) and the more inclined they are (they have a larger depth of water for the subshower initiated by the Michel electron to be fully contained). Also the signal produced by the Michel electron will be delayed in time by ∼(l/c +τ0) with respect to the signal produced by the muon and in fact, SMichel EM might fall outside the time window in which the signal is collected by the detector. This effect is neglected here. The effect of the correction due to µdecay on the dependence of Sµwith θis shown in Fig. 4.8. 4.1.3 Corrections to the electromagnetic signal The electromagnetic signal must be corrected accounting for the following physical effects: 64 Shower containment in the tank The uncorrected approximation to the EM signal assumes that all the energy in the electromagnetic component is fully deposited inside the tank. The fact is that the electromagnetic subshower initiated by secondary electrons, positrons and particularly photons, is not always completely contained in the tank. In fact it is possible for photons to go through the whole tank without producing no signal at all. To implement this effect, we have simulated a large number of photon, electron and positron induced showers in a large volume of water for different energies using the ZHS code to calculate the fraction of EM tracklength that develops inside a given depth. In the bottom panel of Fig. 4.12 we show examples of the resulting curves for photon showers for a wide range of energies. For a fixed depth, the fraction of energy deposited by a photon subshower inside that depth decreases with energy. For instance, on average a 10 MeV (20 MeV) photon deposits ∼90% (∼85%) of its energy in 1.2 m of water. For a fixed photon energy, the fraction increases with depth, or equivalently with the zenith angle of the shower (see Fig. 4.3). We calculate the average depth of water available for the subshower to develop as < d > (θ) = V/Aθ(this is the same equation that gives < Lθ µ>). We parameterize the fraction of tracklength contained inside the tank as a function of the average depth available < d > and of the particle energy for positrons, electrons and photons, and we obtain the following expression: fcont(Ei, θ) = [tanh(A < d >)]B+C(4.38) where the parameters A,Band Cdepend on the particle energy, and the two first depend also on the particle type. Their values can be read in Table 4.1. After correcting the total tracklength used in the first estimation for this effect, the electromagnetic signal becomes: SEM (r) = Aθ Aring(r) k h NEM X i=1 Eiwifcont i(Ei, θ) (4.39) The effect of the correction on the dependence of SEM with zenith angle is shown in Fig. 4.13. This correction produces a reduction on the electromagnetic signal with respect to the uncorrected one, and in fact it is the most important EM correction at all zenith angles. There is a decrease in the uncorrected SEM by about 16.5% at 0◦and about 18.9% at 60◦. This behaviour can be understood looking at Fig. 4.14 where we show the energy spectrum of electrons and positrons (top panel) and photons (bottom panel) at a distance of 1000 m from the shower axis for different zenith angles. For 65 0.5 1 1.5 2 2.5 3 3.5 4 4.5 5 5.5 1 10 100 1000 10000 LEM/EEM [m/GeV] EEM [MeV] e+ showers γ showers eshowers 0 20 40 60 80 100 0.01 0.1 1 10 Fraction of track inside tank [%] d Distance inside water tank [m] d=1.2 m Eγ=1 MeV 2 MeV 4 MeV 10 MeV 20 MeV 40 MeV 100 MeV 200 MeV 400 MeV 1 GeV 2 GeV 4 GeV 10 GeV Figure 4.12: Top panel : kias a function of electromagnetic energy showing the departure from linear scaling with energy. Bottom panel: Fraction of tracklength contained inside the tank in photon induced subshowers of different energies as a function of distance inside the tank. 66 E (MeV) AγBγAe−Be−Ae+Be+C 1. 3.19336 0.80024 0. 0. 2.40062 0.0028627 1. 2. 2.21491 0.805152 0. 0. 1.81648 0.00408368 1. 4. 1.60917 0.823809 231.018 15.1113 97.3032 0.977079 0. 10. 1.16089 0.826219 38.4995 1.35635 32.7209 1.15589 0. 20. 1.0136 0.881936 16.6198 1.2748 13.5751 1.05314 0. 40. 0.92851 0.913472 7.01112 1.09468 5.08748 0.817457 0. 100. 0.9109 1.11474 2.29793 0.91754 1.89732 0.795238 0. 200. 0.86014 1.37832 1.31594 0.994075 1.23131 0.962837 0. 400. 0.76394 1.62843 0.966934 1.24693 0.946752 1.23706 0. 1000. 0.65279 2.05272 0.751189 1.6555 0.74017 1.65012 0. 2000. 0.59384 2.46554 0.652856 2.0446 0.664692 2.03622 0. 4000. 0.54622 2.8765 0.600476 2.38313 0.59809 2.36559 0. 10000. 0.51189 3.59528 0.538004 2.8438 0.532596 2.80372 0. Table 4.1: Results of the fitted parameters for the fraction of total electromagnetic tracklength contained inside the tank given by Eq. 4.38. instance, for θ= 0◦roughly 50% of the photons have energies larger than 40 MeV, and by looking at the curves in the bottom panel of Fig. 4.12 this implies that they are depositing about 20% of their energy outside the tank. It is also important to keep in mind that the electromagnetic component contributes only about 20% to the total signal above 60◦, so the relative contribution of this effect is small, on the order of ∼4% reduction in the total signal at large zenith angles. Effect of δ-rays on the VEM unit As previously mentioned, the effect of the δ-rays emitted by muons must be taken into account in the effective tracklength of a vertical muon, i.e., in the normalization used to express SEM in units of VEM (see Eq. 4.19). This is not a correction due to a process that affects the electromagnetic component, but it is included in this section because it affects the EM signal by increasing the effective tracklength of a calibration muon and producing a constant reduction of the signal by a factor 1.13 (a constant reduction in SEM (1000) of about 10%) as shown in Fig. 4.13. This effect turns out to be the next in importance after the containment of the shower described above. 67 The correction accounting for direct light in the EM component is expected to be less important because at small zenith angles, at which the EM component contributes most to the total signal, it is practically non-existent, and at large zenith angles the total signal is dominated by muons. Although there are not specific studies in which this effect has been addressed, we use the same parameterization for the EM component as in the muonic case as a first order approximation. The effect of direct light in both the EM and muonic signal is the production of an extra signal in the tank. In Fig.4.8 the effect on Sµis an increase of the signal with θso that for very inclined showers the correction can be as large as 10%. The same effect is seen in Fig.4.13 for the electromagnetic signal. Particle deviations from the shower axis So far we have worked under the assumption that all the particles travel parallel to the shower axis, and therefore we have assumed that they all enter the tank with an angle θiequal to the zenith angle of the shower. The fact is that the particles deviate from the shower axis and therefore we must take into account the true zenith angle of the particle in many of the corrections described above: •The area of the detector has to be projected onto the plane perpendicular to the arrival direction of each particle. As a consequence, in Eq. 4.19 the area Aθmust be replaced by Aθi, where the index iruns over all the particles produced in the shower simulation. The sampling area Aring must be also projected onto this plane Aring(r, θi). The latter correction affects the calculation of both muon number densities and electromagnetic energy densities: ρµ(r) = Nµ X i=1 wi Aring(r, θi) ǫEM (r) = NEM X i=1 Eiwi Aring(r, θi)(4.45) •The corrections of Cherenkov efficiency and muons becoming subthreshold depend on the zenith angle of the particle because they depend on the maximum available depth of water: Lθ µ 74 •The correction that accounts for the probability of muon decaying inside the tank also depends on the path traveled by the muon inside the tank Lθ µ, and hence on the angle of the particle. •Particle deviations from the shower axis also affects to the containment of the shower induced by an electron, positron or photon inside the tank. This is due to the fraction of total tracklength contained inside the tank (Eq. 4.38) depending on the available depth of water < d > (θi) = V/Aθifor the shower to develop. •The correction due to direct light depends on the angle of incidence of the particles: fDL(θi). To quantify the effect of this correction we firstly compute the total signal taking into account all the corrections but assuming that particles are parallel to shower axis. The signal can be expressed as: Sµ(r) = Nµ X i=1 wi Aθ Aring (1 −Pdecay) (1 + fDL)1 heff (1 + rhδ +rhpp)Leff µ(1 + rδ+rpp) + Nµ X i=1 wi Aθ Aring SMichel EMi (4.46) SEM (r) = (1 + fDL)Aθ Aring 1 heff (1 + rhδ +rhpp) NEM X i=1 kiEiwifcont i(4.47) We also compute the total signal but this time accounting for the actual direction of the particles obtained in the simulation. The two signals are shown in Fig. 4.17. As one can see, accounting for the directions of the particles modifies very little the S(1000) curve especially at large zenith angles. To gain more insight on the relevance of this correction, we plot in Fig. 4.18 the distribution of incidence zenith angles of the particles at 1000 m from the shower axis for vertical (left panel) and inclined showers (right panel). The fact that this correction is very small at large zenith angles is due to the muonic component being dominant as one can see in Fig. 4.18, and to the fact that muons deviate very little from the shower axis partly because only the more energetic muons reach the ground. There is also some degree of compensation because in a non vertical shower, and in the early region of the shower, the particles will typically enter the tank with zenith angles 75 θsec 1 1.2 1.4 1.6 1.8 2 2.2 2.4 2.6 2.8 3 S(1000) [VEM] 8 10 20 S(1000) (uncorrected) )) i θ=θS(1000) (all corrections + ( )) i θ≠θS(1000) (all corrections + ( Figure 4.17: Attenuation of S(1000) with the zenith angle of the shower including all the corrections. The circles represent the signal assuming that particles travel parallel to the shower axis, θi=θ. The triangles represent the signal accounting for particles not travelling parallel to the shower axis. The squares correspond to the uncorrected signal. The simulations were done for 10 EeV protons with AIRES and the QGSJET model. 76 µ Entries 100 Mean 30.67 RMS 11.64 i θ 0 10 20 30 40 50 60 70 80 90 0 0.01 0.02 0.03 0.04 0.05 0.06 0.07 0.08 µ Entries 100 Mean 30.67 RMS 11.64 EM Entries 100 Mean 35.21 RMS 17.72 EM Entries 100 Mean 35.21 RMS 17.72 µ EM ° = 30θ µ Entries 100 Mean 68.72 RMS 2.213 i θ 0 10 20 30 40 50 60 70 80 90 0 0.05 0.1 0.15 0.2 0.25 0.3 0.35 µ Entries 100 Mean 68.72 RMS 2.213 EM Entries 100 Mean 62.41 RMS 11.46 EM Entries 100 Mean 62.41 RMS 11.46 ° = 70θ Figure 4.18: Distribution of the zenith angle of incidence of a particle in at 30◦(left panel) and 70◦(right panel) showers at 1000 m from the shower axis. The continuous lines correspond to the distribution for muons. The dashed lines correspond to the distribution for electrons, positrons and gammas.The simulations were done for 10 EeV proton with AIRES and the QGSJET model. smaller than the shower zenith angle, while in the late region the opposite behaviour occurs (see Fig. 5.7). However, in showers with small zenith angles the large spread in the arrival directions of muons and electromagnetic particles (left panel of Fig. 4.18) produces an increase of both the muonic and electromagnetic signals with respect to the case in which the direction of the particles was not accounted for. Moreover, for sec θ.1.1 (θ.30◦) the increase is such that the S(1000) curve flattens. This is mainly due to the combination of two effects. On one hand, a large fraction of electromagnetic particles in showers with θ.30◦arrive at ground with θparound 25◦, and therefore their corresponding projected tank areas reach the highest possible value. On the other hand, the EM particles in showers with θ.30◦practically have the same energy spectrum regardless of the zenith angle (see Fig. 4.14). As a consequence and since the effects of all corrections to the EM signal depend only on the energy and θiof the particles (see Eq. 4.47), one expects S(1000) to be rouhgly independent of θfor θ.30◦. 77 4.1.5 Corrected signal The final total signal at r=1000 m from the shower core, after accounting for all the effects described in the previous sections, is shown in Fig. 4.19. The most important correction to the electromagnetic component of the signal is due to the secondary shower not being completely contained inside the tank. For the muons it is the energy loss and inefficiencies in the Cherenkov yield. Most of the corrections tend lower the signal. As a result the corrected S(1000) curve is below the uncorrected one for θ < 64◦. For zenith angles θ > 64◦, the signal is larger then the uncorrected one mainly due to the effect of the direct light. θsec 1 1.2 1.4 1.6 1.8 2 2.2 2.4 2.6 2.8 3 S(1000) [VEM] 1 2 10 20 30 40 50 60 70 S (uncorrected) (uncorrected) EM S (uncorrected) µ S S (all corrections) (all corrections) EM S (all corrections) µ S Figure 4.19: Comparison between uncorrected and corrected electromagnetic, muonic and total signals vs sec θat 1000 m from the shower axis. The simulations were done for 10 EeV proton with AIRES and the QGSJET model. 78 4.2 Comparative study between S(1000) USC code and Geant4 The aim of this section is to compare our code with a well known, accurate and tested simulation of the passage of particles through matter such as Geant4. For this purpose, we compare the response of the tank to single particles instead of the response to particle showers. In this way we avoid having to employ unthinning algorithms that might be different in our S(1000) code and in the Offline framework where Geant4 is implemented which can introduce artificial differences. We have performed a comparative study of the response to vertical and inclined individual particles of the S(1000) USC code and Geant4 code for different particle kinetic energies. We used the fast version of the package Geant4 [72] implemented in the G4FastTankSimulatorPS Offline module. The fast version of Geant4 was developed by the Pierre Auger Collaboration [96] to optimize the performance of the existing Geant4 tank simulator, re-writing the code for fast and efficient tracking of the Cherenkov photons in the Auger tank. The goal was to reach the best compromise between accuracy in the detector response and CPU time. The computing speed is 5 times faster than the original Geant4. The output of Geant4 and the fast Geant4 is statistically indistinguishable. The particles were injected (using the ParticleInjectorOG Offline module) at 1.35 m of height all over the tank surface after projecting it onto the plane tranverse to the particle direction, in order to average over all the impact parameters 3. The signal in Geant4 is given in number of photoelectrons (pe), and it must be converted into VEM units. The simulation of the tank calibration needed for this purpose was performed using as input vertical and centered muons of 1.05 GeV, obtaining that 1 VEM corresponds to 89.53 ± 9.06 pe. The signal in VEM produced by muons, electrons, positrons and photons was obtained for different kinetic energies and zenith angles of the injected particles using Geant4 and the S(1000) USC code. The results of the comparison of the response of both codes to muons, electrons, positrons and gammas are shown in Tables 4.3, 4.4, 4.5 and 4.6 respectively. The relative differences between both codes taking Geant4 as reference are within the RMS of the Geant4 output. For muons, the relative differences between both codes are less than 10% at all energies. At intermediate energies, we expect part of the discrepancy to be due to the fairly simplistic treatment of the muon decay process in the 3Note that in the S(1000) approach the average over impact parameter is done implicitly when calculating the average µ,e−,e+and γtracklengths. 79 S(1000) USC code. For instance at K= 0.4 GeV and θ= 75◦, a muon stops after crossing a mean distance (Ldecay) which is just a bit smaller than the maximum physical distance inside the tank (< Lθ µ>). In this case and in the S(1000) USC code, the muon always decays inside the tank into a Michel electron with < EMichel >= 37 MeV. However in Geant4, the muon might not decay inside the tank because Geant4 takes into account the different tracklengths of the muon inside the tank depending on the sampled impact parameter, evaluates the corresponding decay probability and accounts for the energy distribution of the Michel spectrum. More generally, if the muon decays well outside the tank (Ldecay ≫< Lθ µ>) or well inside the tank (Ldecay ≪< Lθ µ>), we do not expect a significant discrepancy between codes due to the implementation of the muon decay process. Only at those energies and angles at which Ldecay is approximately equal to the available tracklength inside the tank we expect large differences due to the different treatment of the muon decay in Geant4 and S(1000) USC. K(GeV) θi(deg) SG4(VEM) SUSC (VEM) 0.1 45. 0.179 ±0.112 0.196 (9) 0.1 75. 0.171 ±0.118 0.192 (12) 0.4 45. 0.801 ±0.349 0.878 (10) 0.4 75. 0.965 ±0.441 1.260 (30) 1. 45. 1.010 ±0.481 1.016 (0.6) 1. 75. 1.565 ±0.940 1.611 (3) 10. 45. 1.138 ±0.630 1.123 (-1) 10. 75. 1.921 ±1.341 1.814 (-6) Table 4.3: Muon signal in VEM in an Auger tank as obtained in Geant4 and the S(1000) USC code for different kinetic energies and angles of incidence. The results of Geant4 show the average over all impact parameters. The numbers in parenthesis indicate the relative differences ((SG4−SUSC )/SG4in %) using the Geant4 result as reference. For the electromagnetic particles, the relative differences are ∼25% at most, with the electromagnetic signal obtained with the S(1000) USC typically higher. This discrepancy is expected because the total tracklength of an EM subshower in the the ZHS code (the results of which are used in S(1000) USC) is about 10% larger than the track obtained with Geant4 as discussed in [79]. This difference between ZHS and Geant4 is due to the different implementation of the relevant electromagnetic processes. The difference seems to be larger at high θ, however the contribution of the EM component to the total signal in large θshowers is expected to be small (<15%) and hence the 80 impact of these differences in the total signal is expected to be smaller than 3%. In conclusion, the agreement between S(1000) USC code and Geant4 is generally good with differences of less than 10% for muons and 20% for electrons, positrons and gammas, all within the RMS of the Geant4 output. Finally it is important to remark that the distributions of the signal obtained with Geant4 have fairly large RMSs, mainly due to the correct account for the variations in particle tracklengths correlated with the different impact parameter of the particles [97]. K(GeV) θi(deg) SG4(VEM) SUSC (VEM) 0.01 45. 0.022 ±0.024 0.022 (0) 0.01 75. 0.016 ±0.015 0.013 (-19) 0.04 45. 0.138 ±0.063 0.142 (3) 0.04 75. 0.124 ±0.071 0.140 (13) 0.1 45. 0.330 ±0.129 0.376 (14) 0.1 75. 0.351 ±0.172 0.390 (11) 1. 45. 2.115 ±1.143 2.248 (6) 1. 75. 2.681 ±1.472 3.305 (23) Table 4.4: Signal produced by an electron in an Auger tank in Geant4 and the S(1000) USC code for different kinetic energies and angles of incidence. The results of Geant4 show the average over all impact parameters. The numbers in parenthesis indicate the relative differences ((SG4−SUSC )/SG4in %) taking the Geant4 result as reference. 81 K(GeV) θi(deg) SG4(VEM) SUSC (VEM) 0.01 45. 0.021 ±0.018 0.023 (9) 0.01 75. 0.016 ±0.015 0.014 (-12) 0.04 45. 0.137 ±0.077 0.142 (4) 0.04 75. 0.124 ±0.069 0.140 (13) 0.1 45. 0.314 ±0.128 0.373 (19) 0.1 75. 0.352 ±0.160 0.389 (10) 1. 45. 2.056 ±1.133 2.219 (8) 1. 75. 2.797 ±1.391 3.278 (17) Table 4.5: Signal produced by an positron in an Auger tank in Geant4 and the S(1000) USC code for different kinetic energies and angles of incidence. The results of Geant4 show the average over all impact parameters. The numbers in parenthesis indicate the relative differences ((SG4−SUSC )/SG4in %) taking the Geant4 result as reference. K(GeV) θi(deg) SG4(VEM) SUSC (VEM) 0.01 45. 0.028 ±0.018 0.032 (14) 0.01 75. 0.029 ±0.019 0.037 (28) 0.04 45. 0.107 ±0.076 0.128 (20) 0.04 75. 0.126 ±0.087 0.154 (22) 0.1 45. 0.262 ±0.158 0.304 (16) 0.1 75. 0.314 ±0.196 0.380 (21) 1. 45. 1.653 ±1.255 1.632 (-2) 1. 75. 2.359 ±1.527 2.795 (18) Table 4.6: Signal produced by a gamma in an Auger tank as obtained in Geant4 and the S(1000) USC code for different kinetic energies and angles of incidence. The results of Geant4 show the average over all impact parameters. The numbers is parenthesis indicate the relative differences ((SG4−SUSC )/SG4in %) taking the Geant4 results as reference. 82 Chapter 5 Study of the signals in inclined showers: the role of the electromagnetic halo The conventional separation between vertical and horizontal (inclined) showers is based on the zenith angle θof the particle that induces the shower: horizontal showers are defined as those with 60◦< θ < 90◦. The differences between vertical and horizontal showers come from the different atmospheric grammage that the showers have to cross before reaching the ground which increases approximately as sec(θ). For instance, the slant depth of atmosphere for a completely vertical shower θ= 0◦is ∼879.6 g cm−2for the Pierre Auger Observatory altitude, increasing to ∼1760 g cm−2for a shower at 60◦and being about 35 times larger for a completely horizontal shower. Nucleonic cosmic rays initiate air showers at the top of the atmosphere in the first few 100 g cm−2. For instance in Fig. 5.1 we show the typical longitudinal development of a 10 EeV proton shower. The electromagnetic (EM) component of the shower rises as the shower penetrates and reaches a maximum that in this example is at a depth Xmax ∼780 g cm−2. After Xmax the EM component is rapidly absorbed in the atmosphere due to low-energy processes and the photoelectric effect. Meanwhile, non-decaying muons propagate practically unattenuated to the ground, except for energy loss and deflections in the geomagnetic field. Therefore, a 10 EeV energy shower at θ= 0◦reaches the ground level shortly after reaching maximum and the electromagnetic component dominates at ground. However, in horizontal showers muons dominate at the ground level because the electromagnetic component due to cascading processes, i.e. from π0decay is largely absorbed before reaching the ground [89]. However though small there is still an electromagnetic component in inclined showers. This is the so-called 83 5.2 Azimuthal asymmetries without the geomagnetic field In the previous section, we have studied the lateral distribution of the ratio SEM /Sµunder the assumption that the signals are equal at the same distance from the shower axis in the shower plane regardless of the azimuthal angle ζ. This assumption is only an approximation, and in fact there is an azimuthal asymmetry in the signal due to the several effects, the most important being the so-called geometrical effect, the longitudinal development effect and ground screening [90, 91, 92]. Furthermore, the geomagnetic field is another source of asymmetry in inclined showers which for the moment we will neglect and defer its study to the next section. 5.2.1 The geometrical effect In Fig. 5.6 we show a sketch of an inclined shower hitting the ground. This sketck serves us to illustrate that the shower particles do not travel parallel to the shower axis and hence cross different paths before reaching the ground depending on their azimuthal angles. Moreover, particles hit the detectors in the “early” region (before the shower axis hits the ground) “more vertically” than the ones that hit the tanks located in the “late” region. This is essentially the basis of the so-called geometrical effect. To demonstrate this behaviour, we show in Fig. 5.7 the distributions of angles of incidence θiof early and late particles in simulated showers at different zenith angles (θ). One can see that the mean θiof the late muons is always larger than the corresponding mean of the early muons. Also the difference between the mean values decreases with θ. For the electromagnetic particles the same behaviour occurs. As a conclusion by inspecting Fig. 5.7, the asymmetry induced by the difference in the angle of incidence between the particles reaching the early and late regions of a shower, i.e. the geometrical effect, is expected to be more important in showers with small zenith angle. The geometrical effect is expected to affect differently the EM and muonic components of a shower: •The electromagnetic signal is roughly proportional to the area of the tank projected onto the plane perpendicular to the particle direction (see Eq. 4.47). The area decreases with zenith angle θi(see Fig. 4.2). As a consequence we expect the electromagnetic signal to be larger in the early region than in the late one. •The muonic signal is roughly proportional to the track-length in wa- 90 +∆XX X−∆X p θ p θ Core plane θ Ground screening r r Early plane Late plane X shower axis Figure 5.6: Schematic picture of an inclined shower reaching the ground. Three planes are displayed intersecting the ground plane, each one at different depths on the shower development: early (red), late (blue) and core (black) planes. The latter is also called shower plane. ter (see Eq. 4.46), and the mean track-length increases with θi(see Fig. 4.3). On the other hand, the area of the tank projected onto the plane perpendicular to the particle direction decreases with θi. Therefore, there should be a large degree of compensation between both behaviours and we expect the muonic signal to be approximately the same in the early and late regions 2. 5.2.2 The longitudinal development effect The longitudinal development effect can be understood as follows. Particles at the same distance to the core in the shower plane r, but arriving with different azimuthal angles ζtravel along different paths, and they belong to different stages in the evolution of the shower. The importance of this effect depends on the evolution with depth of the lateral distribution and 2In fact the uncorrected muon signal is independent of θ(see Eq. 4.16) 91 Eµ Entries 100 Mean 16.96 RMS 8.654 i θ 0 10 20 30 40 50 60 70 80 90 a.u. 0.05 0.1 0.15 0.2 0.25 0.3 Eµ Entries 100 Mean 16.96 RMS 8.654 Lµ Entries 100 Mean 41.94 RMS 8.918 Lµ Entries 100 Mean 41.94 RMS 8.918 early late µ°30 EM E Entries 100 Mean 19.02 RMS 13.69 i θ 0 10 20 30 40 50 60 70 80 90 a.u. 0 0.02 0.04 0.06 0.08 0.1 0.12 0.14 0.16 0.18 EM E Entries 100 Mean 19.02 RMS 13.69 EM L Entries 100 Mean 53.7 RMS 12.99 EM L Entries 100 Mean 53.7 RMS 12.99 EM Eµ Entries 100 Mean 52.59 RMS 6.552 i θ 0 10 20 30 40 50 60 70 80 90 a.u. 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 Eµ Entries 100 Mean 52.59 RMS 6.552 Lµ Entries 100 Mean 62.67 RMS 2.16 Lµ Entries 100 Mean 62.67 RMS 2.16 early late µ °60 EM E Entries 100 Mean 34.33 RMS 14.66 i θ 0 10 20 30 40 50 60 70 80 90 a.u. 0.05 0.1 0.15 0.2 0.25 0.3 0.35 0.4 EM E Entries 100 Mean 34.33 RMS 14.66 EM L Entries 100 Mean 60.42 RMS 13.74 EM L Entries 100 Mean 60.42 RMS 13.74 EM Eµ Entries 100 Mean 66.61 RMS 2.019 i θ 0 10 20 30 40 50 60 70 80 90 a.u. 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 Eµ Entries 100 Mean 66.61 RMS 2.019 Lµ Entries 100 Mean 70.88 RMS 1.337 Lµ Entries 100 Mean 70.88 RMS 1.337 early late µ°70 EM E Entries 100 Mean 59.05 RMS 12.62 i θ 0 10 20 30 40 50 60 70 80 90 a.u. 0 0.05 0.1 0.15 0.2 0.25 0.3 0.35 EM E Entries 100 Mean 59.05 RMS 12.62 EM L Entries 100 Mean 66.85 RMS 8.728 EM L Entries 100 Mean 66.85 RMS 8.728 EM Figure 5.7: Distribution of the zenith angle of incidence of the shower particles in the early (solid line) and late regions (dashed line) at 1000 m from the shower axis for different shower zenith angles. Top panels: θ= 30◦. Middle panels: θ= 60◦. Bottom panels: θ= 70◦. For each angle we show the distributions for the muonic (left panel) and electromagnetic (right panels) components. The simulations were done for 10 EeV proton showers with AIRES and the QGSJET model. 92 on the attenuation of the total number of particles. In Fig. 5.6 we show two detectors at the same distance from the core and two planes transverse to the shower axis containing their positions. These planes serve to illustrate the different values of atmospheric depth crossed by the particles, and therefore the different stages of shower development. We also draw the transverse plane containing the impact point of the shower axis on the ground (shower plane). From Fig. 5.6 it is evident that the tank in the early region is hit by a younger stage on the evolution of the shower than the tank in the late region (for more details see Section 6.3). For instance, in an event produced by a 10 EeV proton shower at θ= 60◦the depth crossed by an early and a late particle hitting tanks at r= 1000 m in the shower plane differs by ∼370 g cm−2in slant depth. The difference of grammage crossed by the early and late particles increases with the distance from the core. The asymmetry introduced by the effect of the longitudinal development is more important for the electromagnetic component from π0decay. This component is exponentially suppressed after the shower maximum, and as a consequence small changes in the depth crossed induce large differences in the number of EM particles on the ground. However, the muonic component is less attenuated and therefore the asymmetry induced by this effect is smaller. Therefore, we expect that the contribution of this effect to the azimuth asymmetry of the signal is small at large zenith angles (θ > 70◦) at which the electromagnetic component from π0decay on the ground is practically suppressed at all azimuth angles and the electromagnetic halo inherits the behaviour of the muonic component. In Fig. 5.8 we illustrate the effect of the longitudinal development on the electromagnetic and muon components of the signal by plotting the lateral distributions of electromagnetic energy density (top panels) and the lateral distributions of muon number density (bottom panels) in two ranges of the azimuthal angle: early (ζaround 0◦) and late (ζaround 180◦) The size of the azimuthal bins in these plots is ∆ζ= 30◦in order to have an acceptable particle statistics. For showers at θ= 60◦the difference in the EM energy density between the early-late regions is important at even small distances to the core (r∼100 m) because the EM component from π0decay is still significant in the early region, while it is practically absorbed before reaching ground in the late region. However, at θ= 70◦the difference between the densities in the early and late regions is small because the component from π0decay is absorbed for all ζ(see Fig. 5.1) and the EM halo, produced by the decay of muons which are less affected by the longitudinal development effect, dominates. In the bottom panels of Fig. 5.8 we plot the lateral distributions of the muon number density. It can be seen that the difference between the densities in the early and late regions is always small regardless of the zenith 93 angle. The behaviour of the early-late asymmetry with ris shown in Fig. 5.9, where we plot the relative differences between the early and late electromagnetic energy densities (left panel) and muon number densities (right panel). At θ= 60◦the asymmetry in the EM energy density increases rapidly with r. However, at θ= 70◦there is only a slight increase above r= 1000 m which follows the same behaviour seen in the muon number density (see right panel), because at θ= 70◦the electromagnetic component is mostly due to muon decay. 5.2.3 The screening effect Finally, the azimuthal asymmetry in the signal is also induced by the so-called screening effect. This effect is produced by the absorption of the shower core after its impact on the ground. As a consequence, the hadronic core of the shower stops feeding the EM and muonic components in a portion of the late region (see Fig.5.6). The combination of these 3 effects produces an azimuthal asymmetry in the lateral distribution of the electromagnetic and muonic components of the signal. The asymmetry in the signal at a fixed ζcan be quantified defining an asymmetry parameter Asym: Asym(r) = S(ζ)−< S > < S > (5.5) where < S > is the signal averaged over all ζ. In Fig. 5.10 we plot the asymmetry parameter of the muonic (left panel) and electromagnetic (right panel) signals as a function of the distance to the core for showers at θ= 60◦. We have plotted the asymmetry parameter in 4 azimuthal bins of size ∆ζ= 30◦:ζ= 0◦, 90◦, 180◦and 270◦. The behaviour of the asymmetry depends strongly on the type of signal (EM or muonic) and therefore, the asymmetry will not cancel out when the ratio of the EM signal to the muonic signal is calculated. For example, in the left panel of Fig. 5.11 we show the 2-dimensional map of the ratio SEM /Sµin the shower plane for 10 EeV proton showers at θ= 60◦. The arrow shows the shower direction. One can clearly see the azimuthal asymmetry in the ratio SEM /Sµ at a fixed distance to the core (indicated with the head of the arrow). The fraction of electromagnetic signal is larger in the early region. However, at zenith angles greater than 70◦as shown in the right panel of Fig. 5.11, no azimuthal asymmetry in the ratio is observed because there are only muonic 94 [r/m] 10 log 1 1.5 2 2.5 3 3.5 ] -2 [GeV m EM ∈ -4 10 -3 10 -2 10 -1 10 1 10 2 10 3 10 4 10 °60 (early)°0 (late)°180 [r/m] 10 log 1 1.5 2 2.5 3 3.5 ] -2 [GeV m EM ∈ -4 10 -3 10 -2 10 -1 10 1 10 2 10 3 10 4 10 °70 [r/m] 10 log 1 1.5 2 2.5 3 3.5 ] -2 [m µ ρ -3 10 -2 10 -1 10 1 10 2 10 °60 (early)°0 (late)°180 [r/m] 10 log 1 1.5 2 2.5 3 3.5 ] -2 [m µ ρ -3 10 -2 10 -1 10 1 10 2 10 °70 Figure 5.8: Top panels: Lateral distributions of the electromagnetic energy density in the early (ζaround 0◦) and late (ζaround 180◦) regions of the shower plane of showers at θ= 60◦(left panel) and 70◦(right panel). Bottom panels: Lateral distributions of the muon number density in the early (ζaround 0◦) and late (ζaround 180◦) regions of the shower plane of showers at θ= 60◦(left panel) and 70◦(right panel). Each distributions corresponds to the average of 100 proton showers of E= 10 EeV simulated with AIRES. 95 [r/m] 10 log 1 1.5 2 2.5 3 3.5 (late) EM ∈(late)] / EM ∈(early)- EM ∈[ -1 0 1 2 3 4 5 6 7 °60°70 [r/m] 10 log 1 1.5 2 2.5 3 3.5 (late) µ ρ(late)] / µ ρ(early)- µ ρ[ -1 0 1 2 3 4 5 6 7 °60°70 Figure 5.9: Left panel: Early-late asymmetry of the electromagnetic energy density as a function of the distance from the core in the shower plane for showers at θ= 60◦(full circles) and θ= 70◦(empty circles). Right panel: Early-late asymmetry of the muon number density as a function of the distance from the core in the shower plane for showers produced at θ= 60◦(full squares) and θ= 70◦(empty squares). [r/m] 10 log 1 1.5 2 2.5 3 3.5 µ Asym -1.5 -1 -0.5 0 0.5 1 1.5 °60 [r/m] 10 log 1 1.5 2 2.5 3 3.5 µ Asym -1.5 -1 -0.5 0 0.5 1 1.5 (early)°0°90 (late)°180°270 °60 [r/m] 10 log 1 1.5 2 2.5 3 3.5 EM Asym -1.5 -1 -0.5 0 0.5 1 1.5 °60 [r/m] 10 log 1 1.5 2 2.5 3 3.5 EM Asym -1.5 -1 -0.5 0 0.5 1 1.5 (early)°0°90 (late)°180°270 °60 Figure 5.10: Asymmetry of the lateral distribution of the muonic (left panel) and electromagnetic (right panel) signal components with respect to the mean value for different azimuth regions. The distributions correspond to the average of 100 proton showers of 10 EeV simulated at 60◦with AIRES + QGSJET + S(1000) USC code. 96 component and the EM halo on the ground. Since these two components approximately have the same asymmetry, the final asymmetry is practically canceled out when making the ratio of EM and muonic signals. x (m) -3000 -2000 -1000 0 1000 2000 3000 y (m) -3000 -2000 -1000 0 1000 2000 3000 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 °=60θ 10EeV µ /S EM S x (m) -3000 -2000 -1000 0 1000 2000 3000 y (m) -3000 -2000 -1000 0 1000 2000 3000 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 °=70θ 10EeV µ /S EM S Figure 5.11: Left panel: Contour map of the ratio SEM /Sµon the transverse (shower) plane in θ= 60◦showers. Right panel: Contour map of the ratio SEM /Sµon the transverse (shower) plane in θ= 70◦showers. The arrow shows the shower direction and the head of the arrow shows the core position. Each map is obtained using the average of 100 proton showers of E= 10 EeV simulated with AIRES + QGSJET + S(1000) USC code. In Fig. 5.12 we plot the ratio SEM /Sµas a function of the azimuth angle for a fixed distance r= 1000 m. We use a system of polar coordinates (SEM , ζ). At θ= 60◦(blue points) there is a large asymmetry, the signal ratio at ζ= 0◦is more than twice the ratio at ζ= 180◦. Note also the symmetry in the ratio at ζ=±90◦. At θ= 70◦(red squares) there is essentially circular symmetry for the reasons explained before. The fact that the ratio SEM /Sµdepends on the azimuthal angle must be taken into account in the analysis and reconstruction of inclined showers. For this purpose we have performed a parameterization of the ratio SEM /Sµas a function of the distance from the core r, zenith angle θ, and azimuth angle ζ, SEM /Sµ(r, θ, ζ) = < SEM /Sµ(r, θ)>(1 + Aasym(r, θ, ζ)) (5.6) where < SEM /Sµ>is the parameterization of the ratio of signals averaged over all ζangles given in Eq. 5.1 and the parameter Aasym characterizes the azimuthal asymmetry (in absence of geomagnetic field). This correction is more important in the range θ < 70◦. For larger angles the asymmetry is expected to be negligible as shown before (see for instance Fig. 5.11). 97 0.6 0.4 0.2 0 0.2 0.4 0.6 0.6 0.4 0.2 0 0.2 0.4 0.6 0o ξ=90o 180o 270o r=1000 m Figure 5.12: The ratio of the electromagnetic to the muonic signals as a function of the azimuth angle ζfor r= 1000 m in polar coordinates (SEM /Sµ, ζ). The blue circles correspond to the average of 100 proton showers of E= 10 EeV at θ= 60◦. The red squares correspond to showers at θ= 70◦. The showers were simulated with AIRES + S1000 USC code. In Fig. 5.13 we show the ratio SEM /Sµas a function of rin different bins in ζcompared to the mean value (left panel) for showers at 60◦and their corresponding asymmetry parameter Aasym (right panel). We parameterize Aasym using the following equation: Aasym(r, θ, ζ) = D(θ, ζ)r+E(θ, ζ)log10r(5.7) with the distance rin meters and the angles θand ζin degrees. The fit is valid in the range θ∈[60◦,69◦], log10r∈[1., 3.8] m and ζ∈[−180◦,180◦]. D(θ, ζ) in Eq. 5.7 is parameterized as: D= 10−5×[D1+D2sinc(πD3ζ) + sinc2(πD3ζ)] (5.8) where sinc(x) = sin x/x and: D1=−193.143 + 9.454 θ−0.154 θ2+ 0.0008 θ3 D2= 1.6−7.909 ×10−6θ2 1−0.014 θ 98 D3= 21.401 −1.005 θ+ 0.016 θ2−8.208 ×10−5θ3(5.9) E(θ, ζ) in Eq. 5.7 is parameterized as: E=E1cos(E2ζ) (5.10) with: E1= 2.052 −0.053 θ+ 0.0003 θ2 and E2= 1.08 −4.467 ×10−8θ3 1−0.014 θ(5.11) [r/m] 10 log 1 1.5 2 2.5 3 3.5 µ /S EM S -1 10 1 10 (early)°0°90 (late)°180°270 Average °10 EeV 60 [r/m] 10 log 1 1.5 2 2.5 3 3.5 ratio Asym -2 -1 0 1 2 Figure 5.13: Left panel: The ratio of the electromagnetic to muon contributions to the tank signal as a function of the distance from the shower axis in the shower plane in different bins in ζ. Right panel: Asymmetry of the lateral distribution of the ratio SEM /Sµ in different ζbins. The distributions correspond to the average of 100 proton showers with E=10 EeV at θ= 60◦, simulated with AIRES+QGSJET+S1000USC. To test the accuracy of this parameterization, we have compared the simulated ratio and the parameterization at ζ= 0◦and ζ= 180◦in different bins in log10 r. In Fig. 5.14 we show the histogram of the results of this comparison for two zenith angles: θ= 60◦(left panel) and θ= 68◦(right panel). The mean values of the ratio are well reproduced by the parameterization within <10% (the RMS of the distributions are all <25%). 99 value) while at θ= 86◦it is ±1.5 µT (11% of the central value). The central value of the oscillation increases with the zenith angle. By inspecting Fig. 5.20, we expect the horizontal deflection to be minimum at φ= 87◦and maximum at φ= 267◦(the same azimuth angles at which Bg ⊥= 0). [deg]φ 0 50 100 150 200 250 300 350 T]µB [ -25 -20 -15 -10 -5 0 5 10 15 20 25 BperpS ° = 60θ ° = 65θ ° = 70θ ° = 75θ ° = 80θ ° = 85θ ° = 89θ s B g B Figure 5.20: The components Bg ⊥and Bs ⊥of B⊥plotted as a function of the azimuth direction of the shower for different zenith angles ranging from 60◦to 89◦. Note that Bg ⊥ is independent of θ. To illustrate how both types of deflection affect the shape of the signal maps, we show in Fig. 5.21 the muon (on the left-hand side) and electromagnetic (on the right-hand side) signal maps in the shower plane for 10 EeV proton induced showers at θ= 86◦and different shower azimuthal directions (assuming the Auger convention with φ= 0◦corresponding to the geographical East). We choose θ= 86◦because the geomagnetic deviation of muon trajectories is expected to be very important at this angle (see the right panel of Fig. 5.16). In Fig. 5.20 one can also see that the strength of Bs ⊥varies very little with the azimuthal angle for θ= 86◦(the length of the semi-major of the ellipse will not change much) and as a consequence, the change with φ of the shape of the muon maps will mostly depend on the intensity of Bg ⊥. For instance, at φ= 0◦(top panels) and 180◦(bottom panels) the shower is 106 practically arriving perpendicular to the direction of the geomagnetic field. By inspecting Fig. 5.20 one can observe that both the vertical and horizontal deflections are relevant at both azimuthal angles. The “horizontal deflection” is the cause of the lobular structure and the “vertical deflection” produces the asymmetry in the size of the lobes. The sign of Bg ⊥determines which of the two lobes has a higher muon density. At φ= 0◦the component Bg ⊥is <0, so that a fraction of the negative muons are removed from the ground and consequently the µ−lobe becomes smaller (see Fig. 5.21). On the contrary, at φ= 180◦the component Bg ⊥is >0 and a fraction of positive muons are deviated away from the ground, turning the µ+lobe into the one with the smallest muon density. In the middle panel of Fig. 5.21, we show the muon map for a shower at φ= 90◦. In this case, the only relevant deflection is the horizontal one (Bg ⊥= 0) and therefore, the map exhibits a mirror symmetry with respect to the arrival direction of the shower, i.e. both lobes have the same size. In the panels on the right-hand side of the Fig. 5.21, we show the signal maps of the electromagnetic component for the same showers. These maps exhibit a similar behaviour to their corresponding muon maps. This is due to the fact that the electromagnetic particles at large zenith angles come mostly from muon decay and therefore, preserve to some extent the muon spatial distribution. However, the electromagnetic maps have less sharp patterns than the muon maps because the deflection of electrons and positrons is dominated by multiple Coulomb scattering. It is clear that the lateral distribution of the SEM /Sµratio averaged over all ζshould be strongly modified by the presence of the geomagnetic field at very large zenith angles. We can infer from Figs. 5.16 and 5.20 what is the effect of the geomagnetic field on the lateral distribution of SEM /Sµdepend on the shower zenith and azimuth angle. On one hand, the effect of the geomagnetic field is expected to be more relevant the larger the zenith angle. On the other hand, one expects the difference in the lateral distribution of SEM /Sµwith and without geomagnetic field to be minimum in a shower with φ= 90◦and maximum at 180◦. For this reason, we study the effect of the geomagnetic field on the SEM /Sµlateral distribution of showers induced at different zenith angles and at φ= 90◦and 180◦for each θ. In the left panels of Figs. 5.22 and 5.23 we show the lateral distribution of SEM /Sµin the presence of the geomagnetic field for proton induced showers at different θarriving at φ= 90◦and 180◦in each case. In the same panels, we also show the corresponding lateral distribution obtained neglecting the geomagnetic effect. We also plot in the right panels of Figs. 5.22 and 5.23 the relative difference between the distributions with and without the geomagnetic field, taking as reference the case without field (Rt in the figure 107 designates the ratio SEM /Sµ) as a function of the distance to the shower core. In Fig.5.22 one can see that the relative differences are .20% for the two azimuthal angles and for θ≤80◦and log10r > 1.5, and as a consequence the effect of the geomagnetic field remains negligible at θ≤80◦. In Fig. 5.23 one can see that the geomagnetic effect starts to be relevant at θ= 82◦for the case of φ= 180◦(maximum deviation) where the relative difference is >20% for log10r < 2, whereas at φ= 90◦(minimum deviation) the relative difference remains smaller than 20% for log10r > 1.5. At larger angles the situation changes and the geomagnetic field has a strong effect on the SEM /Sµ distribution, even for the azimuth angle of the shower at which the effect is expected to be minimum. For instance, at θ= 86◦the relative difference is much larger than 20% at φ= 180◦in all the range of distances to the shower core, and also at φ= 90◦for log10r < 2.5. SEM /Sµincreases the most near the shower core when the geomagnetic field is included. The reason is that only the highest energy muons are not significantly deflected by the geomagnetic field and there are more likely to suffer hard interactions where an EM shower is produced. As a consequence SEM increases and at the same time Sµdecreases because lower energy muons are being moved away from the core with the overall effect of increasing SEM /Sµ. In conclusion, for the purposes of event reconstruction the effect of the geomagnetic field on the SEM /Sµlateral distribution must be taken into account only when θ&86◦. 108 B µ µ+ − B µ µ− + B µ µ+ − Figure 5.21: Muon (left) and electromagnetic (right) signal maps in the shower plane for 10 EeV proton showers with an incident zenith angle of 86◦as obtained in AIRES simulations for the following azimuthal angles in the Auger convention: 0◦(top), 90◦(middle) and 180◦ (bottom). The reference system has the x-axis pointing out in the shower direction (white arrow). The black arrow indicates the direction of ~ B. Note that ~ Bindicates the total field (and not only the component of it perpendicular to the shower axis that is responsible for the muon and EM deflection). 109 [r/m] 10 log 1 1.5 2 2.5 3 3.5[r/m] 10 log 1 1.5 2 2.5 3 3.5 µ /S EM S -1 10 1 ° = 70θ10 EeV ° = 90φB ° = 180φB No B [r/m] 10 log 1 1.5 2 2.5 3 3.5[r/m] 10 log 1 1.5 2 2.5 3 3.5 NoB ) / Rt NoB -Rt B ( Rt -0.4 -0.2 0 0.2 0.4 ° = 90φ° = 180φ [r/m] 10 log 1 1.5 2 2.5 3 3.5[r/m] 10 log 1 1.5 2 2.5 3 3.5 µ /S EM S -1 10 1 ° = 80θ10 EeV ° = 90φB ° = 180φB No B [r/m] 10 log 1 1.5 2 2.5 3 3.5[r/m] 10 log 1 1.5 2 2.5 3 3.5 NoB ) / Rt NoB -Rt B ( Rt -0.4 -0.2 0 0.2 0.4 0.6 0.8 ° = 90φ° = 180φ Figure 5.22: Left panel: The lateral distribution of the ratio of the electromagnetic to muon signals in the shower plane for showers at φ= 90◦(circles) and φ= 180◦(squares) under the presence of the geomagnetic field, compared with the distribution without the effect of the geomagnetic field. Right panel: The relative differences between the distributions without and with geomagnetic field effect (see text) for showers at φ= 90◦(circles) and φ= 180◦(squares). The distributions correspond to the average of 100 proton showers with 10 EeV simulated with AIRES + S1000 USC code at θ= 70◦(top) and 80◦(bottom). 110 [r/m] 10 log 1 1.5 2 2.5 3 3.5[r/m] 10 log 1 1.5 2 2.5 3 3.5 µ /S EM S -1 10 1 ° = 82θ10 EeV ° = 90φB ° = 180φB No B [r/m] 10 log 1 1.5 2 2.5 3 3.5[r/m] 10 log 1 1.5 2 2.5 3 3.5 NoB ) / Rt NoB -Rt B ( Rt -0.5 0 0.5 1 1.5 2 2.5 3 3.5 ° = 90φ° = 180φ [r/m] 10 log 1 1.5 2 2.5 3 3.5[r/m] 10 log 1 1.5 2 2.5 3 3.5 µ /S EM S -1 10 1 10 ° = 86θ10 EeV ° = 90φB ° = 180φB No B [r/m] 10 log 1 1.5 2 2.5 3 3.5[r/m] 10 log 1 1.5 2 2.5 3 3.5 NoB ) / Rt NoB -Rt B ( Rt -1 0 1 2 3 4 5 6° = 90φ ° = 180φ [r/m] 10 log 1 1.5 2 2.5 3 3.5[r/m] 10 log 1 1.5 2 2.5 3 3.5 µ /S EM S -1 10 1 10 ° = 88θ10 EeV ° = 90φB ° = 180φB No B [r/m] 10 log 1 1.5 2 2.5 3 3.5[r/m] 10 log 1 1.5 2 2.5 3 3.5 NoB ) / Rt NoB -Rt B ( Rt 0 2 4 6 8 10 ° = 90φ° = 180φ Figure 5.23: Left panel: The lateral distribution of the ratio of the electromagnetic to muon signals in the shower plane for showers at φ= 90◦(circles) and φ= 180◦(squares) under the presence of the geomagnetic field, compared with the distribution without the effect of the geomagnetic field. Right panel: The relative differences between the distributions without and with geomagnetic field effect (see text) for showers at φ= 90◦(circles) and φ= 180◦(squares). The distributions correspond to the average of 100 proton showers with 10 EeV simulated with AIRES + S1000 USC code at θ= 82◦(top), 86◦(middle) and 88◦(bottom). 111 Chapter 6 Identification of Neutrino Candidates in surface detector of the Pierre Auger Observatory The main experimental challenge in the detection of neutrino-induced showers with the surface detector of the Pierre Auger Observatory is to identify them in the background of showers initiated by nucleonic cosmic rays. Deeply penetrating highly energetic particles such as neutrinos, can initiate showers very close to the ground level while protons, heavier nuclei and photons interact shortly after entering the atmosphere. As suggested almost 30 years ago, the observation of inclined showers enhances the difference between these two types of showers [53]. Therefore the main signature of down-going neutrino events are inclined showers that interact deep in the atmosphere (Deep Inclined Showers from now on). 6.1 Selection of Inclined Events in the data set recorded by the Surface Detector Down-going neutrino showers are searched among the inclined data set registered with the Surface Detector (SD) of the Pierre Auger Observatory. An efficient selection of inclined events in the data set is crucial for the subsequent identification of neutrino candidates. The SD data set consists of all T3 level triggers acquired by the Central Trigger System (see Section 3.3.2). The higher levels of trigger (T4 and T5 levels) developed for vertical showers to select real events, are not suitable for 113 inclined showers. In the case of inclined showers, the selection of real showers in the data set can be made a posteriori by means of algorithms of station selection and event reconstruction. We have developed an event selection for inclined showers that follows a bottom-up procedure applied to the data set collected at the surface detector and that consists of several steps: preselection of time-constrained configurations of stations; selection of a reconstruction seed; selection of candidate stations studying the space-time compatibility with the seed, and reconstruction of the resulting configuration of candidate stations which is finally accepted as a physical event if the reconstruction succeeds. The selection procedure presented here is aimed at obtaining an optimum efficiency in the selection of inclined showers produced by neutrinos. The advantage with respect to the current algorithms used to select conventional nucleonic showers [71], is that the algorithms presented here have been developed taking into account the topological characteristics of the neutrino events obtained in Monte Carlo simulations, in order to minimize the rejection of potential neutrino candidates. Also the tolerances of the tank selection have been optimized to avoid rejecting the earliest stations in the event, which as we will show below are crucial for neutrino identification. However and despite this fact, the algorithms presented here still share many of the ideas developed by the inclined shower group of the University of Santiago de Compostela [101, 67] to select conventional nucleonic inclined showers. 6.1.1 Station selection Developing criteria for station rejection is essential to select physical events. These criteria have been developed taking into account the particular timing and topological characteristics of inclined showers. The spatial configuration of the stations in an inclined event can be of two types: non-aligned or aligned, which have to be treated with different selection algorithms. In the following, the different steps used to select candidate stations are described in the same order they are applied to real data. The first four algorithms are common to non-aligned and aligned configurations, whereas the remainder are different for each type of configuration. Rejection of Engineering Array stations The Engineering Array (EA) is a small subset of 100 stations originally built for testing the design of the observatory at the beginning of the project. The stations that belong to the EA are different from the ones currently being 114 deployed (components, electronics,...). All the remaining stations have the same characteristics. Due to this, the stations belonging to the EA must be removed from the event. This is easy since each of the stations of the Surface Detector is identified with an unique identifier number (ID), which is less than 100 for the stations belonging to the EA. Treatment of twin stations If both stations of a twin pair 1belong to the same event, the one with the higher ID is removed by default from the event because this station is not part of the central trigger system. Selection of time clusters The main aim here is to preselect configurations of stations that will be tested by the subsequent algorithms until a satisfactory configuration is found. A good criterion to build them up is to look for time-constrained configurations, called time clusters. The procedure to build up a time cluster is as follows. After applying the two previous criteria, the remaining stations are sorted by increasing start-time. We begin with the earliest station and check if the next one in time is closer than 16 µs. If this is the case both stations are grouped into the same cluster. We keep adding stations to this cluster applying the same procedure to the following stations in time. If we find one station farther than 16 µs from its predecessor, a new and different cluster is built up using that station as starting point. If multiple time clusters are found, these are sorted by increasing quality. The quality score is based on compactness in time of the stations belonging to the cluster. The score is defined as the number of stations whose start-time difference is less than 6 µs. At the end of this procedure, we have a set of time clusters sorted by increasing quality. Selection of the best seed The main aim here is to find a good seed of three tanks which is used afterwards as the base for selecting the candidate stations that will be involved in the event reconstruction. The selection of the seed should be robust enough 1A twin pair are two stations that are separated by around 11 meters, which serve the purpose of studying the accuracy of the angular and signal determination. The default separation between tanks in the surface array is 1.5 km. 115 Non-aligned configuration In the case of a non-aligned seed, we study the compatibility of a station iwith the seed using the zenith (θseed) and azimuth angles of the seed (φseed). For every station i, different from the seed stations, Eqs. (6.4) and (6.5) are applied to estimate the shower angles from the triangle formed by station i and the 2 stations in each side of the seed (θil, φil with l= 1,2,3 labeling the sides of the triangle). For each side, we check if the station iis in line with the side lof the seed. If this is not the case, the angles θil and φil are considered to be compatible with the angles of the seed if |sin θil −sin θseed|<0.2 |φil −φseed|<10 deg (6.11) If these conditions are fulfilled for one of the sides of seed, the station gains a score. Finally, station iis accepted as a station of the event if its final score is at least 2. Aligned configuration In the case of configurations in which the seed stations are in line, a different procedure is applied to check the space-time compatibility of a station. It is based on the fact that the apparent transmission speed of the signal between any two stations along the shower direction should not vary significantly for whatever couple of stations in the event. For every station idifferent from the seed stations, we test its compatibility with the apparent transmission speed of the signal in the seed, given by vseed =c/ sin θseed. Eq. (6.8) is applied to calculate the apparent transmission speed of the signal between the station iand each seed station (vil with l= 1,2,3 labeling the stations of the seed). For each station of the seed, we first check if the position of station iis ordered in start-time with the seed station following the procedure described in the previous section (see Eq. 6.6). If this is the case, a speed vil is considered to be compatible with vseed if, |vil −vseed| vseed <0.16 (6.12) The 16% tolerance allows small variations of the apparent transmission speeds of the signal because of the variable curvature of the shower front in the shower direction. For instance, for a pair of stations in the early region of the shower the apparent speed of the signal may be larger than that for a pair of stations in the late region. Station iis accepted as a station of the 122 event if the 3 apparent transmission speeds of the signal between it and the seed stations are compatible with vseed. Finally, the linearity of the full event is tested using as reference value the azimuth angle of the seed φseed and the same tolerance as in the previous cases. If the number of aligned stations is equal the total number of selected stations, the event is labeled as “completely aligned”. Second treatment of isolated stations After selecting the set of candidate stations, the algorithm for rejecting isolated stations is applied again. The aim is to avoid accidental stations that were not removed by the previous algorithms and that can affect the angular reconstruction. One should note that the rejected stations may belong to the event, but it is better to be restrictive rather than to introduce “noisy” stations in the analysis, especially in the case of low multiplicity configurations. Figure 6.2: Footprint of an aligned event where the stations that remain after the selection procedure are too much separated in the array. It is also possible that after applying this algorithm the accepted configuration is not an event. An example is shown in Fig. 6.2 in which two of the selected stations (those at the two ends along the line of tanks) do not have at least 2 stations within 5000 m and they are rejected. As a consequence the final configuration does not pass any of the T3 trigger conditions. 6.1.2 Angular reconstruction After selecting the candidate stations, the next step is the angular reconstruction needed in our work to select inclined events. The reconstruction follows two different methods depending on the event configuration: 123 •Non-aligned configuration: non-aligned event with either an aligned or a non-aligned seed. •Aligned configuration: the event is “completely aligned”. Non-aligned configuration In the case of non-aligned events, the “standard” angular reconstruction can be used to determine the direction of the shower. Plane Fit: non-linear solution The direction of the shower axis is estimated from the start-time of the selected stations under the basic assumption that the shower front is a plane front moving with the speed of light along the shower axis. Thus, the time t(~ri) when the shower plane passes through a given position ~ri= (xi, yi, zi) on the ground is given by, ct(~ri) = ct0−~a ~ri(6.13) where t0is the time at which the impact point of the shower axis reaches ground and ~a = (u, v, w) is a unit vector in the forward direction of the shower axis. Assuming that the positions of the stations are given with no uncertainty and that the only source of uncertainty is that due to the uncertainty σiin the start-time (obtained from [102]), we can obtain the parameters (t0, u, v, w) by minimizing the squares of the differences between the measured start-time and the predicted times, given by Eq. (6.13), χ2= N X i=1 ti−t(~ri) σi2 (6.14) where Nis the number of selected stations. In order to get a good numerical precision, it is better to sum over quantities with small absolute values, so the positions and times of the stations are referred to the signal-weighted barycenter of the candidate stations, which is set as the origin. In reality the ground is not exactly a plane. The curvature of the Earth can be taken into account by projecting the zicoordinate of each station ionto a plane tangential to the ground located at the impact point of the shower axis, assuming again the signal-weighted barycenter of the event as the origin. The coordinate ziis shifted by δzi=−r2 i/2R(6.15) 124 where riis the distance of the station ifrom the origin and Ris the radius of curvature that is assumed constant and equal to radius of the Earth. After this correction, the following non-linear system is obtained, χ2= N X i=1 cti−ct0+xiu+yiv+ziw cσi2 (6.16) with the constraint u2+v2+w2= 1 and therefore 3 independent parameters: t0, u, v. The procedure to solve this system is iterative and converges to an unique solution if the stations are not all along the same straight line. The solution (u, v) corresponds to a physical direction (zenith and azimuth angles) if u2+v2≤1 using Eq. (6.5). The uncertainties are obtained propagating the uncertainties on the directional cosines into the angles (for details see [106]). We use this algorithm to reconstruct the shower direction when searching for neutrino candidates. A more elaborated angular reconstruction requires knowing the position of the shower core and fitting the start-times after correcting them with (for instance) a model of the time delay of muons described in [30]. The aim of our work is however to identify neutrino candidates without a previous knowledge of the core, because current algorithms designed for reconstruction of the core of nucleonic inclined showers are in principle not suitable to determine the core of a neutrino-induced shower. Aligned configuration The angular reconstruction in the case of an aligned configuration is just an estimate of the zenith and azimuthal angles of the shower. The procedure is simple and robust. For every pair of candidate stations (i,i+1) sorted by increasing time, the apparent transmission speed of the signal is calculated along the shower direction given by the azimuthal angle of the seed. The zenith angle corresponding to each pair of stations (θj) is calculated using Eq. 6.9. The corresponding azimuth angle (φj) is calculated as the angle subtended between the vector given by the position of stations (i, i+1) and the x-axis. The mean values of these angles corresponding to different pairs of stations give an estimate of the shower zenith and azimuth, θ=Pθj/σ2 θj P1/σ2 θj φ=Pφj/σ2 φj P1/σ2 φj (6.17) 125 The angular uncertainties σθj and σφj are obtained assuming that the positions of the stations have no uncertainty, and that the only source of uncertainty is that associated to the start-time (obtained from [102]). In general the reconstructed azimuthal angle is different from the azimuthal angle of the seed assumed to be given by the direction of the line of triggered tanks. Due to this the reconstructed zenith angle of an aligned configuration tends to be smaller than the actual zenith angle. 6.2 Characterization and identification of downgoing neutrino showers The first step in the study of the possibility of identifying down-going neutrino events in the background of ordinary proton and nuclei showers, is to characterize the neutrino-induced showers using Monte Carlo (MC) simulations. The interaction of a neutrino with atmospheric nuclei may result in a “purely hadronic” or in a “mixed” shower depending on the neutrino flavour, type of interaction (charged current or neutral current) and on the fraction of the energy of the neutrino carried by the secondary particles in the detector (Section 2.3.2). 6.2.1 Simulation of ν-like events At the present stage of this study, we have assumed that proton primaries interacting deep in the atmosphere produce showers equivalent to the hadronic showers resulting from NC interactions of neutrinos of all flavours, or CC interactions of νµor ντ(neglecting both the possible shower initiated by the µor the τ). The resulting hadronic shower is assumed to carry 20% of the neutrino energy. The validity of this approximation has been studied by comparing proton-induced showers with ν-induced showers in which the interaction of the neutrino is simulated with the Monte Carlo code HERWIG [103], and then the products of such interaction are propagated in CORSIKA [104]. For instance, in Fig. 6.3 we show the comparison of the signal map in the transverse plane of proton- and νµ−induced showers with θ= 80◦, injected at a slant depth measured from the ground ∆X= 910 g cm−2and for a proton energy or an energy carried by the hadronic shower Ep=Esh ≃1018 eV. The agreement between both maps is good at the ∼20% level, although this value depends on the distance to the shower core. This is confirmed in Fig. 6.4 where we show the comparison of the muon and electromagnetic contributions to the tank signal as a function of the distance from the shower axis 126 for the proton-induced and νµ−induced showers. From this plot, no significant difference is appreciable between hadronic showers induced by protons and νµs if the showers carry approximately the same energy. This result is in agreement with the detailed study performed in [105]. x (m) -1500 -1000 -500 0 500 1000 1500 y (m) -1500 -1000 -500 0 500 1000 1500 -1 -0.8 -0.6 -0.4 -0.2 -0 0.2 0.4 0.6 0.8 1 -2 X = 910 g cm∆ ° = 80θ = 1 EeV sh ) E ,cc µ ν +S p ) / 0.5 (S ,cc µ ν -S p (S Figure 6.3: Comparison between the signal map in the transverse plane of proton-induced showers and νµ−induced showers in CC interactions at θ= 80◦,∆X=910 g cm−2and with an energy going into the hadronic shower Esh ≃1018 eV. Each map is obtained as the average of 10 showers simulated with CORSIKA + S1000 USC code. The neutrino interaction is simulated with HERWIG. The coloured scale indicates the relative difference between both maps (Sp−Sνµ)/0.5(Sp+Sνµ). We have generated a library of proton showers using the shower propagation Monte Carlo code AIRES 2.6.0. and the hadronic interaction model QGSJET01. We used a 10−6thinning level that gives a good compromise between CPU time consumption per shower and artificial fluctuations due to the statistical sampling of particles. Showers were generated with energies ranging from 0.1 to 10 EeV, at different incident zenithal angles (from 75◦to 89◦) and injection points3chosen so that the slant atmospheric depth crossed by the shower from the injection point to the ground (∆X) is as large as 5000 g cm−2(measured along the shower axis). The simulations were performed in the conditions of the southern site of the Pierre Auger Observatory, neglecting the effect of the geomagnetic 3The injection point is the vertical depth of the first interaction point 127 1e-04 0.001 0.01 0.1 1 10 100 1000 10000 100000 1e+06 1 1.5 2 2.5 3 3.5 4 S [VEM] log10 [r / m] Esh = 1 EeV, θ = 80o, ∆X =910 g/cm2 SEM Proton SEM νµ Sµ Proton Sµ νµ Figure 6.4: Muon and electromagnetic contributions to the tank signal in VEM as a function of the distance from the shower axis in the shower plane for proton-induced showers and νµ−induced showers at θ= 80◦,∆X=910 g cm−2and with an energy going into the hadronic shower Esh ≃1018 eV. field. Although the magnetic deviations of the muons are very important for ordinary inclined showers, in the case of showers produced deep in the atmosphere the path lengths traveled by muons are in general not large enough to be significantly affected by the geomagnetic field. The simulation of the Surface Detector (SD) of the Pierre Auger Observatory has been performed using the standard modules of the Auger Offline Framework [94] using as input the ground particle files produced by the AIRES code. The response of the tank has been simulated with the Geant4 fast tank simulator [72]. The current Surface Detector central trigger configuration [98] has been applied to select the showers that would trigger the detector using the Central Trigger Simulator module. We have assumed an infinite array, so that the shower is always fully contained inside the array. In those cases in which we only needed the average signal in an Auger tank, we used the S1000 USC code to obtain it. This code provides a fast response to the electromagnetic and muonic components of the shower at the ground, although no time information of the signals can be obtained (see Chapter 4). 128 6.2.2 Signals produced by ν-showers at the ground The current picture of a shower induced by a neutrino interacting deep in the atmosphere close to the ground, is that of a “young” shower sharing many characteristics with vertical showers induced by hadrons. According to this image, the signals in all the triggered surface detectors should exhibit signatures of the presence of a significant electromagnetic component. However, our simulations have shown that this is not quite correct when the resulting shower from the neutrino interaction is an inclined purely hadronic shower or a mixed shower with a relatively large hadronic component. In fact in inclined showers, the azimuthal asymmetries on the time structure and signal size at the ground are very important, and our simulations show that in the case of a νinducing a deep inclined hadronic shower one should expect signals with EM characteristics only in the early part of the shower as we will show below. Azimuthal asymmetries in deep inclined showers In inclined showers, there is an azimuthal asymmetry in the signal due to the combination of several effects, the most important being the geometric effect, the longitudinal development effect and ground screening (for a description of these effects see Section 5.2). The geometric effect can be quite important for inclined nucleonic showers, but for ν-showers the most relevant effects for νidentification are the longitudinal development effect, and the ground screening effect as will be discussed in the following. In Fig. 6.5 we show a sketch of an inclined shower hitting two detectors at the same distance from the core at the ground. We also display three different planes that correspond to three different values of atmospheric depth crossed by the particles in the shower, and therefore three different stages of shower development. From this illustration it is evident that the tank in the early region is hit by a younger stage on the evolution of the shower than the tank in the late region. To quantify the difference between the three depths defined above for a shower of zenith angle θinitiated at ∆Xfrom the ground (in slant depth), we use a simple geometrical approach. The slant depth crossed by the shower core ∆X, and the difference with the corresponding depth for the early plane ∆XE(late plane ∆XL) is denoted by dE(dL), measured along the shower axis (see Fig. 6.5). dEcorresponds to an early tank (ζ= 0◦) at a distance rE from the core and dLcorresponds to a late tank (ζ= 180◦) at a distance rL from the core. They can be obtained as: 129 dE=rEsin θ dL=rLsin θ(6.18) Thus, the slant depth crossed by the early plane and late planes are: ∆XE= ∆X−dE ∆XL= ∆X+dL(6.19) The angles subtended between the shower axis and the paths from the injection point to the early and late tanks are given by: αE=θ−arctan tan θ−rE h(∆X) αL= arctan rLcos θcos θ h(∆XL)(6.20) where h(∆X) and h(∆XL) are the vertical heights as measured in meters corresponding to ∆Xand ∆XL, respectively. Therefore, the depth crossed by the early and late planes along the straight lines from the injection point to the tanks (not parallel to the shower axis) are: ∆X∗ E=∆XE cos αE ∆X∗ L=∆XL cos αL (6.21) It is important to remark that this simple geometrical model is only approximate but helps understanding the difference between the early and late regions of the shower. Their predictions work better for muons than for electrons, because e electrons do not travel in straight lines due to multiple scattering. As a result of the effect described above, for certain ranges of rand ∆X the early region corresponds to a younger stage in the shower development with a significant electromagnetic component, so most of the particles arriving at ground are electrons, positrons and gammas from the cascading processes. However the late region has to cross a much larger atmospheric depth, being in an older stage of evolution where the electromagnetic component becomes more attenuated and only muons and the electromagnetic halo component produced by muon decay and muon interactions arrive at the ground. In Fig. 6.6, we plot the depths ∆X∗ Eand ∆X∗ Lcrossed by the 130 Early plane Late plane Core plane θXL ∆ X∆ XE ∆ αE αL L d rE rL dE Xinj Figure 6.5: Geometrical model of an inclined shower reaching the ground. Three planes are displayed intersecting the ground plane, each one at different depths along the shower development: early (red), late (blue) and core (black) planes. The last one is also called shower plane. shower particles hitting an “early” tank and a “late” tank at r=4.5 km from the core, as a function of ∆X. From this figure, we can see that for instance the particles in a θ= 85◦shower initiated at ∆X∼1500 g cm−2, hit an “early” tank at a distance rE= 4.5 km after crossing ∆X∗ E≃1000 g cm−2 and a “late” tank, at the same distance rL= 4.5 km, after crossing ∆X∗ L≃ 1950 g cm−2. This means that for example if the shower is produced by a 10 EeV proton, we infer from Fig. 5.1 that the particles hitting the “early” tank correspond to a stage in which the electromagnetic component is at its maximum, while those hitting the “late” tank correspond to a stage in which the electromagnetic component is largely attenuated and the muonic component starts to dominate in the total signal. The additional atmosphere crossed by the late region with respect to the early one clearly depends on the distance from the core for a given θand ∆X as shown in Fig 6.7 where we show the difference ∆X∗ L−∆X∗ Eas a function of the distance from the core when the shower axis travels along ∆X= 1500 g cm−2. In this simple geometrical model, ∆X∗ L−∆X∗ Eincreases linearly with r, for example increasing from 640 g cm−2at r= 3 km to 960 g cm−2 131