ABSTRACT Title of Dissertation: A Search for Short Duration Very High Energy Emission from Gamma-Ray Bursts David Noyes, Doctor of Philosophy, 2005 Dissertation directed by: Professor Gregory W. Sullivan Department of Physics Milagro is a water-Cherenkov detector capable of observing air showers produced by gamma rays with primary energies of approximately 100 GeV and higher. The wide eld of view ( 2 sr) and high duty cycle (>90%) of Milagro make it ideal for searching for transient emission from gamma-ray bursts (GRBs). The median energy of photons detected by Milagro is a few TeV, but the e ective area is still relatively large at a few hundred GeV ( 50 m2 at 100 GeV). This results in a gamma- ray uence sensitivity comparable to previous satellite detectors at keV energies. Measurements have been made of GRB spectra up to a few tens of GeV with no sign of a cuto , however much is still unknown about the nature and existence of this Very High Energy (VHE) component. Additionally, gamma/gamma absorption from infrared background photons or from the optically thick region of the burst source complicate observations of this VHE component. However, many models predict VHE emission from GRBs through mechanisms such as synchrotron self-Compton processes. In the absence of a GRB localization provided by another instrument, the Milagro data is searched independently for VHE emission from GRBs. In 2.3 years of searching for bursts with durations ranging from 250 s to 40 s, no signi cant evidence was observed for VHE emission from GRBs. Models for di erent GRB parameters (such as redshift and isotropic energy distributions) are used to constrain the VHE spectrum of GRBs. A Search for Short Duration Very High Energy Emission from Gamma-Ray Bursts by David Noyes Dissertation submitted to the Faculty of the Graduate School of the University of Maryland, College Park in partial ful llment of the requirements for the degree of Doctor of Philosophy 2005 Advisory Committee: Professor Gregory W. Sullivan, Chairman/Advisor Professor Jordan A. Goodman Professor Theodore Jacobson Professor Abolhassan Jawahery Professor M. Coleman Miller Dedicated to: My parents ii ACKNOWLEDGEMENTS Six years seems like an awfully long time to get a few extra letters added to the end of your name. It would not have been possible without the help I received from many people. First of all I would like to thank my parents, Carl and Susan Noyes. The support and encouragement I received from them really meant a lot to me. I would also like to thank my cousin Samantha for always being so excited to see me come home. I would like to thank my advisor Greg Sullivan for his patience, guidance, and for giving me the opportunity to work on Milagro. I would also like to thank all my great collaborators on Milagro, the help I received from them was invaluable. Extra thanks goes to Andy Smith, without his help I don?t know how I would have carried out the work presented in this dissertation. More thanks to Andy and also Pablo for helping get this dissertation into an acceptable form. My time in graduate school was made in nitely more tolerable by all the great friends I had. I cannot summarize in this space all that I owe to iii Marilia for her support over the years. She has been one of the most important people in my life, and I am glad to have spent a part of it with her. Du san Tur can has been a great friend, o cemate, and roommate. He is one of the most generous people I know, and even though we often disagree, I have bene ted a lot from our discussions. Jonathan Ozik has really been a great friend. I have to thank him for lots of good conversation, for all the fun we had going to places like Wonderland, and for not selling out. I would like to thank Sarah Newport for all the tea we drank together while she was here, and for being a really good friend. I feel really lucky to have met Ellen Roche, and in the relatively short time that I have known her she has made me really happy. I would also like to thank Jenni zzle, Brendan, all the Baltimore peeps, Timmmmahh, the University of Maryland Food Coop, Milagro graduate students (especially Liz, Magda, Matt, Aous, Vlasis), and last but not least Jiji and Giogio for their valuable assistance in writing this dissertation. Figure 1: Jiji and Giogio. iv TABLE OF CONTENTS List of Tables viii List of Figures ix 1 Gamma-Ray Bursts 1 1.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1 1.2 Experimental History . . . . . . . . . . . . . . . . . . . . . . . . . . . 2 1.2.1 BATSE . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2 1.2.2 GRB Afterglows . . . . . . . . . . . . . . . . . . . . . . . . . 7 1.2.3 GRB-Supernova Connection . . . . . . . . . . . . . . . . . . . 10 1.2.4 Beaming . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 13 1.3 GRB Theory . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 14 1.3.1 The Compactness Problem . . . . . . . . . . . . . . . . . . . . 15 1.3.2 The Internal/External Shock Scenario . . . . . . . . . . . . . 16 1.3.3 Progenitors . . . . . . . . . . . . . . . . . . . . . . . . . . . . 20 1.4 VHE Emission from GRBs . . . . . . . . . . . . . . . . . . . . . . . . 21 1.4.1 Observations of VHE Emission from GRBs . . . . . . . . . . . 21 1.4.2 Theories Predicting VHE Emission from GRBs . . . . . . . . 23 2 Milagro: A VHE Gamma-Ray Observatory 29 2.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 29 v 2.2 Cosmic Rays, Extensive Air Showers, and Cherenkov Radiation . . . . . . . . . . . . . . . . . . . . . . . . 30 2.3 Detector Description . . . . . . . . . . . . . . . . . . . . . . . . . . . 33 2.3.1 The Pond . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 35 2.3.2 The Outriggers . . . . . . . . . . . . . . . . . . . . . . . . . . 36 2.3.3 Tube Repairs . . . . . . . . . . . . . . . . . . . . . . . . . . . 38 2.4 The Electronics and DAQ System . . . . . . . . . . . . . . . . . . . . 39 2.4.1 Signal Extraction . . . . . . . . . . . . . . . . . . . . . . . . . 39 2.4.2 Triggering . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 42 2.4.3 The Data Aquisition System . . . . . . . . . . . . . . . . . . . 46 2.5 Calibration . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 49 2.5.1 ADC Calibration . . . . . . . . . . . . . . . . . . . . . . . . . 50 2.6 Event Reconstruction . . . . . . . . . . . . . . . . . . . . . . . . . . . 51 2.6.1 Core Fitting . . . . . . . . . . . . . . . . . . . . . . . . . . . . 52 2.6.2 Angle Fitting . . . . . . . . . . . . . . . . . . . . . . . . . . . 54 2.7 Optimal Bin Size . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 56 2.7.1 Point Source Search Techniques . . . . . . . . . . . . . . . . . 59 2.8 Background Rejection . . . . . . . . . . . . . . . . . . . . . . . . . . 60 3 The Gamma-Ray Burst Search 64 3.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 64 3.2 Search Description . . . . . . . . . . . . . . . . . . . . . . . . . . . . 65 3.2.1 Overview of the search . . . . . . . . . . . . . . . . . . . . . . 65 3.2.2 Search Durations . . . . . . . . . . . . . . . . . . . . . . . . . 66 3.2.3 Oversampling . . . . . . . . . . . . . . . . . . . . . . . . . . . 67 3.2.4 Background calculation . . . . . . . . . . . . . . . . . . . . . . 71 3.2.5 Data Sample and Cuts . . . . . . . . . . . . . . . . . . . . . . 73 3.2.6 Monitoring the Results of the Search . . . . . . . . . . . . . . 80 vi 3.2.7 The Estimated Annual Rate . . . . . . . . . . . . . . . . . . . 82 3.2.8 Generating GRB alerts . . . . . . . . . . . . . . . . . . . . . . 86 3.3 Search Results . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 88 3.3.1 Summary of Lowest Probability Locations . . . . . . . . . . . 89 4 Sensitivity to GRBs and Results 93 4.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 93 4.2 De nitions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 94 4.3 Calculating the Number of Photons from a Source: The E ective Area 95 4.4 Modi cations to the Spectrum . . . . . . . . . . . . . . . . . . . . . . 99 4.4.1 Extra-Galactic Background Light . . . . . . . . . . . . . . . . 100 4.4.2 Absorption in the Source Region . . . . . . . . . . . . . . . . 100 4.5 Model Independent Sensitivity Calculations . . . . . . . . . . . . . . 101 4.5.1 Spectral Dependence . . . . . . . . . . . . . . . . . . . . . . . 105 4.5.2 Dependence on EBL Model . . . . . . . . . . . . . . . . . . . 106 4.6 Limits on GRB Properties . . . . . . . . . . . . . . . . . . . . . . . . 106 4.6.1 Poisson Fluctuations . . . . . . . . . . . . . . . . . . . . . . . 107 4.6.2 The Simulations . . . . . . . . . . . . . . . . . . . . . . . . . . 108 4.6.3 Case 1: Measured Redshift and Eiso Distributions . . . . . . . 111 4.6.4 Case 2: Measured Redshift and BATSE Fluence Distributions 116 4.6.5 Case 3: The Lag-Luminosity Relationship . . . . . . . . . . . 117 4.6.6 Case 4: Collapsar/Binary Merger Scenarios . . . . . . . . . . . 119 4.6.7 Case 5: Broken Power Luminosity Distributions . . . . . . . . 122 5 Conclusions 126 5.1 Summary . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 126 5.2 Future Directions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 128 vii LIST OF TABLES 3.1 One block of output from a log le, reformated here for ease of reading. The length of the block is 600 s because it was the rst block the search ran on that day, and this is the minimum amount of time required for building the background map. The columns from left to right are the duration, start time of the block, end time of the block, number of actual searches performed, the smallest probability found in the time window, the signal, background, zenith angle, RA, , time for this probability. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 81 3.2 Simple form of a GRB alert . . . . . . . . . . . . . . . . . . . . . . . 88 3.3 Locations found with estimated annual rate less than 5. . . . . . . . . 90 4.1 Dependence of upper limit on the source spectrum and EBL model. . 115 5.1 Summary of results from di erent models. The column on the left gives the redshift and isotropic energy model and the column on the right gives the 90% CL upper limit on the energy ratio. . . . . . . . . 127 viii LIST OF FIGURES 1 Jiji and Giogio. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . iv 1.1 A collection of BATSE light curves [7]. . . . . . . . . . . . . . . . . . 3 1.2 The bimodal duration distribution observed in BATSE GRBs [8]. . . 4 1.3 The location and uence of BATSE GRBs. There is no clustering in the Galactic plane [10]. . . . . . . . . . . . . . . . . . . . . . . . . . . 6 1.4 Log(N)-Log(S) distribution of BATSE bursts for the 1024 ms trigger time scale [11]. Drawn on the data is a line of slope -3/2. The data clearly deviates from this line. . . . . . . . . . . . . . . . . . . . . . . 7 1.5 GRB970228, the rst GRB afterglow [13]. On the left is the X-ray emission hours a short time after the detection of the prompt emission, on the right is the same region a few days later. . . . . . . . . . . . . 8 1.6 Measured redshift distribution of GRBs. . . . . . . . . . . . . . . . . 9 1.7 Optical afterglow light curve for GRB990510 [20]. . . . . . . . . . . . 11 1.8 Comparison of spectra for two GRB associated with a SN [24]. . . . . 12 1.9 Energy distributions assuming isotropic (top) and beamed emission (bottom) [25]. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 14 1.10 Diagram showing the stages of the reball model [29]. . . . . . . . . . 17 1.11 Progenitor Models: Left: Simulations of the local density structure at di erent stages of the relativistic jet that forms the GRB in the collapsar model [32]. Right: An illustration of a binary merger [33]. . 20 ix 1.12 High energy observations of GRB940217 [35]. The 18 GeV photon was observed by EGRET 90 minutes after the prompt emission. . . . 23 1.13 A distinct high energy component in GRB941017 [37]. E2dN=dE is plotted for 5 di erent time intervals: time since BATSE trigger, from top to bottom, -18-14s, 14-47s, 47-80s, 80-113s, 113-211s. . . . . . . . 24 1.14 Sky map in the region surrounding GRB970417a [41]. . . . . . . . . . 25 1.15 Time averaged theoretical spectra of GRBs in the reball model[46]. The three di erent curves are for di erent fractions of thermal en- ergy carried by the magnetic eld B: B = 0:33 (solid), B = 10 2 (dashed), B = 10 4 (dotted). For all three curves = 600, e = 10 0:5, t = 10 3s. . . . . . . . . . . . . . . . . . . . . . . . . . . . . 26 1.16 Survival probability for high energy photons interacting with IR back- ground photons [49]. Based on the model by Bullock and Primack. . 28 2.1 The cosmic-ray spectrum at the top of the atmosphere [50]. . . . . . . 31 2.2 Diagram of an EAS initiated by a gamma ray. . . . . . . . . . . . . . 32 2.3 An aerial view of Milagro. . . . . . . . . . . . . . . . . . . . . . . . . 34 2.4 Dimensions of the Milagro pond. . . . . . . . . . . . . . . . . . . . . 35 2.5 An inside view of Milagro. The AS layer tubes can be seen attached to the PVC grid crossings, while the MU layer tubes are tied between grid crossings. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 37 2.6 Distribution of shower particles that reach the ground. From the sim- ulation of an air shower induced by a 1.5 TeV gamma ray. The black points indicate the shower particles, the red square indicates the pond dimensions, and the blue point is the location of the shower core. . . 38 2.7 A diver removing a tube to be repaired. . . . . . . . . . . . . . . . . . 40 2.8 Drawing of Milagro?s electronics and DAQ system. . . . . . . . . . . . 41 x 2.9 Trigger rate versus tube trigger threshold for all events and for events for which an angle could not be reconstructed. . . . . . . . . . . . . . 44 2.10 Rise time distributions [52]. Top: rise time of Milagro data events which could be t. Middle: rise time of of Milagro data events which could not be t. Bottom: rise time for simulated gamma-ray induced air showers. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 45 2.11 NAS distribution, showing the relative number of events coming from each trigger. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 47 2.12 Increase in sensitivity between new and old trigger. On the left is shown the number of simulated gamma-ray induced air showers that satisfy the trigger conditions as a function of energy. The new trigger is in red, the old (55 tube threshold) trigger is in blue. The ratio of the two is plotted on the right. . . . . . . . . . . . . . . . . . . . . . . 48 2.13 Distribution of ADC counts for one tube. The peak on the left is the pedestal, the one on the right is the 1 PE peak. . . . . . . . . . . . . 51 2.14 Reconstructed core positions, employing outriggers for the t, for 50,000 events. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 53 2.15 Error in the reconstructed core position. The true core position from the simulation is compared to the tted core position. . . . . . . . . . 54 2.16 Event display showing the shower plane. . . . . . . . . . . . . . . . . 55 2.17 The eo distribution. The median of eo=2 is approximately the an- gular resolution, which is about 0:78 . . . . . . . . . . . . . . . . . . . 56 2.18 The delangle distribution for 100-500 GeV gamma rays. . . . . . . . . 57 xi 2.19 Determination of the optimal bin size. From simulations, the number of gamma rays contained within the bin is determined as a function of bin size. The value of the delangle cut that optimizes the ratio of the number of gamma rays to the square root of the area of the bin determines the optimal bin size. . . . . . . . . . . . . . . . . . . . . . 58 2.20 Signi cance map in the region of the Crab Nebula. The position of the Crab is at RA = 05h35m, Declination = 22h01m. . . . . . . . . . 60 2.21 Di erences between simulated gamma ray (bottom) and proton (top) induced air showers. . . . . . . . . . . . . . . . . . . . . . . . . . . . 61 2.22 X2 distribution for simulated gamma-ray showers (blue), simulated proton showers (black), and data (red). The horizontal axis is the X2 value. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 62 3.1 The number of background events in the search bin divided by the duration as a function of zenith angle. For short durations, the average number of events from any position on the sky is small. . . . . . . . . 67 3.2 Illustration of how a time period from t0 to tend is searched at a du- ration tdur. For each time window of duration tdur the entire sky is searched by comparing the number of signal and background events in each search bin (in grey). A new search bin is centered on each of the 0:2 0:2 ne bins (in black). . . . . . . . . . . . . . . . . . . . 69 3.3 The e ciency to nd a low probability vs the probability threshold for performing the ne over the coarse search. At probability of 10 4 the search e ciency is 99.1%. This is the value used in the search. . 70 3.4 The all-sky rate for MJD 2944. Note that the rate is not constant over the entire day. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 71 xii 3.5 The e ciency map for a 600 s block of data. The background for any bin is obtained by multiplying the e ciency by either the all-sky rate for the time window (duration>10 s), or the all-sky rate in a 10 s interval centered on the start time (duration<10 s). . . . . . . . . . . 72 3.6 The all-sky rate for MJD 2425. Top: The rate is averaged over 1 s. Note the long period where the rate has dropped lower and with many large uctuations, followed by a period where there was no data, then normal running again. Middle: The same plot zoomed in. Time periods containing no data are clearly visible. Bottom: The rate is averaged over 100 s bins. The period containing the time gaps is less visible. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 75 3.7 The distribution of the time between events for MJD 3147. This is for a day with no signi cant time gaps. . . . . . . . . . . . . . . . . . . . 77 3.8 The average all-sky rate vs Julian day for the entire data sample. . . 78 3.9 Fraction of day searched for each day in the data sample. The code x that eliminated the time gap problem was installed on MJD 2808. 79 3.10 Signal, background, and probability distributions from one day of searching (MJD 2702). Durations of 0.001 s, 1 s, and 10 s. . . . . . . 82 3.11 The number of trials for each duration searched. In blue is the number of trials calculated assuming that each search is uncorrelated. In red is the number of trials calculated from the estimated annual rate. In black is the actual number of searches performed, and is only depen- dent on the search algorithm as discussed earlier. . . . . . . . . . . . 85 3.12 The probability at which an annual rate of 365 is expected. Pthresh = 365=Ntrials. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 86 3.13 Probability distributions for entire data set (search durations 0.000251 s - 0.158 s). . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 91 xiii 3.14 Probability distributions for entire data set ( search durations 0.251 s - 39.8 s). . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 92 4.1 Triggered energy distribution of simulated gamma ray induced air showers. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 96 4.2 Median triggered energy versus zenith angle. . . . . . . . . . . . . . . 97 4.3 The e ective area for three di erent zenith angle ranges. . . . . . . . 98 4.4 The e ective area versus zenith angle for an E 2:4 spectrum. . . . . . 99 4.5 Modi cation to an E2:4 spectrum at a redshift of 0.2 due to -pair creation on the EBL. . . . . . . . . . . . . . . . . . . . . . . . . . . . 101 4.6 Pre-trials probability required for a 5 detection. . . . . . . . . . . . 102 4.7 Number of gamma rays dectectable from a GRB as a function of the burst redshift for di erent isotropic energies (curved lines). The hor- izontal lines indicate the minimum number of gamma rays required for a 5 detection. The intersection of the the curved lines with the horizontal lines indicates the maximum observable redshift for that combination of duration and isotropic energy. . . . . . . . . . . . . . 104 4.8 Maximum observable redshift as a function of isotropic energies for di erent durations. . . . . . . . . . . . . . . . . . . . . . . . . . . . . 105 4.9 The e ect of the spectral index on the maximum observable redshift. 106 4.10 The e ect of the EBL model on the maximum observable redshift. . . 107 4.11 BATSE T90 distribution. . . . . . . . . . . . . . . . . . . . . . . . . . 109 4.12 The measured redshift distribution. A t using Eqn.4.14 is drawn on top of the distribution. . . . . . . . . . . . . . . . . . . . . . . . . . . 110 4.13 The implied isotropic energy distribution for GRBs with measured redshifts. No beaming assumed. . . . . . . . . . . . . . . . . . . . . . 111 4.14 Fluence vs. Duration for BATSE GRBs. The line drawn on both gures is an estimate of the minimum uence vs duration for BATSE. 112 xiv 4.15 Cumulative redshift distribution using the t (black) and using the actual distribution (red). The two distributions agree reasonably well. 113 4.16 Results of a simulation based on the measured distributions. . . . . . 114 4.17 BATSE uence distribution for short (T90 < 2s) and long (T90 > 2s) duration bursts. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 117 4.18 Results of a simulation based on the measured redshift distribution and using the measured BATSE uence (which is correlated with du- ration). . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 118 4.19 Redshift and Eiso distributions from the lag-luminosity relationship. . 119 4.20 Results of the simulation using the lag-luminosity relation. . . . . . . 120 4.21 Models of the Star Formation Rate. . . . . . . . . . . . . . . . . . . . 121 4.22 Redshift distribution of NS-NS and NS-BH mergers. . . . . . . . . . . 122 4.23 Results of the simulation based on binary merger and SFR redshift distributions. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 123 4.24 Results of the simulation based on a broken power isotropic energy distribution and the SFR1 redshift distribution. . . . . . . . . . . . . 125 5.1 Upper limit on the fraction of GRBs with a VHE component for the di erent models considered in Chapter 4. . . . . . . . . . . . . . . . . 129 5.2 Maximum observable redshift as a function of isotropic energy for miniHAWC compared to Milagro. The redshift reach of miniHAWC is much greater than that of Milagro. . . . . . . . . . . . . . . . . . . 130 5.3 Results of the simulation based on a t to the measured distributions using the miniHAWC e ective areas. . . . . . . . . . . . . . . . . . . 131 xv Chapter 1 Gamma-Ray Bursts 1.1 Introduction Gamma-ray bursts (GRBs) are one of the most intriguing objects in the known universe. They are intense bursts of photons (typically in the keV to MeV energy range) lasting for as little as a few milliseconds and up to many minutes. Their origin is unknown and they are distributed isotropically across the sky, showing no repetition or clustering in any region of space. The light curves of GRBs show great diversity, with variability timescales as small as a few milliseconds. The typically observed GRB uxes ( 10 6 erg/s/cm2) and redshifts (z 1) imply an enormous amount of energy ( 1053 erg, assuming isotropic emission) being emitted in a short period of time. This is equivalent to converting approximately a tenth of a solar mass entirely into gamma rays in this short time period. The nature of the gamma- ray emission implies ultra-relativistic motion, with bulk Lorentz factors of at least several hundred (the largest bulk Lorentz factors to be observed). While the majority of GRB observations have occurred below a few MeV, no cut-o has been observed in the GRB spectrum and observations have been made which suggest that the spectrum extends to at least tens of GeV. While a great deal has been learned about GRBs since their discovery over three decades ago, much about them still remains 1 a mystery. Particularly, we still do not know the radiation mechanism responsible for the prompt GRB emission, nor do we know what is the underlying source of the GRB. In this chapter, the history and current state of the eld of GRB research will be reviewed. 1.2 Experimental History GRBs were discovered in 1969 by the Vela satellites [1, 2], whose purpose was to monitor for violations of the nuclear test ban treaty in space. The detection was initially classi ed, and the paper describing the observations was not published until 1973 [3]. From this time until the launch of the Compton Gamma Ray Observatory (CGRO) in 1991 there were a number of small experiments (e.g. KONUS [4] and the Gamma Burst Detector aboard the Pioneer Venus Orbiter [5]) which observed GRBs. However, due to poor angular resolution and the small sample of bursts they obtained, no signi cant progress was made in determining with certainty any global GRB properties. 1.2.1 BATSE Much of what we currently know about GRBs has come from the Burst and Transient Source Experiment (BATSE) [6], which ew aboard CGRO from 1991- 2000. BATSE served as an all sky monitor for CGRO and was designed for the detection of a large number of bursts. Consisting of 8 individual 20 inch diameter NaI scintillation crystals (the Large Area Detectors), BATSE was sensitive to energies between 20 keV and 2 MeV and had an angular resolution of about 4a0 . Over the period of its ight, BATSE detected 2704 GRBs [11]. Before BATSE, the largest catalogs of GRBs contained a few hundred bursts. From this large sample of bursts, a number of general characteristics were dis- 2 Figure 1.1: A collection of BATSE light curves [7]. cerned. The light curves of the bursts are vastly diverse (Fig.1.1), showing great variability on time scales ranging from milliseconds to seconds and with durations lasting from several milliseconds to many minutes.. The burst durations, de ned as the time in which 90% of the photons arrive (T90), indicate the existence of two classes of bursts. The duration distribution (shown in Fig. 1.2) is clearly bimodal, with one population above and one below about 2 seconds. It is not yet clear what gives rise to the two populations of bursts, but it has been suggested that they may 3 Figure 1.2: The bimodal duration distribution observed in BATSE GRBs [8]. originate from two di erent classes of progenitor. Spectral properties of GRBs were also investigated by BATSE, with photon ener- gies measured in 4 energy channels (20-50 keV, 50-100 keV, 100-300 keV, and >300 keV). The measured burst spectra are non-thermal, and may be t by a smoothly broken power law (the Band function [9]) parameterized by the low energy spectral index, , the high energy spectral index, , and break energy E0. N(E) = 8 >>< >>: A1E e E=E0 E < ( )E0 A2E E ( )E0; (1.1) where A1 and A2 set the amplitude and are such that N(E) varies smoothly at the break energy. Typical values of the spectral indices are 1 and 2:25. The energy at which a plot of E2N(E) peaks is given by Epeak = (2 )E0. Epeak varies 4 between 100 keV and 1 MeV, with a narrow peak around 200 keV. On average, the short duration bursts have harder spectra (more ux at high energies). The tted GRB spectra are time averaged; they generally evolve in time, from hard to soft, with Epeak shifting to lower energies. Prior to the BATSE measurements there was great debate about the distance to GRBs. If GRBs were of cosmological origin, the uences measured by BATSE would imply isotropic energy releases on the order of 1053erg. This may be compared to a typical supernova which radiates 1051erg over a period of months. Given this extremely high energy released in a much shorter time scale, many favored theories which predicted a galactic origin, which implied a much lower energy release. In this case, however, one would expect the angular distribution of bursts to follow the Galactic plane (i.e. they would cluster about the plane of the Galaxy). Fig.1.3 shows a map, in galactic coordinates, of the BATSE GRBs. It is clear from the gure that there is no clustering in the Galactic plane region (the angular distribution is isotropic), suggesting a non-Galactic origin. Further indirect evidence of a cosmological origin was seen in the number-intensity distribution of the BATSE data. This is done by comparing the number of bursts as a function of intensity, assuming Euclidean geometry (as would be expected for a galactic population). The ux of a GRB at a distance d is S = E4 d2; (1.2) where E is the total energy emitted by the burst. Bursts with S > Smin will be detected out to a distance dmax = s E 4 Smin (1.3) The volume of space occupied by bursts with S > Smin is V = 43 d3max: (1.4) 5 Figure 1.3: The location and uence of BATSE GRBs. There is no clustering in the Galactic plane [10]. Assuming a constant number density (n) and a xed total energy (this assumption may be relaxed somewhat), the number of bursts in this volume is N = nV = n43 ( E4 S min )3=2 S 3=2min : (1.5) Therefore, if we plot, on a log-log plot, the number of bursts with ux greater than Smin, the data should lie on a line with a slope of -3/2. Fig.1.4 shows the peak ux distribution of the BATSE bursts. The data clearly deviate from the -3/2 slope. There appear to be fewer dim bursts than expected. If bursts are of cosomological origin, this deviation can be understood as being due to the expansion of the universe. 6 Figure 1.4: Log(N)-Log(S) distribution of BATSE bursts for the 1024 ms trigger time scale [11]. Drawn on the data is a line of slope -3/2. The data clearly deviates from this line. 1.2.2 GRB Afterglows Although there were a number of signs pointing toward a cosmological origin for GRBs, direct evidence was still lacking. Despite the fact that the BATSE data favored such an origin, there were still arguments proposing, for example, that GRBs were produced by objects in an extended Galactic halo. It was not until the obser- vations by the Italian-Dutch X-ray satellite BeppoSAX [12], launched in 1996, that it was certain that GRBs lie at cosmological distances. The satellite consisted of the Wide Field Camera (WFC), sensitive from 2 to 30 keV, and the Narrow Field Instrument (NFI), sensitive from 0.1 to 10 keV. Both instruments were capable of a few arc minute localizations. GRBs were detected by a crystal scintillator detector system that had a nearly all-sky view in the 100-600 keV band, but with poor angular 7 Figure 1.5: GRB970228, the rst GRB afterglow [13]. On the left is the X-ray emission hours a short time after the detection of the prompt emission, on the right is the same region a few days later. resolution. Once the burst was detected, the position was rst re ned by the WFC, and then the spacecraft was able to slew to the direction of the burst and observe it at lower energies with the NFI. This allowed the detection of the rst GRB afterglow. BeppoSAX detected the prompt gamma-ray emission of GRB970228 [13], and then re-pointed the spacecraft within 8 hours to observe a fading X-ray afterglow at the location of the burst (Fig.1.5). The precise localizations provided by BepposSAX combined with the rapid dis- tribution of the coordinates allowed follow-up observations at other wavelengths by ground-based observers. This led to the indenti cation of GRB host galaxies, and the measurement of burst redshifts. GRB970508 was the rst GRB to have a measured redshift (z = 0:835) [14]. The GRB Coordinates Network (GCN) [16] allows rapid distribution of GRB coordinates from the instruments which detect GRBs to those 8 redshift 0 0.5 1 1.5 2 2.5 3 3.5 4 4.5 5 # of GRBs 0 1 2 3 4 5 6 measured redshift distribution Figure 1.6: Measured redshift distribution of GRBs. interested in performing follow-up measurements. The standard procedure is that the prompt gamma-ray emission is detected by a satellite or potentially a ground- based instrument and a noti cation is sent over GCN. Then emission is searched for at other wavelengths by the same or di erent detectors. With the re ned position obtained from X-ray afterglows, optical follow-ups are possible, potentially allowing the identi cation of the GRB host galaxy and redshift. A number of redshifts have been measured in this manner, although not enough to construct a complete distribution ( 40 as of June 2005 [15]). The redshift distribution is shown in Fig.1.6. Most redshifts are obtained from spectroscopy of the GRB host galaxy, and a few have been obtained from absorption-line systems in the burst spectra. All of the afterglows that have been observed have been located inside a galaxy. Some empirical redshift estimators based on light curve variability or other burst properties have been constructed [17, 18], but are still too uncertain 9 to be used with con dence. Afterglows have been observed primarily at X-ray, radio, and optical frequencies. Most bursts for which a measurement was attempted, had a X-ray afterglow, however not all had an optical afterglow. About 50% of the bursts with an observed X-ray afterglow had no optical afterglows, despite deep and prompt searches for such emission. The reason for these optical non-detections is still unknown, but may be due to bursts at very high redshift, bursts which are much dimmer than average bursts, or dust extinction. One important point to note is that redshifts have only been measured for the long duration bursts, the redshift distribution of short duration bursts is completely unknown. The X-ray and optical uxes are observed to decay as a power laws in time, and may be written as Fo;x t : (1.6) For the optical band 1 and 0:7, while in the X-ray band, 0:9 and 1:4 [19]. The radio ux does not typically follow a power law decay, and has a more complicated time dependence. In addition to this, breaks are often observed in the afterglow light curves, where the power law decay steepens. The reasons for this will be discussed later in this Chapter. An example of the light curve for an optical afterglow is shown in Fig.1.7. 1.2.3 GRB-Supernova Connection While a large number of GRBs have been detected, a lot of information comes from the observations of a few bursts. One example of this is that now there are a number of signs that GRBs may be associated with a peculiar type of supernova. The rst indication of such an association came from GRB980425, for which no optical afterglow was observed, but whose error box contained SN1998bw, a peculiar type Ic supernova [21] [22]. A further sign was the late time re-brightening observed in some afterglow light curves. This could be due to a number of e ects, including the 10 Figure 1.7: Optical afterglow light curve for GRB990510 [20]. existence of a supernova that accompanies the GRB. Binary merger models have a di cult time explaining the existence of this re-brightening since such a system is not thought to be able to produce a supernova, and is expected to occur in a low density environment (where processes that could produce this bump, such as re ection or reprocessing of afterglow light would not be likely to occur). The rst observation of a bump in the afterglow light curve came from GRB980326 [23]. Observations made 22 and 28 days after the burst showed the afterglow to be 60 times brighter than expected from an extrapolation of the light curve. In addition to this, the spectrum was signi cantly more red, which would be expected from a supernova. 11 3000 4000 5000 6000 7000 Rest Wavelength (?) 22 21 20 19 -2.5 log(f ? ) + Constant SN 1998bw at day -7 SN 1998bw at max GRB 030329 Apr 8 UT continuum subtracted Figure 1.8: Comparison of spectra for two GRB associated with a SN [24]. GRB030329 provided rm evidence of a supernova associated with a burst [24]. The burst was very bright, was at a redshift of 0.1685, and had a very well observed afterglow. The spectra started out as a power law typical of burst afterglows, then after 7.67 days the spectra began to show broad peaks typical of a supernova. As time went on the burst afterglow continued to fade, while the supernova spectrum continued to rise. When the spectrum of the supernova is compared to that of SN1998bw, they are seen to be very similar (Fig.1.8). In this case it is clear that the burst is related to the supernova. Whether this is true for all bursts or all long duration bursts remains to be seen. One problem in answering this question is the di culty in measuring SNe spectra at large redshifts (z > 1), while the peak of the GRB redshift distribution is around z = 1. 12 1.2.4 Beaming There are a number of reasons to believe that the burst emission is beamed rather than isotropic. When discussing beaming of the GRB emission, there are two di erent beaming angles to consider. First, the material emitted by the source may be beamed into a cone with a xed opening angle, as is the case in many astrophysical objects. Second, due to the relativistic nature of the burst, any emission will be beamed into a cone of opening angle 1= . This relativistic beaming de nitely exists, and given the large Lorentz factors of 100, initial beaming angles are expected to be quite small. Two immediate consequences of beaming are a reduction in the implied energy release by a factor of 1 cos( jet) 2jet, and an increase in the GRB event rate by the same factor (since this implies we are only seeing a small portion of GRBs). A steepening of the afterglow light curve is an indication of beaming, and the time of such a steepening is predicted by afterglow models. The intrinsic opening angle of the burst is jet, and the beaming angle due to relativistic motion is 0 / 1. If initially 0 < jet, to an external observer there is no indication of beaming. But as decreases, at some point 0 > jet, and the light curve will steepen since we are observing less radiation. Afterglow theory can be used to relate the observed break time to jet. This has been observed to have an interesting implication for the distribution of burst energies. From the observed break times, calculated beaming angles have been used to correct the implied energy releases from the isotropic case. This is shown in Fig.1.9, where it is seen that the isotropic energies have a spread of four orders of magnitude, while the beamed energy distribution is much more narrow, suggesting that GRBs have a standard energy reservoir. The mean beaming corrected energy in Fig.1.9 is 3 1051erg which is not much larger than typical supernova energies, and further suggests a relation between GRBs and SNe. In one study [25], the calculated beaming angles vary from 2.5a0 to 17a0 , with a mean of around 3.6a0 . The mean opening angle implies a GRB rate 500 times the observed rate. 13 Figure 1.9: Energy distributions assuming isotropic (top) and beamed emission (bot- tom) [25]. This result is still quite uncertain however. 1.3 GRB Theory Given that bursts lie at such great distances, the measured uxes imply energies on the order of 1051 1054 erg, assuming isotropic emission. A number of models have been proposed to account for this and other observed burst properties. From these models has emerged a generally accepted picture known as the reball shock model. 14 1.3.1 The Compactness Problem Before even considering how the observed radiation is produced, we immedi- ately run into a problem. We know that the source of the GRB must be very compact. From the variability of burst light curves we can infer a rough estimate of the size of the source region as r c t. Observed values of t of about 1 ms give an estimate of r 3 107cm. We also know that the source region is extremely dense. Given burst luminosities on the order of L 1052erg/s (implied by observed uences of 10 7 10 4ergs=cm2=s and measured distances of z 1), and photon energies on the order E 1 MeV, we see that the number of photons is enormous: N L E t 1055: (1.7) This is a problem because two photons (with energies E1 and E2) may pair produce by ! e+e if s E1E2(1 cos 12) 2 > mec 2: (1.8) This will clearly be satis ed by GRB photons. The optical depth for pair production at the source is given by T fpN4 r2 5 1012: (1.9) where T is the Thomson cross section, and fpN is the number of photons above the pair production threshold. This means that it is very unlikely that the photons will avoid creating pairs (the survival probability is e ). However, radiation produced in this way would result in a thermal spectrum, in contradiction of the observed power law [26, 19, 27]. The only known solution to this problem requires assuming the source is moving relativistically. If we assume that the source is moving toward us with a bulk Lorentz factor , there are two e ects which reduce the opacity. First, the photon energies are reduced by a factor of , how this a ects the opacity depends on the burst spectrum (how many were above the threshold for pair production). If we denote the spectrum 15 as I0E (valid from Emin to Emax), then the number of photons with energy greater than E is found by integrating this spectrum and then multiplying by the variability time and the luminosity distance to the burst [28]: fpN>E = 4 d2I0 tE +1=( 1) (1.10) This is the number that will go into the expression for the opacity. Now we can determine how many photons there are that could pair produce with a photon of some energy E1. The energy of the photon which would annihilate another of energy E1 is given by Ean = ( mec2)2=E1 (1.11) Substituting the value of NEan into the expression for the opacity we can see that it is modi ed from the value calculated assuming non-relativistic motion by a factor of 2( 1). The second e ect is that the size of the source is allowed to be larger by a factor of , since t ! t and so r c t. Combining both e ects ! 2 2 (1.12) Thus if E1 > Epeak then = , the high energy spectral index. And for typical values of , the opacity is reduced by 6:25. If we require < 1 we can use this to set a lower limit on . If, for example, one takes E1 = Emax, the maximum observed photon energies from a burst, then the above may be used to constrain to be greater than a few hundred for typical bursts. This makes high energy observations of GRBs interesting, since they may be used to see what range of Lorentz factors bursts may have. 1.3.2 The Internal/External Shock Scenario The general picture of the reball model is shown in Fig.1.10. There are considered to be four stages in the burst: 16 Figure 1.10: Diagram showing the stages of the reball model [29]. a0 A compact source produces a relativistic reball: a plasma of photons, e+e pairs, and a small quantity of baryons. a0 The reball expands to a distance at which it becomes optically thin. a0 The prompt emission is produced, converting the bulk kinetic energy into the observed radiation. a0 The afterglow emission at various wavelengths is produced as the reball slows down and interacts with its environment. A compact source (r 107cm) releases 1053erg in the form of photons, e pairs, and baryons in a short period of time (0.01-100s). This is the reball, where in general a reball refers to an large concentration of energy in a small volume of space with a relatively small number of baryons. The baryons may be produced by 17 the source, or may be present in the surrounding environment. In order to achieve relativistic velocities, it is assumed that the baryon fraction is small: Mbc2=E 1 (1.13) where E is the total energy of the reball. However, we must also have Mb 10 8M in order to convert most of the radiation into bulk kinetic energy as described below [26]. Initially the reball is at some temperature Ti, which if high enough will allow the photons to form pairs. As the reball expands and cools the pairs would annihilate, and eventually the opacity would be low enough for the photons to escape. However this cannot be the process producing the observed radiation, since this would result in a quasi-thermal spectrum as opposed to the observed power law. Therefore there must be some other process involved in producing the prompt emission. Instead, the following is postulated. If the initial temperature is high enough the radiation cannot escape the reball because of the high opacity due to pair production. In addition to this, the baryons will increase the opacity due to Thomson scattering o electrons associated with the baryons. Because the radiation cannot escape, the reball is assumed to expand adiabatically under its own pressure (p = e=3, where e is the photon energy density). Eqn.1.13 implies that most of the energy is in the form of radiation, and when this is the case, conservation laws imply that: T / r 1 (1.14) (r) / r (1.15) The reball cools as it expands and in the process is accelerated. As the reball accelerates, the baryon kinetic energy increases, and eventually is comparable to the total energy. At this point the above equations no longer hold. Instead we have T / r 2=3 (1.16) 18 (r) = constant (1.17) So the reball now coasts with a constant Lorentz factor. Eventually, the reball will have expanded to radius at which it becomes optically thin, allowing radiation to escape. This is the point where the prompt emission of the GRB is produced. Two mechanism have been considered to produce the prompt emission; internal shocks [30] and external shocks [31]. In external shocks, the radiation is produced in shocks formed when the reball interacts with a surrounding medium such as the interstellar medium (ISM). Internal shocks require a variable source, emitting shells of material with varying Lorentz factors. The radiation is produced in shocks formed when a faster shell emitted at a later time catches up to and interacts with a slower shell released earlier. The variability of burst light curves are accounted for in internal shocks by requiring an intermittent source (each peak in the light curve corresponds to an emission episode from the source). For external shocks, the variability is produced by inhomogeneities in the external medium. A number of arguments have been put forward to show that external shocks would have a di culty explaining the variability in this way. And so internal shocks are believed to be responsible for converting the energy of the reball into gamma-rays. However this is still not resolved, and is interesting since di erent shock models are consistent with di erent progenitor models. The radiation mechanism responsible for the prompt emission is still uncertain. Synchrotron radiation of relativistic shock-accelerated electrons is the currently fa- vored process, but there are a number of observations that are hard to explain if this is the case. For example, synchrotron emission puts a constraint on the al- lowed values of spectral indices which some GRBs seem to violate. A number of possible modi cations to synchrotron emission have been considered [19], but more observations are required to settle the matter. 19 -0.02 -0.01 0.00 0.01 0.02 x (1011cm) 0.00 0.02 0.04 0.06 0.08 y (10 11 cm) log ? -2.4 -0.4 1.6 3.6 5.6 (a) Model JA: 2.1 s -0.2 -0.1 0.0 0.1 0.2 x (1011cm) 0.0 0.2 0.4 0.6 0.8 y (10 11 cm) log ? -5.8 -3.0 -0.2 2.6 5.3 (b) Model JA: 7.2 s Figure 1.11: Progenitor Models: Left: Simulations of the local density structure at di erent stages of the relativistic jet that forms the GRB in the collapsar model [32]. Right: An illustration of a binary merger [33]. The nal stage of the GRB is the afterglow emission. After the prompt emission is produced the reball continues to expand, and eventually encounters some external material. Here, it is believed that external shocks between the reball and the external material produce the afterglow emission at various wavelengths. As opposed to the prompt emission, the afterglow emission is understood quite well. Synchrotron radiation is thought to be responsible for the afterglow emission. Models of the afterglow following these assuptions have been very sucessful in reproducing afterglow light curves. 1.3.3 Progenitors The nature of the burst progenitor is one of the greatest mysteries surrounding GRBs. The favored models take two forms. One associates the burst with the core collapse of a massive star (collapsar or hypernova models), the other with the merger of two compact objects; neutron star-neutron star (NS-NS) or neutron star-black hole (NS-BH) binaries. Both cases result in a few stellar mass BH (the mass of a black hole with a 107cm Schwarzschild radius) surrounded by a torus of material which 20 may then be accreted by the BH. Current observations favor the massive star origin at least for the long duration bursts. Observations of burst host galaxies indicate that GRBs occur in normal star forming galaxies. Additionally, bursts are found to occur near or slightly o set from the center of the galaxy, which is consistent with being located in star forming regions in galaxies. While for binary merger models, the large spatial velocities involved should, in many cases, carry the the binaries outside of the host galaxy, however this is not observed. However, since we have not measured redshifts or observed the host galaxies of short duration bursts, binary mergers are still considered a strong candidate for this class of bursts. Both progenitors are able to satisfy the energy requirements of the burst and power the reball. 1.4 VHE Emission from GRBs 1.4.1 Observations of VHE Emission from GRBs In gamma-ray astronomy, the term Very High Energy (VHE) is typically ap- plied to gamma-ray energies from 30 GeV to 30 TeV (here it is also applied to energies somewhat higher than 30 TeV). For the most part, GRBs have been ob- served at keV to MeV energies. Rather than being inherent to GRBs themselves, this is due primarily to the sensitve energy range of past GRB detectors: BATSE (20 keV - 2 MeV), BeppoSAX (100-600 keV), and HETE (6-400 keV). There have been detections at energies greater than a few MeV by a number of detectors. The Solar Maximum Mission satellite (SMM) was sensitive from 0.3-9 MeV. SMM de- tected 73 bursts over a 3:5 year period, over 60% of which had signi cant emission above 1 MeV [34]. The Energetic Gamma Ray Experiment Telescope (EGRET), one of the instruments on CGRO, was sensitive from 30 MeV to 30 GeV. EGRET detected photons above 30 MeV from 7 GRBs, 4 of which had emission above 100 21 MeV. The most famous of these was GRB940217 [35]. EGRET observed 2 photons with an energy of 3 GeV from this burst during the same time it was observed by BATSE. The earth then passed between EGRET and the burst, but 90 minutes later EGRET detected an 18 GeV photon from the burst (Fig.1.12). The measured EGRET spectra were consistent with being an extension of the measured BATSE spectra to higher energies, implying that there is no cuto in the GRB spectrum up these energies [36]. Recently, BATSE and EGRET data have been combined to obtain GRB spectra from 30 keV to 200 MeV [37]. Out of the 26 bursts analyzed all but one were consistent with an extension of the BATSE spectra to higher energies. However, GRB941017 was found to have a distinct high energy spectral component which appears to evolve independently of the low energy component (Fig.1.13). This component extends from a few MeV to greater than 200 MeV, contains 3 times as much energy as the lower energy component, and the high energy ux remains nearly constant throughout the observation, whereas the low energy ux decays by about 3 orders of magnitude. The fact that the ux is increasing with energy (the tted power law spectral index is 1) means that the peak of this component is at greater than 200 MeV. A number of observations have been attempted at energies of a few hundred GeV and higher. The Tibet air shower array performed searches for 10 TeV burst-like events coincident with BATSE bursts [38]. However, no evidence was observed for this type of emission. Rapid follow-up observations were performed by the Whipple air-Cherenkov telescope [39] and are being performed by Milagro (to be described in the next chapter) [40]. The follow-up observations of individual GRBs at VHE energies can provide useful constraints on VHE spectrum of GRBs, even when no emission is observed. Some of these attempts have claimed evidence of VHE emis- sion, but not with high signi cance. Evidence, at the 3 level, of VHE emission from GRB970417a was observed by Milagrito (the prototype for Milagro) [41]. In a 22 Figure 1.12: High energy observations of GRB940217 [35]. The 18 GeV photon was observed by EGRET 90 minutes after the prompt emission. little over a year of operating, there were 54 BATSE GRBs within the eld of view of Milagrito. The Milagrito data was then searched for an excess of events coincident in space and time with the BATSE GRBs. A signi cant excess was observed for one of these bursts (GRB970417a), which had a post-trials probability of being due to a random uctuation of the background of 1:5 10 3. In the absence of coordinates provided by other instruments, VHE emission from GRBs can be searched for inde- pendently by detectors sensistive to VHE photons. Milagro is capable of performing this type of search. Searches for transient emission on time scales of 40 s to 3 hours (long duration emission from GRBs) [42] and 2 hours and higher ( ares from Active Galactic Nuclei (AGN)) [43] have been carried out using the Milagro data. No evi- dence for long duration VHE emission from GRBs or aring AGN was observed in these searches (although steady emission from known AGN was observed [44]). 1.4.2 Theories Predicting VHE Emission from GRBs VHE photons from GRBs are expected to be produced by a number of processes which could occur in the relativistic reball model [19, 45, 46]. One of the most likely 23 Figure 1.13: A distinct high energy component in GRB941017 [37]. E2dN=dE is plotted for 5 di erent time intervals: time since BATSE trigger, from top to bottom, -18-14s, 14-47s, 47-80s, 80-113s, 113-211s. mechanisms capable of producing these photons is inverse-Compton (IC) scattering. If synchrotron radiation of relativistic electrons in the reball is responsible for the keV/MeV GRB emission, then it is possible for the low energy synchrotron photons to undergo IC scattering o of the relativistic electrons to higher energies. This is referred to as sychrotron self-Compton (SSC) emission. In some cases this can produce more VHE emission than has been observed in the keV band [45, 46]. In [46] detailed numerical calculations of the GRB spectra predicted by the reball model are presented. These calculations take into account synchrotron, SSC as well 24 0 2 4 6 8 10 12 14 16 18 RA (deg) DEC (deg) 45 47.5 50 52.5 55 57.5 60 62.5 65 280 285 290 295 300 305 310 315 Figure 1.14: Sky map in the region surrounding GRB970417a [41]. as other processes. In this model the resulting spectrum depends on 6 parameters, 3 related to the underlying source: the total luminosity (L), the bulk Lorentz factor ( ), and the variability time ( t); and 3 related to the shock physics: the fraction of thermal energy carried by electrons ( e), the fraction of thermal energy carried by the magnetic eld ( B), and the power law index of the accelerated electron?s energy distribution (p). Due to the high density of photons in the reball, if the internal shocks occur at a small enough radius from the source, the reball will be optically thick, and high energy photons will not be able to escape the source. To account for this, models are used to calculate the optical depth for absorption of high energy photons in the source region [48]. Fig.1.15 shows the time averaged spectra for a few di erent parameters. In this gure, it is seen that the amount of VHE emission relative to keV/MeV depends on the ratio of B to e. 25 100 102 104 106 108 1010 1012 10?9 10?8 10?7 10?6 Obs. energy [eV] Flux, ? F ? [erg cm ?2 sec ?1 ] Figure 1.15: Time averaged theoretical spectra of GRBs in the reball model[46]. The three di erent curves are for di erent fractions of thermal energy carried by the magnetic eld B: B = 0:33 (solid), B = 10 2 (dashed), B = 10 4 (dotted). For all three curves = 600, e = 10 0:5, t = 10 3s. VHE photons are also produced in models where GRBs are the source of ultra- high energy cosmic rays (UHECRs). In the reball model, the bulk of the total kinetic energy is eventually carried by hadrons. When the internal or external shocks occur, this kinetic energy is converted into internal energy carried by both electrons and protons. These shocks may be capable of accelerating protons up to very high energies (as high as 1020 eV [47]). The high energy protons may then produce high energy photons via sychrotron radiation. Additionally, the protons could interact with surrounding nucleons (in the external shock model), producing 0 particles 26 which may then produces high energy photons. These are just a few of the models that have been considered, and it is clear that more observations at high energies are required in order to learn more about the radiation processes occuring in burst environments. One di culty with observing VHE emission from GRBs is the attenuation due to pair production ( ! e+e ). As mentioned above, high energy photons produced in the burst environment can interact with the lower energy photons in the burst environment. In addition to this, high energy photons which escape the source may still undergo the same process by interacting with the infrared (IR) photons which are a part of the Extra-Galactic Background Light (EBL). This will be discussed in more detail in Chapter 4. As an example, Fig.1.16 shows the survival probability (e (E;z), where (E; z) is the optical depth as a function of photon energy and redshift) for interacting with IR photons for a particular model of the density of IR photons [49]. It is seen that at a redshift of z = 1, photons at energies greater than 100 GeV are highly attenuated, whereas at smaller redshifts, more high energy photons survive. This implies that, in order to observe VHE emission from a GRB, it must be either relatively nearby or release a large amount of energy. 27 log10(E) (GeV) 1 1.5 2 2.5 3 3.5 4 4.5 5 (E,z) ? - e 0 0.2 0.4 0.6 0.8 1 Survival Probability z=0.02 z=0.03 z=0.04 z=0.1 z=0.2 z=0.3 z=0.5 z=1.0 Figure 1.16: Survival probability for high energy photons interacting with IR back- ground photons [49]. Based on the model by Bullock and Primack. 28 Chapter 2 Milagro: A VHE Gamma-Ray Observatory 2.1 Introduction In order to search for VHE emission from GRBs or other astrophysical sources a detector must take into account a number of experimental facts. The typical ux from a source of VHE gamma rays is rather small. For a VHE source such as the Crab Nebula, the integral ux above 1 TeV is 1:75 10 11=cm2=s (a little more than 1 photon/day in an area of 100 m2). This makes it necessary to have a detector with a large area. For this reason, ground-based detectors have an advantage over satellite-based detectors, since it becomes prohibitively expensive to put a large area detector in space. On the other hand, the atmosphere is optically thick to high energy photons. As a result, ground-based detectors can only observe the secondary particles created in the atmosphere. Finally, there is a large ux of high energy cosmic rays at the earth. Since cosmic rays are charged particles they are bent in external magnetic elds on their way to the earth, and form an isotropic background. It is necessary to accurately characterize, and if possible reject, this background when searching for VHE photon emission. For GRBs in particular, the emission is transient in nature, which makes it useful to have a detector with a wide eld of view (FOV) and high 29 duty cycle, allowing a large region of the sky to be observed a large fraction of the time. Milagro is a ground-based water Cherenkov detector which has a number of these features, making it ideal for searching for VHE emission from GRBs. 2.2 Cosmic Rays, Extensive Air Showers, and Cherenkov Radiation Cosmic rays, consisting of protons and heavier nuclei, are constantly striking the atmosphere. The cosmic-ray rate varies with energy, with a rate of approximately 1 particle per square meter per second above 1 TeV. The cosmic-ray spectrum extends from a few hundred MeV to as high as a few hundred EeV (1020 eV), and is shown in Fig.2.1 [50]. The sources of cosmic rays are still not completely known. This is due to the fact that cosmic rays consist of charged particles, which are bent in interstellar magnetic elds on their way to the earth. Therefore the direction from which they arrive at the earth does not point back to the direction of their source. Sources of cosmic rays include the sun (at very low energies), and various galactic and extra-galactic sources such as Supernova Remnants (SNRs), GRBs, and Active Galactic Nuclei. However, it is not clear whether these sources can account for the entire spectrum of cosmic rays, particularly at very high energies. In addition to cosmic rays, gamma rays from various astrophysical sources are also striking the atmosphere, although at a much lower rate. Depending on the energy of the cosmic ray or gamma ray, it will lose energy in the atmosphere via a number of di erent possible radiation mechanisms [51]. For gamma rays with energies above 80 MeV the dominant energy loss mechanism is pair production. The electron/positron pairs that are created form more high energy photons via bremsstrahlung radiation, which in turn form more pairs, and so on. This cascade of pairs and photons is referred to as an an Extensive Air Shower (EAS). The creation 30 Figure 2.1: The cosmic-ray spectrum at the top of the atmosphere [50]. of new particles continues until the average energy of the shower particles goes below a critical energy ( 80 MeV). Below this energy, ionization losses begin to dominate over bremsstrahlung for the electrons, while Compton scattering and photoelectric absorption begin to dominate over pair production for the gamma rays. This point in the shower development is referred to as shower maximum since this is when the shower contains the greatest number of particles. After shower maximum the cascade continues towards ground level, however the number of particles decreases rapidly due to ionization losses, and there may be very few or no particles that reach the 31 Figure 2.2: Diagram of an EAS initiated by a gamma ray. ground. This makes it important for ground-based detectors to be located at high altitude. A cosmic ray striking the atmosphere follows a similar, but more complicated, pro- cess. This is because the primary cosmic ray produces a hadronic cascade, including neutral and charged pions. The neutral pions decay rather quickly ( = 1:78 10 16 s) to two gamma rays, which then produce electro-magnetic cascades as described above. The charged pions can decay to muons and neutrinos ( = 2:55 10 8 s). The low energy muons can further decay to electrons and positrons, while the high energy muons typically reach ground level. Thus, the cosmic-ray induced EAS con- tains a mixture of electrons, high energy hadrons, photons, and high energy muons. This is in contrast to gamma-ray induced air showers which only contain the electro- magnetic component. One way in which these air shower particles are detected is through the Cherenkov light they produce when passing through di erent media. A particle that enters a medium while travelling faster than the speed of light in that medium emits Cherenkov radiation. This radiation is emitted as a light cone with a xed angle 32 with respect to the direction of the particle?s motion (cos( ) = c=vn, where n is the index of refraction in the medium, v is the particle velocity, a c is the speed of light). Atmospheric-Cherenkov telescopes (ACTs) use telescopes on the ground to observe the ash of Cherenkov light that is produced by shower particles in the atmosphere. As will be described next, Milagro is a water-Cherenkov detector, where the shower particles enter a large water reservoir in which they produce Cherenkov radiation. An advantage of this method comes from the fact that the index of refraction of water is much greater than that of air. This results in a Cherenkov light cone with a large opening angle (41 ), and leads to a greater number of Cherenkov photons being produced. The large Cherenkov angle makes it easy to collect light from almost all of the particles that enter the water. 2.3 Detector Description Milagro is a water-Cherenkov detector consisting of a six million gallon arti cial pond sealed with a light tight cover and instrumented with 723 photomultiplier tubes (PMTs) arranged in two layers. In addition to the pond, Milagro is surrounded by a sparse \outrigger" array of 175 individual Cherenkov counters covering an area of 40,000 m2. Milagro is located in the Jemez Mountains, about 40 miles west of Los Alamos, NM. The detector is at an altitude of 2650 m (8692 ft), latitude 35 52? 45" and longitude 106 40? 37". Fig.2.3 shows an aerial view of Milagro. In this photo, the pond is covered with snow, the black circular objects are the outriggers, the Pond Utility Building (the PUB) and the counting house are the buildings immediately to the left of pond. The PUB contains the water pumps and ltering equipment, fans, patch panels for the high voltage (HV) cables, and the signal cables. Pumps bring water in for ltering and return it to the pond at a recirculation rate of about 200 gallons per 33 Figure 2.3: An aerial view of Milagro. minute. Water taken in from the pond goes through a set of lters including a charcoal lter, a 10 m lter, a 1 m lter, a carbon lter, a UV lamp, and a 0.2 m lter. In order to keep external light out of the detector, it is covered with a 1 mm-thick polypropylene cover, the top of which was painted with highly re ective roo ng paint in order to reduce the temperature inside the pond. While the cover is normally kept in contact with the water by adjustable straps surrounding the pond, for maintenance operations it is necessary to in ate the cover. The vents and fans required for this are housed in the PUB. When fully in ated the internal pressure keeps the cover approximately 20 feet from the surface of the water, and an array of straps serve to keep it stable. Once the cover is in ated, people are able to enter the detector and conduct maintenance operations. A single cable from each PMT carrying both the HV and the PMT signal enter 34 Figure 2.4: Dimensions of the Milagro pond. the PUB and connect to a set of patch panels. The cables then go from the patch panels to run underground to the counting house. The counting house holds all the triggering and data collection electronics, as well as the computers for data reconstruction and storage. Since Milagro is located in an area that has one of the highest rates of lightning strikes in the United States, a Faraday cage was set up in the area surrounding the pond, PUB, and counting house. 2.3.1 The Pond The pond that is used for the central Milagro detector was originally used as a holding pond for geothermal energy experiments, and can hold 6 million gallons of water. The pond is 8m deep and measures 60m 80m at the surface, sloping down to the bottom where it measures 30m 50m (see Fig.2.4). The top layer of tubes, referred to as the air shower (AS) layer, consists of 450 PMTs at 1.5 m below the pond surface. The bottom layer of tubes, referred to as the muon (MU) layer, consists of 273 PMTs at 6 m below the pond surface. The buoyant tubes are tied to a weighted PVC grid that lies at the bottom of the pond. The AS layer tubes are tied to the grid crossings, while the MU layer tubes are tied to the centers of the PVC. This can be seen in Fig.2.5, which was taken with the cover in ated during one of the tube repair operations. The AS layer is used primarily for triggering and event reconstruction. The MU layer is used primarily for background 35 rejection but may also be used in event reconstruction. The PMTs are 20 cm in diameter and are made by Hamamatsu (model R5912SEL). They are encased in a PVC housing to protect the base from the water, and a sin- gle RG-59 coaxial cable carries the high voltage to the tube and the PMT signal from the tube. The connector required careful attention, since any corrosion or leaks would allow water in, causing the PMT base to short circuit. The original connector used when all of the PMTs were installed turned out to have a higher than expected failure rate when in water. This was the cause of many PMT failures and led to the development of an improved connector which had reduced strain and was sealed in heat shrink and glue. The tubes are surrounded by a collar of material (a \ba e") to block out light travelling horizontally and to increase the collection area of each PMT. The ba es are visible in Fig.2.5. They were originally made of anodized aluminum (for increased re ectivity) on the inside with black polypropylene on the outside, but due to corrosion of the aluminum, new ba es are being installed. The new ba es consist of white polypropylene on the inside and black polypropylene on the outside. 2.3.2 The Outriggers The 40,000 m2 area surrounding the Milagro pond is covered by a sparse array of \outriggers". An outrigger is an individual Cherenkov counter which consists of a 1500 gallon tank of water, measuring 2.4 m in diameter and 1 m high. The tanks are lined with Tyvek (to re ect light created in the tank) and are instrumented with a single PMT facing the bottom of the tank. The outrigger array allows a more accurate determination of the location of the shower core, which is important in the angular reconstruction (as will be described later). Since Milagro can trigger on events with cores distant from the pond, it is helpful to sample the shower away from the pond. Fig.2.6 shows the simulated 36 Figure 2.5: An inside view of Milagro. The AS layer tubes can be seen attached to the PVC grid crossings, while the MU layer tubes are tied between grid crossings. distribution of shower particles that reach the ground from an air shower induced by a 1.5 TeV gamma ray. As can be seen from the gure, determination of the shower core would be facilitated by greater sampling of the shower particles away from the pond. In addition to core location, the outriggers can be used in the angular reconstruction, in the estimation of the primary particle energy, and in background rejection. By covering a 40,000 m2 area with 175 individual Cherenkov counters, greater sampling of the shower is achieved. The outrigger installation began in the summer of 2000. Due to a large re in the Jemez Mountains that year, contractors were not available for most of the summer. Instead many people working on Milagro went to Los Alamos to help with the installation. This included digging the ditches which would carry the outrigger cables, laying the cables in protective PVC tubes, untangling massive lengths of RG-59 cables, setting up the area where the outrigger tanks would sit (which included clearing away any trees, rocks, etc in the area to make a level surface), and attaching connectors to the cables. Who knew that there 37 x (m) -500 -400 -300 -200 -100 0 100 200 300 400 500 y (m) -500 -400 -300 -200 -100 0 100 200 300 400 500 Distribution of Shower Particles at the Ground (E = 1.5 TeV) Figure 2.6: Distribution of shower particles that reach the ground. From the simula- tion of an air shower induced by a 1.5 TeV gamma ray. The black points indicate the shower particles, the red square indicates the pond dimensions, and the blue point is the location of the shower core. would be so much physical labor involved in earning a Ph.D in particle astrophysics! The complete installation, calibration, and integration into the data reconstruction of the outriggers was completed in 2003. 2.3.3 Tube Repairs Due to the higher than expected failure rate of the original connectors used on the PMTs (which allowed water into the PMT base), or failure of components in the PMT base, a yearly tube repair has been necessary. After the rst few years of running approximately 70 tubes died each year. After the switch to the new connector, this number dropped by about half. The tube repair requires divers to 38 enter the pond and disconnect the dead tubes from the kevlar string holding it to the PVC grid. The dead tube is removed by lowering a weighted container, which is attached by rope to an empty plastic milk jug, to the area where the dead tube is located. The empty milk jug is buoyant enough to oat at the surface while still being attached to the weighted container. This is then connected to the string which holds the PMT to the PVC grid, allowing the tube to then be disconnected from the grid without having the rope fall to the bottom. A small raft follows the divers to collect the tubes as they oat to the surface. Once the tubes are repaired or replaced, they are returned to the pond and reattached to their ropes by divers. I was fortunate enough to be a diver for the last 2 repair operations, and it was one of the more interesting experiences I had while working on Milagro. Fig.2.7 shows a diver (not me) during one of the tube repairs. 2.4 The Electronics and DAQ System 2.4.1 Signal Extraction Fig.2.8 gives a rough outline of the electronics and data acquisition (DAQ) system of Milagro. Each PMT connects to a custom 16 channel front-end board (FEB). The FEB reads in each PMT signal, and distributes the high voltage to each tube. The PMT signals are processed on the FEB and then go to the \digital boards", where timing and pulse height information is prepared for digitization. For each PMT signal, the arrival time and charge must be determined. The most straightforward manner to do this would be to employ analog-to-digital converters (ADCs), which would directly measure the charge in each tube. However, at the time that Milagro was built, speed issues did not make it possible to use ADCs. Instead, the charge is measured indirectly by employing the time over threshold (TOT) technique. 39 Figure 2.7: A diver removing a tube to be repaired. In the FEB, each PMT signal is split o and sent to high gain ( x 7) and low gain ( x 1) ampli ers. These ampli ed signals are then sent to a pair of discriminators with di erent photo-electron (PE) thresholds. The output from the high-gain ampli er is sent to a low threshold discriminator with a 1/4 PE threshold. The output from the low-gain ampli er is split in two, with one part going to a high threshold discriminator with a 5 PE threshold, and the other going to an output which can optionally be connected to an external ADC for calibration purposes. Whenever the PMT pulse crosses either of the low or high discriminator thresh- olds, an \edge" is generated. This is illustrated in the top left corner of Fig.2.8. For a relatively small PMT pulse, which crosses the low threshold, but not the high threshold, two edges are generated. The rising edge is created when the pulse goes 40 Figure 2.8: Drawing of Milagro?s electronics and DAQ system. 41 above the threshold (a) and the falling edge is created when the pulse goes below the threshold (b). In a similar manner, for a pulse which crosses both the high and the low thresholds, four edges are generated. The time spent over threshold can then be calibrated to the charge. This will be will be described in more detail later. Both sets of discriminator pulses are sent to the digital boards, which combine the high and low threshold pulses, and do additional pre-processing. The edges are then digitized in LeCroy FASTBUS time-to-digital converters (TDCs), which can record up to 16 edges per event with a 0.5 ns resolution. A FASTBUS Latch connected to a GPS clock encodes the common stop time for each event. 2.4.2 Triggering At any moment signals are registered in individual PMTs due to Cherenkov light from air shower particles and from other random particles. In order to select air shower events and trigger the detector, the basic approach is to require some minimum number of tubes to be hit (the tube threshold) within a certain time window. For this purpose, in addition to generating edges the digital boards produce a single 25 mV, 200 ns pulse for each hit AS layer PMT. The trigger pulses are then continuously added together producing the analog sum signal which can be used to determine the number of tubes hit within a xed time window. As the tube threshold is lowered, the rate of detected events increases. From the time Milagro began taking data (January 1999) until March 19, 2002, a simple discriminator threshold proportional to the number of tubes hit in the AS layer was used. The number of tubes required for a trigger varied between 50 and 70 tubes hit within a 200 ns time window. The number of tubes required was set by the maximum data rate that the DAQ system could handle. This was designed to be about 2000 Hz. Lowering the trigger threshold would greatly increase the number of low energy 42 (hundreds of GeV) showers collected, and therefore increase Milagro?s sensitivity at these energies. This in turn improves the chances of observing high energy emission from GRBs. This is due to the fact that GRBs are at cosmological distances, and high energy photons above a few hundred GeV experience absorption on the IR background, with the absorption becoming more signi cant at higher redshifts (see Fig.1.16 in Chapter 1). This makes it desirable to increase the sensitivity to the photons which are not attenuated by the IR background. Therefore, if a way could be found to lower the trigger threshold, while not increasing the rate, it would be very bene cial. It turns out that the increase in the trigger rate that occurs when the trigger threshold is lowered is due mainly to single muon events, which produce enough light to trigger the detector, but cannot be t to a shower plane (how the events are t will be described later). In Fig.2.9 the trigger rate is plotted for a number of di erent trigger thresholds. Shown in the gure is the rate for all events, and the rate for events for which an angle could not be reconstructed (not t). One can see from the gure that as the threshold is lowered, the increase in rate becomes predominately due to events which cannot be t. However, from studies of simulated air showers in Milagro, it was known that gamma-ray events could be reconstructed well with as few as 20 PMTs hit. A high angle muon which travels across the pond nearly horizontally produces light which travels at a speed c/n, while the shower particles are travelling at nearly c. Thus the light produced by high angle muons will arrive over a longer time period compared with that produced by shower particles. It was found that a cut on the rise time of the event allowed the elimination of events with a small number of hit PMTs which could not be t. The rise time of an event is de ned as the time it takes for a prede ned fraction of the tubes to be hit. For our purposes, it is taken to be the di erence in time between when 12.5% of the PMTs in the event have been 43 Threshold (PMTs)10 15 20 25 30 35 40 45 50 55 60 Ungated Data Rate (Hz) 103 104 105 All data Not fit data Ungated Trigger Rate Figure 2.9: Trigger rate versus tube trigger threshold for all events and for events for which an angle could not be reconstructed. hit, and when 88.5% of the PMTs in the event have been hit. This is illustrated in the bottom right of Fig.2.8. Fig.2.10 shows the rise time distributions for data events which were t, data events which were not t, and simulated gamma-ray showers [52]. Random hits were included in the simulated gamma ray showers since these could lengthen the rise time of an event. As can be seen in the gure, the data events which are not t have longer rise times. 44 Risetime(ns)0 50 100 150 200 250 300 350 400 0 500 1000 1500 2000 2500 3000 3500 Nent = 176864 Mean = 116.1 RMS = 46.4 Fit Risetime Risetime(ns)0 50 100 150 200 250 300 350 400 0 1000 2000 3000 4000 5000 6000 Nent = 280762 Mean = 161.9 RMS = 44.98 No Fit Risetime Risetime (ns)0 50 100 150 200 250 300 350 400 0 20 40 60 80 100 120 140 160 180 200 220 240 Nent = 7569 Mean = 80.87 RMS = 37.82 Gamma MC Figure 2.10: Rise time distributions [52]. Top: rise time of Milagro data events which could be t. Middle: rise time of of Milagro data events which could not be t. Bottom: rise time for simulated gamma-ray induced air showers. To take advantage of this fact, a custom VME (Versa Module Europa) trigger card was built that would allow the rise time of the event to play a role in the triggering. The trigger card works by reading in the analog sum, and if a set pre- trigger threshold is crossed, it calculates the rise time. The pre-trigger is set so as to minimize the detector dead time. The dead time will increase with a lower pre- 45 trigger threshold, since while calculating the rise time, the card is in a busy state and the rate of events passing the pre-trigger threshold increases. The card is fully programmable, allowing multiple trigger conditions and external triggering. The fol- lowing set of trigger conditions are the ones used when the card was rst installed, but have been adjusted slightly over time to keep the trigger rate at a desirable level (NAS is the number of tubes hit in the AS layer): a0 NAS > 20 & trise < 50 ns a0 NAS > 53 & trise < 87:5 ns a0 NAS > 74 These values were chosen so as to maximize the number of low energy showers, while still keeping the rate and dead time reasonable as well as retaining an unbiased sample of large showers. The pre-trigger threshold is set to 20 tubes, since this is the minimum number needed to ful ll one of the trigger conditions. Fig.2.11 shows the distribution of the number AS tubes hit in each event, for a large number of events. Fig.2.12 shows the relative increase in sensitivity between the new and old trigger. In the gure, on the left is shown the number of simulated gamma-ray induced air showers that satisfy the trigger conditions as a function of energy. The new trigger is in red, the old (55 tube threshold) trigger is in blue. The ratio of the two is plotted on the right. As can be seen, the simulations predict more than a factor of 4 improvement in sensitivity below 100 GeV. 2.4.3 The Data Aquisition System When any of the trigger conditions are satis ed a common stop is sent to the TDCs and the latch. The digitized data is then read into a dual port memory which sits in the VME crate. This data can then be read from the dual port memory 46 50 100 150 200 250 300 350 4000 500 1000 1500 2000 2500 nhitas nhitas>20 and trise<=50ns and nhitas>53 trise<=87.5ns nhitas>74 Figure 2.11: NAS distribution, showing the relative number of events coming from each trigger. by the DAQ computer, which is a standard PC running Linux. A number of these processes are controlled by the FASTBUS Smart Crate Controller (FSCC) which communicates with the DAQ PC over ethernet. The data is broken up into runs and sub runs. A single subrun contains about 5 minutes of data. The raw data consists of the TOT information for each tube, various event ID bits, and the time of the event from the GPS clock. A single subrun of raw data is approximately 600 MB in size (or 170 GB/day), making it extremely costly to save all of the raw data. Instead, the raw data is calibrated and reconstructed in real time, and only the reconstructed data is saved. The reconstructed data consists of the various properties of the event such as the reconstructed direction 47 Figure 2.12: Increase in sensitivity between new and old trigger. On the left is shown the number of simulated gamma-ray induced air showers that satisfy the trigger conditions as a function of energy. The new trigger is in red, the old (55 tube threshold) trigger is in blue. The ratio of the two is plotted on the right. and core position, the event time, the Modi ed Julian Date (MJD) and timing error information from the gps clock, but no information from individual tubes hit in the event. The size of a single subrun of reconstructed data is about 13 MB, or about 4 GB/day. The reconstructed data is stored on DLT tapes and transfered over the network to large redundant disk arrays in Maryland and Los Alamos. However, since the price of disk has become low relative to that of tape, the DLT tape storage is being replaced by disks. In addition to the reconstructed data, raw data is saved for selected sources, such as the Crab Nebula, the active galaxies Mrk 421 and Mrk 501 if they are in a aring state, the sun and the moon, and a small sample of the raw data as well. The raw data is also saved when there was a GRB in Milagro?s eld of view, which we are noti ed of by other detectors. In order to speed up the reconstruction process, a farm of Linux \worker" PCs are utilized. Each block of data that is read in is sent over a socket connection for 48 one of the workers to reconstruct. Once the data block is reconstructed, it is passed back to the DAQ computer where it builds up the subrun. This process makes it possible to have more complicated reconstruction algorithms by simply adding more computers. 2.5 Calibration As was described above, the timing and charge information is determined using the TOT method. In order to convert the edge times into relative times of hit tubes and the number of PEs, the PMT response must be calibrated. This is done primarily with a laser calibration system. This consists of a pulsed laser, a lter wheel, and optical bers that carry the light to a system of di using balls placed throughout the pond. The laser res a set number of pulses of xed amplitude through a lter wheel. The light from the laser goes through optical ber which connects to an optical switch, allowing the light to be sent to any of the thirty laser balls in the pond. A laser ball is simply the optical ber with a sphere of epoxy at the end which di uses the light out isotropically from the ber. They are kept in place by oating PVC platforms which are tied by string to the PVC grid at the bottom of the pond. The lter wheel allows the intensity of the light sent to the pond to be varied. Laser calibration data is taken periodically to produce new calibration constants. In particular, new calibrations are always performed after repairs, since the repaired tubes may have slightly di erent responses and some tubes may have been replaced altogether. One part of the calibration process involves timing corrections. In general, the start time of a PMT pulse with a large pulse height will be earlier than that of one with a small pulse height. This e ect is referred to as electronic slewing, and is corrected for by measuring the change in start time as a function of TOT using the 49 laser calibration data. The variation of TOT is achieved by varying the position of the lter wheel, and hence the intensity of the light sent to the pond. In addition to this, corrections must be made for the di erences in travel times of the PMT pulses through the signal cables and electronics. In parallel with the timing calibrations, the calibration of the PMT charge is carried out. As discussed above, in the TOT method an edge is generated whenever the PMT signal crosses one of the high or low discriminator thresholds. The time spent above threshold is proportional to the logarithm of the number of photo- electrons. This relation is determined using the occupancy method, which is based on the fact that, at low light levels, the number of photon-electrons created by the PMT obeys a Poisson distribution. This produces a logarithmic relation between TOT and the number of PEs. The exact relation is determined by the laser calibration. 2.5.1 ADC Calibration In addition to the laser calibration, data is taken with an external ADC, allow- ing a direct conversion of TOT to charge. This can then be compared to the results obtained from the laser calibration. A single 16 channel LeCroy 1881M FASTBUS ADC was used for this purpose, and required the ADC to be manually cycled through all the channels. In addition to this, the trigger rate had to be descreased to a level that could be handled by the ADC ( 600 Hz). This is a rather laborious, time consuming procedure and therefore is not done regularly. Fig.2.13 shows the distribution of ADC counts for one PMT. The two peaks clearly visible on the plot are the pedestal (0 PEs) and the 1 PE peak. Measuring the separation of the two peaks gives the gain of the PMT. The TOT values of each event are also measured, and associating a value of TOT to each ADC channel allows conversion from TOT to PEs. PE calibrations using the ADC were found to be in good agreement with those from the laser data, allowing one to have con dence in 50 ADC counts 250 300 350 400 450 210 310 410 ADC: Tube 113 Figure 2.13: Distribution of ADC counts for one tube. The peak on the left is the pedestal, the one on the right is the 1 PE peak. the calibration method. However, at very large PE values the ADC saturates, so the ADC calibration becomes unreliable. This is another reason for having the laser system as the primary means of calibrating the detector. 2.6 Event Reconstruction The times and pulse heights of the tubes which were hit in an event allow the determination of the orientation of the shower plane, and thus the direction of the primary particle. As the EAS traverses the atmosphere it spreads out laterally, forming the shower front. The orientation of the shower front points back to the direction of the particle which initiated the EAS. The time at which the particles reach the ground will vary within the shower front. This relative timing is used to 51 determine the orientation of the shower plane. 2.6.1 Core Fitting To determine the direction of the primary particle which initiated an extensive air shower, the arrival times of the shower particles are t to a plane. However, the shower front does not actually form a plane, it forms a cone whose apex is located at the shower core. The arrival times can be adjusted to form a plane if the shower core is known. The size of this correction is 0.07 ns per meter from the core. This was determined by studying the angular resolution for simulated air showers as a function of the \curvature" correction. This makes it important to know the location of the shower core. Over the course of running Milagro, a number of di erent methods have been used to determine the location of the shower core. The shower core should be associated with the location of the largest number of PMTs hit and with the largest number of PEs. However, the problem is complicated due to the fact that the shower core may be located o the pond. This is in fact the case for most of the events detected by Milagro, which is why the outriggers are so important in determining the location of the core. The rst iteration of the core t was a simple center of mass tter (weighted by the number of PEs) which put all the cores on the Milagro pond. This was improved by using various methods to determine if the shower core was likely to be located o the pond. If the core was on the pond, the center of mass was used, otherwise the core was placed at an arbitrary distance (50 m) from the pond in the direction determined by the center of mass. After the outriggers were installed, it was possible to use them to determine if the core was located on or o the pond. For this purpose, the ratio of the number of outriggers hit to the number of pond tubes hit provides a reasonable measure of whether the core is on or o the pond. If the core is o the 52 X Core Position (cm) -10000 -5000 0 5000 10000 Y Core Position (cm) -10000 -5000 0 5000 10000 0 2 4 6 8 10 12 14 16 18 20 Core Position Figure 2.14: Reconstructed core positions, employing outriggers for the t, for 50,000 events. pond the center of mass of the outriggers is used to determine the core position. If the core is determined to be on the pond, the center of mass of the pond tubes is used. The current core tter performs a least squares t to a 2-D Gaussian using the AS layer tubes and the outriggers. The distribution of core locations found by the core tter is shown in Fig.2.14. The error in the core position, determined by comparison to simulated air showers, with known core positions, is shown in Fig.2.15. Once the core position is known, a correction is applied to adjust for the curvature of the shower plane. In addition to correcting for the curvature of the shower front, a correction must be made for the way in which the shower front is sampled. Since the time of each hit is recorded as the arrival time of the rst Cherenkov photon, the arrival times of tubes hit near the core will be slightly earlier than those hit further from the core. Once both of these corrections have been made, it is possible 53 Core Error (m) 0 50 100 150 200 250 3000 500 1000 1500 2000 2500 3000 3500 Error in Core Location for Simulated Gamma Rays Figure 2.15: Error in the reconstructed core position. The true core position from the simulation is compared to the tted core position. to reconstruct the angle of incidence of the primary particle. 2.6.2 Angle Fitting Fig.2.16 shows the PMT arrival times for a single event. The shower plane can be seen by eye. This plane is t using an iterative least squares t. The tubes are weighted by the number of PEs, and several iterations are made based on the number of PEs. The residuals of small PE hits are very non-Gaussian, so only larger PE hits are used in the rst iteration of the t. The residuals from the t are calculated, and if they are within a pre-set range they are retained for the second iteration, otherwise they are excluded. On the second iteration the large PE requirement is relaxed, allowing more tubes to come into the t. This process is repeated ve times. Using this procedure 97% of simulated gamma-ray showers with NAS > 20 are t. 54 Figure 2.16: Event display showing the shower plane. In the data 80% to 90% of events are able to be successfully t with this method. This reduction is due to the fact that some triggers are not due to air showers (e.g. horizontal muons) and cannot be t to a plane. A measure of the angular resolution is given by the eo variable. eo is de- termined by rst separating the AS tubes into a checkerboard pattern. Then the angular reconstruction is performed twice, once using only the \black" tubes, and then using only the \white" tubes. eo is the space angle di erence between these two reconstructed angles. This is plotted in Fig.2.17 for a sample of Milagro data. eo is expected to be equal to twice the angular resolution of the detector [53]. For the data shown in the gure, the median of eo=2 is about 0:78 . 55 (deg)eo? 0 1 2 3 4 5 6 7 8 9 100 200 400 600 800 1000 1200 eo? Figure 2.17: The eo distribution. The median of eo=2 is approximately the angular resolution, which is about 0:78 . 2.7 Optimal Bin Size Due to the nite angular resolution of the detector, a point source will get smeared out to some extent on the sky. The amount by which it is smeared out is determined by the point spread function of the detector. For Milagro, the point spread function is determined from simulation of the detector response to gamma-ray induced air showers, and is given by the distribution of the space angle di erence between the simulated primary direction and the reconstructed direction (the de- langle distribution). This is similar to eo, but since the true angle is known for simulated showers, delangle is a more appropriate measure. In Fig.2.18 the delangle distribution is shown for simulated 100-500 GeV gamma-ray induced air showers. Given the fact that a point source is smeared, we must decide how large of an 56 delangle (deg) 0 1 2 3 4 5 6 7 8 9 100 0.005 0.01 0.015 0.02 0.025 0.03 0.035 0.04 0.045 Delangle Distribution Figure 2.18: The delangle distribution for 100-500 GeV gamma rays. area on the sky to use to search for a signal. This xed area that is used to search for a signal is called the search bin. The larger one makes the search bin, the greater the fraction of signal contained within the search bin. However, if the search bin is made too large, the background will begin to wash out the signal. For this reason, in a binned search for a point source there will exist an optimal bin size. The optimal bin size will give the greatest signi cance detection of the source. For large numbers of expected background events, the optimal search bin is a circular bin of radius 1:58 . For a Gaussian angular resolution, a search bin of this size will on average contain 72% of the signal events from a point source [53]. For a small number of signal events the search bin should be slightly larger. For large numbers of signal events, the signi cance is simply the ratio of the signal to the square root of the background, and this is maximized to nd the optimal bin size. Since the number of background events is directly proportional to the area 57 delangle (deg) 0 0.5 1 1.5 2 2.5 3 3.5 4 (arbitrary units) bin A / ? N 0 0.2 0.4 0.6 0.8 1 1.2 1.4 1.6 1.8 2 Optimal Bin Size Figure 2.19: Determination of the optimal bin size. From simulations, the number of gamma rays contained within the bin is determined as a function of bin size. The value of the delangle cut that optimizes the ratio of the number of gamma rays to the square root of the area of the bin determines the optimal bin size. of the bin (Abin), it is only necessary to nd the number of gamma-rays (N ) that are t within a bin as a function of the bin size. Then N =pAbin is plotted as a function of bin size, the peak of which is the optimal bin size. As mentioned before, in searching for GRBs it is crucial to maximize the sensitivity for energies below a few hundreds of GeV. For this reason, we should nd the optimal bin size for low energy gamma-ray showers. In Fig.2.19 N =pAbin is plotted as a function of bin size for 100-500 GeV gamma rays. From the plot it is seen that the peak is rather broad and centered between 1:5 and 1:8 . This is for a radial bin, while for simplicity the search uses square bins. In going from a radial bin to a square bin of equal area, there is a small loss in sensitivity. A N detection with an optimal radial bin goes 58 to a 0.985N detection with an optimal square bin [53]. The search uses a 3 square bin, which has roughly the same area as a 1:7 radial bin, and is within the broad peak seen in Fig.2.19. The value of 1.7 was chosen because it is known that for small numbers of signal events the optimal bin is larger than 1.58. 2.7.1 Point Source Search Techniques The basic test that a detector is functioning properly is its ability to detect a known source. The most well studied TeV source is the Crab Nebula. A signi cance map of the region surrounding the Crab Nebula is shown in Fig.2.20. The location of maximum signi cance is within the angular resolution of Milagro, and is at 7.9 . The entire sky binned in Right Ascension (RA) and Declination ( ). The signi - cance map is made by rst adding up the number of events in each 0:1 0:1 =cos( ) ne RA and bin. This is the signal map. The appropriate number of ne bins are added together to get the number of signal events in the 2:1 2:1 =cos( ) bin centered on the Crab. This is compared to an estimate of the background. The number of events expected in an RA, bin is given by Nexp(RA; ) = Z Z E(HA; )R(t) (HA; RA; t)dtd ; (2.1) where E(HA; ) is the acceptance of the detector as a function of local coordinates, R(t) is the event rate of the detector, and (HA; RA; t) is equal to 1 if the HA, RA and sidereal time coincide with the bin for which Nexp is being calculated [54]. It is assumed that E(HA; ) is constant over 2 hours, and is calculated by simply counting the fraction of events that come from each position on the sky over a 2 hour period. The expected number of events is compared to the number observed, and the signi cance of the excess is calculated using the method of Li & Ma [55]. This is slightly di erent from the GRB search (which will be described in the next chapter) since a large number of events are being dealt with in this case. The detection of the 59 -92 -90 -88 -86 -84 -82 -80 -78 -76 -74 Declination (deg) 12 14 16 18 20 22 24 26 28 30 32 Significance -2 0 2 4 6 RA 05h00m05h20m05h40m06h00m Figure 2.20: Signi cance map in the region of the Crab Nebula. The position of the Crab is at RA = 05h35m, Declination = 22h01m. Crab Nebula is important in establishing that Milagro is functioning properly, and is capable of detecting VHE sources. 2.8 Background Rejection For each air shower event that is detected by Milagro, it is not known if the primary particle was a gamma ray or cosmic ray. However, we can attempt to identify 60 Figure 2.21: Di erences between simulated gamma ray (bottom) and proton (top) induced air showers. the primary particle by statistical means. Due to the fact that an EAS induced by a cosmic ray contains muons and hadrons, both of which are more penetrating in water than photons or electrons, the MU layer can be used to identify these kind of events. The muons and hadrons which reach the MU layer illuminate a relatively small number of tubes in a small region. Fig.2.21 shows the distribution of light in the MU layer PMTs for simulated gamma ray (bottom) and proton (top) induced air showers. As can be seen in the gure, the proton induced showers result in a clumpy distribution of light in the MU layer, while for gamma ray induced showers the distribution is relatively smooth. A simple method for distinguishing the two types of showers uses the X2 param- eter: X2 = NMU>2PE max (2.2) where NMU>2 is the number of tubes hit in the MU layer with more than 2 PEs, and PEmax is the maximum PE value in the MU layer. Distributions of this parameter for simulated gamma ray and proton induced air showers and from data events are 61 Figure 2.22: X2 distribution for simulated gamma-ray showers (blue), simulated proton showers (black), and data (red). The horizontal axis is the X2 value. shown in Fig.2.22. As can be seen in the gure, the data and proton simulations agree well and have smaller values of X2 than the gamma-ray simulations. Standard analyses require events to have an X2 value greater than 2.5. For events with greater than 60 tubes hit in the AS layer, more than 20 tubes used in the angular t, and which were t within a 1:2 bin, this cut rejects 90% of the simulated proton showers and 91.5% of the data, while retaining 51% of the gamma-ray induced showers. This results in a relative increase in sensitivity (Q-factor) of 1.6 and 1.7 when comparing the simulated gamma rays to simulated protons and data respectively. The e ectiveness of this cut depends on the energy of the event, and for low energy events, the Q-factor quickly drops below 1. The reason for this is that the electro-magnetic component of hadron induced showers with cores located o the 62 pond trigger the detector, and are indistinguishable from low energy gamma ray induced showers. This is important for the GRB search, since, as mentioned before, it is these low energy events that are the most important in searching for VHE emission from GRBs. Furthermore, when searching for emission at short durations (as is the case for the analysis to be presented in Chapter 3), it is desirable to retain as much of the signal as possible. This is due to the fact that for low background levels (as is the case at short durations), the sensitivity of the search is limited by the number of signal events. For these reasons, no background rejection is used in the GRB search described in the following chapters. 63 Chapter 3 The Gamma-Ray Burst Search 3.1 Introduction In Chapter 1, a number of observations of VHE emission from GRBs by past experiments and theoretical models predicting such emission were discussed. This is the motivation for having a detector capable of independently monitoring the entire sky for VHE emission from GRBs, and was one of the primary reasons for building Milagro. In the absence of GRB coordinates provided by other instruments, it is possible to independently search the Milagro data for VHE emission from GRBs. This is an important search to carry out, especially given the low rate of burst localizations that have been available since BATSE was decommissioned. Given its large eld of view ( 2 sr), a large number of GRBs are observable by Milagro. The all-sky rate of GRBs in the universe implied by BATSE observations is about 700 GRBs/year. This implies approximately 110 GRBs/year in Milagro?s eld of view. An observation at VHE energies would be the rst time that the prompt emission from a GRB was detected at energies other than the typical keV/MeV energies. As described in the previous section, Milagro can be used to gather information from gamma-ray and cosmic-ray induced air showers. For each air shower detected 64 by Milagro the direction and time of the primary particle can be reconstructed. The direction and time of each air shower is used to search the sky for GRBs. In this chapter the procedure carried out to search for VHE emission from GRBs is described, and the results from the search are presented. Note that in the following, the dates are often stated in Modi ed Julian Date (MJD) format, for convience the values given are MJD-50,000. 3.2 Search Description 3.2.1 Overview of the search Due to the large cosmic-ray ux at the earth, the bulk of the events detected by Milagro are air showers induced by cosmic rays. It is therefore necessary to search for emission from a GRB on top of this large cosmic-ray background. However, in contrast to doing a search for a known source, the start time, duration, and coordinates of a potential GRB are not known a priori making it necessary to search over the entire sky, at all times, and over multiple durations. The rst step in this process is to create a background map using a long period of time (compared to the search duration). The background map contains the number of events expected from any location on the sky in a time window speci ed by the start time and the duration. Next, a signal map is created which contains the number of events observed in the time window speci ed by the start time and the duration. Given these maps, the probability that the number of observed signal events was due to a random uctuation of the background is calculated for each point on the sky. The signature of a GRB would be a very low probability, inconsistent with a random uctuation of the background. 65 3.2.2 Search Durations The details of how the search is carried out are dictated by the experimental technique and the nature of the emission that one is hoping to observe. In addition to this, the search is run on the data as soon as it is collected to allow for prompt noti- cation when a GRB is observed. This makes it necessary to have a computationally e cient search algorithm. Gamma-ray bursts have been observed to have durations ranging from a few milliseconds to many minutes. This de nes the range of durations that the search should cover. For technical reasons (explained below), it makes sense to break up the search into separate searches covering di erent ranges of duration. If the amount by which an object on the sky appears to move is small compared to the search bin size, there is no need to track it while it moves. Instead, the sky may be assumed to be stationary, allowing the use of local coordinates speci ed by Hour Angle (HA) and Declination ( ). The earth rotates 0:2 in 48 seconds. This is small compared to the angular resolution of the detector ( 0:78 ) and so if durations shorter than this are being searched, it is reasonable to use local coordinates. Additionally, the < 48 s search may be broken up into two searches. For very short durations, there are typically very few events in the entire sky (see Fig.3.1). This makes it computationally ine cient to search the entire sky at short durations. Instead, computations are carried out only in the neighborhood of regions where at least 2 events were detected. For these reasons a search consisting of 27 logarithmically spaced durations ranging from 250 s to 40 s has been developed. Searches ranging from 40 s to 3 hours and 2 hours and higher have been described elsewhere [42, 43], with the latter focused more on searches for Active Galactic Nuclei (AGN). As will be discussed later, Milagro is more sensitive to emission from GRBs at short durations, so the 250 s to 40 s search is the most likely to observe VHE emission from a GRB. 66 zenith angle (deg) 0 5 10 15 20 25 30 35 40 45 background rate (evts/sec) 0 1 2 3 4 5 6 7 8 Background rate vs zenith angle (MJD 2667) Figure 3.1: The number of background events in the search bin divided by the duration as a function of zenith angle. For short durations, the average number of events from any position on the sky is small. 3.2.3 Oversampling In the absence of a localization by a satellite detector, the location, start time, and duration of a potential GRB are unknown, making it necessary to perform a search of the entire sky, at all times, and over multiple durations. To ensure that a potential signal is not missed, the search is done with overlapping space and time bins at each duration. Without this oversampling, a signal could easily be missed. As an extreme example, assume no oversampling is performed in time. It is conceivable that half of the signal could be contained in one time bin and the other half in the next time bin. Unless the signal is very strong, it may not be signi cant enough to be detected above the background. The spatial search is carried out by binning the sky in Hour Angle (HA) and 67 declination ( ) with 0:2 0:2 ne bins, each of which de nes the center of a 3 3 =cos( ) search bin. The appropriate number of ne bins are added together to form each search bin. In this way the entire sky is searched, with the center of adjacent search bins di ering by 0:2 . The beginning of the time bins shift by 10% of the search duration for each new bin. The time of the rst event de nes the start time of the rst time window, the end point of which is de ned by the duration of the search. The next time window begins at the initial start time plus 0.1 times the duration. For example, for a 1 s search, if the time of the rst event is t0, then the whole sky is searched between t0 and t0 + 1:0, then between t0 + 0:1 to t0 + 1:1, and so on. Fig.3.2 illustrates the operation of the search. For each time window (denoted as tw1, tw2, etc.) of duration tdur the entire sky is searched by comparing the number of signal and background events in each search bin (in grey). A new search bin is centered on each of the 0:2 0:2 ne bins (in black). The maximum number of spatial searches is xed by the number of search bins (the centers of which are o set by 0:2 ) that may be t within the region de ned by the zenith angle cut. With a zenith angle cut of 45 this turns out to be 212,461 spatial searches. However, it is unnecessary to search all of these bins, and in order to increase speed of the search algorithm, only a fraction of them are searched. If there were to be a signi cant excess over the background, it would show up in many neighboring spatial bins at a lower signi cance. Therefore, a coarser search is performed, where every third ne bin de nes the center of a search bin. If the Poisson probability in any bin is below a certain value (10 4), all the bins in that region are examined. This reduces the number of spatial searches by about a factor of 9. This does not decrease the sensitivity of the search to nding a GRB by a signi cant amount. In Fig.3.3 the e ciency of the search is plotted as a function of the probability value (P) used to determine the threshold for searching all the bins. 68 Figure 3.2: Illustration of how a time period from t0 to tend is searched at a duration tdur. For each time window of duration tdur the entire sky is searched by comparing the number of signal and background events in each search bin (in grey). A new search bin is centered on each of the 0:2 0:2 ne bins (in black). 69 log10(probability) -10 -8 -6 -4 -2 0 Efficiency 0 0.2 0.4 0.6 0.8 1 Efficiency of Search vs Coarse/Fine Probability Threshold Figure 3.3: The e ciency to nd a low probability vs the probability threshold for performing the ne over the coarse search. At probability of 10 4 the search e ciency is 99.1%. This is the value used in the search. The e ciency is then calculated by running the search on a small sample of Milagro data for di erent values of the threshold probability. The number of times that a probability less than 10 7 was found is computed for each probability threshold as well as for the case of searching every bin (probability threshold of 1). The e ciency is then computed as (P) = NP=1(< 10 7) NP (< 10 7) (3.1) At a probability of 10 4 the search e ciency is 99.1%. This means that 99.1% of the time, the algorithm will identify a signi cant excess even though only 1/9 of the bins are searched. 70 UT seconds 0 10000 20000 30000 40000 50000 60000 70000 80000 Rate (Hz) 1700 1750 1800 1850 1900 All-Sky Rate (MJD2944) Figure 3.4: The all-sky rate for MJD 2944. Note that the rate is not constant over the entire day. 3.2.4 Background calculation As was mentioned previously, the GRB emission must be searched for on top of a large, isotropic cosmic-ray backgroud. A correct estimation of the background is crucial to the success of the search. If the background is overestimated, a real signal could be missed. If the background is underestimated, a uctuation of the background could appear to be a real signal. The simplest approach to measuring the background would be to average the number of events over a long time period prior to the current time window being searched, but at the same location on the sky. However, the fact that the rate is not stable over the entire day may lead to an incorrect estimate of the background. In Fig.3.4, the all-sky rate (total number events in the entire sky) is plotted for an entire day. The plot shows that the average rate varies over the day, and so a 71 HA (deg) -150 -100 -50 0 50 100 150 (deg) ? -60 -40 -20 0 20 40 60 80 0 0.0005 0.001 0.0015 0.002 0.0025 0.003 0.0035 0.004 Figure 3.5: The e ciency map for a 600 s block of data. The background for any bin is obtained by multiplying the e ciency by either the all-sky rate for the time window (duration>10 s), or the all-sky rate in a 10 s interval centered on the start time (duration<10 s). background estimation method that is insensitive to this is used. This method takes advantage of the fact that, while the all-sky rate varies over time, the fraction of events coming from any particular place on the sky remains constant to a part in 10 4. This is due to the fact that the cosmic-ray rate at the top of the atmosphere is roughly constant, while the variation of the all-sky rate throughout the day is due to variations in atmospheric pressure. The fraction of events coming from any portion of the sky is calculated and stored in a map which we call the e ciency map. In Fig.3.5 the e ciency map is plotted for a 600 s block of data. The map is computed by counting the number of events in each search bin and normalizing by the total number of events in the map. The background for any 72 bin is then obtained by multiplying the e ciency by either the all-sky rate for the time window (duration>10 s), or the all-sky rate in a 10 s interval centered on the start time (duration<10 s). The need for these two methods is that, for durations less than 10 s, there are not enough events in the whole sky to obtain a measurement of the rate with high statistics. For a typical post-cuts event rate of about 1400 Hz, a 10 s interval gives an error on the rate less than 1%. For this reason it is important that the rate be constant at least over a period of 10 s. All of the above may be summarized in the following expression Nexp(HA; ) = (HA; )R(t); (3.2) where (HA; ) is the e ciency map, and R(t) is the rate computed according to the above prescription. 3.2.5 Data Sample and Cuts As was noted earlier, a much improved triggering system was installed on March 19, 2002 (MJD 2352). The new trigger provided an increase in sensitivity at the low energy end of Milagro?s response, and e ectively created two experiments, one before the installation of the new trigger, and one after the installation of the new trigger. Due to the much higher low energy sensitivity, I have chosen to analyze the data covering the period from MJD 2353 to MJD 3372 (March 19, 2002 to January 1, 2005). Although the time period covered by this search is 1020 days, the actual amount of data that was searched is less than this due to a number of cuts. The search algorithm begins by reading the reconstructed data into a large bu er, which is then used to construct the signal and background maps on which the search is carried out. The reconstructed data consists of the event time, the reconstructed primary direction in Right Ascension ( ) and Declination ( ), and a few parameters characterizing the 73 event, such as the number of tubes used in the angular t, the Julian date, and timing error information. From this, the data bu er is lled with the time, HA, , and of each event, provided it passes a set of data quality and event cuts. The rst set of cuts are based on the characteristics of the event. The simplest of these cuts only requires that the angle be reconstructed. Given that the event was t, the reconstructed zenith angle of the event is required to be less than 45 . This cut is used due to the fact that the gamma-ray sensitivity of Milagro (as will be described in the next chapter) greatly decreases with increasing zenith angle, making it very unlikely to detect a GRB at large zenith angles. Typically 85% to 95% of the events pass these cuts, and the average for the entire data set was 89%. As discussed in the previous chapter, no gamma/hadron separation cuts are used in this analysis. This is because of the fact that the currently developed gamma/hadron separation methods do not work e ectively on low energy (hundreds of GeV) showers. However, since the sensitivity of the search is signal limited, un- less a cut could be developed which retained most of the gamma-ray events while rejecting a large fraction of the background, it would not improve the sensitivity. The other set of cuts are to ensure the quality of the data. These cuts are necessary to make sure that the detector is running stably and that the data is not corrupted in any way. For this purpose the code checks for repeated event times (caused by a problem in the electronics where the same event gets read out multiple times), time gaps (caused by lost blocks of data, making the detector unstable), time reversals (caused by the incorrect read out of an event), gps timing errors (caused by too few gps satellites being available to the gps receiver, leading to a possible drift in the event times), and large changes in the all-sky rate (due to a number of causes, including partial power failures and PMT failures). Many of these checks were developed while tracking down the cause of false alerts. Time gaps result in the largest cut on the data. An example of the e ect of time 74 Figure 3.6: The all-sky rate for MJD 2425. Top: The rate is averaged over 1 s. Note the long period where the rate has dropped lower and with many large uctuations, followed by a period where there was no data, then normal running again. Middle: The same plot zoomed in. Time periods containing no data are clearly visible. Bottom: The rate is averaged over 100 s bins. The period containing the time gaps is less visible. 75 gaps is shown in Fig.3.6, where the all-sky rate for MJD 2425 has been plotted. In the top plot, it is seen that at the start of the day the rate is constant within some small range, but then begins to uctuate quite a bit until there is a period where there is no data at all, and then normal running resumes. This behavior is due to time gaps in the data, as can be seen in the middle plot of Fig.3.6. The problem began sometime after the installation of a new DAQ computer system (May 22, 2002), and was not immediately noticed. It is only noticeable here by averaging the rate in 1 s bins. If the rate is averaged over 100 s, the problem is not as clearly visible (see the bottom plot in Fig.3.6). Furthermore, the time gaps would appear randomly and then go away if the run was restarted. And since the data appeared ne otherwise and did not a ect searches for emission on longer time scales, it was a di cult bug to track down. However, many time gaps can lead to errors in the background calculation, since an important assumption in the calculation is that the all-sky rate in the detector is constant over at least a 10 s interval. This is not the case when there are many time gaps in the data. Fig.3.7 shows the distribution of the time between events for an entire day of data. On the day for which this is plotted, the detector was running smoothly, and there were no signi cant periods of time gaps. This is used to set a reasonable cut for an acceptable time between events. For the results presented here, the time between events was required to be less than 0.1 s. Otherwise, the event is thrown out and the background map is rebuilt. If there are too many time gaps close together in time, enough events will not be collected to build up the background map, and that block of data will be discarded. The problem was diagnosed after examining some of the low probability locations identi ed by the search, and noting that the background was much lower than would be expected. After lots of searching, it was found that the time gaps were occuring because of dropped data bu ers by the worker computers. Once it began, the DAQ 76 dt (sec) -410 -310 -210 -110 1 10 210 310 410 510 610 710 810 Distribution of the time between events (MJD 3147) Figure 3.7: The distribution of the time between events for MJD 3147. This is for a day with no signi cant time gaps. system would get stuck in this state, and it would not get out of it until the run was restarted. Since restarting the run xed the problem, and since no other solution was found, whenever too many data bu ers get dropped the run is now automati- cally restarted. This correction was not made until slightly over a year after it rst appeared (June 17, 2003). As a result of the time gaps, 9.4% of the total data set considered in this dissertation was discarded. Another check is made for large changes in the all-sky rate. This could be caused by the loss of a patch of 16 tubes due to a single PMT failing, a loss of a large number of tubes due to a low voltage power failure, or other problems with the electronics. To check for this, the average rate is calculated over consecutive 10 s intervals, and each of these 10 s intervals are then averaged together to obtain a running average. If at any point the average rate over a 10 s interval di ers by more than 50 Hz from 77 MJD 2400 2600 2800 3000 3200 Rate (Hz) 0 200 400 600 800 1000 1200 1400 1600 1800 2000 2200 Average Daily Rate Figure 3.8: The average all-sky rate vs Julian day for the entire data sample. the running average, the bu er is searched and the background is re-evaluated. The value of 50 Hz was chosen by examining histograms of the di erence between the rate averaged over 10 s and the running average rate over the entire data set. As in the case of the time gap cut, if there are many uctutions in the rate in a short period of time, enough data will not be collected to create the background maps, causing that data block to be discarded. This cut does not e ect most of the data. For the entire data set that was searched in this analysis, a little less than 15 hours worth of data were discarded altogether. This is negligible compared to the size of the data set. As stated above, the time period searched by this analysis covers 1020 days. This is reduced by the data cuts described above as well as detector down time. The duty cycle of Milagro was 90% for this period. The average all-sky rate for this time period is shown in Fig.3.8. The gaps where no data was taken are due to various 78 MJD 2400 2600 2800 3000 3200 34000 0.2 0.4 0.6 0.8 1 Total Days Searched: 836.585Fraction of Day Searched Figure 3.9: Fraction of day searched for each day in the data sample. The code x that eliminated the time gap problem was installed on MJD 2808. causes such as down time for PMT repairs, a re, or power outages caused by bad weather. The total number of seconds searched is recorded in the log le for each day, al- lowing one to get the total time searched for the whole data set. For the 1020 day period that was analyzed there was a total live time of 836.585 days searched, giving a duty factor of 82%. Almost all of the loss in exposure is due to the time gap cut. In Fig.3.9 the fraction of the day search (total number of seconds searched divided by 86400 s) is plotted for the entire data sample. 79 3.2.6 Monitoring the Results of the Search For each day of data, the search code produces two pieces of output: a log le and a set of histograms. The log le contains information about the position with the lowest probability at each duration for each data bu er searched. This information includes the , , zenith angle, duration, time, and probability, for the location with the lowest probability in that block. The results from searching a single 600 s bu er is shown in Table 3.1. Histograms are produced for each duration and contain the distribution of the number of signal events, the distribution of the background, and the probability distributions. An entry is made in the histogram each time a probability is calculated. In Fig.3.10 these distributions are plotted for three di erent durations (0.01 s, 1.0 s, 10.0 s) for 1 full day of searching (MJD 2702). There are a number of features worth mentioning about these distributions. Since for durations less than 0.2 s, the sky is searched only if there are at least 2 events in one of the ne bins, entries are made in the histograms only when there were two or more events. This is why, the shorter duration probability distributions extend to a much lower probability, the higher probability values only having 1 or 0 events. Additionally, for the very short durations, the shape of the probability distribution is determined almost entirely by the shape of the background distribution since the range in the number of signal events is rather small. In the plots for the 0.001 s duration search in Fig.3.10, there are only 2, 3, or 4 signal events. In the probability distribution for this duration, the region between a probability of 10 5 to 10 7 is entirely due to those locations with 2 signal events. The remaining piece of the probability distribution is due to 3 and 4 signal events. The feature in the probabilty distribution at a probability of 10 4 comes from the extra oversampling done for locations found below this probabilty (this is only done for > 0.2 s durations). A possible GRB candidate would show 80 up as a number of entries at a small probability away from the main probability distribution. tdur t0 tend Searches Pmin Nsig Nbkg tmin (sec) (UT sec) (UT sec) (deg) (deg) (deg) (sec) 0.00025 190.762 790.762 587741 2:42 10 10 3 0.0011 11.08 339.06 28.00 214.525 0.00040 190.762 790.762 1023430 6:18 10 10 3 0.0015 13.40 339.07 47.40 215.397 0.00063 190.762 790.762 1712557 3:16 10 10 3 0.0012 27.26 9.44 16.20 208.554 0.00100 190.762 790.762 2822888 6:18 10 10 3 0.0015 30.79 335.66 7.40 261.064 0.00158 190.762 790.762 4575823 1:38 10 9 3 0.002 32.46 356.89 4.40 221.196 0.00251 190.762 790.762 7297508 2:21 10 9 3 0.0024 37.08 301.15 49.40 426.493 0.00398 190.762 790.762 11605174 1:15 10 9 4 0.013 19.55 8.41 24.00 728.109 0.00631 190.762 790.762 18402324 2:15 10 9 4 0.015 24.39 359.57 13.00 671.628 0.01000 190.762 790.762 29211009 3:30 10 10 4 0.0095 37.08 301.15 49.40 426.491 0.01585 190.762 790.762 46060092 2:15 10 9 4 0.015 37.08 301.15 49.40 426.485 0.02512 190.762 790.762 71866639 2:59 10 9 6 0.11 10.51 336.74 36.20 567.970 0.03981 190.762 790.762 110600012 8:99 10 10 7 0.18 10.36 336.94 36.00 567.972 0.06310 190.762 790.762 167143959 1:34 10 9 6 0.1 31.60 28.14 51.00 424.502 0.10000 190.762 790.762 245688105 4:39 10 9 6 0.12 35.98 316.60 31.80 500.940 0.15849 190.762 790.762 346816343 1:14 10 8 9 0.58 16.40 342.00 45.00 591.767 0.25119 190.762 790.762 567013375 3:62 10 9 8 0.34 35.98 306.18 28.60 195.199 0.39811 190.762 790.762 357706291 2:13 10 11 11 0.55 35.46 306.18 30.00 195.152 0.63096 190.762 790.762 225599900 1:56 10 10 12 0.86 35.46 306.18 30.00 194.966 1.00000 190.762 790.762 142266901 2:77 10 10 14 1.38 35.46 306.18 30.00 195.200 1.58489 190.762 790.762 89684077 3:80 10 8 24 6.14 12.85 359.17 46.20 383.069 2.51189 190.762 790.762 56505264 2:00 10 8 28 7.85 18.58 327.21 42.80 535.534 3.98107 190.762 790.762 35557305 9:36 10 10 18 2.75 40.64 5.62 -2.20 729.332 6.30957 190.762 790.762 22349920 5:57 10 9 38 12.55 27.23 356.86 62.80 689.636 10.00000 190.762 790.762 14013238 8:52 10 8 48 20.06 27.23 356.85 62.80 686.000 15.84893 190.762 790.762 8740279 9:14 10 7 64 32.79 26.97 347.94 9.00 703.693 25.11886 190.762 790.762 5438192 2:64 10 7 179 119.74 4.74 354.17 39.40 658.114 39.81072 190.762 790.762 3348678 1:44 10 6 260 191.39 4.74 354.14 39.40 644.934 Table 3.1: One block of output from a log le, reformated here for ease of reading. The length of the block is 600 s because it was the rst block the search ran on that day, and this is the minimum amount of time required for building the background map. The columns from left to right are the duration, start time of the block, end time of the block, number of actual searches performed, the smallest probability found in the time window, the signal, background, zenith angle, RA, , time for this probability. 81 sigN -110 1 10 210 210 310 410 510 610 710 810 Signal Distribution (0.001 s) bkgN -6 -5 -4 -3 -2 -1 0 1 2 3 10 210 310 410 510 610 710 Background Distribution (0.001 s) log10(prob) -20 -18 -16 -14 -12 -10 -8 -6 -4 -210 210 310 410 510 610 710 810 Probability Distribution (0.001 s) sigN -110 1 10 2101 10 210 310 410 510 610 710 810 910 1010 Signal Distribution (1 s) )bkglog10(N-6 -5 -4 -3 -2 -1 0 1 2 3 10 210 310 410 510 610 710 810 Background Distribution (1 s) log10(prob) -20 -18 -16 -14 -12 -10 -8 -6 -4 -2 1 10 210 310 410 510 610 710 810 910 Probability Distribution (1 s) sigN -110 1 10 2101 10 210 310 410 510 610 710 810 Signal Distribution (10 s) )bkglog10(N -6 -5 -4 -3 -2 -1 0 1 2 3 10 210 310 410 510 610 710 Background Distribution (10 s) log10(prob) -20 -18 -16 -14 -12 -10 -8 -6 -4 -2 1 10 210 310 410 510 610 710 810 Probability Distribution (10 s) Figure 3.10: Signal, background, and probability distributions from one day of searching (MJD 2702). Durations of 0.001 s, 1 s, and 10 s. 3.2.7 The Estimated Annual Rate Just how small a probability would indicate a GRB candidate is in uenced by the fact that a large number of searches are performed (i.e. trials factor). The greater the number of times a random uctuation is searched for, the more times one is likely to nd it. To take a simple example, consider ipping a coin 5 times. The 82 probability of obtaining 5 heads in a row is rather small. However, if we perform this trial (1 trial being 5 ips of the coin) many times, it becomes very likely to get 5 heads in a row. In a similar manner, by searching many times, we are likely to nd some large uctuations of the background. So the probability that is found must be corrected for the number of trials. The rst idea may be to simply count up the total number of searches performed at each duration over the entire data sample, and multiply the probability by this number to get a trials corrected probability. However, each search is not an indepen- dent trial. This is due to the fact that overlapping spatial and time bins are used, and if an upward uctuation of the background is found in one search bin, it is likely to be found in neighbooring spatial and time bins. Therefore, this simple approach does not give the correct number of trials. The trials factor is approximated by using the results of the search to calculate an estimated annual rate of observing a given probability at some duration of searching. From the probability distributions, the number of times a probability below a certain level was seen is known for each day searched. For example, for MJD 2390, at a duration of 1 s, a probability less than 3:16 10 8 was found 2491 times. Therefore, the annual rate of observing an event with P < 3:16 10 8 may be approximated as 2491 365 = 660115. Then, under the assumption that the distribution is linear, the estimated annual rate for observing an event with P < P0 is P0=3:16 10 8 660115. However, it is better to use more than a single day to calculate this factor since if there were large uctuations observed that day, it could throw o the calculation. For each day, the number of events observed below a probability P was plotted versus P. Then each of these distributions were added together, keeping track of the total amount of time searched for each day. If less than 90% of the day was searched, that day was left out of the average. The answer that we get will depend somewhat on the probability value that is used to calculate the annual rate. 83 If the probability is too small, the result su ers from poor statistics. For durations greater than 0.2 s, the number of searches performed changes at a probability of 10 4, therefore we should not include this probability region in the calculation. The probability value used (P0) was set by requiring there to be at least 1000 events below this probability value. Then the estimated annual rate is determined as Rannual(P) = N(P0)TP 0 P = Ntrials P; (3.3) where N(P0) is the number of times a probability less than P0 was observed and T is the total amount of time used to determine N(P0), measured in years. This is done for each duration. In Fig.3.11 Ntrials is plotted as a function of duration. One year of data was used to calculate this result. In the gure is also shown the number of trials in one year calculated assuming that each search is an independent trial (i.e. by adding up the number of searches performed). Note that at longer durations the number of trials calculated from the annual rate approaches the number calculated assuming independent trials. This is due to the fact that for long durations there are a large number of events in each search bin, making it di cult to contain a low probability uctuation in many neighbooring search bins. While for short durations, where there are very few events in the search bin, it is easy to move the search bin around and still keep those few events in one bin. For further illustration of the search algorithm, the number of searches that were actually carried out is also shown in the gure. Recall from the previous discussion that for durations greater than 0.2 s, only every third spatial bin is examined unless a probability less than 10 4 is found. This results in a factor of 9 reduction (as can be seen in the gure) in the number of searches necessary to actually compute while still knowing that any improbable uctuation will still be found. And while the searches were not ?physically? performed, they still count as trials since it is known that any improbable uctuation will still be found. The estimated annual rate is used to de ne a probability threshold for further 84 log10(duration/sec) -4 -3 -2 -1 0 1 2 1010 1110 1210 1310 1410 1510 1610 1710 1810 Maximum number of trials Trials from estimated annual rate Actual number of searches performed in 1 yeartrialsN Figure 3.11: The number of trials for each duration searched. In blue is the number of trials calculated assuming that each search is uncorrelated. In red is the number of trials calculated from the estimated annual rate. In black is the actual number of searches performed, and is only dependent on the search algorithm as discussed earlier. examining a location found by the search. This threshold is set at an estimated annual rate of once per day (365 times per year), which implies a probability of Pthresh = 365=Ntrials (3.4) Fig. 3.12 shows this probability as a function of duration, along with a linear t. Correlations Between Di erent Durations If an improbable uctuation of the background is found at one duration, it is likely to be found at nearby durations. The correlations between di erent durations is much less than for the spatial and time bins. Assuming that the di erent duration searches were independent, an additional trials factor of 27, for the 27 di erent search durations would be incurred. However, because of the correlations, it should 85 log10(duration/sec) -4 -3 -2 -1 0 1 2 log10(probability) -15 -14 -13 -12 -11 -10 -9 Once per day probability threshold Figure 3.12: The probability at which an annual rate of 365 is expected. Pthresh = 365=Ntrials. be less than this. The amount of correlation is determined by looking at the data. A table is made of GRB candidates with an annual rate less than 10 for the entire data set. From this it can be seen how many times a GRB candidate was found at multiple durations. For example, there were 65 GRB candidates with an annual rate less than 10/year. Of these, 6 were found at two di erent durations and 1 at four di erent durations, and so there were a total of 56 distinct GRB candidates. From this it is concluded that the correct trials factor for the multiple durations is 27 (56=65) = 23:26. 3.2.8 Generating GRB alerts The search software can be run online in near real time or o ine on the stored reconstructed data. Most of the details described so far are relevant to both the 86 online and o ine searches. Now a few issues pertaining only to the online search will be discussed. GRB research has bene ted signi cantly from multi-wavelength studies. If a GRB were detected at GeV-TeV energies, valuable information could possibly be obtained from observations at lower energies, and vice versa. This is why it is important to have a fast search, capable of alerting other instruments in the case of a detection. Wherever there is an oportunity, the code has been written to optimize speed. Some of these optimizations have been discussed already, such as the coarse searching that is done unless a probability less than 10 4 is found. Other speed considerations include performing no unnecessary computations (such as evaluating the e ciency map outside of the region speci ed by the zenith angle cut), using pointers whenever possible, using look-up tables for the poisson probabilities (as opposed to calculating it each time), and not searching empty bins for the shorter duration searches. This makes the code more than fast enough to keep up with the data. Originally the online search read the reconstructed data as soon as the REC les were written to disk. However, since a REC le contains about 5 minutes of data, if a GRB were to occur at the very begining of the le, it could not be detected until at least this long. When the new DAQ system was installed, tools for obtaining quicker access to the data were easily available. The reconstructed data could be accessed via a socket connection to the main DAQ computer, which made each block of data available as soon it was reconstructed by one of the worker computers. This allowed the search to run in near real time, the only limitation on the speed being the amount of time required to read in and reconstruct the data and then to search it. The fastest the search could possibly respond to the alert is set by the longest duration and the amount of oversampling in time. With 10% oversampling in time, for a 40 s duration, the search start time steps forward in 4 s increments. So a minimum of 4 s of new data is needed for the 40 s search to run on, and this is the 87 TITLE: MILAGRO BURST POSITION NOTICE NOTICE_DATE: 04/01/2005 05:45:28 UT NOTICE_TYPE: MILAGRO SHORT TRIGGER_NUM: 314, Seq_Num: 1 GRB_DATE: 13374 TJD; 4 DOY; 04/01/2005 GRB_TIME: 20714.912842 SOD; {05:45:14.91} UT GRB_DURATION: 0.00251 SEC GRB_RUN_INFO: RUN 6093 SUBRUN 83 GRB_MIL_RA: 120.9 DEG {08h 00m 00s} (J2000) GRB_MIL_DEC: 10.0 DEG (J2000) GRB_ERROR: 0.54 [deg, radius] GRB_SIGNAL: 4 EVENTS GRB_BACKGROUND: 0.00104 EVENTS GRB_SIGNIFICANCE: 4.7965e-14 (pre trials) GRB_EST_ANN_RATE: 171.079 GRB_MIL_ZEN: 42.6 DEG Table 3.2: Simple form of a GRB alert minimum response time of the search to a GRB. The estimated annual rate is the parameter upon which the decision to send out an alert or not is based. If a GRB candidate is found that has an estimated annual rate below some preset value, an alert is sent out. The frequency of alerts is determined by this preset value. In Table3.2 an example email alert is shown. It contains all the basic information about the burst candidate such as , , time and date in various formats, the number of signal and background events, the probability and the estimated annual rate. 3.3 Search Results In the entire data set no evidence was found for VHE emission from a gamma- ray burst. Although it certainly would have been much more exciting to have found something, interesting limits may still be placed on GRB properties. These are the subject of the next chapter. In Fig.3.13 the probability histograms for the entire data are shown. No signi cant detection is seen. As mentioned before, this would 88 appear as an number of entries in the histograms at a low probability value, away from the main distribution. This is not seen in any of the histograms. 3.3.1 Summary of Lowest Probability Locations Although no strong evidence of emission from a single GRB was observed, those candidates with the lowest probability were examined in more detail to see if there was anything to distinguish them from a typical background uctuation. In Table3.3.1 all the locations found with an estimated annual rate less than 5 are shown. The lowest of these has an estimated annual rate of 0.56, and has a Poisson probability of 1:64 10 14. Given that the duration of this candidate was 1 sec, the number of trials in one year is 2:87 1013. This should be multiplied by a factor of 23.26 for the correlations between the durations and a factor of 2.8 since this is the number of years covered by the data set. Given the trials, this probability is consistent with being due to a uctuation of the background. In addition to this, for the lowest probability candidates checks are made for any di erence in the distribution of event parameters such as X2, core location, number of hit tubes, etc for background and the events for the GRB candidate. None of these showed any indication of being anything other than background. 89 MJD tdur Nsig Nbkg P Ofact Rate tmin 2381 6.31 94 39.20 8:65 10 14 35 0.70 347.89 37.20 4.53 61038.183 2381 10 123 59.92 6:57 10 13 21 3.70 347.48 37.40 4.90 61034.000 2392 0.00398 6 0.009 5:76 10 16 24 1.40 124.00 34.20 27.43 81697.341 2400 2.51 43 10.98 2:06 10 13 7 3.40 73.89 37.60 14.35 80117.371 2616 1 25 3.51 9:79 10 14 7 3.30 59.84 19.20 16.75 21972.100 2716 6.31 31 6.01 5:04 10 13 15 4.10 221.66 4.40 35.67 32313.849 2763 1.58 20 2.06 1:06 10 13 64 2.50 331.88 69.20 33.43 53111.990 2805 10 45 12.27 5:35 10 13 16 3.00 226.98 2.60 33.62 18187.000 2844 0.251 12 0.39 2:02 10 14 23 2.00 272.99 64.00 28.52 16890.276 2844 0.398 14 0.63 9:78 10 15 59 0.68 272.99 64.00 28.52 16890.174 2868 1 28 4.20 1:64 10 14 35 0.56 136.87 24.40 11.52 66981.900 2934 25.1 151 78.34 2:2 10 13 16 0.61 125.81 28.60 19.77 43249.661 2971 0.398 15 0.86 3:59 10 14 6 2.50 87.81 16.00 27.17 35636.404 2997 0.251 13 0.48 7:68 10 15 46 0.76 201.19 61.60 26.88 48681.288 3016 2.51 24 3.36 2:78 10 13 15 4.60 166.41 39.60 33.43 48986.307 3049 25.1 117 54.80 2:03 10 13 52 0.56 26.83 58.80 23.16 85007.261 3127 25.1 45 12.27 5:35 10 13 27 1.50 358.59 9.00 43.61 67697.851 3170 6.31 29 5.25 5:58 10 13 8 4.50 68.18 36.40 39.92 77144.000 3212 0.00631 5 0.003 2:11 10 15 24 3.70 336.17 -2.00 42.23 28445.878 3251 10 49 14.35 6:52 10 13 7 3.70 198.34 52.80 32.64 610.000 3251 15.8 65 22.94 5:49 10 13 22 2.20 198.34 52.80 32.64 607.014 3252 15.8 57 18.35 4:07 10 13 10 1.60 102.66 8.20 33.46 49150.708 3264 0.01 8 0.047 5:85 10 16 94 0.71 301.72 38.20 4.34 11614.929 3272 3.98 36 7.85 2:19 10 13 16 2.50 355.41 32.00 26.87 15974.846 3333 0.631 13 0.59 9:49 10 14 19 4.60 179.75 41.40 36.97 64756.477 Table 3.3: Locations found with estimated annual rate less than 5. 90 log10(prob) -20 -18 -16 -14 -12 -10 -8 -6 -4 -2 1 10 210 310 410 510 610 710 810 910 1010 Probability Distribution (0.000251 s) log10(prob) -20 -18 -16 -14 -12 -10 -8 -6 -4 -2 10 210 310 410 510 610 710 810 910 1010 Probability Distribution (0.000398 s) log10(prob) -20 -18 -16 -14 -12 -10 -8 -6 -4 -2 1 10 210 310 410 510 610 710 810 910 1010 Probability Distribution (0.000631 s) log10(prob) -20 -18 -16 -14 -12 -10 -8 -6 -4 -2 1 10 210 310 410 510 610 710 810 910 1010 Probability Distribution (0.001 s) log10(prob) -20 -18 -16 -14 -12 -10 -8 -6 -4 -21 10 210 310 410 510 610 710 810 910 1010 Probability Distribution (0.00159 s) log10(prob) -20 -18 -16 -14 -12 -10 -8 -6 -4 -2 1 10 210 310 410 510 610 710 810 910 1010 1110 Probability Distribution (0.00251 s) log10(prob) -20 -18 -16 -14 -12 -10 -8 -6 -4 -2 1 10 210 310 410 510 610 710 810 910 1010 1110 Probability Distribution (0.00398 s) log10(prob) -20 -18 -16 -14 -12 -10 -8 -6 -4 -21 10 210 310 410 510 610 710 810 910 1010 1110 Probability Distribution (0.00631 s) log10(prob) -20 -18 -16 -14 -12 -10 -8 -6 -4 -2 1 10 210 310 410 510 610 710 810 910 1010 1110 Probability Distribution (0.01 s) log10(prob) -20 -18 -16 -14 -12 -10 -8 -6 -4 -2 1 10 210 310 410 510 610 710 810 910 1010 1110 Probability Distribution (0.0158 s) log10(prob) -20 -18 -16 -14 -12 -10 -8 -6 -4 -21 10 210 310 410 510 610 710 810 910 1010 1110 1210 Probability Distribution (0.0251 s) log10(prob) -20 -18 -16 -14 -12 -10 -8 -6 -4 -2 1 10 210 310 410 510 610 710 810 910 1010 1110 1210 Probability Distribution (0.0398 s) log10(prob) -20 -18 -16 -14 -12 -10 -8 -6 -4 -21 10 210 310 410 510 610 710 810 910 1010 1110 1210 Probability Distribution (0.0631 s) log10(prob) -20 -18 -16 -14 -12 -10 -8 -6 -4 -2 1 10 210 310 410 510 610 710 810 910 1010 1110 1210 Probability Distribution (0.1 s) log10(prob) -20 -18 -16 -14 -12 -10 -8 -6 -4 -2 1 10 210 310 410 510 610 710 810 910 1010 1110 1210 Probability Distribution (0.158 s) Figure 3.13: Probability distributions for entire data set (search durations 0.000251 s - 0.158 s). 91 log10(prob) -20 -18 -16 -14 -12 -10 -8 -6 -4 -2 1 10 210 310 410 510 610 710 810 910 1010 1110 1210 Probability Distribution (0.251 s) log10(prob) -20 -18 -16 -14 -12 -10 -8 -6 -4 -2 1 10 210 310 410 510 610 710 810 910 1010 1110 1210 Probability Distribution (0.398 s) log10(prob) -20 -18 -16 -14 -12 -10 -8 -6 -4 -2 1 10 210 310 410 510 610 710 810 910 1010 1110 1210 Probability Distribution (0.631 s) log10(prob) -20 -18 -16 -14 -12 -10 -8 -6 -4 -2 1 10 210 310 410 510 610 710 810 910 1010 1110 1210 Probability Distribution (1 s) log10(prob) -20 -18 -16 -14 -12 -10 -8 -6 -4 -2 1 10 210 310 410 510 610 710 810 910 1010 1110 1210 Probability Distribution (1.58 s) log10(prob) -20 -18 -16 -14 -12 -10 -8 -6 -4 -2 1 10 210 310 410 510 610 710 810 910 1010 1110 1210 Probability Distribution (2.51 s) log10(prob) -20 -18 -16 -14 -12 -10 -8 -6 -4 -2 1 10 210 310 410 510 610 710 810 910 1010 1110 1210 Probability Distribution (3.98 s) log10(prob) -20 -18 -16 -14 -12 -10 -8 -6 -4 -2 1 10 210 310 410 510 610 710 810 910 1010 1110 Probability Distribution (6.31 s) log10(prob) -20 -18 -16 -14 -12 -10 -8 -6 -4 -2 1 10 210 310 410 510 610 710 810 910 1010 1110 Probability Distribution (10 s) log10(prob) -20 -18 -16 -14 -12 -10 -8 -6 -4 -2 1 10 210 310 410 510 610 710 810 910 1010 1110 Probability Distribution (15.8 s) log10(prob) -20 -18 -16 -14 -12 -10 -8 -6 -4 -2 1 10 210 310 410 510 610 710 810 910 1010 1110 Probability Distribution (25.1 s) log10(prob) -20 -18 -16 -14 -12 -10 -8 -6 -4 -2 1 10 210 310 410 510 610 710 810 910 1010 Probability Distribution (39.8 s) Figure 3.14: Probability distributions for entire data set ( search durations 0.251 s - 39.8 s). 92 Chapter 4 Sensitivity to GRBs and Results 4.1 Introduction Given the fact that no evidence of VHE emission from a GRB was seen in this data set, what can be inferred about the nature of GRBs and the existence of a VHE component to the GRB spectrum? In order to answer this question, the response of the detector to a potential GRB must rst be studied using simulations. CORSIKA [56] (a program for simulating extensive air showers induced by cosmic rays and gamma rays) and GEANT [57] (a program for simulating the detector response) are used for this purpose. Using CORSIKA, a large number of 100 GeV to 100 TeV gamma-ray induced air showers are generated. These are then propagated through the detector using GEANT 3. A number of assumptions must be made about the properties of GRBs and the potential VHE component of the GRB spectrum. This is a di cult procedure, involv- ing many unknown factors. Models must be used for the form of the VHE spectrum, the redshift and isotropic energy distributions of GRBs as well as any correlation these may have with the burst duration, the amount of absorption due to the IR background, and the amount of absorption at GRB source. Within the context of these models, limits can be set on the relative amount of energy emitted in the form 93 of VHE photons from GRBs. 4.2 De nitions In order to quantify the sensitivity of Milagro to a GRB, some de nitions are needed. The di erential energy spectrum (dN=dE) describes how many particles per unit of area, per unit of energy, are detected from a source. For example, the spectrum of most GRBs is well modeled in the keV/MeV energy range by a broken power law (the Band function [9]): dN dE = 8 >>> >>< >>> >>: I1( EE0 ) Emin < E < Eb I2( EE0 ) Eb E < Emax; (4.1) where Eb is the break energy and E0 is an arbitrary value where the spectrum is normalized and only a ects the spectrum normalization, I. I has units of number per energy, per area, and the integral of dN=dE over energy gives the total number of photons per unit area, with energies between Emin and Emax. As discussed in Chapter 1, the form of the high energy spectrum of GRBs is very uncertain. EGRET observations suggest that the power law measured by BATSE extends to at least a few GeV [36]. Typical values of the high energy spectral index are 2. However, GRB941017 exhibited a distinct high energy spectral component, with a power law index of 1, up to at least 200 MeV [37]. This must of course cuto at some energy above 200 MeV, but how high in energy it extends is uncertain. Some models predict a signi cant sychrotron self-Compton (SSC) component from GRBs. This would imply a similar spectrum as at keV energies, but shifted to higher energies, and with more or less total energy in the VHE component. For the purposes of calculating Milagro?s sensitivity to GRBs, a power law spectrum is assumed over Milagro?s energy range, and calculations are done for a sample of power law indices 94 and for di erent values of Emin and Emax. Another important quantity is the energy uence of an astrophysical source. Energy uence is de ned as: S[Emin;Emax] = Z Emax Emin EdNdE dE = Z Emax Emin EI0( EE 0 ) dE: (4.2) Energy uence has units of energy per unit area. The uence is related to the total energy released by the source: S[Emin;Emax] = 1 + z4 D2 l Eiso[E min;Emax] ; (4.3) where z is the redshift, Eiso is the total isotropic energy (integrated over 4 ) released in photons between Emin and Emax, and Dl is the luminosity distance. In a at CDM model, the luminosity distance is de ned as: Dl = cH 0 Z z 0 dz0q M(1 + z0)3 + ; (4.4) where M is taken here to be 0.3 and to be 0.7, and c H0 = 9:2516 10 27h 1cm; (4.5) where h = H0=100 [58]. The equations above should be corrected for the redshift, with E ! E=(1 + z), and this is done in the calculations that follow. 4.3 Calculating the Number of Photons from a Source: The E ective Area Milagro has no clearly de ned energy threshold, below which no events are detected. At lower primary particle energies, fewer particles survive to ground level, making the event di cult to detect. While at higher primary particle energies, for a power-law spectrum, there are simply fewer primary particles. 95 log10(energy/GeV) 2 2.5 3 3.5 4 4.5 5 arbitrary scale 2 4 6 8 10 12 610? Triggered Energy Distribution (0-45 deg) Figure 4.1: Triggered energy distribution of simulated gamma ray induced air show- ers. Fig.4.1 shows the triggered energy distribution for Milagro from simulated gamma- ray induced air showers thrown on an E 2:4 spectrum. This is the number of gamma rays that passed the trigger conditions described in Chapter 2, passed any cuts, and were reconstructed within a 1:7 radial bin (equivalent in area to the 3 square search bin used in the GRB search). From the gure, it is seen that most of the photons detected are around a few TeV, but the distribution is fairly broad. The median energy of detected events is a function of zenith angle. This is due to the fact that at higher zenith angles, there is more atmosphere to traverse, and fewer showers particles survive to the ground. This is illustrated in Fig.4.2 which shows that the median triggered energy is roughly constant at small zenith angles, but is an order of magnitude larger at above 40 . Using extensive air showers to measure the direction of the primary particle results 96 zenith angle (deg) 0 5 10 15 20 25 30 35 40 45 log10(energy/GeV) 2 2.5 3 3.5 4 4.5 5 Median Triggered Energy vs. Zenith Angle Figure 4.2: Median triggered energy versus zenith angle. in a very di erent acceptance than that of a detector that directly detects the primary particle. For example, the initial trajectory of the primary need not intersect with the physical area of the detector. In order to determine the sensitivity of this detector, this must be taken into account to compute an e ective area. The e ective area is de ned as: Aeff(E; ) = Athrow Ntrig(E; )N throw(E; ) ; (4.6) where Athrow is the area over which the simulated primary particles are thrown (a little more than 1 km), Ntrig is the number of events that pass the trigger criteria (pass the VME trigger, were able to be t to a plane, and are t within the search bin), and Nthrow is the number of events thrown at that energy and zenith angle. Fig.4.3 shows the e ective area for three ranges in zenith angles ( 0 to 15 , 15 to 30 , 30 to 45 ). At an energy of 1 TeV, for zenith angles between 0 and 15 , the e ective area is 6 107cm2. For a mono-energetic spectrum of 1 TeV photons, 97 log10(energy/GeV) 2 2.5 3 3.5 4 4.5 5 ) 2 area (cm 410 510 610 710 810 910 Effective Area o0-15 o 15-30 o30-45 Figure 4.3: The e ective area for three di erent zenith angle ranges. this would be the e ective area (compared to a physical area of 4:8 107cm2). For an arbitrary spectrum, the e ective area as a function of zenith angle only is given by Aeff( ) = R Emax Emin Aeff(E; ) dN dE dER Emax Emin dN dE dE : (4.7) This is shown for an E 2:4 spectrum in Fig.4.4. Below 15 the e ective area is approximately 1 107cm2, and decreases by order of magnitude at 45 . The e ective area is used to calculate the total number of events seen in the detector. One can think of the e ective area as converting the number of particles arriving at the top of the atmosphere into the number seen in the detector on the ground. The number of events is given by N = Z Emax Emin Aeff(E; )dNdE dE = Z Emax Emin Aeff(E; )Io( EE 0 ) dE: (4.8) Using the above, the number of photons expected in Milagro from a GRB may be calculated given the source spectrum. But before this is done there is another 98 zenith angle (deg) 0 5 10 15 20 25 30 35 40 45 ) 2 area (cm 2 4 6 8 10 12 14 610? Effective Area vs. Zenith Angle Figure 4.4: The e ective area versus zenith angle for an E 2:4 spectrum. important consideration that must be made. 4.4 Modi cations to the Spectrum The GRB source may emit high-energy photons with a power law spectrum, but other processes may occur between emission and detection which would modify the spectrum. Particularly, -pair creation between high energy photons and low energy background photons have a large e ect on the spectrum. The probability for a high energy photon to survive the passage through an optically thick photon background is described by the survival probability e (E;z), where (E; z) is the optical depth and is a function of photon energy and possibly redshift. The power law is then modi ed to be dN dE = I0( E E0 ) e (E;z): (4.9) 99 4.4.1 Extra-Galactic Background Light Photons with energies greater than a few hundred GeV have a high probability to interact with infrared (IR) background photons. These IR photons are a part of the extra-Galactic background light (EBL). The EBL is a di use photon background which lls the space between galaxies, and is produced by the radiation of galaxies at di erent redshifts. The EBL is not well measured in some wavelength bands. Systematic errors due to radiation from our own galaxy complicate the measurements. Upper limits are the best constraints in certain wavelength ranges, particularly in the infrared. Di erent models exist based on di erent methods of computing the EBL. The forward evolution approach uses semi-analytic models of galaxy formation to evolve from theoretical initial conditions to the present [59, 60]. This allows the determination of galaxy luminosity functions and in turn the spectral energy distribution (SED) of the EBL. The backward evolution approach takes observations of present day spectra of galaxies, and assumes some form of luminosity evolution with redshift [61]. This allows one to calculate the SED of the EBL based on empirical models. In addition to these approaches, there are others that use some combination of the two methods. Fig.4.5 shows the modi cation to an E 2:4 spectrum at a redshift of 0.2 due to -pair creation on the EBL for two forward evolution models (Primack) and two backward evolution models (Stecker). 4.4.2 Absorption in the Source Region At the burst source, the reball model predicts an optically thick region of low energy photons. Depending on where in the source region the high energy emission is produced, the high energy photons may pair produce on these lower energy pho- tons [48]. Calculations of the opacity due to low energy photons are based on very 100 log10(energy/GeV) 1 1.5 2 2.5 3 3.5 4 4.5 5 arbitrary scale -510 -410 -310 -210 -110 1 10 210 310 410 Unmodified Spectrum Primack 04 Primack/Bullock 99 Stecker Baseline Stecker Fast Evolution Effect of EBL Models (z = 0.2) Figure 4.5: Modi cation to an E2:4 spectrum at a redshift of 0.2 due to -pair creation on the EBL. uncertain models, and are not considered here. It is likely that the opacity varies from burst to burst, depending on local conditions, and is hard to account for in a general way. 4.5 Model Independent Sensitivity Calculations For a given level of the background, one can ask how many signal events are necessary to make a 5 detection. The pre-trials probability required for 5 will vary with the duration of the search due to the varying number of trials. In Fig.4.6 the pre-trials probability required for a 5 detection is plotted as a function of duration. This is determined by the number of e ective trials calculated from the estimated annual rate. No emission at any duration was detected with a probability below this 101 value. log10(duration/sec) -4 -3 -2 -1 0 1 2 Probability -2310 -2210 -2110 -2010 -1910 ?Pre-Trials Probability for 5 Figure 4.6: Pre-trials probability required for a 5 detection. For a given zenith angle and duration, the background is known from the data, and the minimum number of signal events for a 5 probability is calculated using P = e Nbkg(Nbkg)(Nbkg+N ) (Nbkg + N )! ; (4.10) where N is the number of gamma rays from the GRB and Nbkg is the number of background events. For example, the 5 probability for a duration of 1 second is 2 10 20, so given a background of 3.78 events expected, at least 34 photons must be detected to reach this probability. In order to calculate the number of gamma rays from a GRB, a number of as- sumptions must be made about the properties of VHE emission from GRBs. Since the VHE spectrum of GRBs is very uncertain, a simple power law is assumed be- tween 100 GeV and 10 TeV. The dependence of these results on the spectral index 102 and the energy range will be discussed later. Another assumption is about the total energy output in VHE gamma rays. At keV energies, the redshift and measured uence imply a total energy release of as much as 1054 erg, assuming isotropic emis- sion. Since it is unknown how much energy is output in GeV/TeV gamma rays, the number of photons is calculated as a function of the isotropic energy (Eiso). To calculate the number of photons at the earth for GRBs at di erent redshifts, with di erent isotropic energies, and with a power law spectrum over an energy range from Emin to Emax, it is necessary to rst calculate the intrinsic spectal normalization (i.e. not including absorption on the EBL). Using Eqn.4.2 and Eqn.4.3, I0(z; Eiso) is given by I0(z; Eiso) = (1 + z)4 D2 l EisoR Emax Emin EE dE (4.11) The number of predicted photons detected as function of redshift (z), zenith angle ( ), and Eiso is calculated using this value of I0 in the following: N (z; Eiso; ) = Io(z; Eiso) Z Emax Emin Aeff(E; )E e (E;z)dE: (4.12) In order to compute this integral, the e ective area histograms were t by linerally interpolating between sucessive bins. Then the integral was done numerically using a Romberg integration method. This is done for all other integrals as well. In Fig.4.7 N is plotted as a function of redshift for di erent xed values of Eiso, and for a zenith angle of 20 . The Stecker baseline model was used to account for the EBL. The horizontal lines on the gure indicate the minimum number of gamma rays required for a 5 detection for di erent durations. As can be seen from the gure, the number of photons observed decreases rapidly with redshift. This is due to absorption on the EBL as well as the usual 1=r2 decrease. From Fig.4.7, the maximum redshift to which a GRB could be observed at 5 by Milagro can be determined for any duration, zenith angle, and isotropic energy, assuming a power law spectrum between 100 GeV and 10 TeV. It is only necessary 103 redshift 0 0.2 0.4 0.6 0.8 1 1.2 1.4 1.6 1.8 2 N 1 10 210 310 410 510 )? = 20? spectrum, -2.0 vs. Redshift (E?N 0.001s 0.010s 0.100s 1.000s 10.00s 1e50 erg 1e51 erg 1e52 erg 1e53 erg 1e54 erg Figure 4.7: Number of gamma rays dectectable from a GRB as a function of the burst redshift for di erent isotropic energies (curved lines). The horizontal lines indicate the minimum number of gamma rays required for a 5 detection. The intersection of the the curved lines with the horizontal lines indicates the maximum observable redshift for that combination of duration and isotropic energy. to nd the intersection of the horizontal lines in Fig.4.7 with the curves for di erent values of the isotropic energy. This is shown in Fig.4.8. At extremely high isotropic energies, the redshift reach of Milagro is quite far (greater than a z of 1 when Eiso > 1053erg). While these energies are rather large, some theories do predict this much energy released in VHE photons [62], and at keV energies one burst was observed with Eiso > 1054erg. 104 )isolog10(E 47 48 49 50 51 52 53 54 55 56 redshift -310 -210 -110 1 )? = 20? spectrum, -2.0Maximum Observable Redshift (E 0.001s 0.010s 0.100s 1.000s 10.00s Figure 4.8: Maximum observable redshift as a function of isotropic energies for dif- ferent durations. 4.5.1 Spectral Dependence In general, the above results depend upon the spectrum assumed. In this work, a simple power law has been used in all cases. However, the value of the power-law index will a ect the results. Fig.4.9 shows the maximum observable redshift for a duration of 0.01 s and a zenith angle of 20 for three di erent spectral indices. As can be seen from the gure, the maximum redshift is more sensitive to the spectral index at larger GRB energies. This is expected because at higher energies one can see to larger redshifts, where the e ects of the EBL make the result more sensitive to the number of photons at the low end of the spectrum. So one would expect a steeper spectrum to allow observations to larger redshifts, as is seen in the plot. While at low GRB energies, the fact that there is more e ective area at larger energies makes the less steep spectrum allow one to see to higher redshifts. 105 isoE 47 48 49 50 51 52 53 54 redshift -310 -210 -110 1 )? = 20?Effect of Spectral Index on Maximum Observable Redshift ( -1.5E -2.0E -2.5E Figure 4.9: The e ect of the spectral index on the maximum observable redshift. 4.5.2 Dependence on EBL Model In addition to the di erent spectral indexes, the EBL model in uences the sensitivity. Fig.4.10 shows the di erence in the maximum observable redshift for a duration of 0.01 s and a zenith angle of 20 , with an E 2:0 spectrum. As can be seen from the gure, the di erence grows with isotropic energy. This is as expected since at higher isotropic energies, Milagro is sensitive to GRBs at higher redshifts, where the di erences in the EBL models becomes greater. 4.6 Limits on GRB Properties In the preceding section the main assumption was on the shape of the high energy spectrum of GRBs. However, given a model of the redshift, isotropic energy, and duration distributions of GRBs, constraints may be placed on the VHE compo- 106 )isolog10(E 47 48 49 50 51 52 53 54 redshift -310 -210 -110 1 = 0.01s)d, t? = 20? spectrum, -2.0Maximum Observable Redshift (E Stecker, Baseline Primack04 Figure 4.10: The e ect of the EBL model on the maximum observable redshift. nent of the GRB spectrum. This is done by taking models of these distribution and creating a simulation of a GRB population which draws randomly from them. Ad- ditionally, Eqn.4.12 gives the mean number of photons from a source with isotropic energy Eiso, redshift z, zenith angle , and a power-law spectrum. However, a com- plete analysis should take into account the Poisson uctuations about the mean in both the number of signal events and the number of background events. This is described in the following sections. 4.6.1 Poisson Fluctuations Eqn.4.12 gives the expected number of photons (N ;exp) from a source with isotropic energy Eiso, redshift z, zenith angle , and on a power-law spectrum. In general there will be uctuations in the number of photons emitted. The Poisson uctuations in N ;exp are accounted for by generating a random probabilty (P) and 107 nding what value of N ;obs (the number of gamma rays observed) gives this proba- bility. Where P = X e N ;exp(N ; exp)N ;obs N ; obs! (4.13) One can then solve for N ;obs given P and N ; exp. The same procedure is carried out for calculating the number of background events (i.e, given an Nbkg;exp, an Nbkg;obs is calculated). This allows the calculation of the Poisson probability of observing N ;obs + Nbkg;obs events when expecting Nbkg;exp. 4.6.2 The Simulations Using the above procedure, the Poisson probability of observing N ;obs+Nbkg;obs events when expecting Nbkg;exp may be calculated for any combination of redshift, zenith angle, isotropic energy, duration, and spectral index. If this probability is less than the 5 probability described earlier, a GRB with these parameters should be observed by Milagro. Note that a smaller probability gives a larger signi cance. The observed angular distribution of GRBs is isotropic. For this reason, in all the simulations, the zenith angle distribution is isotropic. To take account of Milagro?s exposure, the burst zenith angle is randomly chosen from a cos( ) distribution. The duration distribution has been well measured by BATSE. This distribution is shown in Fig.4.11, t with the sum of two Gaussians. The t is then used to de ne the distribution from which durations are randomly drawn for the simulation. The GRB redshift distribution is not well measured. At the time of this writing, only 40 GRBs have measured redshifts. The measured redshift distribution is shown in Fig.4.12 [15]. As can be seen from the gure, the distribution peaks around a redshift of 1, but includes four GRBs with z > 3 and a number of bursts below z = 1. The smallest measured redshift is 0.0085 and the largest is 4.5. However, the number of measured redshifts is small, and it is unclear if this re ects the true distribution or is in uenced by sampling biases. Additionally, redshifts have been 108 log10(T90/sec) -2 -1.5 -1 -0.5 0 0.5 1 1.5 2 2.5 30 0.01 0.02 0.03 0.04 0.05 0.06 BATSE T90 distribution Figure 4.11: BATSE T90 distribution. measured only for the long duration bursts (td > 2s), except possibly GRB050509b [63]. Because of the uncertainty in the redshift distribution, models are also used for the distribution, as will be discussed later. The uncertainty in the redshift distribution makes the isotropic energy distribu- tion highly uncertain as well. At the earth, a detector like BATSE simply measures a uence (S). Given this, if the redshift of the burst is known, the implied isotropic en- ergy may be calculated according to Eqn.4.3. It is widely accepted that the emission from GRBs is not isotropic, but is more likely beamed. Studies of this, as discussed in Chapter 1, suggest a standard energy release for GRBs of around 1051 erg with beaming factors of a few to several hundred. The results here are independent of whether the emission is beamed or not, since we are only constraining the implied isotropic energy. One could similarly assume some beaming angle and talk about the beaming corrected energy. In addition to this, one can use models of the uence 109 z 0 0.5 1 1.5 2 2.5 3 3.5 4 # of GRBs 0 1 2 3 4 5 6 measured redshift distribution Figure 4.12: The measured redshift distribution. A t using Eqn.4.14 is drawn on top of the distribution. distribution. The Rate of GRBs in the Universe Simulating a population of GRBs allows the calculation of the number of GRBs per year expected to be observed by Milagro given di erent models. This number will depend on the total rate of GRBs in the universe. BATSE observations imply an all-sky GRB rate of about 700 per year [68]. However, this number is dependent on the uence sensitivity of BATSE. Fig.4.14 shows the uence versus duration for a large number of BATSE GRBs. The line drawn on the gure gives an estimate of the minimum uence versus duration of BATSE GRBs. This line is taken to de ne the minimum uence required for a rate of 700 GRBs/year. In the models that are considered below, the same distribution is computed for the simulated GRBs. Then 110 /erg)isolog10(E 47 48 49 50 51 52 53 54 55 56 # of GRBs 0 2 4 6 8 10 Implied Isotropic Energy for GRBs with Measured Redshifts Figure 4.13: The implied isotropic energy distribution for GRBs with measured redshifts. No beaming assumed. the fraction of GRBs that are above this same line allows the GRB rate implied by the model under consideration to be calculated. If a smaller fraction of bursts are above the line, then a higher GRB rate would be implied. 4.6.3 Case 1: Measured Redshift and Eiso Distributions The simplest case to consider is using the measured redshift and isotropic energy distributions. Fig.4.12 shows the measured redshift distribution. This was t to a function of the form (shown by the line in the gure): N(z) = Az2e (z zo)2= 2 (4.14) This form was used since one would expect the nearby behavior of the distribution to go like z2 (the number should be proportional to the surface area of a sphere 111 log10(duration/seconds) -3 -2 -1 0 1 2 3 )) 2 log10(fluence/(erg/cm -9 -8 -7 -6 -5 -4 -3 -2 Fluence vs. Duration for BATSE GRBs Figure 4.14: Fluence vs. Duration for BATSE GRBs. The line drawn on both gures is an estimate of the minimum uence vs duration for BATSE. with radius z) and the distribution should go to zero at z = 0. Fig.4.15 shows the cumulative redshift distriubutions using both the t and the actual measured values. The two distributions match reasonably well. Fig.4.13 shows the implied isotropic energy of these bursts. For this case, it was t to a gaussian (shown by the line in the gure). The duration distribution is the measured BATSE distribution, t to a sum of two gaussians (Fig.4.11). It is assumed here that both the short duration and long duration bursts have the same redshift distribution. A large number GRBs are then simulated by drawing randomly from these dis- tributions. The results of the simulation are shown in Fig.4.16. The top three plots in the gure are the redshift, duration, and energy distributions of the simulated bursts. A single E 2:0 spectrum and the Stecker baseline EBL model was used for all the bursts. This is a typical value for the high energy spectral index for a Band 112 redshift 0 0.5 1 1.5 2 2.5 3 3.5 4 4.5 50 0.2 0.4 0.6 0.8 1 Cumulative Redshift Distribution Fit Distribution Figure 4.15: Cumulative redshift distribution using the t (black) and using the actual distribution (red). The two distributions agree reasonably well. function spectrum, and it was shown earlier that the results are not very sensitive to the choice of spectral index. The simulations were only carried out to a redshift of 2 since it is unlikely that VHE emission from a GRB would reach Milagro without being absorbed by the EBL from beyond there. The bottom left gure is the redshift distribution of the bursts observed in this simulation (those with a probability less than the 5 threshold). The bursts were thrown out to a zenith angle of 90a0 , with 5720 out of 500000 (1.1%) being detected. The total rate of GRBs assumed for this is 700 GRBs/year. Since bursts were simulated out to 90a0 , this means 350 GRBs/year. But these were only simulated out to z = 2, which contains about 84% of all bursts, given this redshift distribution. This gives 295 bursts, of which Milagro should see 1.1%, or 3.37 per year. Given the total data sample searched here (836.585 days, or 113 Redshift Model: Fit to Measured Distribution: Total GRB Rate: 700/year Duration Model: BATSE T90 Distribution Isotropic Energy Model: Fit to Measured Distribution Redshift 0 0.2 0.4 0.6 0.8 1 1.2 1.4 1.6 1.8 20 1000 2000 3000 4000 5000 6000 7000 Simulated Redshift Distribution log10(duration/sec) -2 -1.5 -1 -0.5 0 0.5 1 1.5 2 2.5 30 2000 4000 6000 8000 Simulated Duration Distribution log10(Eiso/erg) 47 48 49 50 51 52 53 54 55 56 570 5000 Simulated Isotropic Energy Distribution redshift 0 0.2 0.4 0.6 0.8 1 1.2 1.4 1.6 1.8 20 50 100 150 200 250 Redshift Distribution for those Observed in the Simulation Eiso(TeV)/Eiso(keV) -410 -310 -210 -110 1 10 210 310 # GRB/year -410 -310 -210 -110 1 10 210 90% CL on ratio: 0.12 Number of GRB/year expected Figure 4.16: Results of a simulation based on the measured distributions. 2.29 years), 7.7 bursts would be expected above the 5 probability threshold given the assumptions in this simulation. This assumes that the energy emitted in the form of VHE photons is equal to that emitted by keV photons. Since we have not observed a GRB candidate at 5 with Milagro, this analysis may be used to set upper limits on the burst parameters. If we assume that all bursts have a VHE component, and that the total energy in this component is some constant factor times the total energy in the low energy component, a limit may be set on this factor if we assume that it is the same for all bursts. This is shown in the bottom right plot on Fig.4.16. This plot was made by simply scaling the number of photons expected by this factor (since N is just proportional to Eiso) and calculating the number of GRBs expected per year. 114 For a total exposure of 2.29 years, our 90% con dence level (CL) upper limit on the number of events observed is 2:3=2:29 = 1:0 event per year. This implies an energy ratio factor of 0.12, shown by a line drawn on the gure. In other words, the total energy in the VHE component is no greater than 0.12 times that in the low energy component at the 90% CL. Sources of Systematic Error Two sources of systematic error come from the EBL model, and the assumption of a single E 2:0 spectrum for the entire burst simulation. In order to get an estimate of the e ect of these two parameters, the simulation was run four additional times. The four simulations were for all combinations of the Stecker fast evolution EBL, the Primack 04 EBL, an E 1:5 spectrum, and an E 2:5 spectrum. This gives a broad range in the space covered by these parameters, giving a maximum and minimum for the constraint on the energy ratio factor calculated above. The result of these simu- lations are shown in Table 4.1. As can be see in the table, the result is less sensitive to changes in spectra compared to the di erent EBL models. The Primack models in general predict much less absorption than the Stecker models. For this reason, in order to be more conservative, the remaining studies use the Stecker baseline model. Spectral Index EBL Model Upper Limit -1.5 Stecker Fast Evolution 0.52 -2.5 Stecker Fast Evolution 0.10 -1.5 Primack 04 0.0025 -2.5 Primack 04 0.0052 Table 4.1: Dependence of upper limit on the source spectrum and EBL model. 115 4.6.4 Case 2: Measured Redshift and BATSE Fluence Dis- tributions Instead of using the distribution of isotropic energies inferred from the GRBs with measured redshifts, the measured BATSE uence may be used. Then, a redshift is drawn randomly from some distribution (in this case using the t to the measured distribution), and Eiso is calculated using Eqn.4.3. Since the BATSE uence distri- bution is being used, a rate of 700 GRB/year may be used without comparing to the BATSE uence versus duration line. Fig.4.17 shows the measured BATSE uence distribution [11] for long and short burst durations. The BATSE data shows a correlation between uence and duration. Several authors have discussed whether this is intrinsic to the burst or the result of an instrumental bias [64]. For this test, we will simply take the correlation between uence and duration as observed by BATSE to be intrinsic to the GRB population. Fig.4.18 shows the result of a simulation of a population of bursts based on this model. Again all the bursts were assumed to have an E 2:0 spectrum, and the Stecker baseline EBL was used. Out of 500000 simulated bursts, 185 would be detectable by Milagro (0.037%). Assuming equal total energy in VHE photons and keV photons (energy ratio factor of 1), this implies a rate of 0.11 GRBs/year detectable by Milagro. Therefore, none would be expected to have been observed in this data set. The reason that so few GRBs are expected to be observable in this model is due to the nature of the BATSE uence distribution. This distribution is such that shorter duration bursts have smaller uences, and hence lower isotropic energies for a given redshift. Given that this analysis is most sensitive at short durations, the small expected GRB rate from this model makes sense. Again, varying the energy ratio factor described above allows an upper limit to be set. The 90% CL upper limit on the energy ratio factor is 30.20 in this case. While 116 )2log10(fluence/erg/cm -9 -8 -7 -6 -5 -4 -3 # of GRBs 0 20 40 60 80 100 120 BATSE Fluence Distribution (50-300 keV) T90<2 s T90>2 s Figure 4.17: BATSE uence distribution for short (T90 < 2s) and long (T90 > 2s) duration bursts. this is not as hard of a constraint as in the previous model, some models predict more energy in the VHE component than in the keV component [45, 46, 62]. 4.6.5 Case 3: The Lag-Luminosity Relationship In [17] a relationship is found between the lag time ( 0) between the light curves in two di erent energy bands and the GRB luminosity. In [65] this relationship is applied to a large sample of BATSE bursts in order to generate redshift and isotropic energy distributions. The resulting distributions from this study are shown in Fig.4.19. For the pur- pose of the simulation, the Eiso distribution was t to a Gaussian and the redshift distribution was t to a function of the form N(z) = Ae( (z zo) 2 z ) (4.15) 117 Redshift Model: Fit to Measured Distribution: Total GRB Rate: 700/year Duration Model: BATSE T90 Distribution Isotropic Energy Model: From BATSE Fluence Distribution Redshift 0 0.2 0.4 0.6 0.8 1 1.2 1.4 1.6 1.8 20 Simulated Redshift Distribution log10(duration/sec) -2 -1.5 -1 -0.5 0 0.5 1 1.5 2 2.5 30 20 40 60 80 100 120 140 160 310? Simulated Duration Distribution log10(Eiso/erg) 47 48 49 50 51 52 53 54 55 56 570 20 40 60 80 100 120 140 160 180 200 310? Simulated Isotropic Energy Distribution redshift 0 0.2 0.4 0.6 0.8 1 1.2 1.4 1.6 1.8 20 2 4 6 8 10 12 14 16 Redshift Distribution for those Observed in the Simulation Eiso(TeV)/Eiso(keV) -110 1 10 210 310 # GRB/year -410 -310 -210 -110 1 10 90% CL on ratio: 30.20Number of GRB/year expected Figure 4.18: Results of a simulation based on the measured redshift distribution and using the measured BATSE uence (which is correlated with duration). The ts are shown by the lines on the plots. Compared to the t to the measured redshift distribution, this model has fewer GRBs at low redshifts (84% of bursts at z < 2 for the t to the measured distribution, while 66% of the bursts are at z < 2 for this model). On the other hand, the Eiso distributions are nearly identical. Given this, it is expected that slightly fewer GRBs would be detectable by Milagro in this model compared to case 1. In this model as well, a total GRB rate of 700/year is assumed. The results of this simulation are shown in Fig.4.20. For an energy ratio factor of 1, the expected rate of GRBs which would be detectable by Milagro is 1.88/year. The 90% con dence level upper limit on the energy ratio factor is 0.36. 118 redshift 0 1 2 3 4 5 6 7 8 9 10 # of GRBs 0 10 20 30 40 50 60 70 Redshift Distribution /erg)isolog10(E 47 48 49 50 51 52 53 54 55 56 # of GRBs 0 20 40 60 80 100 Isotropic Energy Distribution Figure 4.19: Redshift and Eiso distributions from the lag-luminosity relationship. 4.6.6 Case 4: Collapsar/Binary Merger Scenarios Another manner in which to arrive at a GRB redshift distribution is to consider models for the GRB progenitor. One assumption that is often made is that the GRB rate follows the star formation rate (SFR). This is often done, at least for the long durations, since it is believed that GRBs are related to the collapse of massive stars. Short duration bursts may also follow the SFR or, if they are due to binary mergers, they may trail slighty the SFR. There are a number of parameterizations of the SFR, and following [67] three di erent parameterizations are considered. RSF1(z) = 0:3h65 exp(3:4z)exp(3:8z) + 45M yr 1Mpc 3 (4.16) 119 Redshift Model: Lag-Luminosity Relation: Total GRB Rate: 700/year Duration Model: BATSE T90 Distribution Isotropic Energy Model: Lag-Luminosity Relation Redshift 0 0.2 0.4 0.6 0.8 1 1.2 1.4 1.6 1.8 20 1000 2000 3000 4000 5000 6000 7000 Simulated Redshift Distribution log10(duration/sec) -2 -1.5 -1 -0.5 0 0.5 1 1.5 2 2.5 30 2000 4000 6000 8000 10000 12000 14000 16000 Simulated Duration Distribution log10(Eiso/erg) 47 48 49 50 51 52 53 54 55 56 570 5000 10000 15000 20000 25000 Simulated Isotropic Energy Distribution redshift 0 0.2 0.4 0.6 0.8 1 1.2 1.4 1.6 1.8 20 20 40 60 80 100 120 140 160 180 200 220 Redshift Distribution for those Observed in the Simulation Eiso(TeV)/Eiso(keV) -210 -110 1 10 210 310 # GRB/year -410 -310 -210 -110 1 10 210 90% CL on ratio: 0.36 Number of GRB/year expected Figure 4.20: Results of the simulation using the lag-luminosity relation. RSF2(z) = 0:15h65 exp(3:4z)exp(3:4z) + 22M yr 1Mpc 3 (4.17) RSF3(z) = 0:2h65exp(3:05z 0:4)exp(2:93z) + 15 M yr 1Mpc 3 (4.18) Fig.4.21 shows these di erent parameterizations. The di erence in the three models result from the attempts to take into account uncertainties in the data or to correct for a potential underestimation of the SFR at high-z due to the e ects of dust extinction. In [66], the redshift distribution of binary mergers is calculated for a number of di erent scenarios. There are a large number of uncertainties in these calculations, including di erent initial star formation rates and uncertainties resulting from dif- 120 redshift 0 1 2 3 4 5 6 7 8 9 10 ) -3 Mpc -1 yr SFR (M 0 0.2 0.4 0.6 0.8 1 1.2 Models of the Star Formation Rate SFR1 SFR2 SFR3 Figure 4.21: Models of the Star Formation Rate. ferent merger scenarios. The redshift distribution of BH-NS and NS-NS mergers for one particular model is shown in Fig.4.22. The main feature to note in the gure is that there are more bursts at z < 1 if GRBs follow the binary merger distribution than if they follow the SFR. For the redshift distribution shown in Fig.4.22, 26% of the distribution is at z < 1, while for SFR1 about 9.7% are at z < 1. For the purpose of the simulation, it is assumed that for durations less than 2 s, the GRB redshift distribution follows the binary merger redshift distribution, and for durations longer than 2 s it follows the SFR (using the SFR1 parameterization). Then the simulation is run as described before, but additionally keeping track of how many GRBs would have been detected at short and long durations. For the Eiso distribution the measured BATSE uence distribution (and a GRB rate of 700/year) is used as in Case 2 above. More than the previous models, this model takes into account di erences between short and long duration GRBs. However much of this is 121 redshift 0 2 4 6 8 10 GRB rate (arbitrary scale) 0.2 0.4 0.6 0.8 1 GRB Rate Due to Binary Mergers Figure 4.22: Redshift distribution of NS-NS and NS-BH mergers. still very uncertain. Fig.4.23 shows the results of this simulation. The 90% CL upper limit on the energy ratio is 91.20. This takes into account the di erent total rates for short and long duration GRBs while still assuming a total rate of 700 GRBs/year. From the BATSE duration distribution about 77% of the bursts are long duration while 23% are short duration. Additionally a correction was made for the fact that a di erent fraction of the redshift distribution was simulated for long and short durations by only going out to z = 2. For short durations, z < 2 contains 64.2% of the redshift distribution, while for long durations it only contains 37.4%. 4.6.7 Case 5: Broken Power Luminosity Distributions The BATSE peak ux distribution can be used to determine the GRB red- shift and luminosity distributions [68]. The observed peak ux distribution depends 122 >2s): Total GRB Rate: 700/yeardur<2s)/ SFR (tdurRedshift Model: Binary Merger (t Duration Model: BATSE T90 Distribution Isotropic Energy Model: From BATSE Fluence Distribution Redshift 0 0.2 0.4 0.6 0.8 1 1.2 1.4 1.6 1.8 2 2000 4000 6000 8000 Simulated Redshift Distribution log10(duration/sec) -2 -1.5 -1 -0.5 0 0.5 1 1.5 2 2.5 30 5000 10000 15000 20000 25000 30000 Simulated Duration Distribution log10(Eiso/erg) 47 48 49 50 51 52 53 54 55 56 570 5000 10000 15000 20000 25000 30000 35000 40000 45000 Simulated Isotropic Energy Distribution redshift 0 0.2 0.4 0.6 0.8 1 1.2 1.4 1.6 1.8 20 2 4 6 8 10 12 14 16 18 Redshift Distribution for those Observed in the Simulation Eiso(TeV)/Eiso(keV) -110 1 10 210 310 # GRB/year -410 -310 -210 -110 1 10 210 90% CL on ratio: 91.20 Number of GRB/year expected Figure 4.23: Results of the simulation based on binary merger and SFR redshift distributions. on both the redshift and luminosity distributions, making the problem somewhat complicated. For the luminosity function a broken power law is often used. The local (z = 0) luminosity function for bursts with peak luminosity L may be parameterized as 0(L) = co 8 >< >>: (L=L ) L = 1 < L < L (L=L ) L < L < 2L : (4.19) The peak ux P(L; z) is given by P(L; z) = L4 D2 L C(E1(1 + z); E2(1 + z)) C(E1; E2) (4.20) where C(E1; E2) is the integral of the spectrum from E1 to E2. Then the number of 123 bursts with a peak luminosity greater than P, is then given by N(> P) = Z 0(L)dlogL Z zmax 0 R(z) (1 + z) dV (z) dz F(z; M; )dz (4.21) where zmax is the maximum redshift to which a burst may be detected given a ux limit Plim (which is detector dependent), dV (z)=dz is the comoving volume element, R(z) is given by one of Eqn.4.16 - 4.18, and F(z; M; ) = q M(1 + z)3 + (1 + z)3=2 : (4.22) For a given redshift distribution, Eqn.4.21 is used to determine the parameters of the luminosity function by comparing to the N(> P) distribution measured by BATSE. Once these parameters are found, the luminosity and redshift distributions are speci ed. In most calculations this is done using the peak luminosity (measured in erg=cm2=s), while for the calculations here, the distribution of isotropic energies is needed. In [69] the parameters were appropriately adjusted to apply to the isotropic energy distri- bution. For this simulation, the parameters used were L = 4:61 1051erg, 1 = 30, 2 = 10, = 0:5, = 1:5. For the redshift distributions the SFR1 parameterization was used. In this case, the total GRB rate will be di erent than 700/year. Given this redshift and isotropic energy distribution, 73% of the GRBs were above the minimum uence line shown in Fig.4.14. This implies a total GRB rate of 700/0.73 = 959/year in order to be consistent with BATSE. Fig.4.24 shows the results of this simulation. The 90% CL upper limit on the energy factor in this model is 75.86. 124 Redshift Model: SFR1: Total GRB Rate: 959/year Duration Model: BATSE T90 Distribution Isotropic Energy Model: Broken Power Law Redshift 0 0.2 0.4 0.6 0.8 1 1.2 1.4 1.6 1.8 20 1000 2000 3000 4000 5000 6000 7000 8000 Simulated Redshift Distribution log10(duration/sec) -2 -1.5 -1 -0.5 0 0.5 1 1.5 2 2.5 30 2000 4000 6000 8000 Simulated Duration Distribution log10(Eiso/erg) 47 48 49 50 51 52 53 54 55 56 570 10000 20000 30000 40000 50000 60000 Simulated Isotropic Energy Distribution redshift 0 0.2 0.4 0.6 0.8 1 1.2 1.4 1.6 1.8 20 20 40 60 80 100 120 Redshift Distribution for those Observed in the Simulation Eiso(TeV)/Eiso(keV) -410 -310 -210 -110 1 10 210 310 # GRB/year -410 -310 -210 -110 1 10 210 90% CL on ratio: 75.86 Number of GRB/year expected Figure 4.24: Results of the simulation based on a broken power isotropic energy distribution and the SFR1 redshift distribution. 125 Chapter 5 Conclusions 5.1 Summary In Chapter 1, along with an overview of the history of GRB research, the experimental and theoretical motivations for searching for VHE emission from GRBs were given. In Chapter 2, the Milagro detector, an instrument ideal for searching for VHE emission from GRBs, was described. In Chapter 3, a search of the Milagro data for VHE emission from GRBs was presented. The data was searched from MJD 2353 to 3372 (March 19, 2002 to January 1, 2005), and no evidence for VHE emission from GRBs was observed. In Chapter 4, the sensitivity of Milagro for detecting VHE emission from GRBs was presented. Then given the fact that no evidence for VHE emission from GRBs was observed, constraints were placed on the VHE component of the burst spectrum. This was done by simulating a population of GRBs given di erent models for GRB redshift and isotropic energy distributions, models of the IR background. It was then assumed that the total energy radiated in the form of VHE photons was directly proportional to the total energy radiated at keV/MeV energies. It was further assumed that the ratio of proportionality (the energy ratio) did not vary from burst to burst. While this is not necessarily the case, it provides a handle on a problem with many unknown factors. 126 This allowed the calculation of the number of GRBs/year expected to be observed by Milagro as a function of the energy ratio for di erent models of the burst popula- tion. Given that no evidence of VHE emission was observed by Milagro, the 90% CL upper limit on the number of bursts in this data set is 2.3. Given the 2.3 year period searched in this analysis, the upper limit on the rate of GRBs observed by Milagro is 1 GRB/yr. The value of the energy ratio that gave a GRB rate of 1 GRB/year in Milagro is the upper limit on the ratio. The di erent models considered in Chapter 4 resulted in di erent limits on the energy ratio. Table5.1 summarizes the limits obtained for the di erent models considered in this dissertation. Model 90% CL Upper Limit Measured Redshift and Eiso Distributions 0.12 Measured Redshift and BATSE Fluence Distributions 30.20 Lag-Luminosity Model 0.36 Binary Merger/SFR Model 91.20 SFR Redshift and Broken Power Law Eiso Distribution 75.86 Table 5.1: Summary of results from di erent models. The column on the left gives the redshift and isotropic energy model and the column on the right gives the 90% CL upper limit on the energy ratio. The large spread in values for the upper limit is understandable, given the un- certainty in the di erent models. Redshifts have not been measured at all for the short duration bursts, and are still not well measured for long duration bursts. This makes the isotropic energy distribution uncertain as well. The upper limit is sensi- tive to how many low redshift bursts are bright enough to be detected by Milagro. In general these results can be adapted to constrain any model that predicts a VHE component to the burst spectrum. Which of the models above more accurately re- ects the actual GRB population will be decided with future observations. A new 127 GRB observatory called SWIFT has been launched recently, and will obtain redshifts for a large number of GRBs, including those with short durations [70]. With better knowledge of these distributions more accurate limits will be set. In addition to this, Milagro observations are continuing. If it turns out that VHE emission by GRBs is a rare phenomenon, continued observations could still result in its detection. Another way of interpreting these results is shown in Fig5.1. If it is not assumed that all GRBs have a VHE component, an upper limit on the fraction of GRBs that could have this component can be calculated. It is still assumed that the energy ratio is the same for all bursts that do have a VHE component. As before, the number of GRBs/year observable is determined as a function of the energy ratio. The value of the energy ratio that results in an expected rate of 1 GRB/year is the 90% CL upper limit on the energy ratio. This means that, at values of the energy ratio below the upper limit value, all GRBs could have a VHE component and still be consistent with the Milagro observations. However, at higher values of the energy ratio, a smaller fraction of GRBs could have a VHE component and still be consistent. For example, for the model that used ts to the measured distributions, the upper limit was 0.12, this is shown by the black curve in Fig5.1. Below an energy ratio of 0.12, 100% of GRBs could have a VHE component. At an energy ratio of 1, 30% could have a VHE component. 5.2 Future Directions Milagro is a unique instrument in the eld of VHE astrophysics, and much was learned in its construction and operation. A new instrument is currently being con- sidered which would build on what was learned with Milagro. This new instrument, called miniHAWC, would be located at a much higher altitude (4300 m above sea level, while Milagro is at 2600 m). This would allow the detection of a much greater 128 Eiso(TeV)/Eiso(keV) -210 -110 1 10 210 310 Fraction of GRBs -210 -110 1 Fit to measured distributions BATSE Fluence Lag-Luminosity Relation SFR/Binary Merger Broken Power Law Upper Limit on the Fraction of GRBs with VHE Emission Figure 5.1: Upper limit on the fraction of GRBs with a VHE component for the di erent models considered in Chapter 4. number of shower particles since the detector would be closer to shower maximum. miniHAWC would also be a much larger detector. Currently under consideration is a 150m 150m area detector (Milagro is 60m 80m). This would allow a greater portion of the shower to be contained within the detector. Finally, the PMTs in miniHAWC would be optically isolated from each other. This would prevent light which travels horizontally in the detector from triggering the detector. Simulations of this detector have been created in the same manner as for Milagro. With its much larger physical area and higher altitude, miniHAWC has a larger e ective area than Milagro, particularly at low energies. As done in Chapter 4, the maximum observable redshift is calculated as a function of isotropic energy. This is show in Fig.5.2 for a zenith angle of 20 , a duration of 10 s, and a E 2:0 power-law 129 isoE 47 48 49 50 51 52 53 54 55 redshift -310 -210 -110 1 10 , duration = 10 s)? = 20? spectrum, -2.0Maximum Observable Redshift (E Milagro miniHAWC Figure 5.2: Maximum observable redshift as a function of isotropic energy for mini- HAWC compared to Milagro. The redshift reach of miniHAWC is much greater than that of Milagro. spectrum. This is also plotted for Milagro for comparison. As can be seen from the gure, miniHAWC has a much larger redshift reach than Milagro. At more moderate isotropic energies, miniHAWC is capable of detecting GRBs out to redshifts of 1, and therefore is sensitive to a great fraction of the GRB population. Simulations were run using the same models as described in Chapter 4. Fig.5.3 shows the results from a simulation based on the measured redshift and isotropic energy distributions. In this model, for an energy ratio of 1, 18.5 GRBs/year are expected (compared with 3.4 for Milagro). The 90% CL upper limit on the energy ratio, assuming the same exposure time and no detection of a GRB, is 0.0052 (more than a factor of ten better than the limit set by Milagro). Of course, with the greater sensitivity of an instrument such as miniHAWC, it 130 Redshift Model: Fit to Measured Distribution: Total GRB Rate: 700/year Duration Model: BATSE T90 Distribution Isotropic Energy Model: Fit to Measured Distribution Redshift 0 0.2 0.4 0.6 0.8 1 1.2 1.4 1.6 1.8 20 1000 2000 3000 4000 5000 6000 7000 Simulated Redshift Distribution log10(duration/sec) -2 -1.5 -1 -0.5 0 0.5 1 1.5 2 2.5 30 2000 4000 6000 8000 Simulated Duration Distribution log10(Eiso/erg) 47 48 49 50 51 52 53 54 55 56 570 5000 10000 15000 20000 25000 Simulated Isotropic Energy Distribution redshift 0 0.2 0.4 0.6 0.8 1 1.2 1.4 1.6 1.8 20 200 400 600 800 1000 Redshift Distribution for those Observed in the Simulation Eiso(TeV)/Eiso(keV) -410 -310 -210 -110 1 10 210 310 # GRB/year -410 -310 -210 -110 1 10 210 90% CL on ratio: 0.0052Number of GRB/year expected Figure 5.3: Results of the simulation based on a t to the measured distributions using the miniHAWC e ective areas. is hoped that GRBs will be detected. With an improved sensitivity to gamma ray energies of a few hundred GeV and lower, miniHAWC will be sensitive to GRBs at higher redshifts and therefore to a higher fraction of the total burst population. 131 Bibliography [1] Singer, S. Proc. IEEE, 53, 1935 (1965). [2] Dickinson, H & Tamarkin, P., Proc. IEEE, 53, 1921 (1965). [3] Klebesadel, R., ApJ. 182, L85 (1973). [4] Mazets, E. P., et al., Astrophysics and Space Science, vol. 80, no. 1,(1983). [5] Klebesadel, R., et al., IEEE Trans on Geoscience and Remote Sensing, GE-18, 76 (1980). [6] BATSE http://cossc.gsfc.nasa.gov/batse [7] Bonnell, J. T., http://cossc.gsfc.nasa.gov/images/epo/gallery/grbs/ [8] BATSE http://f64.nsstc.nasa.gov/batse/grb/duration/ [9] Band, D., et al., ApJ., 413, 281 (1993). [10] BATSE http://www.batse.msfc.nasa.gov/batse/grb/skymap/ [11] Paciesas, W. S., et al., ApJ. Supplement Series, 122, 2, 465-495 (1999). [12] BeppoSAX http://www.asdc.asi.it/bepposax/ [13] Costa, E., et al., Nature, 387, 783 (1997). [14] Metzger, M. R., et al., Nature, 387, 878 (1997). 132 [15] Greiner, J., http://www.mpe.mpg.de/ jcg/grb.html [16] Barthelmy, S., http://gcn.gsfc.nasa.gov/gcn/gcn main.html [17] Norris, J. P., Marani, G., & Bonnell, J., ApJ, 534, 248 (2000). [18] Donaghy, T., et al., AIP Conf.Proc. 662, 450-453 (2003). [19] M esz aros, P., International J. of Mod. Phys. A, in press, astro-ph/0311321 (2003) [20] Harrison, F. A., et al., ApJ., 523 L121 (1999). [21] Galama, T. J., et al., Nature, 395, 670 (1998). [22] Kulkarni, S. R., et al., Nature, 395, 663 (1998). [23] Bloom, J. S., et al., Nature, 401, 453 (1999). [24] Stanek, K. Z., et al., ApJ., 591, L17 (2003) [25] Frail D. A., et al., ApJ., 562, L55 (2001). [26] Piran, T., Phys. Rep. 314, 575 (1999). [27] Waxman, E., in Gamma-Ray Bursts: The Underlying Model, in Super- novae and Gamma-Ray bursters, Ed. K. Weiler (springer), 598, 393 (2003). [28] Lithwick, Y. & Sari, R., ApJ., 555, 540 (2001). [29] http://cossc.gsfc.nasa.gov/images/epo/gallery/grbs/ [30] Kobayashi, S., et al., ApJ., 513, 679 (1997). [31] M esz aros, P. & Rees, M. J., ApJ., 405, 405 (1993). [32] Weiqun, Z., Woosley, S. E., & MacFayden, A. I., ApJ., 586, 356 (2003). 133 [33] Leonard P.J.T., http://heasarc.gsfc.nasa.gov/docs/objects/grbs/grbs.html [34] Matz, S. M., et al., ApJ., 288, L37 (1985). [35] Hurley, K., et al., Nature, 372, 652 (1994). [36] Dingus, B. L. in High-Energy Gamma Ray Astronomy (eds. Aharonian, F. A. & Volk, H. J.) 383391 (AIP, New York, 2001). [37] Gonzales, M. M., et al., Nature 424, 847 (2003). [38] Amenomori, M., et al., A&A, 311, 919 (1996). [39] Connaughton, V., et al., ApJ., 479, 859 (1997). [40] Atkins, R., et al., ApJ., in press, astro-ph/0503270 (2005). [41] Atkins, R., et al., ApJ., 533, L119 (2000). [42] Atkins, R., et al., ApJ., 604, L25 (2004). [43] Hays, E. Ph.D Thesis, University of Maryland, College Park (2004). [44] Atkins, R., et al., ApJ., 608, L680 (2004). [45] Dermer, C. D. & Chiang, J., in GeV-TeV Gamma Ray Astrophysics Work- shop : Towards a Major Atmospheric Cherenkov Detector VI, Ed. Dingus, B. L., et al. AIP, (2000). [46] Pe?er, A. & Waxman, E., ApJ., 613, L448 (2004). [47] Waxman, E., Phys. Rev. Lett. 75,386 (1995). [48] Razzaque, S., M esz aros, P., & Bing, Z., ApJ., 613, L1072 (2004). [49] Bullock, J. S., Ph.D Thesis, University of California, Santa Cruz (1999). 134 [50] Swordy, S. P., Space Science Reviews, v. 99, Issue 1/4, p. 85-94 (2001). [51] Gaisser, T.K., Cosmic Rays and Particle Physics, Cambridge University Press, Cambridge England, 1990. [52] Hays, E., Milagro Collaboration Internal Memo (2001). [53] Alexandreas, D. E., et al., Nucl. Inst. Meth., A328, 570 (1993). [54] Atkins, R., et al., ApJ., 595, L803 (2003). [55] Li, T. P. & Ma, Y. Q., ApJ. 272, L317 (1983). [56] D. Heck et al., CORSIKA: A Monte Carlo Code to Simulate Extensive Air Showers, Forschungszentrum Karlsruhe, Wissenschaftliche Berichte FZKA 6019, (1998). [57] CERN Application Software Group, CERN W,5013, Ver- sion 3.21, (1994). [58] Hogg, D. W., astro-ph/9905116 (1999). [59] Primack, J. R., Bullock, J. S., Somerville, R. S. & Maxminn, D. Astropar- ticle Phys., 11, 93 (1999). [60] Primack, J. R., Bullock, J. S., & Somerville, R. S. 2004 in Gamma 2004 Heidelberg, ed. Aharonian, F. A. (2004). [61] de Jager, O. C. & Stecker, ApJ. 566, L738 (2002). [62] Pe?er, A. & Waxman, E., ApJ., 603, L1 (2004). [63] Gehrels, N., et al., astro-ph/0505630 (2005). [64] Hakkila, J., et al., Gamma-ray Bursts, 5th Huntsville Symposium, Huntsville, Alabama, Eds.: Kippen, R. M., et al., AIP Conference Se- 135 ries, Vol. 526, p.48, American Institute of Physics, Melville, New York, (2000). [65] Band, D. L., et al., ApJ. 613, 484 (2004). [66] Fryer, C. L., ApJ., 526, L152 (1999). [67] Porciani, C., & Madau, P., ApJ., 548, L522 (2001). [68] Schmidt, M., ApJ., 552, L36 (2001). [69] Guetta, D. & Piran, T., Astronomy and Astrophysics, Vol. 435, Issue 2 (2005). [70] SWIFT http://swift.gsfc.nasa.gov/docs/swift/swiftsc.html 136