TRANSIENT SEARCHES WITH THE HAWC GAMMA-RAY OBSERVATORY By Alison Peisker A DISSERTATION Submitted to Michigan State University in partial fulfillment of the requirements for the degree of Physics – Doctor of Philosophy 2021 ABSTRACT TRANSIENT SEARCHES WITH THE HAWC GAMMA-RAY OBSERVATORY By Alison Peisker The universe contains many extremely energetic astrophysical sources that emit particles on timescales from fractions of a second to thousands of years. Studying these sources through gamma-ray and multi-messenger astronomy may help to reveal the answers to several remaining fundamental questions in physics. The High Altitude Water Cherenkov (HAWC) gamma-ray observatory, an extensive air shower array, is well-suited to perform transient searches due to its large field of view and high duty cycle. In this work, a catalog search of all very-high-energy gamma-ray sources in HAWC’s field of view is conducted, revealing a total of 65 sources. Then a search for evaporating primordial black holes is performed. None are detected, so upper limits are set on their burst rate density. Finally, the IceCube neutrino observatory is introduced in order to perform a multi-messenger search for flaring gamma rays and neutrinos originating from the same point in the sky. None are detected, so upper limits are placed on the rate density and total energy of sources that produce both gamma rays and neutrinos. Performing these searches contributes to the understanding of the origins of cosmic rays and how astrophysical sources accelerate their particles to the high energies that are observed on Earth. ACKNOWLEDGEMENTS This work would not have been possible without the constant help and support of many, many people. First and foremost, thank you to my advisor, Kirsten Tollefson. The guidance you provided is invaluable. Thank you for teaching me how to be a good scientist, and for always encouraging me to work on whatever makes me happy. I greatly appreciate your complete support in everything that I do. Thank you to my thesis committee, Andrea Albert, Laura Chomiuk, Claudio Kopper, and Artemis Spyrou, for your helpful guidance and suggestions. Thank you to Mehr Un Nisa, who has been super helpful since first arriving at MSU as a postdoc. Thanks for tolerating all my dumb questions and for guiding me through multiple analyses. Thank you to my fellow graduate students and close friends Jessie Micallef, Joe Lundeen, and Sebastian Enrique Sanchez Herrara. You provided a listening ear and many good laughs over the years whenever I needed it. Thank you to the entire HAWC collaboration, whose years of effort made it possible for me to complete the work presented here. Thanks especially to Pat Harding, Kristi Engel, and Henrike Fleischhack for working with me directly on some of the analyses presented herein. Thank you also to Hugo Ayala, Israel Martinez Castellanos, Kelly Malone, and Chad Brisbois for answering my many questions about physics and software. I enjoyed the company of all HAWC members at conferences and collaboration meetings. Thank you to Kim Crosslan and Brenda Wenzlick for helping with navigating paperwork, room reservations, and many other administrative aspects. You were always there to help and made my life a lot easier. To all WaMPS members, thank you for your camaraderie and for letting me serve on your executive board for three years. It was a very rewarding experience. And thank you iii to all the lovely members of the physics choir for providing camaraderie and a break from work each week. Finally, thank you to my parents, Andrew and Natalie Peisker, for providing all the support it took to get to this point. None of this would have been possible without you. iv TABLE OF CONTENTS LIST OF TABLES . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . viii LIST OF FIGURES . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . x CHAPTER 1 INTRODUCTION . . . . . . . . . . . . . . . . . . . . . . . . . . . 1 1.1 Gamma-ray Astronomy . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1 1.1.1 Cosmic Rays . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2 1.1.2 Gamma Rays . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3 1.1.2.1 Pion Decay . . . . . . . . . . . . . . . . . . . . . . . . . . . 4 1.1.2.2 Synchrotron Radiation . . . . . . . . . . . . . . . . . . . . . 5 1.1.2.3 Inverse Compton Scattering . . . . . . . . . . . . . . . . . . 6 1.1.2.4 Bremsstrahlung . . . . . . . . . . . . . . . . . . . . . . . . . 7 1.1.2.5 Extragalactic Background Light . . . . . . . . . . . . . . . . 8 1.2 Sources of Gamma Rays . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 10 1.2.1 Galactic Sources . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 10 1.2.2 Extragalactic Sources . . . . . . . . . . . . . . . . . . . . . . . . . . . 12 1.3 Detection Techniques . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 15 1.3.1 Direct Detection . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 15 1.3.2 Imaging Air Cherenkov Telescopes . . . . . . . . . . . . . . . . . . . 16 1.3.3 Extensive Air Showers . . . . . . . . . . . . . . . . . . . . . . . . . . 18 1.4 Multi-messenger Astronomy . . . . . . . . . . . . . . . . . . . . . . . . . . . 19 1.4.1 Astrophysical Messengers . . . . . . . . . . . . . . . . . . . . . . . . 19 1.4.2 Multi-messenger Sources . . . . . . . . . . . . . . . . . . . . . . . . . 20 CHAPTER 2 THE HAWC OBSERVATORY . . . . . . . . . . . . . . . . . . . . . 24 2.1 Water Cherenkov Detectors . . . . . . . . . . . . . . . . . . . . . . . . . . . 24 2.2 Data Acquisition . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 26 2.3 Monitoring . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 30 2.3.1 ATHENA . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 31 2.3.2 HOMER and Alerts . . . . . . . . . . . . . . . . . . . . . . . . . . . 32 2.4 Event Reconstruction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 35 2.4.1 Angle Reconstruction . . . . . . . . . . . . . . . . . . . . . . . . . . . 35 2.4.1.1 Core Fitting . . . . . . . . . . . . . . . . . . . . . . . . . . . 36 2.4.1.2 Angle Fitting . . . . . . . . . . . . . . . . . . . . . . . . . . 37 2.4.2 Energy Reconstruction . . . . . . . . . . . . . . . . . . . . . . . . . . 39 2.4.3 Background Estimation . . . . . . . . . . . . . . . . . . . . . . . . . . 41 2.4.3.1 Gamma-hadron Separation . . . . . . . . . . . . . . . . . . 41 2.4.3.2 Direct Integration . . . . . . . . . . . . . . . . . . . . . . . 43 CHAPTER 3 HAWC CATALOG SEARCH . . . . . . . . . . . . . . . . . . . . . . 45 3.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 45 v 3.2 Catalog Construction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 46 3.2.1 Data Selection . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 46 3.2.2 Source Search . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 47 3.2.3 Spectral Fitting . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 49 3.3 Results . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 54 3.3.1 3HWC Sources . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 54 3.3.2 Potential Associations with Known Sources . . . . . . . . . . . . . . . 55 3.4 Limitations and Systematic Uncertainties . . . . . . . . . . . . . . . . . . . . 58 3.4.1 Background Fluctuations . . . . . . . . . . . . . . . . . . . . . . . . . 58 3.4.2 Search Methods Limitations . . . . . . . . . . . . . . . . . . . . . . . 59 3.4.3 Systematic Uncertainties . . . . . . . . . . . . . . . . . . . . . . . . . 60 3.5 Conclusions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 63 CHAPTER 4 SEARCH FOR PRIMORDIAL BLACK HOLES . . . . . . . . . . . . 64 4.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 64 4.2 Theory . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 65 4.3 Data Selection . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 68 4.4 Analysis Methods . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 71 4.4.1 Creating the Model . . . . . . . . . . . . . . . . . . . . . . . . . . . . 71 4.4.2 Calculating an Upper Limit . . . . . . . . . . . . . . . . . . . . . . . 73 4.5 Results . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 75 4.6 Conclusions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 77 CHAPTER 5 THE ICECUBE OBSERVATORY . . . . . . . . . . . . . . . . . . . 79 5.1 The IceCube Array . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 79 5.2 Data Acquisition . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 83 5.3 Event Reconstruction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 85 5.3.1 Event Signatures . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 86 5.3.2 Angular Reconstruction . . . . . . . . . . . . . . . . . . . . . . . . . 89 5.3.3 Energy Reconstruction . . . . . . . . . . . . . . . . . . . . . . . . . . 91 5.3.4 Backgrounds to Astrophysical Neutrinos . . . . . . . . . . . . . . . . 92 5.4 Realtime Alert System . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 92 CHAPTER 6 MULTI-MESSENGER SEARCH . . . . . . . . . . . . . . . . . . . . 97 6.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 97 6.2 Analysis Methods . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 98 6.2.1 Data Selection . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 98 6.2.2 Bayesian Blocks . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 100 6.2.2.1 The Algorithm . . . . . . . . . . . . . . . . . . . . . . . . . 101 6.2.2.2 Calibration . . . . . . . . . . . . . . . . . . . . . . . . . . . 104 6.2.2.3 Application to Data . . . . . . . . . . . . . . . . . . . . . . 106 6.2.3 Fake Flare Study . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 107 6.2.4 Simulation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 110 6.2.4.1 Cosmology . . . . . . . . . . . . . . . . . . . . . . . . . . . 110 6.2.4.2 FIRESONG . . . . . . . . . . . . . . . . . . . . . . . . . . . 111 6.2.4.3 Simulating Sources in HAWC . . . . . . . . . . . . . . . . . 113 6.3 Results . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 115 6.4 Conclusions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 120 CHAPTER 7 CONCLUSIONS AND FUTURE OUTLOOK . . . . . . . . . . . . . 123 BIBLIOGRAPHY . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 125 vii LIST OF TABLES Table 2.1: Definitions of the fhit bins, given by the fraction of available PMTs that are triggered during an air shower event. The angular resolution, denoted Ψ68 as the bin containing 68% of the events, improves with larger events [1]. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 40 Table 3.1: List of sources (part 1) found in the 3HWC catalog search and the search in which they were found (point source or source with 0.5◦ , 1.0◦ , or 2.0◦ extensions). Also shown is the nearest TeVCat source if it is within 1◦ of the HAWC source and the distance between them. TeVCat sources are bolded if it is within 0.5◦ of the HAWC source. Secondary sources are indicated by a dagger (†). Table from [2]. . . . . . . . . . . . . . . . 50 Table 3.2: List of sources (part 2) found in the 3HWC catalog search and the search in which they were found (point source or source with 0.5◦ , 1.0◦ , or 2.0◦ extensions). Also shown is the nearest TeVCat source if it is within 1◦ of the HAWC source and the distance between them. TeVCat sources are bolded if it is within 0.5◦ of the HAWC source. Secondary sources are indicated by a dagger (†). Table from [2]. . . . . . . . . . . . . . . . 51 Table 3.3: List of sources (part 1) found in the 3HWC catalog search, the search in which they were found, and their fitted spectra. Each source’s best-fit spectral index and differential flux at 7 TeV is reported. Also shown is the energy range from which 75% of the source’s significance is obtained. The uncertainties displayed are the statistical and systematic uncertainties, respectively. The fit for the source 3HWC J0659+147 did not converge. Table from [2]. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 52 Table 3.4: List of sources (part 2) found in the 3HWC catalog search, the search in which they were found, and their fitted spectra. Each source’s best-fit spectral index and differential flux at 7 TeV is reported. Also shown is the energy range from which 75% of the source’s significance is obtained. The uncertainties displayed are the statistical and systematic uncertainties, respectively. The fit for the source 3HWC J0659+147 did not converge. Table from [2]. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 53 Table 3.5: List of sources found in the 3HWC catalog search that have no known TeV counterpart. For each source, the galactic latitude and longitude are listed along with the sources from the 4FGL [3], ATNF [4], and SNRCat [5] catalogs within 1◦ . The distance from the 3HWC source to the sources from the other catalogs are shown in parentheses. Table from [2]. 57 viii Table 4.1: The 99% upper limits on the PBH burst rate density for the three re- maining lifetimes searched. The uncertainties are systematic only, and are described in Section 4.5. . . . . . . . . . . . . . . . . . . . . . . . . . 76 Table 4.2: The strongest limit on the burst rate density of PBHs set by each of the five detectors most sensitive to PBH burst studies. . . . . . . . . . . . . . 76 Table 6.1: This table summarizes the results from the fake flare study, where flares are simulated and the Bayesian Block algorithm is performed on HAWC’s data with a simulated flare. The results of all combinations of flare intensity in terms of that of the Crab Nebula and flare durations in days are shown here. A check mark (X) represents a detection and an x indicates that the Bayesian Block algorithm did not detect the fake flare. 110 Table 6.2: Summary of the three different types of sources being simulated using the FIRESONG software package. The flux normalization, Φ0 , and spectral index, γ, assuming a power-law energy spectrum are provided for each source. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 112 ix LIST OF FIGURES Figure 1.1: The energy spectrum of cosmic rays follows a power law shape and features a couple of changes in spectral index. Below the knee, the cosmic rays are of galactic origin, and above the knee the cosmic rays originate from extragalactic sources. Figure from [6]. . . . . . . . . . . . 3 Figure 1.2: An illustration of the photon spectral energy distribution and its fea- tures. The shape at low energies is due to synchrotron radiation, while the features at high energies are due to inverse Compton scattering and neutral pion decay. An observation of the pion decay feature indicates a hadronic acceleration mechanism. Figure from [7]. . . . . . . . . . . . 5 Figure 1.3: An illustration of the inverse Compton scattering process. A high- energy electron scatters off an incident photon, which causes the emit- ted photon to increase in energy while the electron loses energy. Figure from [8]. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 7 Figure 1.4: Feynman diagram of the bremsstrahlung process. An electron collides with the electric field of a nucleus and loses energy by emitting a photon. Figure from [9]. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 8 Figure 1.5: This figure shows the attenuation for gamma rays due to the EBL at a variety of redshifts. The attenuation increases as redshift increases, making it more difficult to observe high-energy gamma rays from far extragalactic sources. The different colors represent different models of the EBL. Figure from [10]. . . . . . . . . . . . . . . . . . . . . . . . . . 9 Figure 1.6: Locations of all known TeV gamma-ray sources from TeVCat at the time of writing, in Galactic coordinates. The majority of sources are located along the Galactic plane, which is the horizontal line in the center of the plot. The different colors represent the different types of sources. Figure from [11]. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 10 Figure 1.7: An image of the Crab Nebula taken by the Hubble Telescope. The colors represent the different elements that were expelled during the supernova explosion. Figure from [12]. . . . . . . . . . . . . . . . . . . . . . . . . . 11 Figure 1.8: An illustration of the components of an AGN. AGN have a supermassive black hole that accretes material onto an accretion disk. Some AGN also have two jets pointing in opposite directions along the axis of rotation. Figure from [13]. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 13 x Figure 1.9: An illustration of a space-based gamma-ray detector. The gamma rays are detected via the electron-positron pairs they produce inside the con- verter. The direction and energy of the initial gamma ray are determined using the tracker and calorimeter. Figure from [14]. . . . . . . . . . . . 16 Figure 1.10: A schematic of the geometry of Cherenkov radiation. Cherenkov radi- ation occurs when a particle is traveling faster than the speed of light in a medium. The particle, denoted by the red arrow, will emit blue light along its path, denoted by blue arrows at a Cherenkov angle of θ. Figure from [15]. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 17 Figure 1.11: Illustration of the different astrophysical messenger particles, excluding gravitational waves. Particles are produced in an astrophysical source and encounter obstacles differently on their way to Earth. Figure from [16]. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 21 Figure 1.12: Gamma ray and neutrino emission from a blazar, which is an AGN with its jets oriented towards Earth. Protons are accelerated in the jets via shock acceleration, and gamma rays and neutrinos are among secondary particles produced by pp and pγ interactions. Figure from [17]. . . . . . 22 Figure 2.1: A photograph of the HAWC detector showing the main array and the surrounding outrigger tanks taken in February 2020. The mountain Pico de Orizaba can be seen in the background. . . . . . . . . . . . . . . 25 Figure 2.2: An illustration of an extensive air shower initiated by a gamma-ray inter- acting with Earth’s atmosphere. The gamma-ray produces an electron- positron pair, which in turn produce more gamma rays. Figure from [18]. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 26 Figure 2.3: An illustration of an extensive air shower initiated by a cosmic ray interacting with Earth’s atmosphere. The cosmic ray produces pions, which in turn produce gamma rays, muons, and neutrinos. Figure from [18]. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 27 Figure 2.4: A schematic of a HAWC WCD, with dimensions and the size of a person for scale. The layout of the four PMTs inside the tank can be seen. Also illustrated is the Cherenkov light produced by an air shower particle. Figure from [1]. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 28 Figure 2.5: An illustration of a four-edge event with the signal crossing the high voltage threshold. The time over threshold (ToT) is used later in the reconstruction process. Figure from [18]. . . . . . . . . . . . . . . . . . 29 xi Figure 2.6: A schematic of the HAWC DAQ system. The left side of the figure shows the electronic components, including the PMTs, FBEs, TDCs, and SBCs. The right side shows the software components, including the Reconstruction Client, Analysis Clients, and Experiment Control. Figure from [1]. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 30 Figure 2.7: A schematic of the ATHENA system. The ATHENA server receives measurements from the monitored components and writes it to the on- site database. The ATHENA synchronization copies over the measure- ments to the off-site database. The measurements are displayed in both on-site and off-site versions of HOMER, and are used to send alerts in certain cases. Figure from [19]. . . . . . . . . . . . . . . . . . . . . . . . 32 Figure 2.8: A screenshot of the HOMER main webpage on a day when the detector is operating normally. This page contains essential information to assess the performance of the detector operations. Further monitoring details are located in other pages linked from the main page. . . . . . . . . . . . 33 Figure 2.9: A screenshot of a plot on a HOMER sub-page that contains measure- ments of the temperatures of some of the electronics. In addition to the interactive plot, options to plot a custom time range and to download a printable version are included. . . . . . . . . . . . . . . . . . . . . . . . 34 Figure 2.10: Charge deposited in each PMT for a reconstructed gamma-ray event. Each large circle represents a WCD and each of the 4 smaller circles within represent a PMT. The color scale represents the amount of charge deposited in each PMT. The red star in the center of the dashed circle shows the location of the shower core fit by the SFCF algorithm. . . . . . 37 Figure 2.11: An illustration of the angle reconstruction of the original particle. The secondary particles of an air shower travel in a plane perpendicular to the direction of the original particle, allowing for the reconstruction of the initial angle after corrections due to the curvature of the plane. Figure from [18]. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 38 Figure 2.12: Normalized histogram of the energy distribution of each fhit bin defined in Table 2.1. This figure was made using Monte Carlo simulations as- suming a spectral distribution ∼ E −2.63 at a declination of 20◦ . Figure from [1]. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 40 Figure 2.13: The lateral distribution of charge for a cosmic ray event (left) and a gamma ray event (right). Cosmic ray showers contain more PMT hits with high charges that are far away from the shower core, while the lateral distribution is smoother for gamma-ray events. Image from [1]. . 42 xii Figure 3.1: Significance map of the full sky within HAWC’s field of view with 1523 days of data. The bright band on the left side is the Galactic plane, the two spots in the middle are Markarian 501 and Markarian 421, and the bright spots on the right are the Geminga pulsar halo and the Crab Nebula. This map was made assuming a point source hypothesis. Figure from [2]. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 48 Figure 3.2: Significance map of part of the Galactic plane assuming a point √ source hypothesis. The green lines show√significance contours starting at T S = 26 and increasing in steps of ∆ T S = 2. The top labels indicate the positions of sources from TeVCat and the bottom labels indicate the locations of 3HWC sources. Figure from [2]. . . . . . . . . . . . . . . . 54 Figure 3.3: Significance map of part of the Galactic plane assuming a point source hypothesis. The top labels indicate the positions of sources from TeVCat and the bottom labels indicate the locations of 3HWC sources. Figure from [2]. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 55 Figure 3.4: Significance map of part of the Galactic plane assuming a point √ source hypothesis. The green lines show√significance contours starting at T S = 26 and increasing in steps of ∆ T S = 2. The top labels indicate the positions of sources from TeVCat and the bottom labels indicate the locations of 3HWC sources. Figure from [2]. . . . . . . . . . . . . . . . 56 Figure 3.5: Significance map of part of the Galactic plane assuming a point √ source hypothesis. The green lines show√significance contours starting at T S = 26 and increasing in steps of ∆ T S = 2. The top labels indicate the positions of sources from TeVCat and the bottom labels indicate the locations of 3HWC sources. Figure from [2]. . . . . . . . . . . . . . . . 57 Figure 3.6: All-sky significance map with randomized background. This map was created by Poisson-fluctuating the number of background events in each pixel. The full catalog search was performed on twenty such maps in order to estimate the amount of spurious sources. . . . . . . . . . . . . . 59 Figure 3.7: Significance maps of the region containing the Crab Nebula (bottom right) and the Geminga pulsar (center). On the left is the map produced with a point source hypothesis and on the right is the map produced with a 1◦ extended source hypothesis. Several point sources were found near the Geminga pulsar in the 3HWC search while the map on the right shows a single extended source in the region, illustrating the limitations to the search method. Figure from [2]. . . . . . . . . . . . . . . . . . . . 60 xiii Figure 3.8: Comparisons between the declinations of HAWC’s sources and their counterparts measured by IACT experiments H.E.S.S., MAGIC, and VERITAS. The difference between HAWC’s measured declination and that from the IACTs’ for each source are plotted on the y-axis as a function of HAWC’s declination on the x-axis. HAWC’s measurements of the source locations agree within uncertainties for most of its decli- nation range except for sources whose zenith angles are larger than 30◦ . Error bars are HAWC’s statistical uncertainty. Figure from [2]. . . . . . 61 Figure 4.1: Time-integrated energy spectrum of photons emitted during black hole evaporation over remaining lifetime intervals τ = 0.2, 1, and 10 seconds. Adapted from [20]. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 67 Figure 4.2: Emission time profile of a bursting PBH integrated over energy in the range 50 GeV - 100 TeV. This profile can be described by a power law with an index of -0.5. Figure from [20]. . . . . . . . . . . . . . . . . . . 68 Figure 4.3: Distributions of the p-values of the data, background, PBH, and model histograms for the three remaining lifetimes used in this search. See Sec- tion 4.4.1 for a description of HP BH and Hmodel . The vertical dashed black line indicates pthr , which is the passing threshold for data to be included in this analysis. Note that Hbkg and HP BH were drawn from a much larger sample than Hdata and were therefore scaled in order to match its duration. Figure from [21]. . . . . . . . . . . . . . . . . . . . 70 Figure 4.4: Locations of simulated PBH burst events. The HAWC detector is lo- cated at the center, and axes are in units of parsecs. PBHs were sim- ulated out to 0.5 pc within 50◦ of HAWC’s zenith, resulting in a cone centered around HAWC’s field of view. . . . . . . . . . . . . . . . . . . 72 Figure 4.5: Example signal expected at HAWC (µ) from a simulated PBH burst as a function of zenith angle (θ) using HAWC’s code ZEBRA. This particular simulated PBH is assumed to have a remaining lifetime of 0.2 s and located at a distance of 0.05 pc from HAWC. The signal plateaus above 45◦ because of HAWC’s low sensitivity at those zenith angles. Figure from [21]. . . . . . . . . . . . . . . . . . . . . . . . . . . 73 Figure 4.6: Comparison of HAWC’s upper limits at the 99% confidence level at burst durations 0.2, 1, and 10 s with the upper limits from Whipple [22], CYGNUS [23], the Tibet Air Shower Array [24], H.E.S.S. [25], VERI- TAS [26], Fermi -LAT [27], and Milagro [28]. The expected limits and the 68% and 95% containment bands are also shown. Systematic un- certainty is not shown since it was significantly smaller than the con- tainment bands. Figure from [21]. . . . . . . . . . . . . . . . . . . . . . 77 xiv Figure 4.7: The difference in upper limits obtained by using different detector mod- els corresponding to different sources of systematic uncertainty. The black stars represent the final result for each burst duration searched and the red band represents the total systematic uncertainty obtained by adding each effect in quadrature. Figure from [21]. . . . . . . . . . . 78 Figure 5.1: A schematic of the IceCube detector, showing the in-ice array, Deep- Core, IceTop, and the IceCube Laboratory [29]. The Eiffel Tower is shown for scale. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 80 Figure 5.2: A schematic showing the components of a Digital Optical Module (DOM). Each DOM consists of a PMT, an LED flasher board, and other elec- tronics to readout the events, all contained within a glass sphere. Figure from [30]. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 81 Figure 5.3: A picture of the IceCube Laboratory (ICL) that contains the data acqui- sition system and online computers for the experiment. Figure from [31]. 82 Figure 5.4: Layout of the IceCube Upgrade and the optical and radio components of IceCube-Gen2 compared to the current IceCube detector. Figure from [32]. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 83 Figure 5.5: Diagram of data flow starting from a PMT and ending with information about the signal being sent to the surface. If a PMT signal crosses the discriminator threshold and so does the signal from neighboring DOMs, the waveform is digitized and sent to the surface. Figure from [30]. . . . 84 Figure 5.6: Feynman diagrams for charged current interactions (a/b) and neutral current interactions (c/d) between neutrinos and quarks that make up nucleons in the ice. A W boson is exchanged in charged current inter- actions and produces an outgoing lepton corresponding to the neutrino flavor along with hadronic cascades. A Z boson is exchanged in neutral current interactions and only hadronic cascades are produced. Figure from [33]. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 87 Figure 5.7: A track event in IceCube. A muon neutrino initiates a charged current interaction and produces a muon that creates the track signature in the detector. The colors correspond to the timing of the PMT hits; red is early and green is late. The size of the PMT hits correspond to the amount of charge deposited; more charge is represented by larger spheres. Figure from [33]. . . . . . . . . . . . . . . . . . . . . . . . . . . 88 xv Figure 5.8: A cascade event in IceCube. Cascade events are produced by either neutral current interactions with any neutrino flavor, charged current interactions involving an electron neutrino, or charged current interac- tions involving a tau neutrino below about 1 PeV. Figure from [33]. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 88 Figure 5.9: A double cascade event in IceCube. A tau neutrino above 1 PeV initiates a charged current interaction and produces a hadronic shower and a tau lepton, which then decays and creates a second shower. This produces a double cascade signature in the detector. Figure from [34]. . . . . . . 89 Figure 5.10: This illustration shows neutrinos originating from different angles that interact in the IceCube detector. Down-going cosmic rays come from the Southern hemisphere and create atmospheric muons and neutrinos that are the main background to astrophysical neutrinos. Up-going cosmic rays come from the Northern hemisphere and the muons they produce are filtered out by the Earth. The signal being searched for in this case are the cosmic neutrinos that originate from astrophysical sources. Figure from [33]. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 93 Figure 6.1: IceCube realtime neutrino alert locations on top of HAWC’s significance skymap from the 3HWC catalog search [2]. The IceCube alerts are spaced out evenly across HAWC’s field of view. Eighty-five of the 99 total alerts issued by IceCube are shown here as those are the ones that fall within HAWC’s field of view. . . . . . . . . . . . . . . . . . . . . . . 99 Figure 6.2: The daily flux from data collected by HAWC between March 2015 and June 2019. Any gaps in the data are due to the HAWC detector being down for maintenance. . . . . . . . . . . . . . . . . . . . . . . . . . . . 100 Figure 6.3: Calibration of the ncp prior parameter used in the Bayesian Blocks al- gorithm. Here a false positive rate of 5% is used, which corresponds to an ncp prior value of 6.75. . . . . . . . . . . . . . . . . . . . . . . . . . 105 Figure 6.4: The Bayesian Block algorithm applied to HAWC flux data at the loca- tion of an example IceCube neutrino alert. In this case only one block is generated because the flux is fairly constant over the duration of the data set. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 106 Figure 6.5: A significance map of an example fake flare that was simulated on top of HAWC’s data. This example flare has a duration of one day and a flux normalization equal to that of the Crab Nebula. This flare achieves a detected significance of around 6σ, which is close to HAWC’s expected daily significance from the Crab Nebula of 5σ. . . . . . . . . . . . . . . 108 xvi Figure 6.6: Two example fake flares injected into HAWC’s data following the pro- cedure described above. The black points are HAWC’s data with the fake flare replacing the data during its length and the orange line shows the Bayesian Blocks. The flare is detected in the bottom case (b) but not in the top (a). . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 109 Figure 6.7: Example of a simulated flare injected into HAWC’s data. The Bayesian Block algorithm was run on the data including the injected flare, and it detected the simulated flare in this case. The black points are HAWC’s data with the simulated flare replacing the data during its length and the orange line shows the Bayesian Blocks. This example was a simulated source with the spectrum of TXS 0506+056 whose flare lasted for ten days. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 115 Figure 6.8: Sensitivity, discovery potential, and 90% upper limits for the AGN source type as a function of rate density and total isotropic energy equivalent in neutrinos. Each flare duration is shown separately. Results corresponding to a density value of 10−6 Mpc−3 yr−1 are only calculated for 1-day flares due to computational limitations. . . . . . . . . . . . . . 116 Figure 6.9: Sensitivity, discovery potential, and 90% upper limits for the diffuse source type as a function of rate density and total isotropic energy equivalent in neutrinos. Each flare duration is shown separately. Results corresponding to a density value of 10−6 Mpc−3 yr−1 are only calculated for 1-day flares due to computational limitations. . . . . . . . . . . . . . 117 Figure 6.10: Sensitivity, discovery potential, and 90% upper limits for the TXS 0506+056 source type as a function of rate density and total isotropic energy equiv- alent in neutrinos. Each flare duration is shown separately. Results corresponding to a density value of 10−6 Mpc−3 yr−1 are only calculated for 1-day flares due to computational limitations. . . . . . . . . . . . . . 118 Figure 6.11: Comparison of the sensitivity, discovery potential, and 90% upper limit for a 1-day flare of all source types tested. The vertical lines correspond to rate densities of different source types. [35–38]. . . . . . . . . . . . . 119 Figure 6.12: Comparison of the two procedures of calculating 90% upper limits for a 1-day flare. The limits from the previously published method described in Section 6.3 are shown in the dark green lines and the limits from the Feldman Cousins method described in Section 6.4 are shown in the bright green lines. The sensitivity and discovery potential of all source types tested are also included in red and blue lines, respectively. The vertical lines correspond to rate densities of different source types. [35–38]. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 121 xvii CHAPTER 1 INTRODUCTION The study of particle astrophysics can help to understand the high-energy universe where the most extreme events occur. There are several different mechanisms that produce particles in a variety of astrophysical objects. These particles span a large range of energies and can be emitted on time scales from a fraction of a second to many years. There is a lot that still remains to be understood about the high-energy universe. Al- though cosmic rays were first discovered over 100 years ago, there is still no known single source responsible for their emission. It is also still unknown how the cosmic accelerators can accelerate their particles up to energies as high as have been observed on Earth. There are several different messenger particles which each have their advantages and disadvantages. Studies that use more than one of these messengers are called multi-messenger astronomy. Using a combination of particles can help to answer some of the questions that still remain about our universe. In this thesis, gamma-ray astrophysics will be used to identify astrophysical objects emitting photons in the energy range of 100 GeV to 100 TeV. It will also be used to search for primordial black holes. Finally, multi-messenger astronomy will be used to search for transient emission in both gamma rays and neutrinos. These studies seek to contribute to understanding some of the outstanding questions about the universe. 1.1 Gamma-ray Astronomy Before gamma-ray astronomy came the detection and study of cosmic rays. Cosmic rays are high-energy charged particles that originate from sources in outer space and interact with Earth’s atmosphere. They were first discovered by Victor Hess in 1912 during his balloon experiments. Hess found that the level of radiation increased as his balloon increased in altitude, thus concluding that there were particles reaching Earth from outer space [39]. 1 1.1.1 Cosmic Rays Cosmic rays consist primarily of protons, but also heaver nuclei such as helium, carbon, nitrogen, oxygen, and iron, among others. They also consist of electrons about 2% of the time. Of the 98% of remaining cosmic rays that reach Earth, 87% are protons, 12% are helium, and 1% are heavier nuclei [6]. Different experiments have detected cosmic rays over many orders of magnitude in energy, ranging from 1 GeV to over 10 EeV. The energy spectrum of cosmic rays can be seen in Figure 1.1. The cosmic ray energy spectrum follows a power law shape with an index of -2.7 up to about 100 TeV. It has several features known as the knee, which is a change in spectral index around 107 GeV, and the ankle, which is another change in spectral index that occurs around 1010 GeV. Cosmic rays in energies up to the knee in the spectrum are believed to originate from galactic sources, and cosmic rays above the knee are believed to originate from extragalactic sources [7]. The spectrum cuts off around 1020 eV due to the Greisen-Zatsepin-Kuzmin (GZK) cutoff [40,41], which is caused by collisions between cosmic rays and the cosmic microwave background. To reach such high energies, cosmic rays are accelerated in extreme astrophysical envi- ronments in a process known as diffusive shock acceleration. In this mechanism, which was first proposed by Fermi in 1949, charged particles are deflected multiple times through shock fronts in magnetic clouds [42]. Shock fronts are abrupt changes in pressure, temperature, and density of a material medium. The particles gain energy each time they cross the shock front, which results in a power-law spectrum with index of about 2. Diffusive shock acceler- ation combined with cosmic-ray propagation models yields a spectral index compatible with the observed value of 2.7 [7]. Since cosmic rays are charged particles, they will be deflected by interstellar magnetic fields on the way from their sources to Earth. This makes it difficult to determine the sources of cosmic rays. However, processes that produce cosmic rays are also expected to produce gamma rays. Gamma rays are electrically neutral particles and therefore are not deflected 2 Figure 1.1: The energy spectrum of cosmic rays follows a power law shape and features a couple of changes in spectral index. Below the knee, the cosmic rays are of galactic origin, and above the knee the cosmic rays originate from extragalactic sources. Figure from [6]. by magnetic fields. Since they point back to the sources in which they were created, gamma rays are a good messenger particle with which to study these sources. 1.1.2 Gamma Rays The observation of gamma rays may yield clues about the origins of cosmic rays. If the cosmic rays interact with their environment immediately after being produced, they can generate gamma rays. Studying gamma ray production mechanisms can provide informa- tion about the origin and propagation of cosmic rays. There are two types of gamma-ray production mechanisms, hardonic and leptonic, based on the type of primary interacting par- ticles. Neutral pion decay is the main hadronic process that produces high-energy gamma rays, while there are three different leptonic production mechanisms: synchrotron radiation, inverse Compton scattering, and bremsstrahlung. 3 1.1.2.1 Pion Decay The main hadronic production mechanism is pair production from neutral pion decay. When protons and nucleons collide in an astrophysical environment, they produce neutral and charged pions, among other particles. The neutral pions decay into a pair of gamma rays via the following process, π 0 → γ + γ. (1.1) Using units where the speed of light c = 1, each photon has an energy of Eγ = 67.5 MeV, which is half of the pion’s mass, in the rest frame of the pion. When boosted into the laboratory frame, the photon energy can extend up to several orders of magnitude higher, making it possible to achieve TeV gamma rays from neutral pion decay. The photon energy in the laboratory frame is Eγ = ΓEγ∗ + βΓp ∗γ · cos θ∗ , (1.2) where pγ is the photon momentum, θ is the angle between the photons, and the ∗ represents Eπ |p | the quantities in the pion rest frame. Also, Γ = m π , β = Eπ , where mπ , Eπ , and pπ are π the mass, energy, and momentum of the pion, respectively. Differentiating this expression with respect to cos θ∗ yields the energy spectrum of the photons in the laboratory frame, dN 1 Eγ2 = , (1.3) dEγ 2βΓp∗γ when considering an E −2 proton energy spectrum at the source. The photon energy spectrum from neutral pion decay rises steeply below about 200 MeV and traces the energy distribution of parent protons at energies above a few GeV. This spectral feature is known as the pion- decay bump as it uniquely identifies the presence of gamma rays that originate from neutral pion decays. This feature can be seen on the right side of the plot in Figure 1.2. An observation of this feature would identify the presence of high-energy protons and allow a measurement of the cosmic-ray spectrum directly at the source [7]. 4 Figure 1.2: An illustration of the photon spectral energy distribution and its features. The shape at low energies is due to synchrotron radiation, while the features at high energies are due to inverse Compton scattering and neutral pion decay. An observation of the pion decay feature indicates a hadronic acceleration mechanism. Figure from [7]. 1.1.2.2 Synchrotron Radiation Synchrotron radiation occurs when relativistic electrons are traveling in a magnetic field. An electron with energy E will produce electromagnetic radiation at the critical frequency νc , which can be expressed as  2 E νc = Γ2 νg = νg , (1.4) me c2 where me is the mass of the electron and νg is the electron’s gyrofrequency, eB νg = , (1.5) 2πme c where e is the charge of the electron and B is the strength of the magnetic field [7]. Charged particles emit electromagnetic radiation with their power given by the relativistic Larmor formula, dE 2e2 P =− = 2 a2 , (1.6) dt 3c 5 where a is the particle’s acceleration. Substituting in the relativistic acceleration assuming a circular orbit and using the definition of the Larmor radius, the energy loss rate can be written as dE 4 B2 − = σT cΓ2 , (1.7) dt 3 8π where σT is the Thomson cross-section. Assume that the electron energies follow a power-law distribution with index αe and normalization κ, dN dE = κE −αe dE. (1.8) dE Then using the definition of νg and the relation Eγ = hν where h is Planck’s constant, the flux density can be expressed as dNγ    dE dN Eγ 2 =− dE = (constants) · Eγ−α · B α+1 , (1.9) dEγ dt dE where the gamma ray spectral index is α = (αe − 1)/2. The spectral energy distribution also depends on the magnetic field as B α+1 [7]. Synchrotron radiation produces photons in the infrared and x-ray range, which can be seen on the left side in Figure 1.2. 1.1.2.3 Inverse Compton Scattering The relativistic electrons that produce synchrotron emission can scatter off photons, which causes the photons to gain energy as the electron loses energy. This is known as inverse Compton scattering (IC), which is illustrated in Figure 1.3. Similar to synchrotron radiation, inverse Compton scattering results in a gamma-ray spectral energy distribution with an index α = (αe − 1)/2 when electron energies are in the Thomson regime. When the electrons are ultra-relativistic, the Klein-Nishina cross section must be used instead, π 2 re2   1 σKN = ln 2hν + , (1.10) hν 2 where h is the Planck constant, re is the radius of the electron, and ν is the photon frequency. In this regime the spectral energy distribution becomes dNγ Eγ2 ∝ Eγ−αe ln Eγ . (1.11) dEγ 6 Figure 1.3: An illustration of the inverse Compton scattering process. A high-energy electron scatters off an incident photon, which causes the emitted photon to increase in energy while the electron loses energy. Figure from [8]. In this case, the index of the photon spectral energy distribution is the same as that of the parent electrons [7]. This corresponds to photons in the TeV energy range, as can be seen on the right side of Figure 1.2. 1.1.2.4 Bremsstrahlung The final leptonic production mechanism is bremsstrahlung, or braking radiation. This process occurs when a high-energy electron emits a photon when it loses energy due to colliding with an electric field of another charged particle that is present in the interstellar gas. Figure 1.4 provides an illustration of this process. The energy loss rate of electrons due to bremsstrahlung is cmp n   dE − = E, (1.12) dt X0 where mp is the mass of the proton, n is the number density of the ambient gas, E is the energy of the electron, and X0 is the radiation length, or the distance over which the electron’s energy falls to 1/e of its original value. The lifetime of the electron due to bremsstrahlung is E τ= ≈ 4 × 107 (n/cm−3 )−1 yr, (1.13) −dE/dt 7 Figure 1.4: Feynman diagram of the bremsstrahlung process. An electron collides with the electric field of a nucleus and loses energy by emitting a photon. Figure from [9]. using the radiation length of hydrogen gas X0 = 60g/cm2 . Since the energy loss rate is proportional to the energy of the electron, the spectral shape of the initial power-law distribution of electrons remains the same after bremsstrahlung. This results in a gamma- ray spectrum that also follows a power law shape with the same spectral index [7]. 1.1.2.5 Extragalactic Background Light When observing gamma rays originating from extragalactic sources, attenuation from the Extragalactic Background Light (EBL) must be considered. The EBL is background ra- diation in the ultraviolet to infrared range, and consists of remnant photons emitted by stars, galaxies, and Active Galactic Nuclei (AGN) over the lifetime of the universe. Very- high-energy gamma rays produced in astrophysical sources interact with the EBL via pair production: γE + γ → e+ + e− , (1.14) where γE is the very-high-energy photon and γ is an EBL photon. This process attenuates the flux from very-high-energy gamma rays that originate from extragalactic objects at far distances from Earth. The gamma-ray flux from the source will be attenuated according to 8 Figure 1.5: This figure shows the attenuation for gamma rays due to the EBL at a variety of redshifts. The attenuation increases as redshift increases, making it more difficult to observe high-energy gamma rays from far extragalactic sources. The different colors represent different models of the EBL. Figure from [10]. the following, F (E, z) = F0 e−τγγ (E,z) , (1.15) where F (E, z) is the observed flux, F0 is the flux at the source, and τγγ is the optical depth of the source, which depends on the energy E and redshift z. Since the attenuation increases with redshift, the EBL becomes a larger effect as distance to the source increases. This causes the flux from gamma rays to start being attenuated around 1 TeV and makes it difficult to observe any gamma rays above a few tens of TeV for sources at redshifts of 0.03, and the energy threshold lowers as the redshift of sources increases [7]. Figure 1.5 shows the gamma-ray attenuation due to the EBL for a variety of redshift values. 9 Figure 1.6: Locations of all known TeV gamma-ray sources from TeVCat at the time of writing, in Galactic coordinates. The majority of sources are located along the Galactic plane, which is the horizontal line in the center of the plot. The different colors represent the different types of sources. Figure from [11]. 1.2 Sources of Gamma Rays There are many types of astrophysical objects that can emit gamma rays up to TeV energies. The TeVCat is a catalog of all observed TeV gamma-ray sources, of which there are a total of 276 at the time of writing [11]. Figure 1.6 shows where in the sky these sources are located in Galactic coordinates. These objects can be divided into two categories depending on their location: galactic sources and extragalactic sources. They may also be classified as either steady sources, which are constantly emitting gamma rays, or transient sources, which only emit gamma rays for short amounts of time. 1.2.1 Galactic Sources Supernova remnants (SNR) are the leftover components of stellar explosions. When a massive star (M > 8M ) reaches the end of its evolution lifecycle, it collapses in a powerful explosion. The explosion releases ∼ 1051 erg of kinetic energy that forms a shock wave. Part of the energy of this shock wave is converted into high energy particles, which can accelerate cosmic 10 Figure 1.7: An image of the Crab Nebula taken by the Hubble Telescope. The colors represent the different elements that were expelled during the supernova explosion. Figure from [12]. rays via diffusive shock acceleration up to the knee in the cosmic ray spectrum. The cosmic rays then interact with the surrounding material, which produces gamma rays that follow a power-law energy spectrum. Particles are accelerated while the shock front expands across the interstellar matter, which is typically around 1000 years [7]. The Crab Nebula is one of the most prominent supernova remnants in gamma-ray as- tronomy. Its explosion was first recorded in the year 1054 by Chinese astronomers and it is currently used as a reference source for TeV astronomy as it is the brightest source in the gamma-ray sky. The Crab Nebula and the pulsar at its center are well-studied in almost all wavelength bands of the electromagnetic spectrum, from radio to > 100 TeV gamma rays. Figure 1.7 shows an image of the Crab Nebula from NASA’s Hubble Space Telescope. Pulsars are rotating neutron stars that emit a beam of electromagnetic radiation along 11 their magnetic axis. Neutron stars are created during supernova explosions and have masses of around 1.4M . A pulsar’s radiation can only be observed when its emitted radiation is pointed towards the Earth. The rotation period and interval between observed pulses is very regular, usually occurring between 1.4 ms to 8.5 s. The rotation slows down over time as electromagnetic radiation is emitted [7]. The magnetic fields of pulsars are typically very high around 1012 G, which provides a region where charged particles can be accelerated to very high energies. The acceleration from the curvature of co-rotating electric field lines causes the charged particles to emit synchrotron radiation. Some pulsars have an extended nebula emitting radiation; these are known as pulsar wind nebulae (PWNe). In PWNe, a wind of energetic particles from a pulsar carries most of the rotational power into the surrounding region. Electrons are accelerated during interactions with the interstellar medium. Photons from pulsars and PWNe can reach energies up to the TeV range via synchrotron radiation and inverse Compton scattering [7]. 1.2.2 Extragalactic Sources Active Galactic Nuclei (AGN) are galaxies with a supermassive black hole at their center that accretes material. AGN are typically classified according to their orientation with respect to Earth, the presence of jets, and their radiative efficiency. Figure 1.8 provides an illustration of an AGN showing its black hole, accretion disk, and jets. About 10% of AGN have jets, which are highly collimated beams of relativistic particles. If the jet is pointed towards the observer, the AGN is classified as a blazar; if the jet is angled away from the observer, the AGN is classified as a radio galaxy. AGN without jets are classified as Seyfert galaxies, which are radio-quiet but show emission in the optical range [7]. There are two categories of blazars: BL Lacertae (BL Lacs) and flat spectrum radio quasars (FSRQs). The main distinction is that FSRQs show emission lines in their optical spectrum while BL Lacs do not. Blazars are typically categorized into classes based on their emission. If their emission peaks in the radio wavelengths they are classified as low-energy 12 Figure 1.8: An illustration of the components of an AGN. AGN have a supermassive black hole that accretes material onto an accretion disk. Some AGN also have two jets pointing in opposite directions along the axis of rotation. Figure from [13]. BL Lacs (LBL), and if their emission peaks in the x-ray wavelengths they are classified as high-energy BL Lacs (HBL). Most of the AGN detected in very-high-energy gamma rays are of the HBL type [7]. The low-energy emission from blazars is due to syncrotron emission from relativistic electrons and the high-energy emission is related to the inverse Compton emission of radiation fields. High-energy gamma ray emission may also be explained by a hadronic acceleration model if they are accelerated by shocks in the jets. The emission of AGN can be highly variable. The origin of their variability is still an open question [7]. Gamma-Ray Bursts (GRBs) are extremely intense but short-lived bursts of gamma ray emission. They are the most energetic transient explosions observed in the Universe. GRBs occur isotropically across the sky and are of extragalactic origin. Their emission spectra 13 follow a broken power law shape and range from a few hundred keV up to a few GeV, and there have recently been hints of GRBs emitting up to TeV energies [43]. While GRBs typically only last for a few seconds, they also display an afterglow that can last for days in other wavelengths. This allows for their host galaxies to be identified [7]. GRBs are classified into two categories by their duration: short GRBs last for an average of 0.3s and long GRBs have average durations around 30s. The spectrum of short GRBs is typically harder than that of long GRBs - there are more high-energy photons and the power law index is less negative. Long GRBs are thought to originate in supernovae explosions in active star forming galaxies. Short GRBs likely originate from a merging of two compact objects such as neutron star mergers in galaxies that contain many old stars [7]. While the progenitors of short and long GRBs differ, the mechanism that produces their bursts are the same. This mechanism is known as the fireball model [44]. The progenitor causes an explosion that sends relativistic blast waves through the star. Two opposite jets form where shells collide and accelerate particles through diffusive shock acceleration. The gamma rays are most likely emitted through synchrotron radiation, although inverse Compton scattering may also play a role. After the initial shock, the blast wave pushes the stellar material through the star surface, which collides with external gas and dust. This produces the additional photon emission that is described as the afterglow [7]. The final source type that will be discussed in this thesis is Primordial Black Holes (PBHs). These objects have not yet been discovered, but were proposed by Hawking in 1974 [45]. They are predicted to thermally emit particles with a temperature inversely proportional to their mass, ending with a final burst of particles when they evaporate. PBHs formed with initial masses of ∼ 5 × 1014 g should be expiring today with a burst of gamma rays in the GeV - TeV energy range. These objects will be discussed in more detail in Chapter 4. 14 1.3 Detection Techniques There are two general categories of methods to detect gamma rays: direct and indirect. Direct detection is when a gamma ray directly interacts with a detector, and indirect detec- tion is when the presence of a gamma ray is inferred from secondary particles. Detectors for direct gamma-ray emission must be located above the atmosphere while indirect detectors are located on the ground. 1.3.1 Direct Detection In order to directly detect gamma rays, it is necessary to go above Earth’s atmosphere. This is accomplished by sending detectors to space on board satellites. The detectors are made up of several layers in order to estimate the direction and energy of the gamma rays they detect. First, charged particles are vetoed by a plastic scintillator that rejects the background from cosmic rays. Gamma rays are detected via the electron positron pairs they produce via pair production inside a converter made of material that has a high atomic number. A tracker device is used to record the paths of the electrons and positrons, which are used to reconstruct the arrival direction of the gamma ray. The electrons and positrons also pass through scintillator detectors in order to provide timing information, which also helps with reconstructing the arrival direction. Finally, the electrons and positrons enter a calorimeter, which is used to determine the energy of the initial gamma rays that produced the pairs [7]. Figure 1.9 provides an illustration of the components used to detect and reconstruct gamma rays. An early example of a space-based instrument is the Energetic Gamma Ray Experiment Telescope (EGRET), an instrument on the Compton Gamma-Ray Observatory (CGRO). EGRET operated from 1991 to 2000, observing gamma rays in the energy range of 20 MeV to 30 GeV. During this time, EGRET was able to detect a total of 271 gamma-ray objects [46]. The Large Area Telescope (LAT) on board the Fermi satellite has been operating since 15 Figure 1.9: An illustration of a space-based gamma-ray detector. The gamma rays are detected via the electron-positron pairs they produce inside the converter. The direction and energy of the initial gamma ray are determined using the tracker and calorimeter. Figure from [14]. 2008, detecting gamma rays from 20 MeV to above 300 GeV. Its improved electronics and technology improvements over EGRET provide a larger effective area and field of view, allowing it to detect thousands of gamma-ray sources [3]. The Fermi satellite also carries the Gamma-ray Burst Monitor (GBM) to detect gamma-ray transients. Space-based telescopes are advantageous for the detection of gamma rays because they have a wide field of view and good angular resolution. However, above a few hundred GeV, the gamma-ray fluxes are so small that the effective detection area of space-based instruments is too small to detect enough gamma rays. The study of gamma rays above this energy requires ground-based detectors, such as Imaging Air Cherenkov Telescopes or Extensive Air Shower Arrays. 1.3.2 Imaging Air Cherenkov Telescopes Imaging Air Cherenkov Telescopes (IACTs) are telescopes made of mirrors that focus light to a camera made out of many photomultiplier tubes (PMTs). When a gamma-ray interacts with Earth’s atmosphere, it produces a shower of secondary particles, as will be explained in Section 2.1. If the particles are traveling faster than the speed of light in a medium, in 16 Figure 1.10: A schematic of the geometry of Cherenkov radiation. Cherenkov radiation occurs when a particle is traveling faster than the speed of light in a medium. The particle, denoted by the red arrow, will emit blue light along its path, denoted by blue arrows at a Cherenkov angle of θ. Figure from [15]. this case air, they will emit Cherenkov light. Figure 1.10 shows an illustration of the geometry of Cherenkov radiation. If the particle satisfies the following condition, c vp > , (1.16) n where vp is the velocity of the particle, c is the speed of light, and n is the refractive index of the medium, then it will emit blue light at a Cherenkov angle of θ. Since the emitted light waves travel at speed c/n, they travel a distance of nc t, where t is time. The particle travels a distance of βct, where β = vp /c. Therefore the emission angle is 1 cos θ = . (1.17) nβ In air, the Cherenkov angle is around 0.7◦ , depending on the altitude since n changes with atmospheric depth. This results in light from a gamma-ray shower focusing onto the ground in an ellipse-like region with radius of ∼ 120 m. The long axis of the ellipse corresponds to the vertical extension of the shower and points back towards the position of the source. 17 The energy of the initial gamma ray is roughly proportional to the intensity of the measured Cherenkov light [7]. Some of the currently operating IACTs include H.E.S.S. [39], MAGIC [47], and VERITAS [48]. They all are made up of an array of multiple Cherenkov telescopes, which provides better background rejection and reconstruction due to the stereoscopic view of the air shower. The Chernkov Telescope Array (CTA) is the next generation IACT instrument. It will be composed of more than 100 Cherenkov telescopes and will therefore be about 10 times more sensitive than currently existing instruments [49]. IACTs have great angular resolution reaching less than 0.1◦ , high sensitivity, and achieve energy resolution to about 15%. However, their observation times are limited to dark moon- less, cloudless nights, and their fields of view are limited to ∼ 4◦ . Extensive air shower arrays are needed to increase the field of view and observation time of gamma-ray events [7]. 1.3.3 Extensive Air Showers Extensive Air Shower (EAS) arrays consist of many detectors on the ground that cover a large area. They are typically located at high altitudes in order to detect the shower maximum of extensive air shower particles. The energy range for this type of detector is typically between 100 GeV to 100 TeV. EAS arrays can use a few different techniques to detect the air shower particles. Early experiments such as ARGO-YBJ [50] and Tibet-AS [51] used resistive plate chambers and scintillation detectors, respectively, as their detection technique. Milagro [52] and HAWC [1] detect Cherenkov light emitted by the secondary shower particles, using water as the medium. The recently built LHAASO [53] array uses a combination of techniques to achieve even better sensitivity. The pattern and amount of light on the ground created by the secondary particles allow for the reconstruction of the arrival direction and energy of the original particle. More details on the reconstruction techniques for HAWC, an EAS array, will be given in Section 2.4. While 18 EAS arrays do not have the same angular reconstruction as IACTs, their advantages include a wide field of view and being able to operate nearly continuously. 1.4 Multi-messenger Astronomy The era of multi-messenger astronomy has recently arrived with many observatories op- erating that study the different messenger particles. There has been a large increase in coordinated efforts between the observatories to study the same objects simultaneously, par- ticularly if one of them observes an interesting transient event. There are advantages and disadvantages to studying each messenger, but combining them may help to answer some of the remaining questions about the universe. 1.4.1 Astrophysical Messengers There are four astrophysical messenger particles: cosmic rays, gamma rays, neutrinos, and gravitational waves. Cosmic rays, discussed in Section 1.1.1, are produced in a variety of astrophysical objects. However, since they are charged particles, they are deflected by the interstellar magnetic fields on their way to Earth. This makes it difficult to identify and study the sources in which they are produced. Other messenger particles, such as electrically-neutral gamma rays or neutrinos, point directly back to the sources where they originate. Gamma rays, discussed in Section 1.1.2, are produced when cosmic rays interact either near their sources or during their propagation. They can be produced in either hadronic or leptonic processes. As neutral particles they point directly back to their sources, and they are also easy to detect via the techniques discussed in Section 1.3. However, their interactions with the EBL attenuate their flux, making it difficult to study sources that are far away from Earth. Neutrinos are produced in the same sources as cosmic rays and gamma rays, but they are only produced in hadronic processes. Neutrinos are also neutral particles and therefore 19 point directly back to their sources. Since neutrinos only interact via the weak force, they are not affected by the EBL and therefore can travel long distances from their sources to Earth. However, their weak interactions also make them very difficult to detect, requiring detectors with a large volume. One such neutrino detector, IceCube, will be described in Chapter 5. The final messenger particle is gravitational waves. They are produced by the acceleration of large masses in events such as black hole mergers and neutron star mergers. They are only affected by gravitational interactions, so they travel directly from their sources to Earth. Gravitational waves are detected via km-scale Michelson laser interferometers such as the LIGO [54] and Virgo [55] observatories. While each messenger has different properties that affect their detection at Earth, they may be produced in the same astrophysical objects. Studying them together can help to reveal the sources of cosmic rays and neutrinos along with their acceleration mechanisms, which is difficult to determine by studying a single messenger. Chapter 6 will discuss a multi- messenger search using gamma rays and neutrinos. Figure 1.11 provides an illustration of some of the advantages and disadvantages of the different messenger particles. Gravitational waves are not included in this figure and will not be discussed further in this thesis. 1.4.2 Multi-messenger Sources Neutrinos (ν) can be produced via interactions of protons with either gas (pp interactions) or with photons (pγ interactions). In the hadronic scenario described in Section 1.1.2, the interaction of the protons and nucleons within a source produce charged pions (π) in addition to the neutral pions that decay into gamma rays. These charged pions then decay into leptons and neutrinos via the following processes, π + → µ + + νµ (1.18) 20 Figure 1.11: Illustration of the different astrophysical messenger particles, excluding grav- itational waves. Particles are produced in an astrophysical source and encounter obstacles differently on their way to Earth. Figure from [16]. and µ+ → e+ + νe + ν¯µ , (1.19) where µ+ is a muon, e+ is a positron, and ν̄ is an antineutrino. Therefore both gamma rays and neutrinos are expected to be produced if the source accelerates its particles hadronically. Since gamma rays can be produced via both leptonic and hadronic processes, the observation of neutrinos along with gamma rays can point to a hadronic acceleration mechanism. While a number of sources, including supernovae, pulsar wind nebulae, and gamma ray bursts, are predicted to produce neutrino emission, active galactic nuclei are particularly interesting sources to study. AGN, especially blazars, are some of the most powerful objects and comprise the majority of gamma-ray sources in the extragalactic sky [3]. Protons are ac- 21 Figure 1.12: Gamma ray and neutrino emission from a blazar, which is an AGN with its jets oriented towards Earth. Protons are accelerated in the jets via shock acceleration, and gamma rays and neutrinos are among secondary particles produced by pp and pγ interactions. Figure from [17]. celerated in their relativistic jets via shock acceleration, producing gamma rays and neutrinos in pγ and pp interactions. So far only gamma rays have been discovered from AGN, and there are many different hadronic and leptonic models to describe these objects [56]. Therefore detection of neutrinos from AGN would definitively identify their acceleration mechanism as hadronic. Figure 1.12 provides an illustration of the particle production process from a blazar. The promise of blazars as a target to study multi-messenger astronomy can also be seen from the follow-up from IceCube’s neutrino alert on September 22, 2017. When IceCube reported the detection of a high-energy neutrino, a large campaign ensued to perform follow- up observations. These observations spanned many wavelengths from radio, optical, and x-ray up to very-high-energy gamma rays performed by over a dozen observatories. Gamma- ray emission was found by Fermi -LAT and MAGIC, resulting in the identification of the source of the neutrino, the blazar TXS 0506+056 [57]. With IceCube’s growing list of neutrino alerts and HAWC’s continuous operations collecting more gamma ray data, there is the potential to identify more blazars exhibiting multi-messenger emission. This would 22 help to not only identify the acceleration mechanism in AGN, but also to identify sources of neutrinos. 23 CHAPTER 2 THE HAWC OBSERVATORY The High Altitude Water Cherenkov (HAWC) Gamma-ray Observatory is located in the state of Puebla, Mexico at an altitude of 4100 m above sea level. HAWC detects gamma rays in the energy range from a few hundred GeV to over 100 TeV. It has an instantaneous field of view of 2 sr and can see 2/3 of the sky each day. HAWC is located at a latitude of 19◦ North and its visible sky ranges from −26◦ to 64◦ in declination. It also has a high duty cycle and is running > 95% of the time. These qualities make it an ideal instrument to conduct searches for transient sources, like the ones described in this work. A picture of the HAWC array is shown in Figure 2.1. 2.1 Water Cherenkov Detectors When a gamma ray interacts with Earth’s atmosphere, it produces an extensive air shower of lower-energy secondary particles. These secondary particles are mainly produced through bremsstrahlung and pair production. A high-energy gamma ray interacts with a nucleus in the atmosphere and forms an electron-positron pair via pair production. Then these secondary particles produce more gamma rays through bremsstrahlung. The air shower continues to propagate until the energy of the secondary particles decreases beyond the ionization threshold, at which point fewer particles are produced [58]. An illustration of this process in shown in Figure 2.2. Cosmic rays also produce extensive air showers, and are the main background for HAWC observations. Cosmic rays can be protons or heavier nuclei. When a cosmic ray collides with a nucleus in Earth’s atmosphere, it produces pions and other nuclear fragments. The neutral pions decay into two gamma rays, which will in turn produce electron-photon cascades. The charged pions decay into muons and neutrinos [58]. These additional particles allow for the discrimination between gamma-ray and cosmic-ray showers by detectors. An illustration of 24 Figure 2.1: A photograph of the HAWC detector showing the main array and the surround- ing outrigger tanks taken in February 2020. The mountain Pico de Orizaba can be seen in the background. a cosmic-ray shower is shown in Figure 2.3. If the secondary air shower particles are traveling faster than the speed of light in a medium, they will produce Cherenkov radiation. This blue Cherenkov light is detected by photo-multiplier tubes (PMTs) in HAWC, which uses water as the medium. HAWC consists of 300 Water Cherenkov detectors (WCDs), which are cylindrical tanks filled with water that each have four upward-facing PMTs anchored to the bottom. Each tank is 7.3 m in diameter and 5 m high and contains about 200,000 L of purified water. There are three 8-inch PMTs arranged around one central 10-inch PMT in each WCD. The array of tanks covers a total area of about 22,000 m2 [1]. A schematic of a WCD with its PMTs is shown in Figure 2.4. A building located in the center of the array, the counting house, holds the front-end electronics, computers, and other equipment necessary for the functioning of the experiment. 25 Figure 2.2: An illustration of an extensive air shower initiated by a gamma-ray interacting with Earth’s atmosphere. The gamma-ray produces an electron-positron pair, which in turn produce more gamma rays. Figure from [18]. The HAWC detector was upgraded with a secondary, sparser array known as the outrig- gers that surround the main array. The outrigger array is composed of 345 smaller tanks that are 1.55 m in diameter and 1.55 m tall. They each contain a single upward-facing eight-inch PMT anchored to the bottom. The outrigger tanks can be seen in Figure 2.1. The instrumented area is about 4x bigger than the main array, and is expected to improve HAWC’s sensitivity above 10 TeV by providing an improved reconstruction of showers whose cores are not well-contained within the main array [59]. However, at the time of writing, the outrigger array has not yet been incorporated into the reconstruction, and therefore none of the algorithms or data from the outriggers are included in this work. 2.2 Data Acquisition The PMTs in each tank are connected by RG-59 cables to the Front-End Boards (FEBs), which are located in the counting house. The cables are all equal in length (610 ft) in order to preserve the signal timing. The FEBs process the signals from the PMTs and then the Time 26 Figure 2.3: An illustration of an extensive air shower initiated by a cosmic ray interacting with Earth’s atmosphere. The cosmic ray produces pions, which in turn produce gamma rays, muons, and neutrinos. Figure from [18]. to Digital Converters (TDCs) digitize those signals. They are also responsible for providing the high voltage to the PMTs. There are a total of 10 TDC modules which each have 128 input channels [60]. When the signal from a PMT enters an FEB, it is compared to two different pre-set voltage thresholds. The low threshold corresponds to a charge of about 1/4 photoelectrons deposited in the PMT, and the high threshold corresponds to 4 photoelectrons. This pro- duces a signal with either two or four transitions depending on whether the voltage reaches the high threshold. The time over threshold (ToT) of the PMT signals is measured by a TDC module [60]. This is illustrated in Figure 2.5. The information about the edge times is used to reconstruct the timing and amount of charge deposited by the original PMT signals. 27 Figure 2.4: A schematic of a HAWC WCD, with dimensions and the size of a person for scale. The layout of the four PMTs inside the tank can be seen. Also illustrated is the Cherenkov light produced by an air shower particle. Figure from [1]. The TDCs are synchronized by a GPS Timing and Control (GTC) system that has a timing resolution of 98 ps. The data from the TDCs are transferred to the Single Board Computers (SBCs) via the Versa Module Europa (VME). The SBCs are connected via an Ethernet network to computers that process data in real time and record the filtered data on disk. The online processing system uses ZeroMQ [61] to transfer data between software components at a speed of about 60 MB/s [60]. The online system runs the Reconstruction Client, which reconstructs air shower informa- tion from the raw PMT data. This allows for performing a few online analyses that require prompt information. These analyses are run by Analysis Clients and include searches for Gamma Ray Bursts (GRBs), a flare monitor, and follow-ups of gravitational waves, and 28 Figure 2.5: An illustration of a four-edge event with the signal crossing the high voltage threshold. The time over threshold (ToT) is used later in the reconstruction process. Figure from [18]. neutrinos reported by other experiments. All other analyses are performed using the full reconstruction chain, which is run offsite [60]. The various processes running on the data acquisition system (DAQ) are controlled by a system known as Experiment Control. Experiment Control is responsible for starting and stopping runs, and for obtaining data needed to monitor the experiment. If an error in running the experiment causes it to crash, Experiment Control will attempt to automati- cally restart the run. [60]. The main cause of these errors is lightning, but recovery from these crashes occurs typically on the order of minutes. Figure 2.6 shows a diagram of all components of the HAWC DAQ and online processing system [60]. 29 Figure 2.6: A schematic of the HAWC DAQ system. The left side of the figure shows the electronic components, including the PMTs, FBEs, TDCs, and SBCs. The right side shows the software components, including the Reconstruction Client, Analysis Clients, and Experiment Control. Figure from [1]. 2.3 Monitoring In order to maintain HAWC’s high duty cycle, a monitoring system is needed to oversee the remote operation of the detector’s hardware. This system allows off-site shifters to mon- itor components of the detector’s hardware and to alert experts to problems with operation. It also allows experts in turn to diagnose issues in order to solve them. This is accom- plished using two frameworks: a system to collect data and write it to a database known as Advanced Tracking of HAWC Experiment Notifications and Alerts (ATHENA), and a front-end webpage to display the information known as HAWC Observatory Monitoring for Experiment and Reconstruction (HOMER). 30 2.3.1 ATHENA The main infrastructure of ATHENA can be divided into three major components: the on- site measurements, the ATHENA server, and the ATHENA database synchronization. The on-site measurements contain a variety of systems to be monitored. Examples of components that are monitored include PMT rates, high voltage, the temperature of electronics, the load and storage capacity of on-site computers, the water level in the WCDs, and the main transformer responsible for providing power to HAWC. The system was designed so that more monitored data can be easily added at any point in time [19]. The ATHENA server is a Python library that receives data from on-site measurements, which occurs about once every minute. It is primarily responsible for writing the monitored information into a database using ZeroMQ [61]. Each monitoring stream results in a different table in the MySQL database [62]. Approximately 8 GB of data are written to the database per year. Two copies of the database exist: one on the computers on-site and one on the off-site computers located at the University of Maryland. The database synchronization is responsible for copying all of the data from the on-site computers to the off-site ones. This process occurs about once every minute and must be efficient due to network bandwidth constraints. Therefore, the synchronization script checks both the on-site and off-site versions of the database and only copies over the most recent entries to the off-site database. About 20 Kb of data are sent during each synchronization instance in order to keep the off-site database up to date. Occasionally power outages or network outages occur which cause the synchronization to need to wait before transferring data off-site. When the network comes back up, the synchronization process will then attempt to transfer all of the most recent data off-site. The synchronization script automatically times out and restarts after about one minute to prevent problems with a slow table from interrupting the entire monitoring process. When longer outages occur, a separate script can be used to move data off-site in smaller packets and catch up the monitoring more efficiently. When the off-site database is back up to date, 31 Figure 2.7: A schematic of the ATHENA system. The ATHENA server receives measure- ments from the monitored components and writes it to the on-site database. The ATHENA synchronization copies over the measurements to the off-site database. The measurements are displayed in both on-site and off-site versions of HOMER, and are used to send alerts in certain cases. Figure from [19]. the main synchronization script can be used again. The measurements located in the on-site database are displayed on a simple webpage on-site and are also used to alert the collaboration to detector operation errors. The mea- surements located in the off-site database are displayed on a webpage off-site that are used by shifters and experts to monitor the experiment. These additional components are described in Section 2.3.2. A diagram of the on-site and off-site components of ATHENA, HOMER, and the alerts can be seen in Figure 2.7. 2.3.2 HOMER and Alerts HOMER is the collection of webpages where the measurements from the ATHENA database are displayed in near real-time. It is generated with the PHP language for preprocessing 32 Figure 2.8: A screenshot of the HOMER main webpage on a day when the detector is operating normally. This page contains essential information to assess the performance of the detector operations. Further monitoring details are located in other pages linked from the main page. HTML [63] and uses built-in SQL wrapper libraries to access the ATHENA database. It is used by off-site shifters to monitor the detector daily and also by experts to diagnose any problems with detector operations. The main page contains the most critical information to the operation of the detector. It is displayed in such a way in order to easily determine whether the detector is running smoothly or if it needs intervention by an expert. The information on the main page includes whether the detector is running and for how long, whether the high voltage is on or off, the PMT rates, the temperatures of some of the electronic elements, and the amount of space available on the disks at site. These critical measurements are color-coded according to their values so that the user can easily tell if the detector is operating normally or not. A screenshot of this page can be seen in Figure 2.8. The rest of the monitoring information that contains more details is accessible by tabs from the main page. There is a separate page for each of the monitored components listed in Section 2.3.1. They contain interactive plots of the measurement versus time that are 33 Figure 2.9: A screenshot of a plot on a HOMER sub-page that contains measurements of the temperatures of some of the electronics. In addition to the interactive plot, options to plot a custom time range and to download a printable version are included. created using Google Charts API [64], which allows users to get the exact value of a data point by dragging their cursor over it. Users can also input a custom time range to see how the detector was performing at different points in time. Each page includes information about when it was last updated so that the user can see how recent the measurements were made. An example of one of these monitoring pages can be seen in Figure 2.9. There is also a version of HOMER on the on-site computers that has a few different features from the off-site version. It does not use the Google Charts API but uses Python scripts to make the plots instead. This allows the service to be more light-weight and for it to work even if the network connection is down. It runs on the same computer as the on-site 34 ATHENA database so the data is displayed in real time as it does not need to wait for a database synchronization. It is meant to be used by the experienced site crew to troubleshoot the experiment and contains buttons to turn on or off high voltage channels. This feature does not exist in the off-site version because the high voltage control should not be accessible to non-expert users. The final component of HAWC’s remote monitoring is the alert system. This system polls certain electronics about once a minute to ensure they are operating in a normal state. If a measurement is outside of the set normal bounds, a message is sent to the HAWC Slack workspace, where the off-site monitoring shifters are watching. This allows experts to be notified quickly in case of a problem with the detector operations. Such alerts include a change in the AC voltage or high voltage, a high temperature warning, a file size monitor, a data transfer monitor, and a change in the run status. 2.4 Event Reconstruction In order to make use of the data that HAWC collects, information about the air shower created by the original particle interacting with Earth’s atmosphere needs to be reconstructed from the PMT measurements. HAWC’s reconstruction algorithms use the timing of hits and the amount of charge deposited to reconstruct the direction and energy of the original particle and to determine whether the particle was a gamma ray or a cosmic ray. This information can then be used to perform a wide variety of science analyses, from studying galactic and extra-galactic sources to cosmic ray composition and searches for dark matter and other beyond the Standard Model phenomena. 2.4.1 Angle Reconstruction There are two steps involved to reconstruct the direction of the original particle. The first is identifying the location of the core of the air shower, and the second is determining the direction of the particle that initiated the air shower. 35 2.4.1.1 Core Fitting The first step in reconstructing an air shower is to determine its core, or the location on the ground where the most charge is deposited. This corresponds to the location within the shower where the concentration of secondary particles is the highest, which occurs along the trajectory of the original primary particle. This core location is then used by several successive reconstruction algorithms. The first core reconstruction algorithm uses a simple center-of-mass calculation to obtain a rough estimate of the core location. This algorithm computes the center of mass by taking the average PMT locations and weighting them by their recorded charges. This first guess is used to seed the next fitting algorithm [1]. HAWC’s main core fitting algorithm is known as the Super Fast Core Fitter (SFCF). Figure 2.10 shows the charge deposited in each PMT for a sample event and the fit location of the shower core. The SFCF algorithm fits the lateral charge distribution to a modified Nishimura-Kamata-Greisen (NKG) function:   1 2 N Si (A, r) = A 2 e−ri /2σ + , (2.1) 2πσ (0.5 + ri /Rm )3 where Si is the signal in the ith PMT, σ = 10 m, N = 5 × 10−5 , and Rm = 120 m is the Moliere radius of the atmosphere at HAWC’s altitude. This function has a Gaussian component with width σ to parameterize the charge distribution at the core and a tail component with normalization N , whose values were chosen from HAWC simulations. The Moliere radius, defined as the radius of a cylinder containing 90% of the energy deposition of the air shower, is determined from the radiation length in air, the scattering energy, the critical energy, and the air density at the height of observation. The free parameters to be fit are the core location r and the amplitude A. The full NKG function [65, 66] includes an additional free parameter for the shower age and contains power-law and gamma function evaluation which is computationally intensive. The simplified equation allows the 36 Figure 2.10: Charge deposited in each PMT for a reconstructed gamma-ray event. Each large circle represents a WCD and each of the 4 smaller circles within represent a PMT. The color scale represents the amount of charge deposited in each PMT. The red star in the center of the dashed circle shows the location of the shower core fit by the SFCF algorithm. minimization to converge faster because it has no pole at the center and the derivatives can be computed analytically [1]. The median error of the SFCF is approximately 2 m for large events and 4 m for smaller events for cores which land on the HAWC array. The error increases for events whose cores land off the array. For example, the error is about 35 m for a shower with a core that is located 50 m from the edge of the array [1]. 2.4.1.2 Angle Fitting After the location of the core of the air shower is fitted, the direction of the original particle can be determined. This is accomplished by using the relative timing information of all the 37 Figure 2.11: An illustration of the angle reconstruction of the original particle. The secondary particles of an air shower travel in a plane perpendicular to the direction of the original particle, allowing for the reconstruction of the initial angle after corrections due to the curvature of the plane. Figure from [18]. PMT hits involved in the detection of the air shower. However, there are two effects that need to be accounted for first. One of the effects that needs to be corrected is the shower curvature. As can be seen in Figure 2.11, the secondary particles of the shower travel to the ground at the speed of light on a curved plane. However, not all of the secondary particles arrive on the ground at the same time. The particles whose paths are along the shower axis have less distance to travel so they will reach the ground first, and the particles further away from the core will arrive later. The second effect, known as sampling, is due to variation of arrival times at ground level. This is caused by the scattering of shower particles. The time delay between the earliest and 38 latest particles in the shower is about 5 to 10 ns near the axis and larger further away from the core. The HAWC electronics only measure the arrival time of the earliest particle, so a position-dependent timing correction is needed to overcome this effect [1]. These two effects are corrected together in the reconstruction as they are hard to disen- tangle. This timing correction is based on a combination of simulation and data since the simulation cannot successfully reproduce the subtleties of the PMT timing on its own. After the corrections are applied, a χ2 function is used to fit a plane to the times of the PMT hits. The angle can be reconstructed to a resolution of between 0.1◦ and 1◦ , depending on the energy and zenith angle of the primary particle [1]. 2.4.2 Energy Reconstruction HAWC uses a variable known as fhit as a simple estimator of the energy of the primary particle that produced the air shower. This variable is defined as the ratio of PMTs partici- pating in the reconstruction (nhit ) of the air shower to the total number of available PMTs in the reconstruction (NChAvail ). Typically around 1000 PMTs are available in any given air shower event as some may be temporarily undergoing maintenance. This variable can be used as an estimator since higher energy particles produce a larger shower footprint on the ground and therefore trigger more PMTs. The angular resolution also improves with energy since there are more PMTs available to use in the reconstruction. The events are divided into nine bins; the bin definitions and their corresponding angular resolutions can be seen in Table 2.1 [1]. This method of estimating energy has a few limitations. The variable fhit is dependent on the spectral assumption of the observed source as well as the zenith angle of the primary particle. Additionally, the detector becomes saturated above about 10 TeV, making it diffi- cult to distinguish the energy of particles above this energy. Further, the energy distribution of each bin is wide, as can be seen in Figure 2.12. To address some of the limitations of the fhit energy estimator, two independent energy 39 Bin Lower Edge (%) Upper Edge (%) Ψ68 (deg) 1 6.7 10.5 1.03 2 10.5 16.2 0.69 3 16.2 24.7 0.50 4 24.7 35.6 0.39 5 35.6 48.5 0.30 6 48.5 61.8 0.28 7 61.8 74.0 0.22 8 74.0 84.0 0.20 9 84.0 100.0 0.17 Table 2.1: Definitions of the fhit bins, given by the fraction of available PMTs that are trig- gered during an air shower event. The angular resolution, denoted Ψ68 as the bin containing 68% of the events, improves with larger events [1]. B=1 1.0 B=2 B=3 B=4 B=5 B=6 0.8 B=7 B=8 N Events - Scaled B=9 0.6 0.4 0.2 0.0 -2 10 10-1 100 101 102 103 log10(E/TeV) Figure 2.12: Normalized histogram of the energy distribution of each fhit bin defined in Table 2.1. This figure was made using Monte Carlo simulations assuming a spectral distribution ∼ E −2.63 at a declination of 20◦ . Figure from [1]. 40 estimation algorithms have been developed by the HAWC collaboration [67]. These esti- mators use the reconstructed zenith angle direction and the charge distribution of the air shower footprint on the ground around the shower core, allowing for a more precise estima- tion of the energy of the primary particle, particularly above 10 TeV. However, these energy reconstruction methods were not used in the analyses presented in this work and therefore will not be discussed further. 2.4.3 Background Estimation Over 99% of the events that HAWC records are due to hadronic cosmic rays [1]. In order to perform gamma-ray astronomy, it is first necessary to remove this overwhelming background from the data. This is accomplished using gamma-hadron separation algorithms. The re- maining background is then estimated using a method known as direct integration. These two methods are described in Sections 2.4.3.1 and 2.4.3.2, respectively. 2.4.3.1 Gamma-hadron Separation As explained in Section 2.1, air showers caused by cosmic rays and gamma rays contain different characteristics, which cause their footprint on the ground to look different. Since cosmic ray showers contain muons that have high transverse momentum, they will cause PMTs far from the shower core to contain high amounts of charge. In contrast, gamma-ray showers have smoother and more compact distributions in the detector. Two variables have been created to take advantage of these features in order to distinguish between the two types of particles, compactness and PINCness (Parameter for Identifying Nuclear Cosmic rays). The first parameter, compactness (C), relates to the number of PMTs hit far from the shower core. It is defined as Nhit C= , (2.2) CxP E40 41 4 SFCF Fit 4 SFCF Fit PINC Movi g Average <ζ > PINC Movi g Average <ζ > Qeff Qeff 3 3 log10(Qeff) log10(Qeff) 2 2 1 1 0 0 −1 −1 0 20 40 60 80 100 120 140 0 20 40 60 80 100 120 140 PMT Dista ce to Reco structed Core [m] PMT Dista ce to Reco structed Core [m] Figure 2.13: The lateral distribution of charge for a cosmic ray event (left) and a gamma ray event (right). Cosmic ray showers contain more PMT hits with high charges that are far away from the shower core, while the lateral distribution is smoother for gamma-ray events. Image from [1]. where Nhit is the number of PMTs hit during the air shower and CxP E40 is the highest charge observed in a PMT more than 40 m from the reconstructed shower core. Since the muons in a cosmic ray shower will deposit charge far from the shower core, cosmic rays have smaller values of compactness than gamma rays [1]. The second variable, PINCness (P), aims to quantify the smoothness of the distribution of charge around the shower axis. It is defined as N 1 X (Qi − hQi i)2 P= 2 , (2.3) N σQ i=0 i where N is the number of hit PMTs, Qi is the logarithm of the effective charge measured by the ith PMT, hQi i is the average charge of PMTs within a 5 m radius of the shower core, and σQi is the error on the logarithm of the charge as measured from a sample of likely gamma-ray events from the Crab Nebula. As can be seen in Figure 2.13, there are more PMT hits with high amounts of charge further away from the shower core in cosmic ray showers than in gamma ray showers. Therefore, cosmic rays will have larger values of PINCness than gamma rays [1]. The optimal cuts on PINCness and compactness are determined separately for each fhit 42 bin since they depend on the shower size. The goal is to filter out as many cosmic rays as possible while also trying to leave most of the gamma rays. The cuts are determined using Monte Carlo simulations and observations from the Crab Nebula. With these optimized cuts, 99% of hadrons are rejected from the data at the highest energies [1]. 2.4.3.2 Direct Integration Even though the gamma-hadron separation cuts remove about 99% of hadrons from the data, the majority of remaining events are still cosmic rays. This is because there are about 104 times as many cosmic rays than gamma rays arriving at HAWC, depending on the energy. Therefore, it is necessary to estimate the amount of remaining background in the data. This is accomplished using a method known as direct integration, first developed by the Milagro experiment [68]. This method relies on the isotropic arrival directions of cosmic rays, which is due to their deflection by magnetic fields on their way to Earth. It also assumes that the detector is stable over the 2 hour integration time period. To estimate the background at a given location in right ascension and declination, Nexp (ra, δ), the all-sky rate is convolved with the normalized local direction probability distribution for a given location in the sky, as follows: Z Z Nexp (ra, δ) = E(ha, δ)R(t)dt, (2.4) where E(ha, δ) is the detector efficiency in local coordinates (hour angle and declination) and R(t) is the event rate of the detector as a function of time. This calculation is performed separately for each fhit bin and each declination band, since HAWC’s sensitivity depends on declination. An additional smoothing of the background over a radius of 0.5◦ is necessary since some of the bins in which the integration is performed may be sparse. This may occur at high fhit bins, near the edges of the field of view, or in the case of transients when only a short time period of data is being analyzed. To avoid bias 43 in the calculation, events from known gamma-ray sources like the Galactic plane, the Crab Nebula, and Mrk 421 and Mrk 501, are not included in the background estimation [1]. 44 CHAPTER 3 HAWC CATALOG SEARCH HAWC’s high duty cycle and wide field of view as described in Chapter 2 make it an ideal instrument to perform a catalog search. This involves conducting a blind search through every pixel in HAWC’s visible sky and compiling a list of gamma-ray sources in its energy range that are significantly detected. Then we can compare the list of sources that HAWC sees to lists from other gamma-ray detectors in order to learn more about those sources. This analysis is a follow-up of the one presented in [69], and the results were published in [2]. 3.1 Introduction The HAWC detector was completed in 2015 and has been surveying the Northern gamma- ray sky since then with only small interruptions in observations. HAWC’s first all-sky catalog of gamma-ray sources using 507 days of full detector operations, 2HWC, was published in 2017 [69]. The first catalog looking at the Galactic plane and using data from the partially completed detector was published in 2016 [70]. This analysis is HAWC’s third catalog of very-high-energy gamma-ray sources, 3HWC, and it uses 1523 days of data, which is three times the exposure of the previous catalog [2]. The goal is to identify all significantly detected gamma-ray sources within HAWC’s field of view, fit each source’s energy spectrum, and look for potential associations with other catalogs of gamma-ray sources. Finding new astrophysical objects contributes to answering a few fundamental physics questions. It can help us to study the origin of cosmic rays, understand the mechanisms with which the objects accelerate their particles, and constrain physics beyond the standard model such as dark matter. It also allows for follow-up from other detectors to perform multi-wavelength and multi-messenger studies of the sources in order to learn more about their behavior. 45 3.2 Catalog Construction 3.2.1 Data Selection In this analysis, HAWC data collected between November 2014 and June 2019 is used. Before the catalog search can be performed, it is first necessary to transform the raw data from HAWC’s PMTs into a usable format. First, the methods in Section 2.4 are used to reconstruct the direction and energy of the primary particle that created the air shower. The events are then distributed into the nine analysis bins shown in Table 2.1, which allows us to fit the energy spectrum of the sources. Gamma-hadron separation cuts are applied in each bin to reduce the background due to cosmic rays, as described in Section 2.4.3.1. We also use the method of direct integration, described in Section 2.4.3.2, in order to estimate the amount of remaining background events in the data. The events are then binned further into spatial bins using the Hierarchical Equal Area isoLatitude Pixelization (HEALPix) binning scheme [71]. This scheme divides a spherical surface into pixels of equal areas that are centered at a given right ascension and declina- tion. They are characterized by the parameter Nside , where the total number of pixels is Npix = 12 × Nside 2 . Here we use an N side value of 1024, which corresponds to a pixel size of approximately 0.05◦ . This is chosen to be smaller than HAWC’s angular resolution. The angular resolution is the uncertainty on the direction of the primary particle, which is known as the point spread function (PSF). In HAWC the PSF is modeled by a sum of two Gaussian fits. The PSF depends on the zenith angle of the primary particle and the fhit bin. It improves with smaller zenith angles because the secondary particles travel through less atmosphere if coming from directly overhead, and with larger fhit bins since there are more PMTs hit ant therefore more information when reconstructing event direction. 46 3.2.2 Source Search A blind search is performed over the entire sky that falls within HAWC’s field of view. A likelihood value is calculated for each pixel in the visible sky. This calculation assumes that the counts in each pixel are Poisson-distributed with the mean equal to the estimated background counts plus the expected number of gamma-ray signal counts. The log-likelihood value for an assumed signal is given by Nbins X log(L) = log(Li ), (3.1) i (Si + Bi )Ni e−(Si +Bi ) log Li = , (3.2) Ni ! where Si is the expected signal counts in bin i, Bi is the estimated background counts in bin i, and Ni is the total number of observed counts in bin i. Here there are both the spatial bins and the fhit bins. To compare signal and background hypotheses in each pixel, a log likelihood ratio test is used. A test statistic is defined as follows: Ls+b   T S = 2 log , (3.3) Lb where Lb corresponds to the background only hypothesis and Ls+b corresponds to the source hypothesis. To estimate the signal counts, HAWC uses a forward folding method. Simulation is used to create an expected detector response, which is then used to transform an assumed true energy spectrum into a reconstructed spectrum. Here a simple power law energy spectrum is assumed as follows:  α dN E = F0 , (3.4) dE 7 TeV where F0 is the flux normalization, E is the energy, 7 TeV is the pivot energy, and α is the spectral index. For a first pass where we are simply searching for sources, the spectral index 47 All-sky view; 0.0 ; 1523 days 360 0 -3 -0 3 6 9 12 15 18 21 24 TS Figure 3.1: Significance map of the full sky within HAWC’s field of view with 1523 days of data. The bright band on the left side is the Galactic plane, the two spots in the middle are Markarian 501 and Markarian 421, and the bright spots on the right are the Geminga pulsar halo and the Crab Nebula. This map was made assuming a point source hypothesis. Figure from [2]. α is fixed to −2.5 and the flux normalization is the only free parameter. This parameter is fit by maximizing the TS value, using the likelihood framework described in [72]. Since the source model and null model are nested, the TS follows a χ2 distribution √ with one degree of freedom (since we only have one free parameter). Therefore, T S can be interpreted as the significance of the rejection of the null model, according to Wilks’ theorem since we have sufficiently large statistics [73]. The significance is allowed to be negative when the best-fit flux normalization is negative, or when the background fluctuates √ above the signal. Sources are considered to be significantly detected when T S > 5. The resulting significance map assuming a point source after fitting the TS in each pixel across the sky can be seen in Figure 3.1. This likelihood fitting process is repeated for point sources and extended disk-like sources with radii of 0.5◦ , 1.0◦ , and 2.0◦ . To compile the catalog, we first make a list of all the point 48 sources whose significance values are > 5σ. Then sources with a 0.5◦ extension that do not appear in the point source list are added. Then the sources with 1.0◦ and 2.0◦ extensions are added in the same manner to complete the list of sources. Extended sources are only included if they are more than 1.0◦ from any source that is already in the list. In some cases, several local TS maxima may occur very near to each other due to Poisson fluctuations of the signal. In order to avoid double counting these sources, a TS valley criterion is imposed on the catalog search. A source is considered a primary source if it is √ separated from all other sources within 5◦ by a valley of ∆( T S) > 2. A source is considered √ a secondary source if it is separated by 1 < ∆( T S) < 2 from all other sources within 5◦ . Secondary sources are denoted by a dagger (†) in Tables 3.1 and 3.2. Tables 3.1 and 3.2 show the resulting list of sources and the corresponding search in which they were found that makes up the catalog. The nearest known TeV source from the TeVCat catalog of known TeV sources [11] is included if it is within 1.0◦ from the HAWC source. Also shown is the location of each source in celestial coordinates, the statistical position uncertainty, the test statistic, and the distance from the nearest TeVCat source. 3.2.3 Spectral Fitting After all of the sources are identified, their energy spectra can be fitted. Again a simple power law spectrum is assumed, but this time both the spectral index and flux normaliza- tion are fit for each source. To do this, a piece of HAWC software known as ZEnith BAnd Response Analysis (ZEBRA) is used [74]. ZEBRA uses simulation to characterize the de- tector response and angular resolution as a function of zenith angle. A separate instrument response is calculated for each combination of fhit bin and zenith band, and for each tested spectral assumption during the fitting procedure. This process allows for a more complete characterization of the PSF since HAWC’s response depends strongly on zenith angle. ZEBRA can be used to estimate the number of counts observed from a source during an arbitrary period of time. Here it is used to estimate the counts from the catalog sources over 49 Nearest TeVCat Source Name Radius TS R.A. Decl. 1σ Stat. Unc. Dist. Name (deg) (deg) (deg) (deg) (deg) 3HWC J0534+220 0.0 35736.5 83.63 22.02 0.06 0.01 Crab 3HWC J0540+228† 0.0 28.8 85.17 22.87 0.11 0.77 HAWC J0543+233 3HWC J0543+231 0.0 34.2 85.78 23.11 0.10 0.29 HAWC J0543+233 3HWC J0617+224 0.0 32.3 94.39 22.47 0.14 0.17 IC 443 3HWC J0621+382 0.5 28.0 95.32 38.21 0.30 14.56 ... 3HWC J0630+186 0.0 38.9 97.69 18.68 0.10 1.18 ... 3HWC J0631+107 0.0 26.5 97.78 10.73 0.09 3.84 ... 3HWC J0631+169† 0.0 39.1 97.95 16.96 0.17 0.44 Geminga 3HWC J0633+191 0.0 27.5 98.44 19.12 0.28 1.35 ... 3HWC J0634+067 0.5 36.2 98.66 6.73 0.22 0.28 HAWC J0635+070 3HWC J0634+165† 0.0 30.0 98.53 16.53 0.11 0.92 Geminga 3HWC J0634+180 0.0 60.0 98.75 18.05 0.09 0.38 Geminga Pulsar 3HWC J0659+147† 0.0 28.8 104.85 14.75 0.12 0.50 2HWC J0700+143 3HWC J0702+147 0.0 28.9 105.56 14.75 0.19 0.60 2HWC J0700+143 3HWC J1104+381 0.0 3025.3 166.11 38.16 0.06 0.04 Markarian 421 3HWC J1654+397 0.0 227.8 253.52 39.74 0.09 0.05 Markarian 501 3HWC J1739+099 0.0 28.2 264.99 9.93 0.06 2.87 ... 3HWC J1743+149 0.0 25.9 265.82 14.94 0.10 4.61 ... 3HWC J1757-240 1.0 28.6 269.30 -24.09 0.48 0.74 HESS J1800-240B 3HWC J1803-211† 0.0 38.4 270.97 -21.18 0.27 0.54 HESS J1804-216 3HWC J1809-190 0.0 264.8 272.46 -19.04 0.08 0.31 HESS J1809-193 3HWC J1813-125 0.0 51.9 273.34 -12.52 0.14 0.17 HESS J1813-126 3HWC J1813-174 0.0 416.0 273.43 -17.47 0.06 0.18 2HWC J1814-173 3HWC J1819-150† 0.0 93.8 274.79 -15.09 0.14 0.05 2HWC J1819-150* 3HWC J1825-134 0.0 2212.5 276.46 -13.40 0.06 0.00 2HWC J1825-134 3HWC J1831-095 0.0 237.7 277.87 -9.59 0.14 0.31 HESS J1831-098 3HWC J1837-066 0.0 1542.7 279.40 -6.62 0.06 0.06 2HWC J1837-065 3HWC J1843-034 0.0 876.6 280.99 -3.47 0.06 0.24 2HWC J1844-032 3HWC J1844-001† 0.0 33.2 281.07 -0.19 0.15 1.20 ... 3HWC J1847-017 0.0 338.1 281.95 -1.75 0.09 0.17 HESS J1848-018 3HWC J1849+001 0.0 427.5 282.35 0.15 0.06 0.16 IGR J18490-0000 3HWC J1852+013† 0.0 126.7 283.05 1.34 0.09 0.06 2HWC J1852+013* 3HWC J1857+027 0.0 763.5 284.33 2.80 0.06 0.14 HESS J1857+026 3HWC J1857+051† 0.0 35.3 284.33 5.19 0.11 1.23 ... Table 3.1: List of sources (part 1) found in the 3HWC catalog search and the search in which they were found (point source or source with 0.5◦ , 1.0◦ , or 2.0◦ extensions). Also shown is the nearest TeVCat source if it is within 1◦ of the HAWC source and the distance between them. TeVCat sources are bolded if it is within 0.5◦ of the HAWC source. Secondary sources are indicated by a dagger (†). Table from [2]. 50 Nearest TeVCat Source Name Radius TS R.A. Decl. 1σ Stat. Unc. Dist. Name (deg) (deg) (deg) (deg) (deg) 3HWC J1907+085 0.0 75.5 286.79 8.57 0.09 0.07 2HWC J1907+084* 3HWC J1908+063 0.0 1320.9 287.05 6.39 0.06 0.14 MGRO J1908+06 3HWC J1912+103 0.0 198.2 288.06 10.35 0.09 0.24 HESS J1912+101 3HWC J1913+048† 0.0 44.7 288.32 4.86 0.13 0.11 SS 433 e1 3HWC J1914+118 0.0 102.9 288.68 11.87 0.09 0.15 2HWC J1914+117* 3HWC J1915+164 0.0 27.5 288.76 16.41 0.10 2.92 ... 3HWC J1918+159 0.0 31.6 289.69 15.91 0.11 1.99 ... 3HWC J1920+147† 0.0 55.4 290.17 14.79 0.10 0.80 W 51 3HWC J1922+140 0.0 176.6 290.70 14.09 0.06 0.10 W 51 3HWC J1923+169† 0.0 46.1 290.79 16.96 0.10 1.54 ... 3HWC J1928+178 0.0 216.7 292.10 17.82 0.08 0.06 2HWC J1928+177 3HWC J1930+188 0.0 115.6 292.54 18.84 0.09 0.09 SNR G054.1+00.3 3HWC J1935+213† 0.0 26.2 293.95 21.38 0.12 1.89 ... 3HWC J1936+223 0.0 28.2 294.08 22.31 0.11 1.62 ... 3HWC J1937+193 0.0 25.2 294.39 19.31 0.13 1.72 ... 3HWC J1940+237 0.0 27.2 295.05 23.77 0.36 0.28 2HWC J1938+238 3HWC J1950+242 0.0 25.1 297.69 24.26 0.11 0.32 2HWC J1949+244 3HWC J1951+266† 0.5 35.6 297.90 26.61 0.61 2.14 ... 3HWC J1951+293 0.0 68.7 297.99 29.40 0.14 0.25 2HWC J1953+294 3HWC J1954+286† 0.0 48.3 298.70 28.63 0.14 0.12 2HWC J1955+285 3HWC J1957+291 0.0 41.0 299.36 29.18 0.10 0.75 2HWC J1955+285 3HWC J2005+311 0.0 33.1 301.46 31.17 0.12 3.01 ... 3HWC J2006+340 0.0 67.4 301.73 34.00 0.14 0.23 2HWC J2006+341 3HWC J2010+345† 0.0 27.6 302.69 34.55 0.17 1.01 ... 3HWC J2019+367 0.0 1227.5 304.94 36.80 0.06 0.07 VER J2019+368 3HWC J2020+403 0.0 93.9 305.16 40.37 0.09 0.40 VER J2019+407 3HWC J2022+431 0.0 29.0 305.52 43.16 0.10 1.80 ... 3HWC J2023+324 1.0 30.7 305.81 32.44 0.73 3.96 ... 3HWC J2031+415 0.0 556.9 307.93 41.51 0.06 0.07 TeV J2032+4130 3HWC J2043+443† 0.5 28.6 310.89 44.30 0.24 2.88 ... 3HWC J2227+610 0.0 52.5 336.96 61.05 0.19 0.16 Boomerang Table 3.2: List of sources (part 2) found in the 3HWC catalog search and the search in which they were found (point source or source with 0.5◦ , 1.0◦ , or 2.0◦ extensions). Also shown is the nearest TeVCat source if it is within 1◦ of the HAWC source and the distance between them. TeVCat sources are bolded if it is within 0.5◦ of the HAWC source. Secondary sources are indicated by a dagger (†). Table from [2]. the full 1523 days of data. In Chapter 4, it is used to estimate the counts from a bursting source over a time scale of seconds. In each case, the estimated counts are used to maximize the likelihood; here the fitted spectral index is obtained. The results of the spectral fitting can be seen in Tables 3.3 and 3.4. The spectral index and flux normalization at 7 TeV are displayed along with the search in which the source was found (point source or source with 0.5◦ , 1.0◦ , or 2.0◦ extension). Also shown is the energy range expected to yield 75% of each source’s significance. 51 Name Radius Index Flux Energy Range (deg) (10−15 TeV−1 cm−2 s−1 ) (TeV) 3HWC J0534+220 0.0 −2.579 ±0.005 +0.067 −0.018 234.2 ±1.4 +55.2 −34.0 1.6 – 37.4 3HWC J0540+228 0.0 −2.84 ±0.14 +0.14 −0.03 4.8 +0.7 +1.1 −0.8 −0.8 1.1 – 28.5 3HWC J0543+231 0.0 −2.13 +0.15 +0.16 −0.16 −0.01 4.2 ±0.9 +1.9 −0.6 10.5 – 205.9 3HWC J0617+224 0.0 −3.05 ±0.11 +0.06 −0.02 4.5 ±0.8 +1.2 −0.8 0.5 – 12.0 3HWC J0621+382 0.5 −2.41 +0.12 +0.08 −0.13 −0.01 8.9 +1.4 +2.6 −1.5 −1.2 5.7 – 138.3 3HWC J0630+186 0.0 −2.21 +0.13 +0.14 −0.14 −0.02 5.1 +0.8 +2.1 −0.9 −0.8 8.4 – 183.8 +0.17 +0.10 +0.9 +1.4 3HWC J0631+107 0.0 −2.23 −0.19 −0.01 4.0 −1.0 −0.6 8.6 – 186.9 3HWC J0631+169 0.0 −2.51 +0.13 +0.13 −0.14 −0.03 5.9 ±0.7 +2.0 −0.9 3.3 – 95.0 +0.14 +0.10 3HWC J0633+191 0.0 −2.64 −0.15 −0.02 4.8 ±0.7 +1.3 −0.7 2.1 – 61.9 3HWC J0634+067 0.5 −2.27 ±0.10 +0.09 −0.01 9.0 +1.3 +2.7 −1.4 −1.2 7.6 – 171.1 3HWC J0634+165 0.0 −2.52 +0.19 +0.13 −0.22 −0.03 +0.7 +1.6 5.0 −0.8 −0.8 3.2 – 93.0 3HWC J0634+180 0.0 −2.47 +0.10 +0.11 −0.11 −0.02 7.3 ±0.7 +2.2 −1.1 3.7 – 102.0 3HWC J0659+147 0.0 ... ... ... 3HWC J0702+147 0.0 −1.99 +0.19 +0.18 −0.24 −0.02 +1.1 +1.9 3.6 −1.3 −0.5 14.1 – 261.8 3HWC J1104+381 0.0 −3.04 ±0.01 +0.06 −0.02 69.4 ±1.3 +16.6 −11.4 0.7 – 15.0 3HWC J1654+397 0.0 −2.91 ±0.04 +0.06 −0.02 20.0 ±1.0 +5.3 −3.2 1.3 – 29.4 3HWC J1739+099 0.0 −1.98 +0.13 +0.12 −0.14 −0.02 3.3 ±0.8 +1.4 −0.6 14.8 – 274.0 3HWC J1743+149 0.0 −2.37 +0.15 +0.08 −0.16 −0.02 4.0 +0.7 +1.3 −0.8 −0.6 5.6 – 139.3 3HWC J1757-240 1.0 −2.80 ±0.11 +0.07 −0.03 +10.9 +17.6 72.0 −11.3 −9.2 5.6 – 158.4 3HWC J1803-211 0.0 −2.59 +0.13 +0.09 −0.14 −0.04 27.9 +5.4 +10.3 −5.7 −4.4 9.7 – 206.8 3HWC J1809-190 0.0 −2.59 ±0.05 +0.09 −0.03 68.0 +4.4 +24.5 −4.5 −9.7 7.7 – 177.3 3HWC J1813-125 0.0 −2.81 ±0.10 +0.09 −0.03 19.9 ±2.0 +5.9 −3.1 2.6 – 69.2 3HWC J1813-174 0.0 −2.54 ±0.04 +0.09 −0.02 74.0 ±4.0 +26.7 −10.7 7.7 – 174.8 3HWC J1819-150 0.0 −2.90 +0.08 +0.10 −0.07 −0.03 33.9 ±2.5 +10.3 −5.4 2.2 – 62.5 3HWC J1825-134 0.0 −2.35 ±0.02 +0.11 −0.02 125.4 ±3.4 +53.1 −18.3 9.2 – 183.4 3HWC J1831-095 0.0 −2.61 ±0.06 +0.12 −0.03 33.8 +1.8 +12.2 −1.9 −5.3 4.2 – 106.7 3HWC J1837-066 0.0 −2.73 ±0.02 +0.11 −0.03 81.9 ±1.7 +25.2 −13.2 2.2 – 57.3 3HWC J1843-034 0.0 −2.36 ±0.03 +0.14 −0.02 47.2 ±1.6 +19.6 −7.2 6.2 – 142.6 3HWC J1844-001 0.0 −2.76 ±0.12 +0.12 −0.03 7.4 ±1.0 +2.4 −1.3 1.9 – 51.2 3HWC J1847-017 0.0 −2.87 ±0.04 +0.10 −0.04 27.7 ±1.3 +8.4 −5.1 1.4 – 33.1 3HWC J1849+001 0.0 −2.17 ±0.04 +0.16 −0.01 23.5 ±1.4 +11.6 −3.9 9.9 – 195.3 3HWC J1852+013 0.0 −2.79 ±0.08 +0.13 −0.04 14.5 +1.0 +4.4 −1.1 −2.6 1.7 – 40.2 3HWC J1857+027 0.0 −2.83 ±0.03 +0.10 −0.03 +10.2 35.9 ±1.2 −6.5 1.3 – 31.8 3HWC J1857+051 0.0 −3.03 ±0.12 +0.09 −0.04 5.9 ±1.0 +1.8 −1.2 0.7 – 15.5 Table 3.3: List of sources (part 1) found in the 3HWC catalog search, the search in which they were found, and their fitted spectra. Each source’s best-fit spectral index and differen- tial flux at 7 TeV is reported. Also shown is the energy range from which 75% of the source’s significance is obtained. The uncertainties displayed are the statistical and systematic un- certainties, respectively. The fit for the source 3HWC J0659+147 did not converge. Table from [2]. 52 Name Radius Index Flux Energy Range (deg) (10−15 TeV−1 cm−2 s−1 ) (TeV) 3HWC J1907+085 0.0 −2.95 ±0.09 +0.09 −0.04 8.4 +0.9 +2.5 −1.0 −1.6 0.8 – 19.7 3HWC J1908+063 0.0 −2.12 ±0.03 +0.18 −0.02 44.7 ±1.4 +20.7 −7.1 8.9 – 182.7 3HWC J1912+103 0.0 −2.85 ±0.05 +0.09 −0.03 14.5 ±0.9 +4.4 −2.6 1.1 – 27.1 3HWC J1913+048 0.0 −2.37 ±0.13 +0.14 −0.03 +0.9 +2.8 6.8 −1.0 −1.1 5.9 – 144.1 3HWC J1914+118 0.0 −2.9 ±0.1 ±0.1 9.5 +0.9 +2.9 −1.0 −1.9 1.0 – 23.5 3HWC J1915+164 0.0 −2.60 ±0.13 +0.07 −0.01 4.6 ±0.7 +1.4 −0.7 2.5 – 69.5 3HWC J1918+159 0.0 −2.49 +0.15 +0.13 −0.17 −0.02 5.1 ±0.7 +1.8 −0.8 3.6 – 101.9 3HWC J1920+147 0.0 −2.98 ±0.13 +0.11 −0.05 6.0 +0.9 +2.0 −1.0 −1.3 0.7 – 16.8 3HWC J1922+140 0.0 −2.62 ±0.06 +0.09 −0.02 13.0 ±0.8 +3.7 −2.0 2.1 – 60.0 3HWC J1923+169 0.0 −3.07 ±0.10 +0.07 −0.03 5.2 ±0.8 +1.4 −0.9 0.5 – 11.1 3HWC J1928+178 0.0 −2.30 ±0.07 +0.17 −0.02 13.6 ±0.9 +5.3 −2.1 5.9 – 140.5 3HWC J1930+188 0.0 −2.76 ±0.07 +0.09 −0.03 10.2 ±0.8 +2.9 −1.7 1.4 – 36.6 3HWC J1935+213 0.0 −2.73 ±0.16 +0.08 −0.03 4.4 ±0.7 +1.3 −0.7 1.6 – 43.0 3HWC J1936+223 0.0 −2.9 ±0.3 +0.2 −0.1 4.1 +0.9 +1.5 −1.2 −1.2 0.9 – 20.9 3HWC J1937+193 0.0 −2.90 ±0.16 +0.09 −0.04 4.2 +0.7 +1.2 −0.8 −0.8 0.9 – 22.2 3HWC J1940+237 0.0 −3.14 +0.12 +0.07 −0.11 −0.03 4.0 ±0.8 +1.1 −0.8 0.4 – 7.9 3HWC J1950+242 0.0 −2.49 +0.18 +0.15 −0.19 −0.02 4.3 ±0.7 +1.5 −0.7 3.7 – 102.2 3HWC J1951+266 0.5 −2.36 +0.12 +0.11 −0.13 −0.02 +1.2 +2.6 8.5 −1.3 −1.3 5.5 – 133.7 +0.09 +0.10 3HWC J1951+293 0.0 −2.47 −0.10 −0.02 7.9 ±0.8 +2.7 −1.2 4.0 – 108.6 3HWC J1954+286 0.0 −2.42 ±0.12 +0.15 −0.02 6.4 ±0.8 +2.4 −1.0 4.8 – 119.7 3HWC J1957+291 0.0 −2.54 +0.13 +0.11 −0.14 −0.02 6.2 ±0.7 +2.1 −1.0 3.3 – 93.9 3HWC J2005+311 0.0 −2.58 ±0.13 +0.11 −0.02 5.6 ±0.7 +1.8 −0.9 3.0 – 83.6 3HWC J2006+340 0.0 −2.56 +0.12 +0.15 −0.13 −0.03 8.5 +0.8 +3.1 −0.9 −1.5 3.3 – 83.5 3HWC J2010+345 0.0 −2.91 ±0.13 +0.10 −0.03 5.4 ±0.9 +1.5 −0.9 1.1 – 24.3 3HWC J2019+367 0.0 −2.04 ±0.02 +0.16 −0.02 34.7 ±1.3 +17.8 −5.7 11.7 – 211.7 3HWC J2020+403 0.0 −3.11 ±0.07 +0.08 −0.04 11.4 ±1.2 +3.4 −2.2 0.7 – 14.8 3HWC J2022+431 0.0 −2.34 ±0.12 +0.10 −0.02 6.0 ±1.1 +2.2 −0.9 8.3 – 181.9 3HWC J2023+324 1.0 −2.70 +0.13 +0.06 −0.12 −0.03 13.8 +1.9 +3.2 −2.0 −2.1 2.0 – 52.9 3HWC J2031+415 0.0 −2.36 ±0.04 +0.14 −0.02 30.7 +1.3 +13.1 −1.4 −4.9 6.3 – 147.8 3HWC J2043+443 0.5 −2.33 +0.14 +0.11 −0.15 −0.02 9.7 +2.1 +3.9 −2.2 −1.5 8.5 – 185.7 3HWC J2227+610 0.0 −2.43 ±0.10 +0.10 −0.02 +5.9 +13.1 30.8 −5.8 −4.1 14.3 – 292.7 Table 3.4: List of sources (part 2) found in the 3HWC catalog search, the search in which they were found, and their fitted spectra. Each source’s best-fit spectral index and differen- tial flux at 7 TeV is reported. Also shown is the energy range from which 75% of the source’s significance is obtained. The uncertainties displayed are the statistical and systematic un- certainties, respectively. The fit for the source 3HWC J0659+147 did not converge. Table from [2]. 53 Galactic plane I; 0.0 ; 1523 days 3 1 b[ ] 1 3 84 82 80 78 76 74 72 70 68 66 64 l[ ] 3 0 3 6 9 12 15 18 21 24 TS Figure 3.2: Significance map of part of the Galactic plane √assuming a point source hypothesis. The √ green lines show significance contours starting at T S = 26 and increasing in steps of ∆ T S = 2. The top labels indicate the positions of sources from TeVCat and the bottom labels indicate the locations of 3HWC sources. Figure from [2]. 3.3 Results 3.3.1 3HWC Sources The catalog search found a total of 65 very-high-energy gamma-ray sources. Seventeen of these sources are considered to be secondary sources since they do not have large separation from their neighbors. The majority of these sources fall within the Galactic plane and are supernova remnants, pulsar wind nebulae, or pulsar halos. Markarians 421 and 501 are the only extragalactic sources found in this search and are active galactic nuclei. Figures 3.2, 3.3, 3.4, and 3.5 show closer views of the significance maps made assuming a point source along the Galactic plane. The locations of the 3HWC sources and TeVCat sources are labeled in the maps. Twenty-eight of the sources found in the 3HWC catalog were not present in the previous 54 Galactic plane II; 0.0 ; 1523 days 3 1 b[ ] 1 3 64 62 60 58 56 54 52 50 48 46 44 l[ ] 3 0 3 6 9 12 15 18 21 24 TS Figure 3.3: Significance map of part of the Galactic plane assuming a point source hypothesis. The top labels indicate the positions of sources from TeVCat and the bottom labels indicate the locations of 3HWC sources. Figure from [2]. 2HWC catalog. These new sources can mostly be explained by the increased live time over which the search was performed, which caused the significance of the sources to pass the detection threshold. In Section 3.3.2, comparisons with other catalogs of known gamma-ray sources are performed to determine what these new sources may be. 3.3.2 Potential Associations with Known Sources To learn more about HAWC’s 28 newly detected sources, several catalogs of known gamma- ray sources were scanned. Sources from these catalogs are considered to be potential associa- tions if they are within 1◦ of the nearest 3HWC source. Eight of these sources are potentially associated with TeV sources from the TeVCat [75]. These sources are 3HWC J0540+228, 3HWC J0543+231, 3HWC J0617+224, 3HWC J0634+067, 3HWC J2227+610, 3HWC J1913+048, 3HWC J1757-240, and 3HWC J1803-211, and their closest TeVCat coun- 55 Galactic plane III; 0.0 ; 1523 days 3 1 b[ ] 1 3 44 42 40 38 36 34 32 30 28 26 24 l[ ] 3 0 3 6 9 12 15 18 21 24 TS Figure 3.4: Significance map of part of the Galactic plane √assuming a point source hypothesis. The √ green lines show significance contours starting at T S = 26 and increasing in steps of ∆ T S = 2. The top labels indicate the positions of sources from TeVCat and the bottom labels indicate the locations of 3HWC sources. Figure from [2]. terparts can be seen in Tables 3.1 and 3.2. Twenty of the new sources have no potential association with previously known TeV gamma-ray sources. Three other source catalogs were searched: the fourth Fermi -LAT source catalog (4FGL) [3], the Australia Telescope National Facility (ATNF) pulsar catalog [4, 76], and the Galactic supernova remnant catalog (SNRCat) [5, 77]. Fourteen of HAWC’s new sources are potentially associated with a previously observed GeV source in the 4FGL catalog. Four of these sources also have potential associations with a supernova remnant from the SNRCat. Two of HAWC’s new sources have no GeV counterpart in 4FGL but are within 1◦ of pulsars from the ATNF catalog. Four of HAWC’s new sources have no counterpart in any of the catalogs that were scanned. Table 3.5 lists the 20 3HWC sources with no nearby TeVCat source, their nearest counterpart in the 4FGL, ATNF, and SNRCat catalogs, and the distance to their nearest counterparts. 56 Galactic plane IV; 0.0 ; 1523 days 3 1 b[ ] 1 3 24 22 20 18 16 14 12 10 8 6 4 l[ ] 3 0 3 6 9 12 15 18 21 24 TS Figure 3.5: Significance map of part of the Galactic plane √assuming a point source hypothesis. The √ green lines show significance contours starting at T S = 26 and increasing in steps of ∆ T S = 2. The top labels indicate the positions of sources from TeVCat and the bottom labels indicate the locations of 3HWC sources. Figure from [2]. HAWC l(◦ ) b(◦ ) 4FGL (◦ ) ATNF (◦ ) SNRCat (◦ ) 3HWC J0621+382 175.44 10.97 4FGL J0620.3+3804 (0.22) J0622+3749 (0.42) ... 3HWC J0630+186 193.98 4.02 ... J0630+19 (0.94) ... 3HWC J0631+107 201.08 0.43 4FGL J0631.5+1036 (0.15) J0631+1036 (0.14) ... 3HWC J0633+191 193.92 4.85 ... ... ... 3HWC J1739+099 33.89 20.34 4FGL J1740.5+1005 (0.22) J1740+1000 (0.13) G034.0+20.3 (0.13) 3HWC J1743+149 39.13 21.68 ... ... ... 3HWC J1844−001 31.95 1.50 4FGL J1848.2−0016 (0.99) J1843−0000 (0.27) ... 3HWC J1857+051 38.22 1.06 4FGL J1855.2+0456 (0.56) J1857+0526 (0.24) ... 3HWC J1915+164 50.19 2.35 4FGL J1912.0+1612 (0.74) B1913+16 (0.32) ... 3HWC J1918+159 50.16 1.33 ... J1918+1541 (0.26) ... 3HWC J1923+169 51.58 0.89 4FGL J1925.1+1707 (0.50) B1921+17 (0.14) ... 3HWC J1935+213 56.90 0.39 4FGL J1935.2+2029 (0.89) J1936+21 (0.24) G057.2+00.8 (0.59) 3HWC J1936+223 57.76 0.73 4FGL J1932.2+2221 (0.94) J1938+2213 (0.44) G057.2+00.8 (0.47) 3HWC J1937+193 55.29 −0.98 4FGL J1936.6+1921 (0.21) J1936+20 (0.77) ... 3HWC J1951+266 63.23 −0.13 4FGL J1951.6+2621 (0.25) J1952+2630 (0.24) ... 3HWC J2005+311 68.74 −0.40 4FGL J2006.2+3102 (0.15) J2006+3102 (0.15) G068.6−01.2 (0.81) 3HWC J2010+345 72.14 0.56 ... ... ... 3HWC J2022+431 80.52 3.54 ... ... ... 3HWC J2023+324 71.85 −2.77 4FGL J2024.0+3202 (0.43) ... ... 3HWC J2043+443 83.74 1.10 4FGL J2047.5+4356 (0.79) ... ... Table 3.5: List of sources found in the 3HWC catalog search that have no known TeV counterpart. For each source, the galactic latitude and longitude are listed along with the sources from the 4FGL [3], ATNF [4], and SNRCat [5] catalogs within 1◦ . The distance from the 3HWC source to the sources from the other catalogs are shown in parentheses. Table from [2]. 57 Four of the new 3HWC sources have no counterpart in any of the previously mentioned catalogs that were searched. Three of these sources (3HWC J0633+191, 3HWC J2010+345, and 3HWC J2022+431) are not well-isolated from other extended sources in the 3HWC catalog. More detailed studies of each of these regions are required to determine whether or not these sources are part of existing sources. The fourth source, 3HWC J1743+149, is located off the Galactic plane and is just over the detection threshold. More data is needed to determine whether this source is a new gamma-ray source or an upwards fluctuation of the background. 3.4 Limitations and Systematic Uncertainties There are several caveats with the catalog search methods and other sources of uncer- tainty that must be considered when interpreting the results. The uncertainties mainly stem from statistical fluctuations and the modeling of the detector. There are other limitations due to the fitting of the sources in the catalog, which may be possible to resolve with studies of each source that go more in depth than can be performed here. 3.4.1 Background Fluctuations It is possible for the background to fluctuate upwards enough to pass the detection thresh- old and imitate a source. To estimate the frequency of such false positive sources, twenty simulated significance maps were created using the background counts from the original source search. The simulated number of signal counts in each pixel is obtained by Poisson- fluctuating the number of background counts in the corresponding pixel. An example of one of the randomized background maps is shown in Figure 3.6. Each of these maps are then run through the full source search as described in Section 3.2.2, including point and extended searches. In the twenty maps, a total of fifteen local maxima are found with a T S > 25. Therefore, the estimated number of false positive sources in the search is 15/20 = 0.75. These fluctuations are distributed evenly across the 58 Figure 3.6: All-sky significance map with randomized background. This map was created by Poisson-fluctuating the number of background events in each pixel. The full catalog search was performed on twenty such maps in order to estimate the amount of spurious sources. sky and typically occur just above the threshold value of T S = 25. 3.4.2 Search Methods Limitations As described in Section 3.2.2, the 3HWC catalog is compiled by first listing all point sources and then adding the sources with larger extensions. The point source list is added first because the extended source searches may merge nearby sources with lower extension to- gether. On the other hand, bright extended sources may be found as multiple nearby point sources. This effect is particularly prominent in the region around the Geminga pulsar, which is shown in Figure 3.7. Consequently, the search in which the source first appears cannot be interpreted as a measurement of its extension. More detailed morphology studies are required in order to determine the extension of each source, but are beyond the scope of this analysis. 59 Galactic plane V; 0.0 ; 1523 days Galactic plane V; 1.0 ; 1523 days 15 15 9 9 b[ ] 3 b[ ] 3 3 3 9 9 -143 -149 -155 -161 -167 -173 -143 -149 -155 -161 -167 -173 l[ ] l[ ] 3 0 3 6 9 12 15 18 21 24 3 0 3 6 9 12 15 18 21 24 TS TS Figure 3.7: Significance maps of the region containing the Crab Nebula (bottom right) and the Geminga pulsar (center). On the left is the map produced with a point source hypothesis and on the right is the map produced with a 1◦ extended source hypothesis. Several point sources were found near the Geminga pulsar in the 3HWC search while the map on the right shows a single extended source in the region, illustrating the limitations to the search method. Figure from [2]. 3.4.3 Systematic Uncertainties The systematic uncertainty on HAWC’s pointing increases as a function of zenith angle. HAWC’s absolute pointing uncertainty is 0.1◦ for sources that transit within 30◦ of directly overhead, but it increases to 0.15◦ for sources at declinations of −10◦ or +50◦ and may be up to 0.3◦ at declinations of −20◦ or +60◦ . Figure 3.8 shows the differences between the declinations of sources measured by HAWC compared to the declinations of the same sources measured by other experiments. Imaging air Cherenkov telescopes like the ones in Figure 3.8 typically have better angular resolution than extensive air shower arrays, as discussed in Chapter 1. New algorithms in HAWC are currently being produced in order to reduce the effects that cause this systematic pointing uncertainty on sources near the edges of the field of view. There are a few other systematic uncertainties related to the modeling of the HAWC 60 VERITAS 0.5 MAGIC H.E.S.S. 0.4 0.3 δHAWC − δIACT [◦ ] 0.2 0.1 0.0 −0.1 −0.2 −20 −10 0 10 20 30 40 50 60 HAWC declination [◦ ] Figure 3.8: Comparisons between the declinations of HAWC’s sources and their counterparts measured by IACT experiments H.E.S.S., MAGIC, and VERITAS. The difference between HAWC’s measured declination and that from the IACTs’ for each source are plotted on the y- axis as a function of HAWC’s declination on the x-axis. HAWC’s measurements of the source locations agree within uncertainties for most of its declination range except for sources whose zenith angles are larger than 30◦ . Error bars are HAWC’s statistical uncertainty. Figure from [2]. 61 detector. These effects include PMT efficiency, late light simulation, PMT thresholds, and charge uncertainty. The PMT efficiency cannot be determined using the calibration system, so an event selection based on charge and timing cuts is used instead to identify incident vertical muons. Since vertical muons provide a monoenergetic source of light, they can be used to measure the relative efficiency of each PMT by matching the muon peak position to the position expected from simulations. The late light uncertainty is due to a mismodeling of late light in the air shower and stems from a discrepancy between the time width of the laser pulse used for calibration and the time structure of the air shower. The PMT threshold uncertainty is based on observations of the cosmic-ray rate. The charge uncertainty considers how much a PMT measurement varies for a fixed amount of light and the relative differences in photon detection efficiency from PMT to PMT [67]. To compute the uncertainties, the spectral fits were repeated with different detector models. As there were no correlations found between these effects, they have all been added in quadrature. An extra uncertainty of 10% is assigned to the flux normalizations to account for additional effects that have not been considered. These systematic uncertainties on the spectral indices and flux normalizations are reported along with the statistical uncertainties in Tables 3.3 and 3.4. Finally, there are also systematic uncertainties that affect the spectral fits. The spectral fits assume a simple power-law spectrum and do not consider any other spectral models. For the two extragalactic sources found in the catalog search, absorption due to the extra- galactic background light has not been taken into account. Galactic diffuse emission is also not accounted for here, and may affect the flux normalizations of Galactic sources. These effects must be characterized in more detail and each source also needs to be studied independently in order to provide a more accurate depiction. Several such studies have been completed [78–82] or are ongoing but are beyond the scope of this analysis. 62 3.5 Conclusions The 3HWC catalog search yielded 65 sources emitting gamma rays in the TeV energy range. This search used 1523 days of data and provided an unbiased view of the Northern sky. HAWC is an ideal instrument to conduct such a survey due to its large field of view and high uptime. Improvements on the analysis presented here are already in the works in order to more ac- curately characterize the catalog of sources. These include improving the systematic pointing uncertainty, incorporating effects from the Galactic diffuse emission, performing multi-source fits, and studying the morphology of each source. Using the energy estimators mentioned in Section 2.4.2 can provide more accurate spectral fits and using the outriggers mentioned in Section 2.1 can help extend this analysis to even higher energies. This catalog provides many new targets for other observatories to follow up. These may include confirmation and more detailed studies by IACTs, extending the energy spectrum to other wavelengths, and multimessenger searches. Future gamma-ray observatories that are more sensitive such as the Southern Wide-field Gamma-ray Observatory [83] and the Cherenkov Telescope Array [49] may be able to add more sources to this catalog. Improve- ments on HAWC’s analysis combined with information from other observatories will help to solve several of the open problems in high-energy astrophysics. 63 CHAPTER 4 SEARCH FOR PRIMORDIAL BLACK HOLES In addition to searching its visible sky for steady gamma-ray sources, HAWC can also be used to search for transient sources such as Primordial Black Holes (PBHs). Its large field of view and continuous operations allow it to detect bursts of gamma rays with short emission durations. This analysis uses data from a transient burst search in HAWC to set upper limits on the local PBH burst rate density. The results from this study were published in [21]. 4.1 Introduction Primordial Black Holes were created by density fluctuations in the early Universe. They may have been formed at times from the Planck time (5.391 × 10−44 s) to one or more seconds after the Big Bang. Depending on when they were formed, their initial masses could be as small as the Planck mass (2.176 × 10−8 kg) or as large as supermassive black holes (105 M ) [84]. Hawking predicted that a black hole will thermally emit fundamental particles with a temperature inversely proportional to its mass in a calculation that convolved quantum field theory, thermodynamics, and general relativity [45]. This process is known as Hawking Radiation and results in a final burst of particles at the end of the PBH’s lifetime. PBHs formed with initial masses of ∼ 5 × 1014 g should be expiring today with a burst of gamma rays in the GeV - TeV energy range. This allows HAWC to observe their evaporation as it is sensitive to energies between a few hundred GeV to over 100 TeV. Confirmed detection of a PBH burst would prove their existence and the theory of Hawk- ing radiation, which has yet to be observed. It would also allow for a measurement of their relic density and rate-density of evaporation. Additionally, PBHs in certain mass ranges may make up a non-negligible fraction of dark matter. Current limits constrain the frac- tion of dark matter that is made up of PBHs to . 10% [85]. Even non-detection of PBH 64 evaporations would provide important insight into the early Universe and the physics of inflation. Searching for PBHs can constrain the cosmological density fluctuation spectrum in the early Universe on scales smaller than those constrained by the cosmic microwave background [20]. Detection or upper limits on the number density of PBHs can also inform inflationary models [84]. Various observatories have searched for PBH bursts using both direct and indirect meth- ods, from cosmological scales down to parsec distance scales. At cosmological scales, the PBH density can be probed using the 100 MeV extragalactic gamma-ray background [84, 86, 87]. On the galactic scale, the local PBH density and anisotropy can be detected using 100 MeV gamma-ray measurements [88]. On the kiloparsec scale, the PBH burst limit can be set using the Galactic antiproton background [89], and on the parsec scale, the PBH burst limit is set directly by searching for individual evaporating PBHs [25–28]. In this work, the search for PBH bursts is performed on the parsec scale. 4.2 Theory The properties of the PBH’s final burst depend on the particle physics at the GeV - TeV energy scale. The predicted energy spectrum is determined by the high-energy particle physics model used. Here the Standard Evaporation Model (SEM) [90] is assumed as the emission model. According to this model, the PBH will emit fundamental Standard Model particles whose wavelengths are comparable to the size of the black hole. As the temperature of the black hole approaches the Quantum Chromodynamics confinement scale (∼ 200 − 300 MeV), quarks and gluons will be radiated which then decay into photons, neutrinos, electrons, positrons, protons, and anti-protons [90]. The emission rate of fundamental particles from a black hole as predicted by Hawking [45, 91] can be written as d2 N 27 = n ψs (x), (4.1) dEdt 2π~(8π)2 dof where ~ is the reduced Planck constant, ndof is the number of degrees of freedom of the 65 particle species, and ψs (x) is a dimensionless emission rate function defined by the quantity x and the particle spin s. The quantity x is defined as E x= , (4.2) kTBH where E is the energy of the radiated particle, k is Boltzmann’s constant, and TBH is the temperature of the black hole. The distribution ψs (x) peaks at different values of x depending on the spin of the particle [20]. The temperature of a black hole is inversely proportional to its mass following the relation,  10  10 kg kTBH = 1.058 GeV, (4.3) MBH where MBH is the mass of the black hole. The black hole loses mass at the following rate, dMBH α(MBH ) =− 2 , (4.4) dt MBH where the function α(MBH ) incorporates all directly emitted particle species and their de- grees of freedom. For black holes whose temperatures are > 100 GeV (MBH < 108 kg), the parameter α = 8.0 × 1017 kg3 s−1 [20]. The black hole mass can be written as a function of the remaining evaporation lifetime (τ ) of the black hole as  τ 1/3 MBH = 1.3 × 106 kg, (4.5) 1s where units are used such that the Boltzmann constant k = 1. Combining Equations 4.3 and 4.5, the black hole temperature (TBH ) can be expressed in terms of its remaining lifetime as  τ −1/3 TBH = 7.8 × 103 GeV, (4.6) 1s for TBH < Tp , where Tp is the Planck temperature (1.22 × 1019 GeV) [20]. The emission rate increases as the black hole loses mass. To obtain the energy spectrum of the black hole, the emission rate (Equation 4.1) is integrated over the remaining lifetime. dN The total time-integrated photon energy spectrum, dEγ , can be parameterized as γ    3/2   3/2  1 GeV 1 GeV GeV−1  dNγ  TBH Eγ for Eγ < TBH ≈ 9 × 10 35 , (4.7) dEγ  1 GeV  3 GeV−1    E for Eγ ≥ TBH γ 66 Time-integrated photon energy spectrum 1031 0.2s 1029 1s 10s 1027 dN/dE [GeV−1 ] 1025 1023 1021 1019 100 101 102 103 104 105 106 E [GeV] Figure 4.1: Time-integrated energy spectrum of photons emitted during black hole evapora- tion over remaining lifetime intervals τ = 0.2, 1, and 10 seconds. Adapted from [20]. for gamma-ray energies above 1 GeV [20]. This energy spectrum includes contributions from photons that were directly emitted from the black hole as well as gamma rays that decayed from other particles species that were directly emitted. Figure 4.1 shows the photon energy spectrum for the PBH remaining lifetimes that will be studied in this work. The light curve of a bursting black hole can be found by integrating the differential emission rate in Equation 4.1 over photon energy, where the energy range is determined by the energy range of the detector being used. This yields the following relation for the time dependence of the black hole burst [20], dNγ ∝ τ −1/2 . (4.8) dt Figure 4.2 shows the gamma-ray time profile of the emission from the black hole burst in the 50 GeV to 100 TeV energy range. 67 Figure 4.2: Emission time profile of a bursting PBH integrated over energy in the range 50 GeV - 100 TeV. This profile can be described by a power law with an index of -0.5. Figure from [20]. 4.3 Data Selection The dataset used for this search consists of 959 days of HAWC data collected between March 2015 and May 2018. The data was collected from an all-sky self-triggered transient search originally designed to detect gamma-ray bursts (GRBs), which are short-lived bursts of intense gamma-ray light. Since PBHs also emit bright gamma-ray light on a similar timescale to short GRBs, this data can be easily optimized to search for PBHs instead. The data were reconstructed according to the procedure described in Section 2.4. How- ever, unlike in the 3HWC search described in Section 3.2.1, this analysis is performed using only a single energy bin. Events were required to have hit at least 70 PMTs, and a loose gamma-hadron separation cut was applied. The remaining background coming from the cosmic rays that survived the cut is estimated using a 1.75 hour direct integration. Gamma 68 rays from steady sources may also contribute to the background, but for the timescales rel- evant to this search these gamma rays have a negligible contribution to the overall rate of events [92]. The self-triggered transient search continuously looks for transient events occurring through- out HAWC’s field of view at energies above a few hundred GeV. It uses sliding time windows and spatial bins to search the entire dataset for clusters of events above the cosmic-ray rate and compares each cluster to a background-only hypothesis. The pre-trials probability for each of the candidate events is stored if it passes a reporting threshold of p ≤ 2.87 × 10−7 , corresponding to a significance 5σ [92]. The search reports p-values as opposed to number of counts due to the zenith angle dependence of the background. A higher background is expected near zenith since the particles have less atmosphere to travel through than at the edges of the field of view, so more signal counts would be required to produce a significant event. Thus the p-value is used to compare candidate burst events across the sky. The spatial search is performed in a grid of right ascension and declination within a zenith angle of 50◦ from the detector using square bins of size 2.1◦ × 2.1◦ . The spatial bin is then advanced by an amount smaller than the bin size in order to allow for desired overlap [92]. The size of the overlap is different for each search duration, and is optimized to balance fine tuning on the spatial position of the burst and avoiding strong correlations between pixels. Points outside of the 50◦ zenith angle are excluded since most photons at the energies expected from a transient burst signal do not have sufficient energy to reach HAWC due to the larger amount of atmosphere they would have to traverse. After the spatial search is performed, the time window is advanced and the spatial search is repeated. Three different time window lengths are used in this analysis – 0.2, 1, and 10 seconds. The fixed-width time windows were used instead of the more sensitive procedure of fitting a light curve profile in order to maintain the computational efficiency necessary to search the entire HAWC dataset [92]. After both the spatial and temporal searches are completed, the p-values of each event are corrected for trials and added to a histogram. 69 Burst Duration 0.2 s Burst Duration 1 s 104 104 HPBH HPBH frequency Hmodel frequency Hmodel 3 3 10 10 Hbkg Hbkg 102 Hdata 102 Hdata 10 10 Excluded Excluded 1 from limit 1 from limit calculation calculation 10−1 10−1 10−2 −17 −16 −15 −14 −13 −12 −11 −10 10−2 −17 −16 −15 −14 −13 −12 −11 −10 10 10 10 10 10 10 10 10 10−9 10−8 10−7 10−6 10 10 10 10 10 10 10 10 10−9 10−8 10−7 10−6 p-value p-value (a) τ = 0.2 s (b) τ = 1 s Burst Duration 10 s 104 HPBH frequency 103 Hmodel Hbkg 102 Hdata 10 Excluded 1 from limit calculation 10−1 10−2 −17 −16 −15 −14 −13 −12 −11 −10 10 10 10 10 10 10 10 10 10−9 10−8 10−7 10−6 p-value (c) τ = 10 s Figure 4.3: Distributions of the p-values of the data, background, PBH, and model his- tograms for the three remaining lifetimes used in this search. See Section 4.4.1 for a de- scription of HP BH and Hmodel . The vertical dashed black line indicates pthr , which is the passing threshold for data to be included in this analysis. Note that Hbkg and HP BH were drawn from a much larger sample than Hdata and were therefore scaled in order to match its duration. Figure from [21]. Figure 4.3 shows the distribution of probability values for both the data (Hdata ) and the background (Hbkg ) for each time duration searched. The background distribution was drawn from a large sample in order to fill out the histogram to low p-values and was then scaled to match the duration of the data collection. The sharp cutoff on the right side of the data distribution is due to the reporting threshold that each candidate event must pass. 70 4.4 Analysis Methods 4.4.1 Creating the Model In order to interpret the data from the transient search to look for PBHs, it is first necessary to create a model of PBH bursts. This is accomplished using a Monte Carlo simulation to generate a random set of PBH burst events in HAWC’s field of view. Burst source points are randomly generated assuming a burst rate density of 104 pc−3 yr−1 . This burst rate density was chosen because it is close to HAWC’s predicted sensitivity to PBH bursts [93], but it could have been set to any positive non-zero value. Events are excluded if their zenith angle is larger than 50◦ or if they are farther than 0.5 pc away. Beyond 0.5 pc, even PBH bursts simulated at HAWC’s zenith for 10 s did not produce an observable signal within HAWC. This procedure produces a cone centered around HAWC’s field of view, as can be seen in Figure 4.4. The simulation was run for thousands of years in order to produce enough statistics but then scaled down to match the duration of the time during which the data was collected. To determine the expected number of photons, µ, observable at HAWC for each of the simulated PBH bursts as a function of remaining lifetime τ , distance to the PBH r, and zenith angle θ, the following expression is used: 1 − f E2 dN (τ ) Z µ(r, θ, τ ) = A(E, θ)dE, (4.9) 4πr2 E1 dE where dN/dE is the PBH gamma-ray spectrum from Equation 4.7 integrated from the re- maining lifetime τ to 0. The energies E1 and E2 correspond to the energy range of the detector, A(E, θ) is the effective area of the detector as a function of energy and zenith angle, and f is the dead-time fraction of the detector. A forward-folding approach is used to perform the calculation of Equation 4.9 with HAWC’s software ZEBRA [74], which was first introduced in Section 3.2.3. In this case the emission is assumed to all come from a single zenith angle since the timescale of the emission is short. ZEBRA was used to calculate the expected number of photons from each 71 Figure 4.4: Locations of simulated PBH burst events. The HAWC detector is located at the center, and axes are in units of parsecs. PBHs were simulated out to 0.5 pc within 50◦ of HAWC’s zenith, resulting in a cone centered around HAWC’s field of view. simulated PBH burst for the three burst durations of 0.2, 1, and 10 s being studied in this analysis. Figure 4.5 shows the number of estimated counts from a simulated PBH at a distance of 0.05 pc with a remaining lifetime of 0.2 s. Now that the expected number of counts from each simulation has been determined, the last step is to calculate the p-value of each simulated burst event in order to be able to compare the simulation with the data. The Poisson probability of obtaining N or more counts given the background B for each burst event is ∞ X B i e−B γ(N, B) p = P(n ≥ N ) = = , (4.10) i! Γ(N ) i=N where N = µ+B and γ is the lower incomplete gamma function. The resulting p-values were 72 6 5 4 µ 3 2 1 0 10 20 30 40 50 ◦ θ[] Figure 4.5: Example signal expected at HAWC (µ) from a simulated PBH burst as a function of zenith angle (θ) using HAWC’s code ZEBRA. This particular simulated PBH is assumed to have a remaining lifetime of 0.2 s and located at a distance of 0.05 pc from HAWC. The signal plateaus above 45◦ because of HAWC’s low sensitivity at those zenith angles. Figure from [21]. compiled into a histogram for each of the three remaining lifetimes studied. This histogram is called HP BH and can be seen in Figure 4.3. To produce a complete model of the data, the histogram of p-values from simulation was added to the background distribution. This new histogram, Hmodel = Hbkg + HP BH , can also be seen in Figure 4.3. 4.4.2 Calculating an Upper Limit Since no statistically significant signal from PBHs was found in the data, upper limits are computed on the local burst rate density of PBHs. The upper limits are computed separately for each remaining lifetime studied. 73 First an analysis threshold, pthr , is defined to be the p-value which corresponds to 10 events in the background histogram, Hbkg . This threshold is chosen in order to ensure that the limit placed is conservative and not sensitive to fluctuations in the data. Above pthr , a fiducial region exists where it can be verified that the background and data distributions are statistically equivalent in order to ensure that the overall rate normalization seen in data matches the background estimate from simulation. Since it was confirmed that the data are well-described by a background-only distribution above the threshold, the extrapolation to the low-statistics region of the background distribution is satisfactory. This threshold is indicated by a vertical dashed black line in Figure 4.3, and all p-values > pthr are excluded from the limit calculation. To compute an upper limit on the local burst rate density of PBHs at the 99% confidence level, a test statistic (TS) is defined as follows, T S = [2 ln L1 − ln L0 ], (4.11) where L1 and L0 are the model and background Poisson likelihoods, respectively. For binned data, the log-likelihood of a Poissonian variable x and mean λ for a sample that is indepen- dently and identically distributed with size n, is generally given by, n n ! ! X Y ln[L(λtot |x1 , x2 , ...xn )] = xi ln λi − nλtot − ln xi ! . (4.12) i=1 i=1 Then the background log-likelihood can be written as X ln L0 = [Hdata (p) ln(Hbkg (p)) − Hbkg (p)], (4.13) p where Hdata is the histogram of p-values produced from HAWC’s transient burst search and Hbkg is the histogram of p-values generated by Monte Carlo using the background distributions from the transient burst search. Similarly, the model log-likelihood can be written as X ln L1 = [Hdata (p) ln(Hmodel (p)) − Hmodel (p)], (4.14) p 74 where Hmodel (p) = Hbkg (p) + HP BH (p) is the histogram of p-values created from the model of PBH bursts in HAWC. Here the sum is over the bins p instead of the sample size. Both of these expressions neglect the factorial term in the likelihood from Equation 4.12 since it will cancel when evaluating Equation 4.11. Now that the likelihoods and test statistic are defined, the next step is to scan over the PBH burst rate density, ρ̇, until the value is found that maximizes the TS in Equation 4.11, T Smax . The 99% upper limit is computed by again scanning over ρ̇ until T Smax is decreased by ∆T S. To find the amount ∆T S by which to decrease T Smax , the equation below is used, Z √∆T S/σ 1 2 2 √ e−x /2σ dx = CL, (4.15) −∞ 2πσ 2 where CL is 0.99, the desired confidence level. Wilks’ theorem states that the TS follows a χ2 distribution with one degree of freedom in the limit of a large data sample [73], which allows the use of this equation. Solving for ∆T S, the value 5.41 is obtained, meaning the value of ρ̇ that corresponds to a T S99 = T Smax − 5.41 gives the upper limit at a 99% confidence level. The resulting upper limits for each remaining lifetime searched are shown in Section 4.5. 4.5 Results The 99% upper limits on the PBH burst rate density for each of the three remaining lifetimes searched are listed in Table 4.1. In order to be conservative, the limit corresponding to HAWC’s strongest predicted sensitivity, which occurred for the remaining lifetime of 10 s, is quoted as the official result of this analysis. This corresponds to an upper limit of 3400+400 −3 −1 −200 pc yr , which is the strongest limit placed on the local PBH burst rate density at the time of writing. The long tails in the model distribution (Hmodel , see Figure 4.3) provided the statistical power that allowed for the setting of strong constraints. This upper limit is compared to the upper limits set on PBH burst rate density by other gamma-ray experiments in Table 4.2. 75 Burst duration Burst Rate Upper Limit 0.2 s 3300 +300 −3 −1 −100 pc yr 1s 3500 +400 −3 −1 −200 pc yr 10 s 3400 +400 −3 −1 −100 pc yr Table 4.1: The 99% upper limits on the PBH burst rate density for the three remaining lifetimes searched. The uncertainties are systematic only, and are described in Section 4.5. Experiment Burst Rate Upper Limit Search Duration Reference Milagro 36000 pc−3 yr−1 1s [28] VERITAS 22200 pc−3 yr−1 30 s [26] H.E.S.S. 14000 pc−3 yr−1 30 s [25] Fermi-LAT 7200 pc−3 yr−1 1.26 × 108 s [27] HAWC 3 yr. 3400 pc−3 yr−1 10 s This Work Table 4.2: The strongest limit on the burst rate density of PBHs set by each of the five detectors most sensitive to PBH burst studies. Figure 4.6 shows the resulting upper limits compared to other gamma-ray experiments for all of the burst durations studied by the different instruments. Note that the burst duration is a search parameter and not a physical parameter. This means that any differences between the limits placed for different durations are only due to variations in the signal-to-background ratio of the detector. HAWC’s expected limits and the 68% and 95% containment bands for the null hypothesis are also shown in Figure 4.6. The expected limits and containment bands were created using 1000 simulations containing no PBHs. The containment bands are purely statistical and indicate that the upper limits are compatible with the expected sensitivity of HAWC assuming only background events. The systematic uncertainties mainly stem from discrepancies in the modeling of the detector. The dominant sources of these discrepancies are PMT efficiency and late light simulation. The systematic uncertainties that contributed less are PMT thresholds, angular resolution discrepancy, and charge uncertainty. These effects are the same as in Section 3.4.3 and are described in more detail in [67]. To estimate the effect of these uncertainties in this analysis, a new model histogram 76 PBH Burst Rate Density Upper Limits Burst Rate Density [pc−3 yr−1 ] Whipple (2006) Milagro (2014) 108 CYGNUS (1993) HAWC 3 Year Limit (this work) Tibet Air Shower Array (1995) HAWC Expected Limit HESS (2013) 68% containment VERITAS (2017) 95% containment 7 10 Fermi LAT (2018) 106 105 104 103 10−4 10−3 10−2 10−1 100 101 102 108 109 Burst Duration [s] Figure 4.6: Comparison of HAWC’s upper limits at the 99% confidence level at burst du- rations 0.2, 1, and 10 s with the upper limits from Whipple [22], CYGNUS [23], the Tibet Air Shower Array [24], H.E.S.S. [25], VERITAS [26], Fermi -LAT [27], and Milagro [28]. The expected limits and the 68% and 95% containment bands are also shown. Systematic un- certainty is not shown since it was significantly smaller than the containment bands. Figure from [21]. Hmodel was created for each source of uncertainty based on different detector models and then the analysis was repeated. Since these effects are independent, they were all added in quadrature. This resulted in an overall systematic uncertainty of +10.6%, -4.3%. The contributions from each source of uncertainty as well as the total systematic uncertainty can be seen in Figure 4.7. 4.6 Conclusions A search for PBH bursts in three years of HAWC’s data has been conducted. No such burst events were found, so upper limits were placed on the local PBH burst rate density at the 99% confidence level. The limit of 3400 pc−3 yr−1 is currently the strongest constraint 77 5 × 103 Nominal PMT efficiency Late light Burst Rate Density [pc−3 yr−1] Charge uncertainty PMT threshold Total systematic uncertainty 4 × 103 3 × 103 10−1 100 101 Burst Duration [s] Figure 4.7: The difference in upper limits obtained by using different detector models corre- sponding to different sources of systematic uncertainty. The black stars represent the final result for each burst duration searched and the red band represents the total systematic uncertainty obtained by adding each effect in quadrature. Figure from [21]. on this quantity placed for nearby PBH bursts. To improve this analysis in the future, a dedicated PBH search optimized on HAWC’s raw data rather than adapting results from a previous transient search can be performed. Such a search would incorporate more days of data, data from the outrigger array, and improvements to reconstruction algorithms in HAWC that are currently in the works. It could also include more burst durations to search and information about the light curve shown in Figure 4.2 that was not used in this analysis. In order to improve the analysis statistically, a future analysis would use full likelihood profiles in both time and space and include likelihood stacking. Future gamma-ray instruments may also be more sensitive to bursting PBHs. For example, preliminary studies demonstrate that the Southern Wide field of view Gamma-ray Observaotry (SWGO) will be able to place limits that are more constraining than the ones presented here by an order of magnitude [94]. 78 CHAPTER 5 THE ICECUBE OBSERVATORY The analyses described in this work so far used a single messenger, gamma rays, to study astrophysical objects. In order to perform multi-messenger astronomy, another messenger particle is needed. Here neutrinos are chosen to execute such a study, and this chapter describes the IceCube Neutrino Observatory, which is used to detect the neutrinos. Similar to HAWC, IceCube also has a high up-time and wide field of view of the sky, making it an ideal instrument to study transient events. 5.1 The IceCube Array The IceCube Neutrino Observatory is a cubic kilometer of instrumented ice located at the South Pole. IceCube detects neutrinos from 10 GeV up to PeV energies. Its field of view covers the entire sky and it achieves an uptime of 99%. IceCube consists of several components: an in-ice array of Digital Optical Modules (DOMs), a more densely instru- mented array known as DeepCore, and an array at the surface called IceTop. The IceCube Laboratory (ICL) hosts the electronics for the experiment [30]. A schematic of IceCube is shown in Figure 5.1. IceCube’s in-ice array is made up of 5160 DOMs on 86 vertical strings located between depths of 1450 m and 2450 m below the surface of the ice. Each string contains 60 DOMs along with a cable that connects in the ICL. The vertical separation of DOMs on 78 of the strings is 17 m and the horizontal spacing between them is 125 m. The strings are arranged in a hexagonal pattern, optimized to detect astrophysical neutrinos from TeV to PeV energies. The main science goals of IceCube are to identify and characterize astrophysical neutrino sources, which includes participating in multi-messenger collaborations [30]. Each DOM contains a 10-inch downward-facing PMT along with readout electronics and an LED flasher board. These components are all contained within a glass sphere designed 79 Figure 5.1: A schematic of the IceCube detector, showing the in-ice array, DeepCore, IceTop, and the IceCube Laboratory [29]. The Eiffel Tower is shown for scale. to withstand the high pressure in the deep ice. Neutrinos interact with the nuclei in the ice and produce secondary particles which are detected by the PMTs inside the DOMs via the Cherenkov radiation the particles produce by traveling faster than the speed of light in the ice. The LED flasher board is included for calibration purposes. It helps to verify the timing response of the DOMs in the analysis software, measure the position of the DOMs in the ice, measure the optical properties of the ice, and verify the performance of shower reconstruction algorithms [30]. A schematic of a DOM can be seen in Figure 5.2. The DeepCore array is a subset of DOMs with denser spacing located in the center of the main in-ice array. It is composed of 8 special strings and 7 nearby standard IceCube strings. The special strings have a smaller vertical spacing of 7 m and horizontal spacing of 72 m between DOMs. They are located at depths between 2100 m and 2450 m below the surface. 80 Figure 5.2: A schematic showing the components of a Digital Optical Module (DOM). Each DOM consists of a PMT, an LED flasher board, and other electronics to readout the events, all contained within a glass sphere. Figure from [30]. Six of the DeepCore strings use PMTs that have 35% higher quantum efficiency than the standard IceCube DOMs. DeepCore lowers the energy threshold to about 10 GeV, allow- ing it to perform searches for neutrino oscillations, dark matter annihilation, and Galactic supernovae [30]. The IceTop array is located on the surface of the ice and is designed to detect cosmic rays in the PeV to EeV energy range. It consists of 162 tanks arranged in 81 stations at locations that follow approximately the same grid as the in-ice array. Each tank is filled with ice up to a height of 0.9 m and is instrumented with two DOMs. IceTop has a denser infill array made up of eight stations in the center separated by 10 m that corresponds to the spacing of DeepCore, which lowers the detection threshold to 100 TeV. It is designed to study PeV gamma rays, transient events, and acts as a partial veto for detection of down-going neutrinos [30]. The final component of the IceCube detector is the ICL, located on the surface of the 81 Figure 5.3: A picture of the IceCube Laboratory (ICL) that contains the data acquisition system and online computers for the experiment. Figure from [31]. ice at the center of the array. It is the central operations building for the experiment. The cables from the array are routed up two towers on either side of the lab and into a server room. The server room houses the data acquisition and online filtering computers [30]. A picture of the ICL can be seen in Figure 5.3. The IceCube detector has a couple of planned upgrades to improve its sensitivity both at very high energies and at lower energies in order to better detect astrophysical neutrinos as well as study neutrino oscillations. The IceCube Upgrade will add seven new densely spaced strings containing a total of 700 DOMs near the bottom center of the in-ice array [95]. IceCube-Gen2 will add 120 new strings each containing 80 DOMs with vertical spacing of 16 m and average horizontal spacing of 240 m, which increases the overall instrumented volume to 7.9 km3 . It will also add a radio array consisting of 200 antenna stations covering an area of 500 km2 [32]. Figure 5.4 shows the design and scale of the IceCube Upgrade and both the optical and radio components of IceCube-Gen2. These upgrades will improve almost all neutrino analyses with IceCube, but are not included in this analysis as they have 82 Figure 5.4: Layout of the IceCube Upgrade and the optical and radio components of IceCube- Gen2 compared to the current IceCube detector. Figure from [32]. not been completed at the time of writing. 5.2 Data Acquisition The data flow starts from the individual PMTs. The signal from a PMT is compared to a discriminator threshold of 0.25 photoelectrons. When this threshold is crossed, the process of capturing and digitizing the waveform is started. The amount of information about the signal that is sent to the surface computers via the cables depends on whether any neighboring DOMs also detected a signal within a time period of 1 µs. If only an isolated signal occurred, then a timestamp and brief charge summary are included. This is known as a Soft Local Coincidence (SLC). If neighboring DOMs also detected a signal, then the full waveform is compressed and included along with the timestamp and charge summary. This is called a Hard Local Coincidence (HLC) [30]. A diagram of the flow of data from a PMT to the surface is shown in Figure 5.5. Each DOM has its own free-running clock, so timestamps must be translated to the syncrhonized clocks in the ICL, which are then translated to UTC. This is accomplished using a calibration procedure that continuously runs by sending pulses to and from each DOM. The waveforms are calibrated by subtracting the common baseline from each sample and multiplying by the gain. The baseline is the mean of the pedestals of each sample in a waveform, which consists of 128 effectively independent digitizers, and the gain is the input 83 PMT Hit records sent to surface DAQ Trigger Flag to Neighbor DOMs SPE Trigger Clock Capture Time Stamp, Discriminator Charge Summary Waveform Data if Coincidence Summary SPE Trigger Flags "Launch" Coincidence? bytes from Neighbor DOMs 40 Msps 256 Samples ADC Digital 300 Msps Low Gain > Threshold ? Waveform Analog Compression Capture & Med. Gain > Threshold ? (lossless) Delay Digitizers High Gain 128 Samples Preamplifiers (each gain) Figure 5.5: Diagram of data flow starting from a PMT and ending with information about the signal being sent to the surface. If a PMT signal crosses the discriminator threshold and so does the signal from neighboring DOMs, the waveform is digitized and sent to the surface. Figure from [30]. voltage change that is needed to increase the readout value by one count. This is calculated by the DAQ software at the beginning of each run [30]. Most DOM hits originate from dark noise, which is PMT signals resulting from effects not due to an external photon source. Therefore several different trigger algorithms and filters are required to help determine which hits are due to neutrino events. The DAQ is responsible for detecting patterns among the DOMs that are likely to be caused by particle interactions and recording those hits as an event. The Simple Multiplicity Trigger (SMT) requires a certain number of HLC hits to occur within a time window of a few microseconds, and it continues recording hits until there is a period of time of a few microseconds without any HLC hits. The multiplicity requirement is a different value for the IceCube in-ice array, DeepCore, and IceTop. There are other triggers designed for certain physics analyses or for 84 calibration [30]. The online Processing and Filtering system (PnF) is responsible for reducing the data volume from all triggered events that were collected by the DAQ. Each filter is designed to select events that are useful for a particular physics analysis. The main filters are the muon track filter, the shower event filter, and the high-charge filter. The muon track filter selects high-quality track events, which is useful for point source and transient neutrino searches. The shower event filter selects events that deposit large amounts of energy in the instru- mented volume, which is used to search for high-energy shower events from atmospheric and astrophysical neutrinos. The high-charge filter selects events that deposit at least 1000 photoelectrons in the instrumented volume, which is used to search for high-energy astro- physical neutrinos. There are also other filters that are used to select cosmic rays in IceTop, lower-energy neutrino events contained within DeepCore, and others for more specialized searches [30]. In order to control and monitor IceCube’s operations, a system called LiveControl is used. This system is responsible for controlling the state of the DAQ and the PnF, as well as starting and stopping data-taking runs and recording the parameters of the runs. It also maintains a database of monitoring quantities and sends alerts if any value is outside a given range or hasn’t been reported in a given amount of time. The monitoring data is displayed on the IceCube Live website, which has one copy on the IceCube network at the South Pole and one copy in the Northern Hemisphere. These pages display information about active alerts, plots of event rates, and the start and stop time of each run, among other quantities [30]. 5.3 Event Reconstruction In order to make use of the data that IceCube collects, information about the original neutrino that interacted within the IceCube detector volume needs to be reconstructed from the PMT measurements. IceCube’s reconstruction algorithms use the timing of hits and the amount of charge deposited to reconstruct the direction and energy of the original neutrino, 85 to determine the flavor of neutrino, and to determine whether the event was a background or signal event. This information can then be used to perform a wide variety of science analyses, from searching for astrophysical sources of neutrinos to studying neutrino oscillations and searching for dark matter. 5.3.1 Event Signatures The two main ways that neutrinos interact with nucleons in the ice are via charged current interactions or neutral current interactions. In charged current interactions, the neutrino exchanges a W boson with a quark that makes up the nucleon to produce a lepton that is the same flavor as the neutrino. In neutral current interactions, the neutrino exchanges a Z boson with a quark and does not produce a lepton [33]. The Feynman diagrams for both charged and neutral interactions are shown in Figure 5.6. Different flavors of incoming neutrinos as well as different interaction types result in different event signatures within the detector [96]. A charged current interaction initiated by a muon neutrino results in an outgoing muon that travels in a straight line through the detector. The muon loses energy as it travels through the ice via bremsstrahlung and pair production and is detected via the Cherenkov light these secondary particles produce. This creates a track event topology within IceCube, shown in Figure 5.7 [96]. A charged current interaction initiated by an electron neutrino produces a free electron and creates an electromagnetic shower. The electron loses energy through bremsstrahlung and pair production and is detected via Cherenkov radiation. Although the energy loss mechanisms are the same for electrons as they are for muons, electrons lose their energy faster because of their lower mass. This results in a cascade event topology within IceCube, as can be seen in Figure 5.8 [96]. A charged current interaction initiated by a tau neutrino produces a tau lepton. If the tau has an energy above about a PeV, it will travel a small distance before it decays, producing 86 Figure 5.6: Feynman diagrams for charged current interactions (a/b) and neutral current interactions (c/d) between neutrinos and quarks that make up nucleons in the ice. A W boson is exchanged in charged current interactions and produces an outgoing lepton corresponding to the neutrino flavor along with hadronic cascades. A Z boson is exchanged in neutral current interactions and only hadronic cascades are produced. Figure from [33]. a double cascade event topology, shown in Figure 5.9. If it has a lower energy, it will produce a cascade event that is indistinguishable from cascades produced by charged current electron neutrino interactions [96]. Neutral current interactions initiated by any flavor of neutrino produce cascade event topologies in the IceCube detector as in Figure 5.8. The outgoing neutrino in this type of interaction is not detected, and the outgoing quark or antiquark creates a cascade of hadronic particles [96]. Track events from muon neutrinos result in better angular resolution that cascade events due to their longer lever arm for the directional reconstruction. This makes them the ideal event type with which to perform neutrino astronomy. As the subject of the analysis in 87 Figure 5.7: A track event in IceCube. A muon neutrino initiates a charged current interaction and produces a muon that creates the track signature in the detector. The colors correspond to the timing of the PMT hits; red is early and green is late. The size of the PMT hits correspond to the amount of charge deposited; more charge is represented by larger spheres. Figure from [33]. Figure 5.8: A cascade event in IceCube. Cascade events are produced by either neutral cur- rent interactions with any neutrino flavor, charged current interactions involving an electron neutrino, or charged current interactions involving a tau neutrino below about 1 PeV. Figure from [33]. 88 Figure 5.9: A double cascade event in IceCube. A tau neutrino above 1 PeV initiates a charged current interaction and produces a hadronic shower and a tau lepton, which then decays and creates a second shower. This produces a double cascade signature in the detector. Figure from [34]. Chapter 6 is multi-messenger astronomy, only track events will be used. Therefore the following sections only describe the angular and energy reconstruction of track events and not cascades. 5.3.2 Angular Reconstruction IceCube’s reconstruction process starts with simpler algorithms, which are then used to seed more sophisticated and accurate algorithms. The first and simplest of these algorithms is known as LineFit, which reconstructs the position and direction of the muon track. This algorithm assumes that the muon will propagate through the detector at a constant velocity and performs a least-square fit to the locations and times of each hit DOM. To partially improve the robustness of LineFit to outlier hits, a Huber penalty function is included in the fit [97]. This results in a median angular resolution of a few degrees, and is used to seed the next reconstruction algorithm [98]. 89 The next algorithms are likelihood-based algorithms for angular reconstruction known as SPE (single photoelectron) and MPE (multi photoelectron). SPE is the standard unbinned likelihood defined as NDOM Nhit [pj (ti )]qi , Y Y L(Θ) = (5.1) j=1 i=1 where NDOM and Nhit are the total number of DOMs and hits respectively, p(t) is the photon arrival probability distribution function (PDF), and qi is the observed charge in each DOM. A numerical minimizer such as Minuit is used to minimize − log L(Θ) [99]. However, since the above likelihood only uses the arrival time of single photons in each DOM, the likelihood can be improved by adding information from other hits [98]. The MPE likelihood is a modified version of the SPE likelihood that accounts for the total charge observed on each DOM by using data that is weighted with the observed charge per hit and the total observed charge per DOM. This likelihood is defined as NDOM Q −q [pj (t1 )]q1 · (1 − Pj (t1 )) j 1 , Y L(Θ) = (5.2) j=1 P where qi is the observed charge per hit and Q = i qi is the total observed charge per DOM, and t1 and q1 are the time and charge corresponding to the first DOM hit. This likelihood uses the photon arrival PDF of the first photon that hits the DOM. Since the first photon that arrives to a DOM is less likely to have undergone significant scattering compared to the average single photon, this likelihood results in a more accurate reconstruction of the six track positional and orientation parameters, Θ = (x, y, z, t, θ, φ) [98]. A spline parameterization is used to determine the PDF in the MPE likelihood in Equa- tion 5.2. This reconstruction is known as SplineMPE, and it includes a more accurate representation of the ice properties. This is accomplished by first inferring the optical prop- erties of the ice from fitting an ice model to flasher data from the LEDs in the DOMs. Then muon tracks are simulated with a given ice model and the resulting photons are recorded 90 in histograms. These histograms are fit with a multi-dimensional spline surface in order to represent the photon arrival PDF. This reconstruction algorithm improves the median angular resolution to a few tenths of a degree [98]. 5.3.3 Energy Reconstruction After the direction of a muon track is reconstructed, IceCube’s energy reconstruction algo- rithms use a Poisson likelihood model to estimate the energy that the muon deposits inside the detector along the track. The observed number of photoelectrons k in a PMT are com- pared to the expected light output of a simulated event Λ to estimate the deposited energy E as follows, ln L = k ln EΛ + ρ − (EΛ + ρ) − ln k!, (5.3) where ρ is the number of expected noise photons. The expected light output takes into account the optical properties of the ice as well as the expected detector response. Contri- butions from all DOMs are included, and the likelihood is then maximized with respect to the one free parameter E [96]. This energy reconstruction technique is called the MuEX algorithm. The Millipede reconstruction algorithm performs an unfolding in order to determine the energy E of the muon. It uses an adapted form of the Poisson likelihood in Equation 5.3. This algorithm takes into account that the muon undergoes stochastic losses due to Brehsstrahlung and pair production. It models the muon more accurately by treating it as a linear combination of energy losses in the form of electromagnetic cascades. The unfolding is performed by rewriting the likelihood in Equation 5.3 in terms of vector operations as follows, ~k − ρ~ = Λ · E, ~ (5.4) 91 where ~k are the observed photons on all DOMs, ρ~ are the noise photons, E ~ describes the energy losses, and the matrix Λ contains the predicted light yield at every point in the detector from every source position. Millipede unfolds the energy loss distribution for a given track position and direction [96]. Millipede can also fit for the position and direction of the muon by repeating the energy unfolding for different track hypotheses. This fits the 5 + N parameters (x, y, z, θ, φ, E),~ where N is the length of E,~ or the number of energy losses that occur along the muon track. This fit is performed on an angular grid using the HEALPix package to apply Millipede at each pixel [71]. This results in a resolution on total deposited energy along the muon track of ∼ 10 − 15% [96]. 5.3.4 Backgrounds to Astrophysical Neutrinos The main background to detecting astrophysical neutrinos comes from atmospheric neutri- nos. Atmospheric neutrinos are created by cosmic ray air showers. As described in Section 2.1, cosmic rays that interact with Earth’s atmosphere create a shower of secondary particles that include muons and neutrinos. These secondary muons and neutrinos can create signa- tures in the IceCube detector that are similar to those created by astrophysical neutrinos. As will be discussed in the following section, the event selections have more stringent cuts on tracks that travel downwards through the detector in order to filter out the events caused by cosmic rays. The tracks that travel upwards through the detector comprise a purer sample of astrophysical neutrinos because the Earth acts as a shield of cosmic-ray muons. Figure 5.10 shows an illustration of upgoing and downgoing cosmic rays and neutrinos in IceCube. 5.4 Realtime Alert System IceCube has developed a realtime alert system to enable multi-messenger astronomy. The main goal of the system is to alert the community when it detects high-energy neutrinos that are likely to be of astrophysical origin, to allow other observatories to conduct follow- 92 Figure 5.10: This illustration shows neutrinos originating from different angles that interact in the IceCube detector. Down-going cosmic rays come from the Southern hemisphere and create atmospheric muons and neutrinos that are the main background to astrophysical neutrinos. Up-going cosmic rays come from the Northern hemisphere and the muons they produce are filtered out by the Earth. The signal being searched for in this case are the cosmic neutrinos that originate from astrophysical sources. Figure from [33]. 93 up observations. This can help with understanding the origin of cosmic rays, understanding how astrophysical sources accelerate their particles, and identifying sources of neutrinos. The realtime alert system started operating in April 2016 and introduced new neutrino candidate selections in June 2019 [100]. As described in Section 5.2, the IceCube DAQ records events when a trigger condition has been satisfied. For neutrino alerts, eight DOMs must receive photon signals within a 5 µs time window. The DAQ records all DOM hits inside a window of +6/-4 µs around the trigger time into a single event. Then the event is processed through the online processing and filtering system, which includes reconstruction algorithms and event selections [101]. The reconstruction algorithms use the patterns of light deposited in the detector during each event to determine the direction, position, and energy of each event. From these algorithms, about 1% of the events are selected to be potentially of astrophysical neutrino origin. The events are chosen if they are well-constructed track-like events that pass a certain charge threshold. This threshold is more stringent for events coming from the Southern sky because of the background from atmospheric muons and neutrinos in that region. The online filtering system then selects events from this set to be sent as alerts [101]. There are three event selections responsible for identifying neutrinos that are likely to have originated from astrophysical sources. The Gamma-ray Follow Up (GFU) selection is designed to select high-quality track events that are similar to those used in point-source analyses. For events from the Northern sky, the energy is estimated using the reconstructed muon energy and for events from the Southern sky, the energy is estimated using the to- tal charge of photoelectrons measured by the detector. The High Energy Starting Event (HESE) selection picks out neutrino events whose interaction vertex occurs within the de- tector volume. It contains cuts on measured track length and declination in order to ensure the event is track-like and more likely to be of astrophysical origin. The Extreme High Energy (EHE) selection searches for track-like neutrino events whose energies are between 500 TeV and 10 PeV. It accomplished this by choosing well-constructed track-like events 94 where the number of photoelectrons is greater than 4000 [100]. A cascade event selection has been introduced in 2021 but will not be described here as events from this selection are not included in this analysis. If the neutrino is selected as a candidate, it is sent as an alert through either a ”Gold” or a ”Bronze” stream, depending on its signalness. Signalness is defined as follows, NS (E, δ) Signalness(E, δ) = , (5.5) NS (E, δ) + NB (E, δ) where NS and NB are the number of signal and background events, respectively, above the estimated neutrino energy E at a declination δ. Events are sent out as Gold alerts if the candidate’s signalness is greater than 50% and they are sent out as Bronze alerts if the signalness is between 30% and 50%. This choice results in about 10 Gold alerts and 20 Bronze alerts per year [100]. Before the alerts are sent, it is important to first make sure that the detector is operating in a stable manner. This is done by checking quantities such as the rate of primary DAQ triggers, the rate of well-defined muon track events selected by the online filter system, the rate of events selected by the selection criteria, and the number of active DOMs. These values are recorded in ten minute intervals in order to look for deviations from the long-term average. If there are no such deviations, the alert will be sent [101]. The alerts are distributed as Gamma-ray Coordinates Network (GCN) Notices [102]. These notices contain information such as the date and time the neutrino event occurred, whether it was a Gold or Bronze event, the location in right ascension and declination, the statistical error on position, the estimated energy, the signalness, and the false alarm rate of each event. This information allows other instruments to perform follow-up observations on the location of the neutrino alert. Information about the event, including the light arrival information from the DOMs, are also sent to computers located at the University of Wisconsin Madison in order to run more sophisticated reconstructions. These reconstructions are more computationally intensive and 95 provide a better angular estimation as well as improved distinction between track-like and shower-like events. This is accomplished by evaluating the likelihood of different arrival directions using increasingly finer directional grids on the sky in order to determine the best- fit direction. These algorithms improve the angular resolution by about 50% to 1◦ for tracks and a revision to the original alert is sent out with this updated information, generally a few hours after the initial alert [101]. The alerts sent out by IceCube are followed up by a wide range of instruments including those that observe in the optical, radio, X-ray, high-energy and very-high-energy gamma ray regions as well as other neutrino observatories. Some of the instruments that have conducted follow-up observations of IceCube’s neutrino alerts include Fermi -LAT [3], the Zwicky Transient Facility [103], INTEGRAL [104], HAWC [1], MAGIC [47], VERITAS [48], ANTARES [105], Swift [106], and VLASS [107]. IceCube also conducts follow-up observa- tions in the location of the alert event to search for any other neutrinos that originated from the same source. An example of an alert that resulted in detections from other instruments occurred on September 22, 2017, when MAGIC and Fermi -LAT follow-up observations re- vealed high-energy gamma rays originating from the same location [57]. 96 CHAPTER 6 MULTI-MESSENGER SEARCH Now that two instruments have been introduced, HAWC and IceCube, that detect two different messenger particles, gamma rays and neutrinos respectively, it is possible to perform a multi-messenger search. The fields of view of these two observatories have a large amount of overlap as IceCube’s view of the full sky encompasses HAWC’s entire field of view. Further, IceCube has higher sensitivity in the Northern sky due to the better background reduction for up-going neutrinos as described in Section 5.3.4. Both detectors also have a high up-time, making it more likely that both are taking data when an interesting event occurs. These factors make them both good instruments to study transient events and to compare the data between them. The analysis described in this chapter consists of a search for flares in HAWC’s gamma-ray data that originate from the same location in the sky as IceCube’s neutrino alerts that were described in Section 5.4. 6.1 Introduction This search seeks to provide insight into some of the open questions in particle astro- physics that were described in Chapter 1. First, emission of both gamma rays and neutrinos from the same source indicates a hadronic production mechanism, which would help to ex- plain how the source accelerates its particles. Additionally, IceCube has detected a diffuse flux of astrophysical neutrinos of unknown origin [108]. By spatially correlating gamma rays with the astrophysical neutrinos that IceCube detects, there is a higher potential for the sources of neutrinos to be identified. The most notable example to date of such a corre- lation is the detection of high-energy gamma rays by MAGIC and Fermi -LAT following a neutrino alert from IceCube from the blazar TXS 0506+056 [57]. This motivates a search for high-energy gamma-ray flares originating from the locations of IceCube’s high-energy astrophysical neutrinos that prompt the issuing of an alert. 97 Previous analyses that have searched for coincident emission of neutrinos from IceCube and gamma rays from HAWC consider only simultaneous emission of both messengers. HAWC performs a real-time follow-up search in its nearest transit within the time that an IceCube alert was issued and reports its results in GCN notices. Another analysis searches for simultaneous emission of sub-threshold events that currently runs in real time via the Astrophysical Multimessenger Observatory Network (AMON) [109]. The analysis described in this chapter follows many of the steps from the coincident sub-threshold event analysis described in [110]. The neutrino flare from TXS 0506+056 that prompted the alert from IceCube to be sent occurred on September 22, 2017 [57]. However, this was not the only flare from this blazar in neutrinos. An archival search performed by IceCube for other neutrinos from the direction of TXS 0506+056 found another flare that occurred in 2014 [111]. The analysis presented in this chapter aims to search for a similar occurrence of a flare in gamma rays. Here an archival search for gamma-ray flares that are spatially correlated with IceCube’s neutrino alerts but occurring at any point in time is presented. 6.2 Analysis Methods 6.2.1 Data Selection Data from both IceCube and HAWC are used in this analysis. From IceCube, there are about ten years of neutrino alerts as described in Section 5.4. Even though the realtime alert system did not start running until 2016, the IceCube collaboration applied the realtime alert selection criteria on archival data going back to 2010 to compile a list of neutrino events that would have triggered an alert had the realtime system been in place then. Therefore the set of IceCube alerts being used here begins in October 2010 and goes through October 2020. Information about the location of the neutrino events in right ascension and declination along with the statistical uncertainty region of the location are used in this analysis. Of the 98 Figure 6.1: IceCube realtime neutrino alert locations on top of HAWC’s significance skymap from the 3HWC catalog search [2]. The IceCube alerts are spaced out evenly across HAWC’s field of view. Eighty-five of the 99 total alerts issued by IceCube are shown here as those are the ones that fall within HAWC’s field of view. 99 total alerts issued during the ten-year time frame, 85 of them fall within HAWC’s field of view. Figure 6.1 shows the location of the neutrino alerts displayed on top of HAWC’s significance skymap from the 3HWC catalog search in Section 3.2.2. The data used from HAWC consists of daily flux information from each day HAWC was running between March 2015 and June 2019. This data was reconstructed using the methods described in Section 2.4. The flux from each pixel in HAWC’s field of view is calculated and stored in a HEALPix map as described in Sections 3.2.1 and 3.2.2. The only difference in this case is that there is a separate map for each day of data taken by HAWC as opposed to all the data being combined into one map as is the case for the 3HWC catalog. An example of HAWC’s daily flux from the ∼ 4 years of data in one location on the sky is shown in Figure 6.2. 99 Figure 6.2: The daily flux from data collected by HAWC between March 2015 and June 2019. Any gaps in the data are due to the HAWC detector being down for maintenance. 6.2.2 Bayesian Blocks To detect any flares in HAWC’s data, a method known as Bayesian Blocks is used. This is a dynamic histogramming method which optimizes a fitness function in order to identify an optimal binning scheme for a given set of data, where the bins can have differing widths. The aim is for this algorithm to identify blocks of time where the daily flux during a flare rises above the average flux over the rest of the data. It is able to identify any number of flares that last for any amount of time throughout the four years of flux data. This algorithm is advantageous because it can handle data that is not evenly spaced or contains gaps that may occur due to detector downtime. Other HAWC analyses that use this algorithm include a study of Mrks 421 and 501 [112] and a real-time flare monitor [113]. Other experiments like Fermi -LAT, Swift, and others have used Bayesian Blocks to study the variability of astrophysical sources like blazars, quasars, and magnetars, and perform a search for GRBs [114–116], for example. 100 6.2.2.1 The Algorithm The Bayesian Block algorithm is a method used to identify and characterize local variability that is statistically significant within sequential time series data. It finds an optimal seg- mentation of the full observational interval by maximizing a fitness function to the data in each segment, or block. The goal is to find light curve features that are distinct from global signals that are present for most of the observational interval [117]. In this case the goal is to identify flares, or an elevated flux state that only occurs for a short amount of time compared to the length of the entire data set. There are two modes and three different types of data to which the algorithm can be applied. The two modes are a real-time trigger mode where it triggers on the first statistically significant rate change as data is being collected, and a retrospective mode where it is applied to a set of data after it has been collected. The three types of data are event data, binned event data, and point measurements. Event data, or time-tagged event data, is where light collected by a telescope is recorded as a set of detection times of individual photons. Binned event data is very similar to time-tagged event data except the events are collected into bins. Point measurements is where measurements of an observable are collected at a sequence of points in time. This occurs when photons are not detected individually but instead in flux measurements [117]. In this analysis, the retrospective mode is used as the goal is to perform an archival search, and the data type is point measurements because it is collected in flux per day. In order to find the optimal block configuration, a piecewise-constant model, or step function, is fit to the data, which is found by maximizing a fitness function. All model parameters are constant within each block but undergo discrete jumps at the change points, which mark the edges of the blocks. The model of the observation contains three parameters: cp Ncp , the number of change points, tk , the change-point starting block k, and Xk , the signal amplitude in block k, where there are a total of Nblocks = Ncp +1 blocks. The model for each block has two parameters: the signal amplitude Xk and the length of the block. The signal 101 amplitude is determined after the change points have been located by taking the average of the signal inside each corresponding block [117]. The fitness function provides a measure of how well the data inside a given block can be represented by a constant signal. It depends only on the data inside the block and not on any outside of it. The fitness of a block configuration is the sum of the fitnesses of all the individual blocks that make up the configuration: Nblocks X F (C) = f (Bk ), (6.1) k=1 where F (C) is the total fitness of the configuration C, Nblocks is the total number of blocks, and f (Bk ) is the fitness of block k. The best model is found by maximizing F over all possible configurations C [117]. In the point measurements data type, the time series data in the form of a signal s(t) measured at a sequence of times tn with observational errors zn can be modeled as xn ≡ x(tn ) = s(tn ) + zn n = 1, 2, ..., N. (6.2) The measurements do not have to be evenly spaced, and the measurements along with their errors are assumed to be independent of each other. Then assuming that the errors obey a Gaussian probability distribution with a mean of zero and given variance, and the model signal is the constant s = λ, the likelihood of measurement n is 1 xn − λ 2   ! 1 Ln = √ exp . (6.3) σn 2π 2 σn Then the likelihood of block k is the product of the likelihoods of each measurement within that block. To remove the dependence of the likelihood product on the parameter λ, find the value of λ which maximizes the block likelihood. The quantity that needs to be maximized is xn − λ 2   1X − . (6.4) 2 n σn 102 Then the maximum value of λ, λmax , is X xn X 1 λmax = 2 / 2 . (6.5) n σn 0 σ 0 n n Define the quantities ak and bk as 1X 1 ak = (6.6) 2 n σn2 and X xn bk = − s. (6.7) n σ n Then inserting λmax into the product of the likelihood in Equation 6.3 and taking the logarithm results in the fitness function of b2k log Lkmax = . (6.8) 4ak As can be seen from Equation 6.8, the fitness function only depends on the data xn and the errors σn inside each block. The times of measurements tn do not appear in the fitness function as they are only used to determine which data points belong in a block [117]. The total number of possible block configurations is 2N , where N is the number of data points. The Bayesian Block algorithm uses a dynamic programming approach so that the computations are completed in order N 2 time. This approach is a similar concept to mathematical induction, where it starts with the first data point and adds one more at each step until it has gone through the entire data set [117]. The algorithm begins with the first data point where R = 1, where R can be considered as the step number that the algorithm is on. In this case there is only one possible configuration, so it is already the optimal configuration. Now assume that step R has been completed, where the optimal configuration has been determined. At each step, the value of optimal fitness is stored in the array best, and the location of the last change point of the optimal configuration is stored in the array last. The next step is to obtain the optimal configuration for R + 1 data points [117]. 103 For some r consider the set of all configurations of the first R + 1 data points whose last block starts with data point r and ends at R + 1. The fitness of this last block is F (r). The only member of this set that could possibly be optimal is the optimal configuration of r − 1 followed by this last block. Using the fact that the total fitness of the configuration is the sum of the fitnesses in each block as in Equation 6.1, the fitness of this configuration is the sum of F (r) and the fitness of the configuration of r − 1 that was saved from the previous step,    0  r=1 A(r) = F (r) + (6.9)   best(r − 1)  r = 2, 3, ..., R + 1 Here A(r = 1) is the special case where there is only one block that covers the entire range of data. The value of r that results in the optimal configuration of R + 1 data points is the value that maximizes A(r) [117]. After repeating this computation until R = N , the last remaining step is to find the locations of the change points of the optimal configuration. This information is found in the array last, in reverse order [117]. 6.2.2.2 Calibration In order to determine the total number of blocks that the algorithm will find, it is necessary to define a prior distribution, which will indicate the relative probabilities of obtaining a smaller or larger number of blocks. Since in most cases the number of blocks will be much smaller than the number of measurements in an observation period, the prior assigns a smaller probability to a large number of blocks, 1−γ P (Nblocks ) = N +1 γ Nblocks . (6.10) 1−γ Based on this probability distribution, finding k + 1 blocks is less likely than finding k blocks by the constant factor γ. Therefore, to incorporate the prior for the number of blocks, the constant − log γ is subtracted from the fitness of each block. This prior is denoted as 104 Figure 6.3: Calibration of the ncp prior parameter used in the Bayesian Blocks algorithm. Here a false positive rate of 5% is used, which corresponds to an ncp prior value of 6.75. ncp prior, and it is a function of the number of data points and the false positive rate of detecting a change point [117]. To find this value of ncp prior, a light curve is simulated that only contains statistical fluctuations and has no signal present. The simulated light curve contains 1500 days of data since that is about the amount of data that will be used in this analysis. Then the Bayesian Block algorithm is run on this simulated data for a range of ncp prior values. This is repeated 1000 times and the false positive rate is plotted against the values of ncp prior that were tested. The value of ncp prior which yields the chosen false positive rate is the value that will be used with the algorithm in the analysis. Figure 6.3 illustrates this procedure, which was repeated for a couple of different values of the mean and width of a Gaussian distribution that were chosen to represent the values that will be seen in this analysis. Here a false positive rate of 5% is chosen, which results in an ncp prior value of 6.75. 105 Figure 6.4: The Bayesian Block algorithm applied to HAWC flux data at the location of an example IceCube neutrino alert. In this case only one block is generated because the flux is fairly constant over the duration of the data set. 6.2.2.3 Application to Data Now that the Bayesian Block algorithm is calibrated, it is ready to be applied to HAWC data. However, it is first necessary to obtain the correct data to feed the algorithm. This analysis is performing a search for flares in HAWC’s data at the same locations of IceCube’s neutrino alerts. However, IceCube’s alert locations come with an error region of varying size. To account for this error region, the pixel in the 3HWC map within IceCube’s error region that recorded the highest flux is used. Then HAWC’s flux for each day corresponding to that pixel’s location is used as input to the algorithm, along with the flux errors. Any negative fluxes are zeroed out so that the Bayesian Blocks do not trigger on underfluctuations of the signal with respect to the background. Figure 6.4 shows the Bayesian Block algorithm applied to the daily flux data in the example in Figure 6.2. In this example, the flux is approximately constant over the entire observation period, within errors, so the Bayesian Block algorithm only generated a single block that covers the entire data set. In the next section, fake flares will be simulated and the Bayesian Block algorithm will be run on HAWC flux data including a simulated flare. 106 6.2.3 Fake Flare Study To determine the performance of the Bayesian Block algorithm on HAWC data, fake flares are generated and the algorithm is run on the data with the fake flare included. Fake flares are generated by injecting a source lasting for a short duration into HAWC data. This is accomplished using the ZEBRA software that was first described in Section 3.2.3, which simulates a source on top of an existing map. Here this was chosen to be a daily map in the middle of HAWC’s range of data while the detector was running. A few different values of flare duration are studied here: 1, 3, 5, and 10 days. The energy spectrum of the injected source that produces the fake flare is assumed to follow a simple power law. The spectral index 2.0 was chosen since it is similar to that expected from the types of sources that will be studied in this analysis. Since the Crab Nebula is a source that has been well-studied within HAWC, the flux normalizations of the fake source were chosen to be in multiples of the flux normalization of the Crab. This value, 2.53 × 10−13 TeV−1 cm−2 s−1 , is taken from HAWC’s most recent study of the Crab Nebula [67]. A few different values of flux normalization are studied here: 1, 2, 3, and 5 times that of the Crab Nebula. The location in the sky for this test is chosen to be a location of one of IceCube’s realtime alerts which occurs at a declination that is towards the middle of HAWC’s field of view. The Bronze alert that occurred on May 12, 2020 was chosen, which occurred at a right ascension of 295.18◦ and a declination of 15.79◦ . A significance map of HAWC’s data with a fake flare on top of it is shown in Figure 6.5. In this example, a flare lasting one day with a flux normalization equal to that of the Crab Nebula is shown. The next step is to determine the daily flux from the location that the fake flares were injected for all of the other days besides the fake flare between March 2015 and June 2019. This data is the same as that described in Section 6.2.1. The flux values from the fake flares observed at HAWC are found in the same way, but the maps including the fake flares are used instead of the normal data maps for the duration of the flare. Then the Bayesian Block 107 Figure 6.5: A significance map of an example fake flare that was simulated on top of HAWC’s data. This example flare has a duration of one day and a flux normalization equal to that of the Crab Nebula. This flare achieves a detected significance of around 6σ, which is close to HAWC’s expected daily significance from the Crab Nebula of 5σ. algorithm is run on this flux data for all of the combinations of flare durations and flux normalizations mentioned above. Table 6.1 summarizes the results of this study. It has been determined that a flare that lasts for only one day would need to be at least three times as bright as the Crab Nebula in order for it to be detected by the Bayesian Blocks. A flare that is just as bright as the Crab would need to last for at least five days in order for it to be detected. A few example plots showing the data with the injected flare and the Bayesian Blocks can be seen in Figure 6.6. 108 (a) Fake flare lasting one day that is the same brightness as the Crab Nebula. (b) Fake flare lasting ten days that is three times as bright as the Crab Nebula. Figure 6.6: Two example fake flares injected into HAWC’s data following the procedure described above. The black points are HAWC’s data with the fake flare replacing the data during its length and the orange line shows the Bayesian Blocks. The flare is detected in the bottom case (b) but not in the top (a). 109 Flare Intensity Flare Duration 1× Crab 2× Crab 3× Crab 5× Crab 1 day x x X X 3 days x X X X 5 days X X X X 10 days X X X X Table 6.1: This table summarizes the results from the fake flare study, where flares are simulated and the Bayesian Block algorithm is performed on HAWC’s data with a simulated flare. The results of all combinations of flare intensity in terms of that of the Crab Nebula and flare durations in days are shown here. A check mark (X) represents a detection and an x indicates that the Bayesian Block algorithm did not detect the fake flare. 6.2.4 Simulation 6.2.4.1 Cosmology Since the goal of this analysis is to find gamma-ray flares in HAWC from neutrino sources that generated an alert in IceCube, it is useful to simulate neutrinos sources that also produce gamma rays in order to determine HAWC’s sensitivity to such sources. Here a study of transient neutrino source populations is performed to put the archival results into context. The sources being considered in this search are of extragalactic origin, so it is necessary to take cosmological evolution into account. The simulation software package that will be described in the next section assumes the ΛCDM model of a flat universe [118]. According to this model, the observed flux density at Earth can be described as Lν Sobs = , (6.11) 4πd2L (z) where Lν is the neutrino luminosity of the source, z is its redshift, and dL (z) is its luminosity distance [119]. The luminosity distance is given by dz 0 Z z c dL (z) = (1 + z) p , (6.12) H0 0 ΩM (1 + z)3 + ΩΛ where c is the speed of light, H0 is the Hubble parameter, ΩM is the matter density, and ΩΛ is the energy density [120], whose values have been adopted from [121]. 110 The flux density can also be expressed in terms of the observed flux, Z Emax Sobs = E · Φ(E)dE, (6.13) Emin where Φ(E) is the flux from the source and Emin and Emax are the energy bounds of the detector. Using a simple power law assumption, the differential neutrino flux can be expressed as E −γ   dΦ Φ0 (z) = (1 + z)2−γ , (6.14) dE 4πd2L (z) E0 where γ is the spectral index, Φ0 is the flux normalization, and E0 is the pivot energy. In addition to the flux, the other important quantity to determine is the total number of neutrino sources in the universe. This is given by 4π(c/H0 )d2L (z) Z zmax N (zmax ) = ρ(z) p dz, (6.15) 0 (1 + z)2 ΩM (1 + z)3 + ΩΛ where ρ(z) is the local comoving source number density and zmax is the largest redshift value that contains the neutrino sources [122]. The density is a function of the source evolution, ρ(z) = A × P (z), where P (z) is the source evolution function. For transient sources, the above calculations need to be modified [123]. Assume that all transient sources have the same duration ∆t and the spectral index does not change with time. For a source at redshift z, the duration observed at Earth will be (1 + z)∆t. Then the differential neutrino fluence, or the differential flux integrated over the time of the flare, is given by E −γ   dΨ Φ0 (z) = (1 + z)3−γ ∆t. (6.16) dE 4πd2L (z) E0 6.2.4.2 FIRESONG The IceCube software package FIRESONG (FIRst Extragalactic Simulation of Neutrinos and Gamma-rays) is used to generate neutrino sources [124]. FIRESONG populates the universe with neutrino sources according to a cosmological model and generates the neutrino 111 Source type Φ0 [GeVcm−2 s−1 ] γ Reference Diffuse neutrino flux 1.44e-8 2.28 [126] TXS 0506+056 2014 flare 1.6e-8 2.2 [111] AGN HBL model 6.3e-10 1.1 [56] Table 6.2: Summary of the three different types of sources being simulated using the FIRESONG software package. The flux normalization, Φ0 , and spectral index, γ, assuming a power-law energy spectrum are provided for each source. fluxes measured at Earth of the simulated sources, using the relations in Section 6.2.4.1. It includes several different choices for source evolution models and luminosity functions. Here, the Candels and Clash model for redshift evolution is used [125] and a standard candle luminosity function is assumed. FIRESONG can be used in two different modes, one for steady sources and one for transient sources. Since this analysis is searching for flares, the transient source mode is used. Here flares with four different durations are studied: those lasting one, three, five, and ten days. FIRESONG also takes an energy spectrum as an input, making it possible to study different types of sources. If a luminosity value is not provided to the simulation, it will calculate it so that the input flux normalization is fully saturated. Here three different types of source spectra are studied. First is the assumption that each simulated neutrino source will reproduce the diffuse neutrino flux that was measured by IceCube in [126]. Second is the assumption that each simulated neutrino source behaves like the blazar TXS 0506+056 during its 2014 flare [111]. Third is the assumption that the simulated neutrinos sources behave like an HBL type AGN, whose emission peaks in the high-energy range [56]. Table 6.2 provides a summary of the spectral parameters used for each of the three source types. The final two parameters in FIRESONG are the local density and luminosity of neutrino sources, or the rate density and total isotropic equivalent energy in neutrinos in the transient mode. The simulation is run for a variety of rate density and total energy values in order to explore this parameter space, which is useful when studying populations of sources. A higher rate density means the sources can be less intense and a lower rate density means 112 the sources have to be brighter in order for them to be detectable. The simulation is run ten times for each combination of input rate density and total energy values. Note that when a luminosity or total energy value is given to FIRESONG, it ignores the spectral flux normalization parameter. The output of the simulation is a list of locations of the generated neutrino sources in declination, right ascension, and redshift, along with their expected neutrino fluxes at Earth. Since the goal is to search for gamma rays that originate from these same neutrino sources, it is necessary to convert the generated neutrino flux into a gamma-ray flux so that these sources can be studied with HAWC. To calculate the expected gamma-ray flux from the neutrino flux, the following relation is used, − d 2 X Eγ Fγ (Eγ ) = e λγγ Eν Fνα (Eν ), (6.17) 3K ν α where Fγ and Fνα are the gamma-ray and neutrino fluxes, respectively; Eγ = 2Eν are the gamma-ray and neutrino energies; α corresponds to the neutrino flavor; K is the ratio of charged to neutral pions; d is the distance to the source; and λγγ accounts for the attenuation of gamma rays due to the extra-galactic background light [127]. The constant K = 1 for proton-γ interactions and K = 2 for proton-proton interactions in the source; proton-γ interactions are assumed in this case. 6.2.4.3 Simulating Sources in HAWC Now that the neutrino sources were generated and their fluxes in gamma rays were calculated, it is possible to determine HAWC’s sensitivity to detecting these sources. The goal is to determine if HAWC would be able to detect sources that would produce a realtime neutrino alert in IceCube, assuming that the source produces gamma rays in addition to neutrinos. This is accomplished by injecting the simulated sources into HAWC’s data and running the Bayesian Block algorithm to determine how many can be detected. 113 A few cuts must be imposed on the simulation results generated in Section 6.2.4.2 before the sources can be studied with HAWC. The first is a cut on declination. FIRESONG generates sources at all declinations while HAWC’s visible declination ranges from −26◦ to 64◦ . The second is a cut on redshift, as HAWC can only see sources out to 0.3 in redshift due to attenuation by the EBL [128]. The final cut eliminates any sources that will not produce a realtime neutrino alert in IceCube. In order to determine if a simulated neutrino source would have generated a realtime alert in IceCube, it must first be determined how many neutrinos the simulated source is expected to produce. Since, as seen in Section 5.4, IceCube’s alerts are generated when a single neutrino is detected that is likely to be of astrophysical origin, the simulated neutrino source only needs to produce one neutrino that is detectable in IceCube. To calculate how many neutrinos IceCube expects to detect from each simulated source, the following equation is used, Z N= Aalert ef f (E, δ) × Φ(E) dE, (6.18) where N is the number of neutrinos, Aalertef f (E, δ) is the effective area of IceCube’s real time alerts as a function of energy and declination [100], and Φ(E) is the energy spectrum of each simulated neutrino source obtained from FIRESONG. If N > 1, then the simulated source passes the final cut. Then the simulated sources can be injected into HAWC’s data, which will help to deter- mine its sensitivity to such flares. This procedure is similar to the one described in Section 6.2.3, where the ZEBRA software is used to inject the simulated sources on top of an existing map that was chosen to be a daily map in the middle of HAWC’s range of data while the de- tector was running. Each injected source is assumed to follow a simple power-law spectrum with normalization equal to the gamma-ray flux given by FIRESONG’s output and index from Table 6.2 depending on the source type. It also takes into account EBL attenuation using the model from [129]. This results in a HEALPix map that contains the injected simulated source, from which 114 Figure 6.7: Example of a simulated flare injected into HAWC’s data. The Bayesian Block algorithm was run on the data including the injected flare, and it detected the simulated flare in this case. The black points are HAWC’s data with the simulated flare replacing the data during its length and the orange line shows the Bayesian Blocks. This example was a simulated source with the spectrum of TXS 0506+056 whose flare lasted for ten days. the expected flux at HAWC can be obtained following the procedure described in Section 3.2.2. Then HAWC’s daily flux data at the location of each simulated source is obtained. The simulated flare replaces the real data for the duration of the flare and then the Bayesian Blocks are run. Figure 6.7 shows an example of the Bayesian Blocks result tested on an example of a simulated flare. 6.3 Results The results from running the simulations are used to compute the sensitivity and discov- ery potential of this analysis. The number of expected flares from background fluctuations is assumed to be about one per year since the false positive rate was calibrated to be 5% and there are a total of 85 alerts. This results in a total background expectation of four flares since the simulation was run for four years to match the amount of HAWC’s data. To determine the sensitivity for a given rate density, the total energy value is found that corresponds to a distribution that crosses the median of the Poisson background distribution 90% of the time. To determine the 5σ discovery potential for a given rate density, the total 115 (a) 1-day flare (b) 3-day flare (c) 5-day flare (d) 10-day flare Figure 6.8: Sensitivity, discovery potential, and 90% upper limits for the AGN source type as a function of rate density and total isotropic energy equivalent in neutrinos. Each flare duration is shown separately. Results corresponding to a density value of 10−6 Mpc−3 yr−1 are only calculated for 1-day flares due to computational limitations. energy value is found that produces a p-value smaller than 2.87 × 10−7 with respect to the Poisson background distribution 50% of the time. The red dashed lines and blue dotted lines on the plot in Figures 6.8, 6.9, and 6.10 show the sensitivity and discovery potential, respectively, for each of the three sources and four flare durations tested. 116 (a) 1-day flare (b) 3-day flare (c) 5-day flare (d) 10-day flare Figure 6.9: Sensitivity, discovery potential, and 90% upper limits for the diffuse source type as a function of rate density and total isotropic energy equivalent in neutrinos. Each flare duration is shown separately. Results corresponding to a density value of 10−6 Mpc−3 yr−1 are only calculated for 1-day flares due to computational limitations. 117 (a) 1-day flare (b) 3-day flare (c) 5-day flare (d) 10-day flare Figure 6.10: Sensitivity, discovery potential, and 90% upper limits for the TXS 0506+056 source type as a function of rate density and total isotropic energy equivalent in neutri- nos. Each flare duration is shown separately. Results corresponding to a density value of 10−6 Mpc−3 yr−1 are only calculated for 1-day flares due to computational limitations. The search described in Section 6.2.2.3 was performed on the data described in Section 6.2.1. No flares were found by the Bayesian Block algorithm so a limit is set at the 90% confidence level. To calculate the upper limit, with νb = 4 background events, nobs = 0 up observed signal events, the following equation is solved for the upper limit νs , −ν up +ν Pnobs (νsup +νb )n e s b n=0 n! β= n , (6.19) −ν n obs ν b P e b n=0 n! 118 Figure 6.11: Comparison of the sensitivity, discovery potential, and 90% upper limit for a 1-day flare of all source types tested. The vertical lines correspond to rate densities of different source types. [35–38]. where β = 0.1, which is the confidence level subtracted from 1 [130]. This results in an up upper limit of nus = 2.3. In order to interpret this result in terms of the parameter space in Figures 6.8, 6.9, and 6.10, a similar procedure to the one used to determine the sensitivity and discovery potential is used. The 90% upper limit is the total energy value that results up in more than νs + νb = 6.3 detections 90% of the time. These limits are also shown in the plots in Figures 6.8, 6.9, and 6.10. Figure 6.11 compares the sensitivity, discovery potential, and upper limits of all three source models assuming a flare with a duration of one day. This is effectively a study of the systematic effects of changing the spectral index in the simulation, since all other parameters between source types are the same. It can be concluded that HAWC’s sensitivity to such sources decreases by a few orders of magnitude for harder spectral indices. 119 6.4 Conclusions A search for gamma-ray flares in four years of HAWC data at the locations of neutrino alerts from ten years of IceCube’s realtime system was conducted using Bayesian Blocks. A simulation was performed to determine the sensitivity and 5σ discovery potential. No flares were found so upper limits were placed at the 90% confidence level on the rate density of flaring sources that produce both neutrinos and gamma rays. The method used to calculate the upper limits were adapted from a previously pub- lished paper that also used both HAWC and IceCube data to study multi-messenger astro- physics [110]. However, this procedure overcounts fluctuations and results in too conservative of a limit. Instead, the methods described by Feldman and Cousins in [131] are a more ap- propriate procedure to calculate the upper limits, since it provides the correct coverage of confidence intervals. The tables provided in [131] are used to determine the upper limit assuming b = 4 background events and n0 = 0 signal events. The resulting 90% upper limit is 1.01, and to interpret this in the parameter space of total energy and rate density, a similar procedure to the one described above is used. The 90% upper limit is the total energy value that results in a mean value of 1.01 detection over the 10 trials run. The resulting upper limits using this procedure result in an improved limit and are shown in the bright green lines in Figure 6.12. Future studies are needed to evaluate these differences. There may also be further improvement in the estimation of background events. The fact that no detections were observed in the data when a background of four was expected indicates that the estimated background may be too high. To study this effect further, 85 random locations at the same declinations of the 85 IceCube alerts were chosen. HAWC sources were masked out and the Bayesian Block algorithm was run on the flux values at the chosen locations. This procedure was repeated 100 times and no detections were found by the Bayesian Blocks. This suggests that further investigation is needed to fully understand the expected background and false positive rate. 120 Figure 6.12: Comparison of the two procedures of calculating 90% upper limits for a 1-day flare. The limits from the previously published method described in Section 6.3 are shown in the dark green lines and the limits from the Feldman Cousins method described in Section 6.4 are shown in the bright green lines. The sensitivity and discovery potential of all source types tested are also included in red and blue lines, respectively. The vertical lines correspond to rate densities of different source types. [35–38]. There are several additional possible ways to improve this analysis in the future. One way would be to study different methods besides the Bayesian Blocks to detect the flares in HAWC’s data. One such possibility is to train a neural network to identify flares from the data. Another simple way to improve the sensitivity of this analysis is to add more data from both HAWC and IceCube. On the IceCube side, there is about another year of alerts that can be searched for in the HAWC data. IceCube has also recently added a cascade alert stream to the realtime system, so those alerts can be used as well. In addition, HAWC collects data every day so there is about another year and a half of data at the time of this writing that can be added to the analysis. There are improvements to the reconstruction algorithms that are currently in progress and data from the outrigger array that can be included as well. Other gamma-ray instruments may be able to improve upon the sensitivity 121 of this analysis as well. For example, SWGO [83] will be able to study the neutrino alerts that occurred in the Southern hemisphere, and LHAASO, the Large High Altitude Air Shower Observatory [53], has a higher sensitivity in a similar declination range as HAWC. Finally, this analysis can be modified in order to run it in realtime. After HAWC receives an alert from IceCube, it can automatically collect flux data for each day in HAWC’s data set from the location of the neutrino alert in the sky and run the Bayesian Block algorithm to try to detect any gamma-ray flares that occurred at any time at that location. If a flare was detected, it would allow other observatories such as the ones mentioned in Section 5.4 to either search through their archival data or to perform follow-up observations on any sources revealed by HAWC, thus enabling multi-messenger astronomy. 122 CHAPTER 7 CONCLUSIONS AND FUTURE OUTLOOK In this work, three analyses were performed with data from the HAWC and IceCube obser- vatories in order to explore some of the fundamental questions in particle astrophysics. The goal was to contribute to the understanding of the sources of cosmic rays, their acceleration mechanisms, and the nature of dark matter. The detection techniques and reconstruction methods for both observatories were described, along with the properties that make them ideal instruments to perform such searches. In Chapter 3, a catalog was composed of all sources in HAWC’s field of view, which allows for the identification of some of the sources that produce cosmic rays. A total of 65 sources were found using 1523 days of data and their locations, extensions, distance to the nearest TeVCat source, and spectral fit parameters were listed. Potential associations with other catalogs including the 4FGL catalog of MeV gamma-ray sources, the ATNF pulsar catalog, and the SNRCat catalog of Galactic supernova remnants were investigated as well. Current studies are underway to improve the catalog search methods in order to more accurately fit sources in crowded regions of the Galactic plane and to account for the galactic diffuse emission. Chapter 4 describes a search for evaporating primordial black holes in 959 days of HAWC data. A detection of a PBH burst would prove their existence and the theory of Hawking radiation and would provide a target with which to search for dark matter. Even a non- detection would place important constraints on the early universe and the physics of inflation. To conduct this analysis, a model of PBHs was created from a simulation and compared to data from a previous transient search in HAWC. No burst events were found, so upper limits were placed on the burst rate density. These limits were the strongest constraints set on this quantity at the time of publication [21]. Future studies include a search optimized to PBHs on HAWC’s raw data, testing more burst durations, and improved likelihood methods. 123 Some of the processes that produce gamma rays are also expected to produce neutrinos. In Chapter 6, a multi-messenger search for flares in gamma rays from known locations of probable high-energy astrophysical neutrinos was conducted using data from both HAWC and IceCube. Finding such emission would indicate a hadronic acceleration mechanism and could help to start revealing some neutrino sources. Three different spectral models were studied: the diffuse flux, the blazar TXS 0506+056, and an AGN of HBL type. No such gamma-ray flares were found, so upper limits were set on the burst density of sources that could produce both gamma-ray and neutrino emission. Studying new methods for detecting flares in addition to a new catalog of neutrino alert events from IceCube would help to improve this analysis in the future. There are even more improvements possible on the work presented in this thesis due to HAWC’s progress in improving its reconstruction algorithms. The newest set of recon- struction algorithms will be available shortly, which will improve the angular and energy reconstruction. The next set of reconstruction algorithms are also currently under develop- ment, which will incorporate the outriggers mentioned in Section 2.1. This will improve the angular reconstruction even further and also extend the high end of HAWC’s energy range beyond 100 TeV. Additionally, both HAWC and IceCube are continuing to add more data every day due to their constant operations. This means that their sensitivity to sources will keep increas- ing over time. There are many other currently operating and future planned detectors that observe different messenger particles or different energy ranges. Increasing communications between the different instruments to combine observations increases the potential to dis- cover more sources and study their acceleration mechanisms. The era of multi-messenger astronomy is just beginning with many more exciting results to come. 124 BIBLIOGRAPHY 125 BIBLIOGRAPHY [1] A. U. Abeysekara et al., “Observation of the Crab Nebula with the HAWC gamma-ray observatory,” The Astrophysical Journal, vol. 843, p. 39, jun 2017. [2] A. Albert et al., “3HWC: the third HAWC catalog of very-high-energy gamma-ray sources,” The Astrophysical Journal, vol. 905, p. 76, dec 2020. [3] S. Abdollahi et al., “Fermi large area telescope fourth source catalog,” The Astrophys- ical Journal Supplement Series, vol. 247, p. 33, mar 2020. [4] R. N. Manchester, G. B. Hobbs, A. Teoh, and M. Hobbs, “The australia telescope national facility pulsar catalogue,” The Astronomical Journal, vol. 129, pp. 1993–2006, apr 2005. [5] G. Ferrand and S. Safi-Harb, “A census of high-energy observations of Galactic super- nova remnants,” Advances in Space Research, vol. 49, pp. 1313–1319, may 2012. [6] P. Zyla et al., “Review of Particle Physics,” PTEP, vol. 2020, no. 8, p. 083C01, 2020. [7] M. Spurio, Particles and Astrophysics: A Multi-Messenger Approach. Springer, 2015. [8] A. Bennun, “High energy dimensioning the quantum space-time of the electron,” ””, 06 2020. [9] B. Carlson, Terrestrial gamma-ray flash production by lightning. PhD thesis, Stanford University, 01 2009. [10] J. R. Primack, A. Domı́nguez, R. C. Gilmore, and R. S. Somerville, “Extragalactic background light and gamma-ray attenuation,” AIP Conference Proceedings, vol. 1381, no. 1, pp. 72–83, 2011. [11] “TeVCat homepage.” http://tevcat.uchicago.edu/. [12] J. Hester and A. Loll, 2000. [13] C. Urry and P. Padovani, 2003. [14] D. J. Thompson, “Space detectors for gamma rays (100 mev–100 gev): From egret to fermi lat,” Comptes Rendus Physique, vol. 16, no. 6, pp. 600–609, 2015. Gamma-ray astronomy / Astronomie des rayons gamma. [15] A. Horvath, “Chernkov radiation,” March 2006. [16] “The IceCube Collaboration, IceCube Internal Gallery..” https://gallery.icecube. wisc.edu/. Accessed: 10/29/21. [17] “The IceCube Collaboration, IceCube Internal Gallery..” https://gallery.icecube. wisc.edu/. Accessed: 9/1/21. 126 [18] Z. Hampel-Arias, Cosmic Ray Observations at the TeV Scale with the HAWC Obser- vatory. PhD thesis, University of Wisconsin-Madison, 2017. [19] J. Pretz, “A.T.H.E.N.A. advanced tracking of HAWC experiment notifications and alerts.” HAWC Internal Technical Note, 2015. [20] T.N.Ukwatta, D.R.Stump, J.T.Linnemann, J.H.MacGibbon, S.S.Marinelli, T.Yapici, and K.Tollefson, “Primordial black holes: Observational characteristics of the final evaporation,” Astroparticle Physics, vol. 80, pp. 90–114, jul 2016. [21] A. Albert et al., “Constraining the local burst rate density of primordial black holes with HAWC,” Journal of Cosmology and Astroparticle Physics, vol. 2020, pp. 026–026, apr 2020. [22] E. T. Linton and et al., “A new search for primordial black hole evaporations using the Whipple gamma-ray telescope,” Journal of Cosmology and Astroparticle Physics, vol. 01, p. 013, 2006. [23] D. E. Alexandreas and et al., “New limit on the rate-density of evaporating black holes,” Physical Review Letters, vol. 71, no. 16, p. 2524, 1993. [24] M. Amenomori et al., “Search for 10 TeV Gamma Bursts from Evaporating Primordial Black Holes with the Tibet Air Shower Array,” in Proceedings, 24th International Cosmic Ray Conference (ICRC1995): Rome, Italy, August 28 - September 8, 1995, p. 112, 1995. [25] J.-F. Glicenstein, A. Barnacka, M. Vivier, and T. Herr, “Limits on Primordial Black Hole evaporation with the H.E.S.S. array of Cherenkov telescopes,” in Proceedings, 33rd International Cosmic Ray Conference (ICRC2013): Rio de Janeiro, Brazil, July 2-9, 2013, p. 0930, 2013. [26] S. Archambault, “Search for Primordial Black Hole Evaporation with VERITAS,” PoS, vol. ICRC2017, p. 691, 2018. [27] M. Ackermann et al., “Search for Gamma-Ray Emission from Local Primordial Black Holes with the Fermi Large Area Telescope,” Astrophys. J., vol. 857, no. 1, p. 49, 2018. [28] A. A. Abdo et al., “Milagro Limits and HAWC Sensitivity for the Rate-Density of Evaporating Primordial Black Holes,” Astropart. Phys., vol. 64, pp. 4–12, 2015. [29] “The IceCube Collaboration, IceCube Internal Gallery..” https://gallery.icecube. wisc.edu/. Accessed: 9/1/21. [30] M. Aartsen et al., “The IceCube neutrino observatory: instrumentation and online systems,” Journal of Instrumentation, vol. 12, pp. P03012–P03012, mar 2017. [31] “The IceCube Collaboration, IceCube Internal Gallery..” https://gallery.icecube. wisc.edu/. Accessed: 9/1/21. 127 [32] M. G. Aartsen et al., “IceCube-gen2: the window to the extreme universe,” Journal of Physics G: Nuclear and Particle Physics, vol. 48, p. 060501, apr 2021. [33] M. Ahlers, K. Helbin, and C. Pérez de los Heros, “Probing particle physics with Ice- Cube,” The European Physical Journal C, vol. 78, 2018. [34] “The IceCube Collaboration, IceCube Masterclass..” https://masterclass.icecube. wisc.edu/en/learn/detecting-neutrinos. Accessed: 9/23/21. [35] D. Wanderman and T. Piran, “The rate, luminosity function and time delay of non- Collapsar short GRBs,” Monthly Notices of the Royal Astronomical Society, vol. 448, pp. 3026–3037, 03 2015. [36] G. R. Farrar and A. Gruzinov, “GIANT AGN FLARES AND COSMIC RAY BURSTS,” The Astrophysical Journal, vol. 693, pp. 329–332, mar 2009. [37] L.-G. Strolger, T. Dahlen, S. A. Rodney, O. Graur, A. G. Riess, C. McCully, S. Ravin- dranath, B. Mobasher, and A. K. Shahady, “THE RATE OF CORE COLLAPSE SUPERNOVAE TO REDSHIFT 2.5 FROM THE CANDELS AND CLASH SUPER- NOVA SURVEYS,” The Astrophysical Journal, vol. 813, p. 93, nov 2015. [38] K. Murase and M. Fukugita, “Energetics of high-energy cosmic radiations,” Phys. Rev. D, vol. 99, p. 063012, Mar 2019. [39] V. Hess, “On the observations of the penetrating radiation during seven balloon flights,” 2018. [40] K. Greisen, “End to the cosmic-ray spectrum?,” Phys. Rev. Lett., vol. 16, pp. 748–750, Apr 1966. [41] G. T. Zatsepin and V. A. Kuz’min, “Upper Limit of the Spectrum of Cosmic Rays,” Soviet Journal of Experimental and Theoretical Physics Letters, vol. 4, p. 78, Aug. 1966. [42] E. Fermi, “On the origin of the cosmic radiation,” Phys. Rev., vol. 75, pp. 1169–1174, Apr 1949. [43] R. Mirzoyan, “Major change in understanding of grbs at tev,” 2020. [44] P. Mészáros, “The fireball model of gamma-ray bursts,” Progress of Theoretical Physics Supplement, no. 143, pp. 33–49, 2001. [45] S. W. Hawking, “Black hole explosions,” Nature, vol. 248, pp. 30–31, 1974. [46] R. C. Hartman et al., “The third EGRET catalog of high-energy gamma-ray sources,” The Astrophysical Journal Supplement Series, vol. 123, pp. 79–202, jul 1999. [47] J. Aleksić et al., “The major upgrade of the magic telescopes, part ii: A performance study using observations of the crab nebula,” Astroparticle Physics, vol. 72, pp. 76–94, 2016. 128 [48] J. Holder and othersr, “The first veritas telescope,” Astroparticle Physics, vol. 25, no. 6, pp. 391–401, 2006. [49] B. Acharya et al., Science with the Cherenkov Telescope Array. WSP, 11 2018. [50] P. Bernardini and the ARGO-YBJ Collaboration, “ARGO-YBJ experiment in tibet,” Journal of Physics: Conference Series, vol. 120, p. 062022, jul 2008. [51] M. Amenomori et al., “Search for steady emission of 10-tev gamma rays from the crab nebula, cygnus x-3, and hercules x-1 using the tibet air shower array,” Phys. Rev. Lett., vol. 69, pp. 2468–2471, Oct 1992. [52] G. B. Yodh, “The milagro gamma ray observatory,” Nuclear Physics B - Proceedings Supplements, vol. 52, no. 3, pp. 264–268, 1997. [53] Z. Cao et al., “Ultrahigh-energy photons up to 1.4 petaelectronvolts from 12 γ-ray galactic sources,” Nature, vol. 594, pp. 33–36, jun 2021. [54] B. P. Abbott et al., “LIGO: the laser interferometer gravitational-wave observatory,” Reports on Progress in Physics, vol. 72, p. 076901, jun 2009. [55] V. Collaboration, “Advanced virgo baseline design,” 2009. [56] P. Padovani, M. Petropoulou, P. Giommi, and E. Resconi, “A simplified view of blazars: the neutrino background,” Monthly Notices of the Royal Astronomical Society, vol. 452, pp. 1877–1887, 07 2015. [57] M. G. Aartsen et al., “Multimessenger observations of a flaring blazar coincident with high-energy neutrino icecube-170922a,” Science, vol. 361, no. 6398, p. eaat1378, 2018. [58] F. Aharonian, J. Buckley, T. Kifune, and G. Sinnis, “High energy astrophysics with ground-based gamma ray detectors,” Reports on Progress in Physics, vol. 71, p. 096901, aug 2008. [59] V. Marandon, A. Jardin-Blicq, and H. Schoorlemmer, “Latest news from the HAWC outrigger array,” 2019. [60] A. Abeysekara et al., “Data acquisition architecture and online processing system for the HAWC gamma-ray observatory,” Nuclear Instruments and Methods in Physics Research Section A: Accelerators, Spectrometers, Detectors and Associated Equipment, vol. 888, pp. 138–146, 2018. [61] “ZeroMQ homepage.” https://zeromq.org/. [62] “MySQL homepage.” https://www.mysql.com/. [63] “PHP homepage.” https://www.php.net/. [64] “Google Charts homepage.” https://developers.google.com/chart. 129 [65] K. Greisen, “Cosmic ray showers,” Annual Review of Nuclear Science, vol. 10, no. 1, pp. 63–108, 1960. [66] K. Kamata and J. Nishimura, “The lateral and the angular structure functions of electron showers,” Progress of Theoretical Physics Supplement, vol. 6, pp. 93–155, feb 1958. [67] A. U. Abeysekara et al., “Measurement of the Crab Nebula spectrum past 100 TeV with HAWC,” The Astrophysical Journal, vol. 881, p. 134, aug 2019. [68] R. Atkins et al., “Observation of TeV gamma rays from the Crab Nebula with Milagro using a new background rejection technique,” The Astrophysical Journal, vol. 595, pp. 803–311, oct 2003. [69] A. U. Abeysekara et al., “The 2HWC HAWC observatory gamma-ray catalog,” The Astrophysical Journal, vol. 843, p. 40, jun 2017. [70] A. U. Abeysekara et al., “Search for TeV gamma-ray emission from point-like sources in the inner GALACTIC plane with a partial configuration of the HAWC OBSERVA- TORY,” The Astrophysical Journal, vol. 817, p. 3, jan 2016. [71] K. M. Gorski, E. Hivon, A. J. Banday, B. D. Wandelt, F. K. Hansen, M. Reinecke, and M. Bartelmann, “HEALPix: A framework for high-resolution discretization and fast analysis of data distributed on the sphere,” The Astrophysical Journal, vol. 622, pp. 759–771, apr 2005. [72] P. W. Younk, R. J. Lauer, G. Vianello, J. P. Harding, H. A. A. Solares, H. Zhou, and M. Hui, “A high-level analysis framework for hawc,” 2015. [73] S. S. Wilks, “The large-sample distribution of the likelihood ratio for testing composite hypotheses,” The Annals of Mathematical Statistics, vol. 9, no. 1, pp. 60–62, 1938. [74] I. M. Castellanos, Search for Gamma-ray Counterparts of Gravitational Wave Events and Other Transient Signals with HAWC. PhD thesis, University of Maryland, 2019. [75] S. P. Wakely and D. Horan, “TeVCat: An online catalog for Very High Energy Gamma- Ray Astronomy,” in International Cosmic Ray Conference, vol. 3 of International Cosmic Ray Conference, pp. 1341–1344, Jan. 2008. [76] “ATNF pulsar catalog homepage.” https://www.atnf.csiro.au/research/pulsar/ psrcat/. [77] “SNRcat homepage.” http://snrcat.physics.umanitoba.ca/. [78] A. Albert et al., “Spectrum and morphology of the very-high-energy source HAWC J2019+368,” The Astrophysical Journal, vol. 911, p. 143, apr 2021. [79] A. Abeysekara et al., “HAWC observations of the acceleration of very-high-energy cosmic rays in the Cygnus Cocoon,” Nature Astronomy, vol. 5, pp. 465–471, may 2021. 130 [80] A. Albert et al., “Evidence of 200 TeV photons from HAWC J1825-134,” The Astro- physical Journal, vol. 907, p. L30, jan 2021. [81] S. C. de León, A. C. Alonso, D. Rosa-González, and A. L. Longinotti, “Spectral anal- ysis of the blazars Markarian 421 and Markarian 501 with the HAWC Gamma-Ray Observatory,” 2019. [82] A. Nayerhoda, F. Salesa Greusa, S. Casanova, D. Gaggeroc, D. Grasso, and A. Marinelli, “Gamma Ray Diffuse Emission from the Galactic Plane with HAWC Data,” PoS, vol. ICRC2019, p. 750, 2019. [83] A. Albert et al., “Science Case for a Wide Field-of-View Very-High-Energy Gamma- Ray Observatory in the Southern Hemisphere,” 2019. [84] B. J. Carr, K. Kohri, Y. Sendouda, and J. Yokoyama, “New cosmological constraints on primordial black holes,” Phys. Rev. D, vol. 81, p. 104019, May 2010. [85] B. Carr, F. Kuhnel, and M. Sandstad, “Primordial Black Holes as Dark Matter,” Phys. Rev. D, vol. 94, no. 8, p. 083504, 2016. [86] D. N. Page and S. W. Hawking, “Gamma rays from primordial black holes,” ApJ, vol. 206, pp. 1 – 7, 1976. [87] F. Halzen, E. Zas, J. H. MacGibbon, and T. C. Weekes, “Gamma-rays and energetic particles from primordial black holes,” Nature, vol. 353, pp. 807–815, 1991. [88] E. L. Wright, “On the Density of Primordial Black Holes in the Galactic Halo,” ApJ, vol. 459, p. 487, Mar. 1996. [89] K. Abe, H. Fuke, S. Haino, T. Hams, M. Hasegawa, A. Horikoshi, K. C. Kim, A. Kusumoto, M. H. Lee, Y. Makida, S. Matsuda, Y. Matsukawa, J. W. Mitchell, J. Nishimura, M. Nozaki, R. Orito, J. F. Ormes, K. Sakai, M. Sasaki, E. S. Seo, R. Shinoda, R. E. Streitmatter, J. Suzuki, K. Tanaka, N. Thakur, T. Yamagami, A. Ya- mamoto, T. Yoshida, and K. Yoshimura, “Measurement of the cosmic-ray antiproton spectrum at solar minimum with a long-duration balloon flight over antarctica,” Phys. Rev. Lett., vol. 108, p. 051102, Jan 2012. [90] J. H. MacGibbon and B. R. Webber, “Quark- and gluon-jet emission from primordial black holes: The instantaneous spectra,” Phys. Rev. D, vol. 41, pp. 3052 – 3079, 1990. [91] S. W. Hawking, “Particle creation by black holes,” Communications in Mathematical Physics, vol. 43, no. 3, pp. 199 – 220, 1975. [92] J. Wood, An All-sky Search for Bursts of Very High Energy Gamma Rays with HAWC. PhD thesis, University of Maryland, 2016. [93] A. Abdo et al., “Milagro limits and hawc sensitivity for the rate-density of evaporating primordial black holes,” Astroparticle Physics, vol. 64, pp. 4–12, 2015. 131 [94] R. López-Coto, M. Doro, A. de Angelis, M. Mariotti, and J. P. Harding, “Prospects for the observation of primordial black hole evaporation with the southern wide field of view gamma-ray observatory,” 2021. [95] A. Ishihara, “The IceCube Upgrade - Design and Science Goals,” PoS, vol. ICRC2019, p. 1031, 2019. [96] M. G. Aartsen et al., “Energy reconstruction methods in the IceCube neutrino tele- scope,” Journal of Instrumentation, vol. 9, pp. P03009–P03009, mar 2014. [97] P. J. Huber, “Robust Estimation of a Location Parameter,” The Annals of Mathemat- ical Statistics, vol. 35, no. 1, pp. 73 – 101, 1964. [98] R. Abbasi et al., “A muon-track reconstruction exploiting stochastic losses for large- scale cherenkov detectors,” Journal of Instrumentation, vol. 16, p. P08034, aug 2021. [99] F. James and M. Roos, “Minuit - a system for function minimization and analysis of the parameter errors and correlations,” Computer Physics Communications, vol. 10, no. 6, pp. 343–367, 1975. [100] E. Blaufuss, T. Kintscher, L. Lu, and C. F. Tung, “The Next Generation of IceCube Real-time Neutrino Alerts,” in Proceedings of 36th International Cosmic Ray Confer- ence — PoS(ICRC2019), vol. 358, p. 1021, 2019. [101] M. Aartsen et al., “The icecube realtime alert system,” Astroparticle Physics, vol. 92, pp. 30–41, 2017. [102] “Gcn: The gamma-ray coordinates network.” https://gcn.gsfc.nasa.gov/. Ac- cessed: 10/1/2021. [103] E. C. Bellm, “The zwicky transient facility,” 2014. [104] O. Pace, D. Pawlak, and C. Winkler, “INTEGRAL: Mission and Satellite,” The As- trophysical Journal Supplement Series, vol. 92, p. 339, June 1994. [105] M. Ageron et al., “Antares: The first undersea neutrino telescope,” Nuclear Instru- ments and Methods in Physics Research Section A: Accelerators, Spectrometers, De- tectors and Associated Equipment, vol. 656, no. 1, pp. 11–38, 2011. [106] D. N. Burrows et al., “The swift x-ray telescope,” Space Science Reviews, vol. 120, pp. 165–195, 2005. [107] M. Lacy et al., “The karl g. jansky very large array sky survey (VLASS). science case and survey design,” Publications of the Astronomical Society of the Pacific, vol. 132, p. 035001, jan 2020. [108] M. Aartsen et al., “Evidence for high-energy extraterrestrial neutrinos at the icecube detector,” Science, vol. 342, no. 6161, p. 1242856, 2013. [109] “AMON homepage.” https://www.amon.psu.edu/. 132 [110] H. A. A. S. andothers, “Multimessenger gamma-ray and neutrino coincidence alerts using HAWC and IceCube subthreshold data,” The Astrophysical Journal, vol. 906, p. 63, jan 2021. [111] M. Aartsen et al., “Neutrino emission from the direction of the blazar txs 0506+056 prior to the icecube-170922a alert,” Science, vol. 361, no. 6398, pp. 147–151, 2018. [112] A. U. Abeysekara et al., “Daily monitoring of TeV gamma-ray emission from mrk 421, mrk 501, and the crab nebula with HAWC,” The Astrophysical Journal, vol. 841, p. 100, may 2017. [113] A. U. Abeysekara et al., “The HAWC real-time flare monitor for rapid detection of transient events,” The Astrophysical Journal, vol. 843, p. 116, jul 2017. [114] M. Meyer, J. D. Scargle, and R. D. Blandford, “Characterizing the gamma-ray vari- ability of the brightest flat spectrum radio quasars observed with the fermi LAT,” The Astrophysical Journal, vol. 877, p. 39, may 2019. [115] L. Lin, E. Göǧüş, Y. Kaneko, and C. Kouveliotou, “DETAILED INVESTIGATIONS OF THE DIMMEST BURSTS FROM TWO MAGNETARS, SGR j0501+4516 AND SGR j1550–5418,” The Astrophysical Journal, vol. 778, p. 105, nov 2013. [116] A. von Kienlin, P. Veres, O. J. Roberts, R. Hamburg, E. Bissaldi, M. S. Briggs, E. Burns, A. Goldstein, D. Kocevski, R. D. Preece, C. A. Wilson-Hodge, C. M. Hui, B. Mailyan, and C. Malacaria, “Fermi-GBM GRBs with characteristics similar to GRB 170817a,” The Astrophysical Journal, vol. 876, p. 89, may 2019. [117] J. D. Scargle, J. P. Norris, B. Jackson, and J. Chiang, “STUDIES IN ASTRONOM- ICAL TIME SERIES ANALYSIS. VI. BAYESIAN BLOCK REPRESENTATIONS,” The Astrophysical Journal, vol. 764, p. 167, feb 2013. [118] I. Taboada, C. F. Tung, and J. Wood, “Constrains on the extragalactic origin of IceCube’s neutrinos using HAWC,” in Proceedings of 35th International Cosmic Ray Conference — PoS(ICRC2017), vol. 301, p. 663, 2017. [119] M. Kowalski, “Status of high-energy neutrino astronomy,” Journal of Physics: Con- ference Series, vol. 632, p. 012039, aug 2015. [120] D. W. Hogg, “Distance measures in cosmology,” 2000. [121] P. A. R. Ade et al., “Planck 2015 results - xiii. cosmological parameters,” A&A, vol. 594, p. A13, 2016. [122] K. Murase and E. Waxman, “Constraining high-energy cosmic neutrino sources: Im- plications and prospects,” Phys. Rev. D, vol. 94, p. 103006, Nov 2016. [123] M. Ahlers and F. Halzen, “Pinpointing extragalactic neutrino sources in light of recent icecube observations,” Phys. Rev. D, vol. 90, p. 043005, Aug 2014. 133 [124] C. F. Tung, T. Glauch, M. Larson, A. Pizzuto, R. Reimann, and I. Taboada, “Firesong: A python package to simulate populations of extragalactic neutrino sources,” Journal of Open Source Software, vol. 6, no. 61, p. 3194, 2021. [125] L.-G. Strolger, T. Dahlen, S. A. Rodney, O. Graur, A. G. Riess, C. McCully, S. Ravin- dranath, B. Mobasher, and A. K. Shahady, “THE RATE OF CORE COLLAPSE SUPERNOVAE TO REDSHIFT 2.5 FROM THE CANDELS AND CLASH SUPER- NOVA SURVEYS,” The Astrophysical Journal, vol. 813, p. 93, nov 2015. [126] J. Stettner, “Measurement of the diffuse astrophysical muon-neutrino spectrum with ten years of IceCube data,” in Proceedings of 36th International Cosmic Ray Confer- ence — PoS(ICRC2019), vol. 358, p. 1017, 2019. [127] M. Ahlers and K. Murase, “Probing the galactic origin of the icecube excess with gamma rays,” Phys. Rev. D, vol. 90, p. 023010, Jul 2014. [128] A. Albert et al., “A survey of active galaxies at TeV photon energies with the HAWC gamma-ray observatory,” The Astrophysical Journal, vol. 907, p. 67, jan 2021. [129] A. Domı́nguez, J. R. Primack, D. J. Rosario, F. Prada, R. C. Gilmore, S. M. Faber, D. C. Koo, R. S. Somerville, M. A. Pérez-Torres, P. Pérez-González, J.-S. Huang, M. Davis, P. Guhathakurta, P. Barmby, C. J. Conselice, M. Lozano, J. A. Newman, and M. C. Cooper, “Extragalactic background light inferred from AEGIS galaxy-SED-type fractions,” Monthly Notices of the Royal Astronomical Society, vol. 410, pp. 2556–2578, 01 2011. [130] G. Cowan, Statistical Data Analysis. Oxford, 1998. [131] G. J. Feldman and R. D. Cousins, “Unified approach to the classical statistical analysis of small signals,” Phys. Rev. D, vol. 57, pp. 3873–3889, Apr 1998. 134