Detection of Pulsed Very High Energy Gamma-Rays from the Crab Pulsar with the MAGIC ...

Diss. ETH No. 18322

Detection of Pulsed Very High Energy Gamma-Rays from the Crab Pulsar with the MAGIC telescope using an Analog Sum Trigger

A dissertation submitted to the

Swiss Federal Institute of Technology Zurich for the degree of

Doctor of Natural Sciences

presented by

Michael Thomas Rissi Dipl. Phys. ETH born August 31, 1979 citizen of Wartau (SG)

accepted on the recommendation of Prof. Dr. Felicitas Pauss, examiner Prof. Dr. Christoph Grab, co-examiner Dr. Adrian Biland, co-examiner 2009

Für Regina

Je mehr man kennt, je mehr man weiss, erkennt man: alles dreht im Kreis. (Goethe, Zahme Xenien VI)

Abstract Pulsars (PULSating stARS) are strongly magnetized and fast spinning neutron stars emerging from the massive core of a star after the supernova at the end of its lifetime. Their rotational frequency ranges from 0.1 Hz to nearly 1 kHz. In the 1960s, Pulsars were discovered at Radio frequencies. Until the mid 1990s, pulsating electromagnetic radiation was measured over the whole spectrum (radio, infrared, optical, ultraviolet, X-rays) from a few 1000s pulsars. Seven pulsars are even known to radiate pulsating γ-rays, being the photons with the shortest yet measured wavelengths and the highest energies ever observed. Even though pulsars have been known for the last 40 years, the mechanisms leading to the emission of electromagnetic radiation is hardly known. The main models – the so-called polar cap and outer gap models – have in common that the emission of electromagnetic radiation is related to relativistic electrons and positrons. These particles are accelerated by the electric field induced by the strong rotating magnetic field. Photons are emitted by those particles, e.g. by means of the synchrotron mechanism. The difference between the two models lies in the region, where the acceleration of the charged particles takes place. The polar cap model assumes the acceleration in a region close to the neutron star above the magnetic poles, while in the outer gap model, the acceleration takes place in an area far out in the pulsar’s magnetosphere. The magnetosphere is the space in the vicinity of a pulsar, where a strong magnetic dipole field co-rotates with the pulsar. It is filled with a plasma. Due to different absorption mechanism for γ-rays, the two models predict different energy spectra for the pulsed radiation at the highest energies. Gamma rays with an energy in the MeV to a few GeV range are usually called High Energy (HE) γ-rays and are measured by satellite experiments. Photons with energies above a few tens of GeV are called Very High Energy (VHE) γ-rays and can be detected from ground with Imaging Air Cherenkov Telescopes (IACT) such as the MAGIC (Major Atmospheric Gamma-ray Imaging Cherenkov) Telescope. These kind of telescopes detect the Cherenkov flashes from electromagnetic atmospheric air showers initiated by a γ-ray hitting molecules or atoms in the upper atmosphere. Among all IACTs, MAGIC has the lowest energy threshold of ≈ 55 GeV. In the last 20 years since the emergence of VHE γ-ray astrophysics with ground based instruments, several unsuccessful attempts to detect VHE γ-rays from pulsars have been made. The reasons for the non-detection is the sharp cut-off of the energy spectra predicted by the two pulsar models and the too high energy threshold of IACTs which did not allow to collect enough pulsating γ-rays. Following a hint for pulsed VHE γ-ray emission above 60 GeV from the Crab Pulsar from data taken with the MAGIC telescope in 2005, we developed a new trigger system for the MAGIC telescope: the sum trigger. The sum trigger allowed to reduce the energy threshold from 55 GeV to 25 GeV. Between October 2007 and February 2008, we observed the Crab Pulsar with the MAGIC telescope using the sum trigger for about 22 hours under good weather conditions. The V

analysis of this data reveals pulsating VHE γ-radiation on a confidence level of 6.2σ. In case of the Crab Pulsar, the two pulsar models discussed above both predict a cut-off of the pulsed γ-ray spectrum of exponential (outer gap) or super exponential (polar cap) shape in the range between a few GeV and some tens of GeV. For both cases, we calculated the cut-off energy to Ecut = (17.7 ± 2.8stat ± 5.0syst ) GeV (exponential), and Ecut = (23.2 ± 2.9stat ± 6.6syst ) GeV (super exponential), respectively. Our measurements contradict the polar cap’s model predictions. In a next step, the pulsed spectrum was determined in the cut-off region. I found that the pulsar still emits pulsed γ-rays above 60 GeV which is higher than predictions from the outer gap model indicate. The pulsed spectrum above 25 GeV rather follows a power law (f (E) ∝ E −α ) with spectral index α = 3.1 ± 0.6stat ± 0.5syst , than the steep exponential cut-off.

VI

Zusammenfassung Pulsare sind stark magnetisierte, rotierende Neutronensterne, welche aus den massiven Kernen von Sternen nach der Supernova am Ende derer Lebenszeit entstehen. Ihre Rotationsfrequenz ist zwischen 0.1 Hz bis fast 1 kHz. Gepulste Radioemission wurde vom ersten Pulsar in den 1960er Jahren entdeckt. Gepulste elektromagnetische Strahlung wurde bis Mitte der 1990er Jahre über das gesamte Spektrum (von Radio über Infrarot, optischen Wellenlängen, Ultraviolett bis Röntgenstrahlen) von ein paar 1000 Pulsaren gemessen. Von sieben Pulsaren ist sogar bekannt, dass sie γ-Strahlung emittieren. Diese Photonen haben die kürzeste je gemessene Wellenlänge und die höchsten bisher beobachteten Energien. Obwohl Pulsare seit 40 Jahren bekannt sind, wissen wir heute immer noch relativ wenig über deren Emissionsmechanismus. Die beiden Hauptmodelle (das “polar cap” und das “outer gap” Modell) haben gemein, dass die Emission der elektromagnetischen Strahlung mit relativistischen Elektronen und Positronen zusammenhängt. Diese Teilchen werden im elektrischen Feld, welches vom rotierenden Magnetfeld induziert wird, beschleunigt. Die Photonen werden dann zum Beispiel über den Synchrotron Mechanismus ausgestrahlt. Die beiden Modelle unterscheiden sich in der angenommenen Region, in der die Beschleunigung der geladenen Teilchen stattfindet. Das “polar cap” Modell nimmt an, dass die Beschleunigung oberhalb der magnetischen Polkappen nahe am Pulsar stattfindet. Im “outer gap” Modell dagegen findet die Beschleunigung weit aussen in der Magnetosphäre statt. Wegen unterschiedlicher Absorptionsmechanismen der hochenergetischen Gammastrahlung sagen die beiden Modelle unterschiedliche Energiespektren für die höchsten Energien voraus. Gammastrahlen mit einer Energie zwischen MeV und wenigen GeV werden mit Satellitenexperimenten gemessen. Die Photonen mit Energien oberhalb ein paar Dutzend GeV können mit Experimenten auf der Erde nachgewiesen werden; zum Beispiel mit dem MAGIC (Major Atmospheric Gamma-ray Imaging Cherenkov) Teleskop. Dieses Teleskop registriert die Cherenkov-Lichtblitze von elektromagnetischen atmosphärischen Luftschauern, die von einem Gammateilchen initiiert werden, welches Atome oder Moleküle der Atmosphäre trifft. Von allen existierenden Cherenkovteleskopen hat MAGIC die tiefste Energieschwelle von 55 GeV. Seit den Anfängen der Gammastrahlen-Astrophysik mit erdgebundenen Instrumenten vor 20 Jahren wurde erfolglos versucht, gepulste hochenergetische Gammastrahlung von Pulsaren zu messen. Ein Grund für diese Nichtdetektion ist der von den Pulsarmodellen vorausgesagte starke “cut-off” (Abbruch) des Energiespektrums und die zu hohe Energieschwelle der Cherenkovteleskope, welche es verunmöglichte, eine genügende Anzahl gepulster Gammateilchen zu messen. Nachdem in Daten, die 2005 genommen wurden, mit dem MAGIC Teleskop ein Hinweis für gepulste hochenergetische γ-Strahlung oberhalb 60 GeV vom Krebspulsar gefunden wurde, haben wir ein neues Trigger-System1 für das MAGIC Teleskop entwickelt (der sogenannte Summentrigger, oder sum trigger). Dieser neue Trig1

Ein Trigger-System ist ein Auslösemechanismus, welcher bewirkt, dass die gemessenen Daten gespeichert werden.

VII

ger reduzierte die Energieschwelle von 55 GeV auf 25 GeV. Zwischen Oktober 2007 und Februar 2008 beobachteten wir den Krebspulsar mit dem Summentrigger während gut 22 Stunden bei guten Wetterbedingungen. Die Analyse dieser Daten führte zur Entdeckung von gepulster hochenergetischer γ-Strahlung mit einer Signifikanz von 6.2σ. Die beiden obengenannten Pulsarmodelle sagen einen steilen cut-off des gepulsten Gammaspektrums mit exponentiellem (outer gap) respektive super-exponentiellem Verlauf (polar cap) voraus. Für beide Fälle habe ich die charakteristische Energie, bei welchem das Spektrum abbricht, berechnet: Ecut = (17.7 ± 2.8stat ± 5.0syst ) GeV (exponentiell), respektive Ecut = (23.2 ± 2.9stat ± 6.6syst ) GeV (super-exponentiell). Die hier vorgestellte Messung widerspricht den Voraussagen des polar cap Modelles und schliesst dieses daher aus. In einem weiteren Schritte habe ich das gepulste Energiespektrum bestimmt und herausgefunden, dass der Pulsar Gammastrahlen sogar oberhalb 60 GeV aussendet; d.h. höher als das outer gap Modell voraussagt. Das gepulste Spektrum oberhalb 25 GeV folgt eher einem Potenzgesetz (f (E) ∝ E −α ) mit spektralem Index α = 3.1 ± 0.6stat ± 0.5syst , als dem steilen exponentiellen Abbruch.

VIII

Contents Abstract

V

Zusammenfassung

VII

Introduction

12

1 Cosmic Rays

17

2 The γ–Ray Universe 2.1 Sources of γ-Rays . . . . . . . . . . . . . . . . . . 2.1.1 VHE γ-Ray Sources . . . . . . . . . . . . 2.1.1.1 Extragalactic Sources . . . . . . 2.1.1.2 Galactic Sources . . . . . . . . . 2.2 Production and Absorption Mechanism of γ-Rays

. . . . .

. . . . .

. . . . .

. . . . .

. . . . .

. . . . .

. . . . .

. . . . .

. . . . .

. . . . .

. . . . .

21 22 22 24 24 24

3 Pulsars 3.1 Introduction: The Crab Pulsar as a Canonical Pulsar . . . 3.2 Measurements of Pulsed γ-Rays from Pulsars before 2007 3.3 Emission of Electromagnetic Radiation from a Pulsar . . . 3.3.1 Polar Cap and Slot Gap Model . . . . . . . . . . . 3.3.2 Outer Gap Model . . . . . . . . . . . . . . . . . . . 3.4 Conclusion . . . . . . . . . . . . . . . . . . . . . . . . . . .

. . . . . .

. . . . . .

. . . . . .

. . . . . .

. . . . . .

. . . . . .

. . . . . .

. . . . . .

. . . . . .

. . . . . .

29 30 30 35 38 40 41

. . . . .

. . . . .

. . . . .

. . . . .

4 Observation Technique and Data Analysis 4.1 Direct Measurement of γ-Rays with Satellites . . . . . . . . . . . . . . . . . 4.2 Imaging Air Cherenkov Telescopes (IACTs) . . . . . . . . . . . . . . . . . . 4.3 The MAGIC Telescope . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4.3.1 Mirrors, Frame and Drive . . . . . . . . . . . . . . . . . . . . . . . . 4.3.2 Camera . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4.3.3 The MAGIC Standard Trigger . . . . . . . . . . . . . . . . . . . . . 4.3.4 The Data Acquisition System . . . . . . . . . . . . . . . . . . . . . . 4.3.5 The Calibration System . . . . . . . . . . . . . . . . . . . . . . . . . 4.3.6 The Central Pixel . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4.3.7 Pyroscope . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4.3.8 The Rubidium and the GPS Clock . . . . . . . . . . . . . . . . . . . 4.4 Data Analysis . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4.4.1 Extraction of the FADC Signals . . . . . . . . . . . . . . . . . . . . . 4.4.2 Image Cleaning . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4.4.2.1 Standard Image Cleaning . . . . . . . . . . . . . . . . . . . 4.4.2.2 Sum Image Cleaning . . . . . . . . . . . . . . . . . . . . . . 4.4.2.3 Calibration Image Cleaning . . . . . . . . . . . . . . . . . . 4.4.3 Image Parameters . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4.4.4 Event Classification with Random Forest and Estimation of the Energy

45 46 46 51 51 53 54 55 55 58 58 58 58 60 60 60 61 63 66 67 IX

4.4.5

4.5

5 The 5.1 5.2 5.3 5.4

5.5 5.6 5.7

5.8 5.9

Determination of the Differential Flux . . . . . . . 4.4.5.1 Unfolding Methods . . . . . . . . . . . . . Analysis of Periodic Signals from Pulsars . . . . . . . . . . 4.5.1 Determination of the Cut-Off Energy of the Pulsed

. . . . . . . . . . . . . . . . . . . . . . . . . . . γ-Ray Emission

. . . .

68 70 71 72

New Sum Trigger 75 Reaching the Design Goal of the MAGIC Telescope . . . . . . . . . . . . . . 76 The Sum Trigger: The Basic Idea . . . . . . . . . . . . . . . . . . . . . . . . 77 Photomultiplier Tubes Afterpulses . . . . . . . . . . . . . . . . . . . . . . . 79 Implementation of the Sum Trigger and Afterpulses into the Detector Simulation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 80 5.4.1 Afterpulse Simulations . . . . . . . . . . . . . . . . . . . . . . . . . . 80 5.4.2 Sum Trigger Simulations: Topology, Thresholds and Shape . . . . . 81 Electronics and Installation on Site . . . . . . . . . . . . . . . . . . . . . . . 83 The Calibration of the Sum Trigger . . . . . . . . . . . . . . . . . . . . . . . 86 Comparison with the MAGIC Standard Trigger . . . . . . . . . . . . . . . . 86 5.7.1 Homogeneity . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 86 5.7.2 Energy Threshold . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 86 5.7.3 Collection Area and Sensitivity for Low Energies . . . . . . . . . . . 88 Future Improvements . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 88 Conclusion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 91

6 Observation of the Crab Pulsar using the Analog Sum Trigger 6.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 6.2 Data Sample . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 6.3 Data Selection . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 6.4 Monte Carlo Simulations considering the Geomagnetic Field . . . . . . . . . 6.5 Data Analysis . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 6.6 The Optical Lightcurve . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 6.7 The Detection of Pulsed High Energy γ–Rays . . . . . . . . . . . . . . . . . 6.7.1 Number of Excess Events and Significance . . . . . . . . . . . . . . . 6.7.2 MC - Data Comparison . . . . . . . . . . . . . . . . . . . . . . . . . 6.7.3 Do P1 and P2 Follow a Different Spectrum? . . . . . . . . . . . . . . 6.8 Cut-Off Energy . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 6.8.1 Exponential and Super Exponential Cut-Off . . . . . . . . . . . . . . 6.8.2 Power Law Cut-Off . . . . . . . . . . . . . . . . . . . . . . . . . . . . 6.8.3 Single Pulse Cut-Off Energy . . . . . . . . . . . . . . . . . . . . . . . 6.9 The Pulsed Spectrum above 25 GeV . . . . . . . . . . . . . . . . . . . . . . 6.9.1 The Total Pulsed Spectrum in P1 and P2 . . . . . . . . . . . . . . . 6.9.2 Single Pulse Spectra . . . . . . . . . . . . . . . . . . . . . . . . . . . 6.9.3 Systematic Errors Influencing the Spectral Measurements . . . . . . 6.10 Phase Diagram . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 6.10.1 Arrival Time of Each Pulse . . . . . . . . . . . . . . . . . . . . . . . 6.10.2 Pulse Profile Morphology . . . . . . . . . . . . . . . . . . . . . . . . 6.10.3 P2/P1 Ratio as a Function of Energy . . . . . . . . . . . . . . . . . . 6.11 Measurement of High Energy γ-rays from the Crab Nebula at Low Energies

93 93 94 100 100 100 102 105 105 106 110 110 110 114 115 118 119 122 124 127 127 131 133 133

7 Summary, Conclusion and Outlook 137 7.1 Implications of our Measurements for Pulsar Models . . . . . . . . . . . . . 138 X

7.1.1

7.2

Conclusion derived from the Total and Phase Resolved Flux Measurements . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 7.1.2 Conclusion Derived from the Pulse Profiles at Different Energies . . 7.1.2.1 The Phase Separation between P1 and P2 . . . . . . . . . . 7.1.2.2 An Upper Limit for the Distance between the Production Regions of Optical and VHE γ-Ray Photons . . . . . . . . 7.1.2.3 A Lower Limit on the Mass Scale of Quantum Gravity . . . 7.1.2.4 The P2/P1 Ratio: Comparison with the Predictions of a Multicomponent Model . . . . . . . . . . . . . . . . . . . . Outlook . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

A Gamma/Hadron Separation Using Shower Autocorrelations A.1 General Idea . . . . . . . . . . . . . . . . . . . . . . . . . . . A.2 Conclusion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . A.3 Future Development: Implementation of the Algorithm using CUDA . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

138 143 143 143 144 145 146

149 . . . . . . . . 149 . . . . . . . . 151 a GPU and . . . . . . . . 152

B List of Acronyms and Abbreviations

153

List of Figures

154

Bibliography

158

Acknowlegdement

169

XI

Introduction Since human beings looked into the sky at night for the first time, they have been fascinated by the luminary’s illustrious equability. The majestic imperturbability of their movement in contrast to the unpredictable and chaotic life on Earth was interpreted by prehistoric peoples as divine signs. Consequently, they started to research the constellation’s periodic emergence and set the cornerstone of one of the eldest sciences: Astronomy. Observations of celestial objects were done at different time epochs in independent cultures, e.g. in ancient Egypt, where 2.000 BC the high correlation between the fruitful inundation of the Nile and the heliacal rising of Sirius was discovered, or in prehispanic Mexico, where the Maya predicted stellar constellations and calculated the duration of a tropical year to the very precise number of 365.2421 days. The methods to observe celestial objects and to predict their apparent movement were improved over the centuries. Until the beginning of the 20th century, all observations of celestial objects were performed by eye, using optical instruments, thus remaining constrained to visible wavelengths. The predictions of electromagnetic waves by Maxwell in the middle of the 19th century, their experimental affirmation by H. Hertz in the 1880s and the understanding of Hertz that light behaves like the electromagnetic waves he measured (refraction, polarization and reflection), led to the conclusion that light is an electromagnetic wave as well. In 1895, an important discovery was made by Willhelm Conrad Röntgen, when he discovered a new kind of radiation which he called X-rays. He found that a running, but closed cathode ray tube was able to make a barium platinocyanide screen shimmer placed more than a meter away from the tube. In order to extent astronomical observations to the whole electromagnetic spectrum, huge technical advances had to be made. The first radio signals from a celestial body were measured by the engineer Karl Guthe Jansky in 1930, studying short wave transatlantic transmission. He noticed static noise on the transmission, correlated to the siderial rotation of the Earth. He concluded that the source had to be a celestial object, not the Sun, and could localize the source at the center of the Milky Way in the Constellation of Sagittarius. Unfortunately, the Earth’s atmosphere is impervious to most electromagnetic wavelengths. Rayleigh scattering prevents e.g. the ultraviolet radiation to reach the surface of the Earth. To observe the radiation before being absorbed by the atmosphere, people have used balloon borne detectors and from the 1960’s on also experiments on satellites. The first cosmic X-ray source was detected by Riccardo Giacconi in 1962, with a X-ray detector on board the research rocket Aerobee. They failed to produce a X-ray image of the moon, but discovered a source of X-rays from the direction of Scorpius X-1. The photons with the highest energies and the shortest wavelengths are called γ-rays. They were discovered in 1900 by P. Villard who studied radioactive elements. He realized that one component of the radioactive radiation is not deflected by magnetic fields in contrast to the charged α and β radiation and called them γ-rays. It was later found that γ-rays belong to the electromagnetic spectrum. In the interactions of γ-rays with matter, however, particle reactions like the photo effect dominates over coherent wave behavior. At the beginning of the 20th century, there was a mystery about an ubiquitous ionizing

14

radiation which caused electrometers to discharge without an apparent reason. It was first thought Earth itself emits this “Earth Radiation”, believed to be radiation from radioactive elements. If the Earth were the source, however, the radiation would become weaker with increasing height. Victor Hess tested this hypothesis in 1911, using electroscopes on a balloon. He found that the radiation became stronger the higher he flew. He concluded that the ionizing radiation has an origin in outer space. Robert Millikan confirmed this measurements and dubbed the radiation cosmic rays. He believed the cosmic rays to consist of γ-rays. Arthur Compton on the other hand was of the opinion that the cosmic rays are mostly charged particles like protons. Bruno Rhossi found in 1934 that the cosmic ray rate depends on the Latitude. He concluded that cosmic rays are affected by the Earth’s magnetic field and must thus be charged. It took until 1948 when Melvin Gottlieb and James Van Allen found with balloon-borne emulsions chambers that the cosmic rays mostly consist of protons. Today it is known, that cosmic γ-rays contribute only to a small fraction of the cosmic rays of around 10−4 − 10−5 . In the 1960s, the first cosmic γ-ray source was discovered by a satellite borne gamma detector. The first ground based detection of cosmic γ-rays using the air shower technique (explained later) succeeded as late as 1989 when the Whipple collaboration discovered γrays from the direction of the Crab Nebula [1]. In figure 0.1, images of the Crab Nebula at Radio wavelengths to γ-rays are displayed. In this thesis, the Crab Pulsar, the highly magnetized and fast spinning neutron star in the center of the Crab Nebula, is established as an emitter of electromagnetic radiation, in particular at VHE γ-rays2 . The Crab Pulsar spins with a high frequency of ≈ 33 Hz. Its strong co-rotating magnetic field of ≈ 8 · 1012 G induces an electric field which accelerates particles. It is believed that electromagnetic radiation is emitted by relativistic electrons and positrons. Measurements of the pulsed γ-ray emission up to an energy of ≈ 10 GeV with satellite borne experiments revealed that the energy spectrum follows a power law f (E) ∝ E −α , with the spectral index α = 2.022 ± 0.014 [2]. Until 2007, no instrument was sensitive enough to measure γ-rays between 10 GeV (the former flux sensitivity limit of space borne experiments) and ≈ 50 GeV (the lower energy threshold of ground based experiments like MAGIC). Above 50 GeV, no pulsed VHE γ-rays were detected before 2007. This thesis is organized in the following way: Chapter 1: The measured cosmic rays on Earth are described, including their energy spectrum, possible relation to cosmic γ-rays and their believed origin. Chapter 2: An overview of HE and VHE γ-ray sources is given and possible emission mechanisms of VHE γ-rays are discussed. Chapter 3: Pulsars are considered as emitters of electromagnetic radiation, in particular the HE and VHE γ-ray emission from young pulsars like Crab. The measurements of pulsed γ-rays from the Crab Pulsar before 2007 is discussed and the two current main theoretical models, the polar cap and the outer gap model is reviewed. Chapter 4: In this chapter, I will explain the observation techniques used to measure HE and VHE γ-rays. The Imaging Air Cherenkov Telescope (IACT) MAGIC is illustrated and 2

Note that “Crab” consists of two potential VHE γ-ray sources: the Crab Pulsar in the center and the Crab Nebula surrounding. They are separated in space by about 7 light years.

15

(a) Radio wavelengths

(d) X-rays

(b) Infrared

(c) Optical

(e) γ-rays

Figure 0.1: The Crab Nebula as seen in radio, infrared and optical wavelengths, in Xand γ-rays. Note that the Nebula in X-rays is about 40% as large as the optical Crab Nebula, which is in turn about 80% as large as the Nebula at Radio wavelengths. The angular resolution of the γ-ray image (taken with the MAGIC telescope) is limited to some 0.1◦ , while the optical image of the Crab Nebula has about a diameter of 0.01◦ . Thus, the image in γ-rays cannot resolve any structure of the Crab Nebula. Image credit: NASA/CXC/SAO (X-ray), Palomar Obs (Optical), 2MASS/UMass/IPAC- Caltech/NASA/NSF (Infrared), NRAO/AUI/NSF: (Radio), MAGIC collaboration (γ-rays) [3, 4].

16

the analysis of data taken with MAGIC is explained, where special emphasis will be put on the analysis of data close to the trigger threshold and to the analysis of data from sources, that are expected to show pulsation. Chapter 5: The standard trigger of the MAGIC telescope provides an energy threshold of ≈ 55 GeV. To measure pulsed γ-rays from the Crab Pulsar, we realized that a lower threshold is mandatory. In 2006/07, we therefore designed a new trigger system, the sum trigger, which allowed to lower the energy threshold to 25 GeV. In this chapter, the new trigger is described. Chapter 6: The new trigger system was used for the observations of the Crab Pulsar in the winter 2007/08. In this chapter, the analysis of this data is described, leading to the discovery of pulsed emission from the Crab Pulsar above an energy threshold of 25 GeV. The spectrum of the pulsed emission and the pulse profile morphology are studied. Chapter 7: The results obtained in chapter 6 are reviewed and compared with theoretical predictions. An outlook on future observations of the Crab Pulsar is given.

1 Cosmic Rays After their discovery by V. Hess in 1912 and their identification as mostly protons1 in the 1940s, cosmic rays were studied in great detail. The charged cosmic ray spectrum was measured over many orders of magnitudes in energy, from a few GeV energy to up to 1020 eV, a macroscopic energy of about 20 joules, corresponding to the kinetic energy of a tennis ball flying with 100 km/h. The cosmic ray spectrum measured by different experiments is shown in figure 1.1. For low energies < 1013 eV, charged cosmic rays can be directly measured with satellite- and balloon-borne experiments. For particle energies above 1013 eV, the small flux, however, makes the direct detection with the small space and air borne experiments difficult. In 1938, Pierre Auger found that the cosmic radiation arrives in short correlated bunches and concluded that a cosmic ray particle with a high energy (he estimated the energy to be above 1015 eV) hitting molecules of the Earth’s atmosphere initiates an extended air shower, where the energy of the incident cosmic ray particle is distributed over millions of secondary particles. To measure charged cosmic ray particles with energies above 1013 eV, ground based experiments are used which sample the secondary particles created in the shower. Based on these measurements, the energy, the incoming direction and the kind of the incident particle are inferred (figure 1.2). The small flux of cosmic rays at the highest energies of only 1 particle/km2 in a few hundred years requires experiments covering huge areas. The currently largest experiment is the AUGER array in Argentina covering an area of 3000 km2 and observing charged cosmic rays for energies above ≈ 1017 eV. Like the photons in optical astronomy, cosmic rays are messengers from distant objects. During their propagation through space, however, charged cosmic rays are randomly deflected by galactic and extragalactic magnetic fields. Therefore, their incoming direction on Earth is randomly distributed and cannot be related with a known source observed e.g. at optical wavelengths. An exception might be the cosmic rays at the highest energies, for which the larmor radius2 extends the radius of our galaxy. In 2007, the AUGER collaboration found a first hint for a correlation of these events with nearby active galactic nuclei (AGN) [6]. Nevertheless, the origin of the cosmic rays is still not completely solved to this date. It is believed that the location where γ-rays are emitted and cosmic hadrons are accelerated are related, where in particular supernova remnants have been considered as accelerators of cosmic ray particles (e.g. [7, 8] and references therein). Thus, the observation of γ-rays can reveal details about the acceleration mechanisms and locations of the charged cosmic rays. A mechanism explaining the high energies of the cosmic rays and their power law spectrum was proposed in 1949 by Enrico Fermi who postulated that charged particles which cross back and forth through moving and magnetized shock fronts would statistically gain energy. The higher the energy of a particle, the shorter the time while it is confined within 1

The composition of cosmic rays at energies of a few GeVs is: 98% nuclei, where 87% are protons, 12% helium nuclei and 1% heavier nuclei, < 2% electrons and < 10−4 γ-rays. p 2 with p the particles momentum, q its charge and B the magnetic field. also called gyroradius: rL = |q|B

18

1. Cosmic Rays

Figure 1.1: The cosmic ray spectrum, measured by different experiments. Note that the flux at energies of ≈ 1019 eV is below 1 particle per km2 per year. Below 3 · 1015 eV (the so-called “knee”), the spectrum follows a steep power law dN/dE ∝ E −α , with α ≈ 2.7, above it steepens to α ≈ 3. Image credit from [5].

1. Cosmic Rays

19

Figure 1.2: Illustration of an extended air shower initiated by a cosmic ray particle. The secondary particles created in the air shower are sampled with water filled tanks measuring the Cherenkov radiation of the relativistic particles passing through the water. The roundish detector in the middle registers the fluorescent light emitted by atmospheric molecules that were excited by the ionizing radiation within the air shower. (Image credit: AUGER collaboration) the front and the higher its escape probability. This naturally leads to a power law distribution of the particle’s energies which leave the acceleration region [9]. The γ-rays are then emitted by the charged particles, either by electrons or positrons, by means of the synchrotron (in magnetic fields) and self compton mechanism, or by protons which may hit an atom from the denser matter surrounding the acceleration region and, in turn, produce neutral pions among others, which will then decay into γ-rays. More details are given in chapter 2. Due to their strong electric field induced by the spinning magnetic field, pulsars have also been considered as electrostatic particle accelerators contributing to the hadronic cosmic hadrons (e.g. [10]), to cosmic electrons (e.g. [11]) or to the cosmic positrons (e.g. [12]). This will be covered in chapter 3 in more detail. There are still many open questions about the details of the cosmic accelerators, the types of cosmic ray sources, the propagation through the interstellar medium or the cosmic ray composition at the highest energies. It is hoped that current cosmic ray experiments like AUGER, together with γ-ray experiments like MAGIC, H.E.S.S., VERITAS or FGST (see chapter 2), will help to shed light on the mystery of cosmic rays.

2 The γ–Ray Universe

Figure 2.1: A selection of scientific targets of very high energy γ-ray astronomy. Gamma-rays are produced by the most energetic processes that occur in the universe, e.g. in jets emitted by AGNs, in supernova remnants or from the vicinity of pulsars with magnetic fields that are 100 Million times stronger than the ones we can produce in laboratories on Earth. Scientific topics of γ-ray astrophysics include the understanding of the production mechanisms of γ-rays in different sources or the correlation of the emission of γ-rays with the acceleration of cosmic rays, as well as fundamental physics questions, including the nature of Dark Matter, where γ-rays are potential annihilation products detectable on Earth or a possible vacuum dispersion due to quantum gravitational (QG) effects, causing time of flight delays of photons with high energies. A selection of scientific targets is displayed in figure 2.1. In this chapter, potential γ-ray sources and production mechanisms of high energy γ-rays

22

2. The γ–Ray Universe

Figure 2.2: The locations of the sources emitting γ-rays above 100 MeV detected with EGRET [13].

will be discussed. In chapter 3, pulsars as γ-ray sources will be discussed in detail.

2.1 Sources of γ-Rays In the 20 years after the first detection of cosmic γ-rays with the Vela satellite in the 1960s, the satellite COS-B discovered 25 point-like sources in total. Due to the rather poor angular resolution, many sources could not be identified with known objects observed in radio and optical wavelengths or in X-rays. A large boost of the number of HE γ-ray sources was achieved with the Compton γ-ray satellite (CGRS), which started in 1991. The EGRET (Energetic Gamma Ray Experiment Telescope) experiment on board of this satellite performed a full scan of the sky. In total, γ-rays were detected from more than 200 sources. Amongst the galactic sources, EGRET detected HE γ-rays from pulsars, from a solar flare or from the large Magellanic cloud. From the extragalactic ones, HE γ-rays from AGNs and from one radio galaxy (Cen A) were discovered. Additionally, γ-rays were detected from several unidentified sources. The so called third EGRET catalog is a list of all those sources with their positions, fluxes and possible counterparts in other wavelengths (see figure 2.2) [13].

2.1.1 VHE γ-Ray Sources In the VHE range, the first discovered source was the Crab Nebula, detected with the Whipple telescope in 1989 [1]. Until 2008, in total 25 extragalactic and 54 galactic sources were discovered by Imaging Air Cherenkov Telescopes (IACTs), see figure 2.3. In figure 2.4, the number of X-ray and γ-ray sources are displayed as a function of time. Among the individual VHE γ-ray sources are:

2. The γ–Ray Universe

VHE γ -ray Sky Map

23

+90

o

(Eγ>100 GeV)

VHE γ -ray sources Blazar (HBL) Blazar (LBL) Flat Spectrum Radio Quasar Radio Galaxy Pulsar Wind Nebula Supernova Remnant Binary System Open Cluster Unidentified

+180

o

o

-180

-90o 2008-10-29 - Up-to-date plot available at http://www.mppmu.mpg.de/~rwagner/sources/

Figure 2.3: The VHE γ-ray sources for energies above 100 GeV. From [14]

Figure 2.4: The known number of X-ray, HE and VHE γ-ray sources as a function of time (before 2007). [15]

24

2. The γ–Ray Universe

2.1.1.1 Extragalactic Sources AGN: Galaxies, which produce more (electromagnetic) emission than deduced from their stellar content, stellar remnants and interstellar medium are called Active Galactic Nuclei (AGN). It is assumed that they contain in their center a supermassive black hole with around 106 to 1010 M⊙ .1 If this black hole accretes matter from a surrounding accretion disk, a jet can form perpendicular to this disk with an extension of up to 1024 cm. Within these jets, more radiation can be produced than in the whole galaxy. It is assumed that in jets charged particles can be accelerated up to the highest known energies. VHE γ-ray emission is measured from the direction of several AGNs. It is still unclear how VHE γ-rays are produced. There are leptonic or hadronic emission models (e.g. [16]), explained in section 2.2. Gamma Ray Bursts: These are very bright and violent short term phenomena. The nature of these bursts is still not fully understood. In a few years, the γ-ray satellite Swift [17], dedicated to the measurement of gamma ray bursts, has detected several 100s of these most energetic objects. MAGIC has observed several Gamma Ray Bursts so far, but no VHE γ-ray emission has been detected. 2.1.1.2 Galactic Sources Supernova Remnants: A supernova may occur when a star at the end of its lifetime runs out of the fuel needed to provide the fusion. Stars that are massive enough will collapse and eject their outer shell. Depending on the mass of the progenitor star, a black hole, a neutron star or a white dwarf is formed in the center. The ejected material builds up a spherically shaped nebula, called supernova remnant. A shock is formed when the ejected material flows into the interstellar medium, where Fermi acceleration can take place. The large amount of relativistic particles may give rise to VHE γ-rays. Pulsars: If the neutron star from a supernova explosion spins with a periodicity below a few seconds, it is called a pulsar. Due to their strong co-rotating magnetic field of > 1012 G, an electric field is induced, accelerating charged particles which will emit γ-rays. See chapter 3. Pulsar Wind Nebulae: Most of the rotational energy lost by the pulsar is emitted in a relativistic particle wind or in a Poynting (radiation) flux. At some distance from the pulsar, a stationary shock is formed from where the VHE γ-ray emission originates. γ-Ray Binary Systems: These objects consist of a compact object (i.e. a neutron star or a black hole) and a massive companion star. The compact objects accrete matter from the companion. In the binary systems called microquasars, a relativistic jet can be formed where particles are accelerated similar to AGNs [18, 19].

2.2 Production and Absorption Mechanism of γ-Rays Photons with energies above 100 keV (wavelengths shorter than 10 · 10−12 m = 10 pm) are usually called γ-rays. Their interactions are dominated by particle like behavior, e.g. 1

M⊙ = 1.9891 · 1030 kg, the mass of the sun.

2. The γ–Ray Universe

25

inelastic Compton scattering, the Photoelectric Effect or Pair Production2 and surmount wave behavior like coherent scattering by a large factor. Up to an energy of several MeV, γ-rays can be emitted by decaying radioactive nucleii. This emission comes from the relaxation of excited nuclei3 and from the annihilation of electron-positron pairs, where the positrons are produced in β + decays. The energy features of this kind of emission is a sharp line at the relaxation energy or at 511 keV, respectively. In a few sources (e.g. Cas A: 1.157 MeV [20], Galactic Center: 511 keV [21]), line features were discovered. However, for most γ-ray sources the energy spectrum is continuous. For VHE γ-rays the spectrum often follows a power law f (E) ∝ E −α with the spectral index α. Thermal emission from a hot body can be excluded as a main production mechanism of HE and VHE γ-rays. The hottest measured structures in the universe are accretion discs around compact objects. Those objects have their thermal peak emission in X-rays and emit photons up to tens of keV. Furthermore, the black body spectrum does not follow a power law as it is observed from VHE γ-ray sources. In general, there are two concurring models trying to explain the HE and the VHE γradiation: Hadronic and leptonic γ-ray production mechanisms. In the following sections, these processes will be briefly discussed. Emphasis will be laid on the particle acceleration and VHE γ-ray emission mechanisms important in the vicinity of a pulsar. Hadronic γ-Ray Emission The charged cosmic rays are dominated by hadronic particles (mostly protons). When these high energy particles hit cosmic matter, (neutral) pions can be produced. The main decay mode of π 0 is into two γ-particles: π 0 → γγ (98.798%)

(2.1)

Assuming that a cosmic particle accelerator is surrounded by dense matter, a large amount of π 0 are produced by the highly energetic protons. These pions will then decay into γ-rays. An overview over VHE γ-ray emission processes initiated by cosmic hadrons can be found in [22]. Leptonic γ-Ray Interactions Relativistic electrons and positrons in strong magnetic fields efficiently lose energy through synchrotron and curvature emission and through inverse Compton scattering that transmits the energy from high energy electrons to lower energy γ-rays, thus producing HE and VHE γ-rays. These mechanism are displayed in the following sections. Pulsars, as described in chapter 3, provide a strong electrostatic field needed to accelerate charged particles and a strong magnetic field needed for the observed synchrotron and curvature radiation and are thus considered as candidates of HE and VHE γ-ray emission. In figure 2.5, the leptonic interaction processes important for the production and absorption of HE and VHE γ-rays in the vicinity of a pulsar are shown. • Synchrotron Emission: Electrons and positrons in magnetic fields efficiently lose energy through synchrotron radiation. For small magnetic fields and electron energies Ee B ≪ 1, with the electron’s energy Ee , its mass me , and natural constant me c2 Bcr Bcr =

m2e c3 e¯ h

≈ 4.4 · 1013 G, the synchrotron emission can be treated ignoring quantum

2

above 1022 keV

3

For example: 60 Co → 60 Ni∗ + e− + ν¯e . γ(1.332 MeV) + γ(1.173 MeV)

β−

60

Ni∗ relaxates under emission of two γ-rays:

60

γ

Ni∗ →

60

Ni +

26

2. The γ–Ray Universe

Synchrotron Radiation

Curvature Radiation

γ

γ e

B

Magnetic Pair Production

γ

B

e

Photon Splitting

e+

γ

strong magn. field

γ2 strong magn. field

e-

Inverse Compton Scattering

Pair Production γLE

γHE

γ1

e+

γ

e-

eHE

LE

γHE

eLE

Figure 2.5: Different interaction mechanism between γ-rays and leptons important in a pulsar’s magnetosphere. The index “LE” and “HE” indicates low energy and high energy, respectively. Note that in case of photon splitting and magnetic pair production, the recoil is absorbed by the strong (quantized) magnetic field.

2. The γ–Ray Universe

27

effects. In this regime, a particle with energy E and mass m emits synchrotron radiation with the power [22]: P =

2Ke2 Γ4 v 4 , 3c3 r 2

(2.2)

where K is the number of unit charges of the radiating particle, e the unit charge, Γ = mEeec2 , v ≈ c the electron’s velocity and r the radius of its trajectory. Depending on the magnetic field B⊥ perpendicular to the flight path of the electron, the peak energy of the emitted photons lies at Epeak (eV) = 5 · 10−9 B⊥ /1GΓ2 . In case of pulsars, magnetic fields can reach 1012 G. Thus, the classical treatment of synchrotron radiation is only possible for electron energies below 0.01 GeV. Above, the synchrotron emission has to be treated within quantum mechanics. This includes photon splitting [23, 24, 25]: Ee B me c2 Bcr

≥1

γ −−−−−−−→ γγ

(2.3)

and magnetic pair creation [24, 26, 27, 28]: Ee B m c2 Bcr

≥1

γ −−e−−−−−→ e+ e− .

(2.4)

In strong magnetic fields, electrons will radiate γ-rays through synchrotron (and curvature) radiation. Those γ-rays will create e+ e− pairs, either by magnetic pair production or by γγ → e+ e− . By repeating this procedure, an electromagnetic cascade forms [22]. • Curvature Radiation: In the strong and bent magnetic fields in the vicinity of pulsars (e.g. [29]), an energy loss similar to synchrotron emission occurs to charged relativistic particles. The particles moving along the field lines rapidly lose their perpendicular energy due to synchrotron radiation. The “spiraling” trajectory around the field line is treated within quantum mechanics - similar as a 2D harmonic oscillator. Due to synchrotron radiation the particle ends up in the ground state. Being in the ground state, the particle does not emit further synchrotron radiation. The trajectory of these particles can thus be treated quasi 1D and emission of γ-rays will only occur due to curvature radiation [30]. The energy loss of a particle emitting curvature radiation is similar to the loss by synchrotron radiation: P =

2Ke2 Γ4 c , 2 3rC

(2.5)

with rC the radius of the bent magnetic field line [31]. In common pulsar models (outer gap model and polar cap model), this is the most important VHE γ-ray emission mechanism [32, 33]. • Inverse Compton Scattering: The upscattering of low energy photons with high energy electrons (or positrons) is called Inverse Compton Scattering: γLE + eHE → γHE + eLE ,

(2.6)

an important process to transfer energy from particles to photons. Depending on the product of γ ray and electron energy, one distinguishes the following cases. In case

28

2. The γ–Ray Universe

of small energies ( Ee Eγ ≪ m2e c4 ), i.e. within the Thomson limit, the process can be described classically and the cross section becomes:   Eγ σγe = 8/3πre2 1 − . (2.7) me c2 For large energies (i.e. Ee Eγ ≫ m2e c4 ), i.e. the Klein-Nishina limit, quantum mechanics effects have to be taken into account. The cross section then decreases with rising γ-ray energy:   2 2Eγ 1 2 me c + σγe = πre ln . (2.8) Eγ me c2 2 Inverse Compton Scattering of photons emitted through synchrotron radiation by the same electrons (the so-called synchrotron self Compton process), is one of the mechanisms used to explain the VHE γ-rays from several source types. In pulsars, photons produced through synchrotron and curvature radiation are Compton scattered to TeV energies too. Due to efficient absorption mechanism in the vicinity of the pulsar, however, it is believed that these TeV γ-rays will not leave the pulsar’s magnetosphere. • Absorption: In the strong magnetic fields in the vicinity of a pulsar, magnetic pair production and photon splitting play an important role in the absorption of γ-rays. Another important absorption mechanism is pair creation: γγ → e+ e−

(2.9)

This process has the highest probability to happen if the mean photon energy lies at 511 keV. A detailed study of all emission and absorption mechanisms in the vicinity of a pulsar can be found in [34].

3 Pulsars Since their discovery at Radio wavelengths in 1967 by Jocelyne Bell [35] (see figure 3.1), pulsars have been observed in nearly the whole electromagnetic spectrum, from Radio waves to γ-rays. In 1934 already, Zwicky and Baade theoretically deduced that after a supernova explosion, the core of a star may remain as a dense neutron star. The discovery of pulsed radio emission from the Crab Pulsar in the center of the Crab Nebula in 1969 [36], which is the remnant of the known historic supernova in 1054, confirmed this theory. In the years after the discovery of pulsars, the observed periodic radio signal was explained due to the rotation [37] or the oscillation [38] of the neutron star. Because of emission of gravitational waves by the oscillation of the pulsar, however, the vibrational energy would be lost quickly [39]. Thus, the oscillation theory was dropped. Electromagnetic radiation is emitted in a (co-rotating) beam. The periodic crossing of the line of sight of the observer with this beam explains the observed pulsation. Today it is known that a pulsar is a fast rotating neutron star being the compact remnant of a star in a certain mass class1 after its supernova explosion at the end of its lifetime (e.g. [40, 41]). The size of the star is dramatically reduced during the collapse (from a radius of ≈ 1011 cm to ≈ 106 cm), reducing its moment of inertia by ten orders of magnitude, resulting inR small rotational periods of down to fractions of a second. The star’s magnetic flux Φ = B · n da is conserved during the collapse. To compute the pulsar’s surface magnetic field, we assume that its progenitor had a surface (dipole) magnetic field strength of B ≈ 102 G. During the collapse, the star’s surface area a is reduced by a factor ≈ 1010 . Thus, the surface magnetic field strength B is increased by this factor to B ≈ 1012 G.2 Even though pulsars have been known for more than 40 years, the precise emission mechanism are still scarcely known. While the radio emission is believed to originate from above the magnetic polar caps emitted by coherent highly energetic electrons, there exist several theories about the location within the pulsar’s magnetosphere where the emission of the remaining electromagnetic emission occurs, displayed in figure 3.2. All theories assume that the emission is connected to electrons and positrons which are accelerated in the strong induced electric field. An overview in more details will be given later in this chapter. Before 2007, only seven pulsars (Crab Pulsar, PSR B1509-58, Vela, PSR B1951+20, Geminga, PSR B1055-52, PSR B1706-44) were known to radiate HE γ-rays above 1 GeV, and no pulsed VHE γ-radiation was detected. In the following sections, I will first expound the general phenomenology of young pulsars, where Crab is considered in particular. Then, an overview of the γ-ray measurements from pulsars before 2007 is given. In section 3.3, the mechanisms for the electromagnetic emission from pulsars is discussed. In the last sections, I will review the current main models explaining the γ-ray emission from pulsars: The polar cap (which was extended to the slot gap model) and the outer gap model. 1

The core mass of the progenitor star must be between 1.44 M⊙ and 3 M⊙ in order to form a neutron star. 2 The energy density of such a magnetic field is ≈ 4 · 1020 J mm−3 . So in 10 mm3 , the energy mankind uses in one year is stored.

30

3. Pulsars

Figure 3.1: The discovery plot of the first pulsar in 1967 by J. Bell [35]. The small peaks indicate the periodicity of the signal. J. Bell realized that the signal was correlated with the siderial day (i.e. a day relative to the fix stars), indicating an extra-solar origin of the pulsating signal.

3.1 Introduction: The Crab Pulsar as a Canonical Pulsar The Crab Pulsar is the most famous representative of a canonical pulsar, characterized by a mass of M ≈ 1.4 M⊙ , a diameter around 20 km (see figure 3.3) and a surface magnetic field of B ≈ 1012 G. Canonical pulsars have rotation periods P between 0.01 - 8 seconds. −18 - 10−12 s s−1 , shown in figure 3.4. In The rotation period derivative P˙ = dP dt is within 10 this diagram, pulsars start in the upper left corner, having a small rotational period and a 3 large P˙ . They lose rotational energy by means pof magnetic dipole radiation and thus slow down. While the surface magnetic field B ∝ P P˙ (see formula 3.5 below) stays constant, they “wander” during their lifetime diagonally down towards the lower right corner where they end up in the “graveyard”, where no more pulsed electromagnetic emission takes place. It is not known yet why the emission stops, but it is assumed that either the induced electric field is too weak to pull charges out of the pulsar’s surface or the acceleration of electrons to high enough energies to initiate pair creation (E > 511 keV) is prevented [42, 43, 44]. The Crab Pulsar is one of the youngest known pulsars and appears thus on the upper left corner.

3.2 Measurements of Pulsed γ-Rays from Pulsars before 2007 Today, more than 2000 radio pulsars are known. The first detection of pulsed HE γ-rays was in the 1990s, when the EGRET experiment onboard CGRS discovered pulsating HE γ-rays from 7 pulsars [45]. In figure 3.5, the pulse profile (also dubbed light curve, phaseogram, phase4 distribution) of those pulsars is shown. Four objects show distinct features: Two rather broad pulses are visible, with a phase separation ∆φ between 0.35 < ∆φ < 0.45. The first pulse is called main pulse (P1), the second one inter pulse (P2). The region in between the two pulses is called bridge. 3 4

Other possible radiation mechanisms are gravitational waves or magnetic radiation of higher order. The phase φ corresponds to the rotational angle of the pulsar, divided by 2π. Within one rotation of the pulsar, φ goes from 0 to 1.

31

open field lines closed field lines α

Light Cylinder

polar cap +

+

+

+

- -

nu su ll c r fa h a ce rge

Rotation axis

ma sym gne me tic fi try eld ax is

3. Pulsars

+

- -

ou

ter

Pulsar

ga

p

slot gap -

+

+

-

+

-

+

-

-

-

+

Figure 3.2: The polar cap, slot gap and outer gap regions are possible candidates for locations within the pulsar’s magnetosphere where the acceleration of particles and emission of electromagnetic radiation might take place. In case of the Crab Pulsar, another distinct feature can be observed: The pulses P1 and P2 appear at the same phase for radio waves, optical emission, X-rays and γ-rays. For energies up to 10 GeV, EGRET measured the pulsed γ-ray spectra for the 8 phase regions5 : Component Leading Wing Peak 1 Trailing Wing Bridge Leading Wing Peak 2 Trailing Wing Offpulse

1 1 2 2

Name LW1 P1 TW1 Bridge LW2 P2 TW2 OP

Phase Interval φ 0.88 − 0.94 0.94 − 0.04 0.04 − 0.14 0.14 − 0.25 0.25 − 0.32 0.32 − 0.43 0.43 − 0.52 0.52 − 0.88

The measured lightcurve for energies above 100 MeV and the phase regions are displayed in figure 3.6. While the measured pulse profiles contain the information about the location of the emission region and the geometry of the emission beam, the energy spectrum also allows 5

Measurements of the spectrum in phase regions are called “phase resolved” spectrum measurements.

32

3. Pulsars

Figure 3.3: A canonical pulsar like the Crab Pulsar has about the size of Zurich, but the mass of 1.4M⊙ and spins 33 times per second around its axis. (Map: Google maps).

to draw conclusions on the mechanisms within the pulsar’s magnetosphere that lead to the emission of the measured photons. The photons in different phase regions are believed to be produced in different locations, or the emission directions of the photons differ. Thus, the measurement of the spectrum in phase resolved bins (i.e. the spectrum of the photons whose phase lies within a certain range) yields more information on the emission and absorption mechanisms within the pulsar’s magnetosphere. In figure 3.7, the phase resolved spectrum measurements of several γ-ray instruments are displayed. It was found that above an energy of ≈ 100 MeV, the spectrum follows a power law: dN = f0 · f (E) = dE dA dt



E Enorm

−α

(3.1)

with the flux normalization f0 , a fixed normalization energy Enorm and the power law’s spectral index α. For γ-ray energies between 100 MeV and 5 GeV, the spectral index of the γ-ray spectrum in each phase interval has a value of α = (2.022 ± 0.014) [2]. Furthermore, the P1 spectrum at the highest energies shows a hint for a spectral turnover. Despite various attempts in the last 20 years, observations with ground based telescopes did not allow to detect pulsed γ-rays from any pulsar before 2007 [46, 47, 48, 49, 50, 51]. Studying the flux sensitivities of these instruments reveals that they would have detected pulsed emission from Crab if the spectrum followed the extrapolated EGRET power law up to the highest energies. Therefore, the spectrum must show a turn over at energies between 10 GeV and < 100 GeV. To account for this spectral change, the power law measured with

3. Pulsars

33

sars

ul gp

n you

canonical pulsars

Figure 3.4: The pulsar period P and its derivative P˙ for the known radio pulsars. The dashed lines denote lines of constant characteristic age τC = 2PP˙ , see formula 3.7. The hatched area in the upper left corner denotes pulsars with a characteristic age below 10-100 kyrs, dubbed “young pulsars”. The dot dashed lines are lines p ˙ with constant surface magnetic field B ∝ P P (depicted in formula 3.5). The reddish area denotes the canonical pulsars, described in the text. During their lifetime, canonical pulsars decelerate. Given that their surface magnetic field is conserved, they follow the direction of the green arrow until they end up in the “graveyard”. For pulsars within this area, no more electromagnetic emission is observed. Image adapted from [40].

34

3. Pulsars

Figure 3.5: The pulse profiles of the 7 pulsars detected with EGRET [45], at radio and optical wavelengths and for X- and γ-rays. On the x-axis of each pulse profile the rotational phase is given. The Crab, Vela, B1951+32 and Geminga pulsars show two distinct peaks in γ-rays with a phase separation of 0.4. In case of Crab, the two peaks are at the same phase for all energies.

350

LW2

TW1 Bridge

TW2

LW1

P2

OP

TW1

P1

300

LW2 Bridge

TW2

LW1

OP

P2

P1

counts [a.u.]

250 200 150 100 50 -1

-0.8

-0.6

-0.4

-0.2

-0 Phase φ

0.2

0.4

0.6

0.8

1

Figure 3.6: The pulse profile of the Crab Pulsar measured with the EGRET experiments for energies above 100 MeV with the phase regions defined in [2]. For better visibility, two identical cycles are shown.

3. Pulsars

35

Figure 3.7: The phase resolved spectrum measurements of the Crab Pulsar [2] for γ-ray energies between 0.1 keV and ≈ 6 GeV. For energies above 100 MeV, the spectra in all narrow phase intervals follow a power law with spectral index α = 2.022, measured up to an energy of ≈ 6 GeV. This power law component is indicated by the dashed line. EGRET is adjusted with a cut-off: dN f (E) = = f0 · dE dA dt



E Enorm

−α

“

·e

−E Ecut

”β

.

(3.2)

Ecut denotes the cut-off energy and β the superexponential cut-off index. This steep turn over led to enhanced efforts to lower the energy threshold of the MAGIC telescope. This is described in the chapters 5 and 6.

3.3 Emission of Electromagnetic Radiation from a Pulsar In a first approximation, the rotating neutron star – the pulsar – consists of a conducting rotating sphere with a co-rotating dipole magnetic field, inclined to the axis of rotation by the angle θ. The magnetosphere surrounding the pulsar is filled with a conducting plasma. A rotating conducting body in a magnetic field is called unipolar inductor, because an

36

3. Pulsars

electric field is induced between the rotation axis and a point on its surface. This was first measured by Faraday in 1831. He found that a spinning disk in a magnetic field induces an electric potential between a contact at the center of the disk and at its perimeter. This effect is described by Faraday’s induction law and included in Maxwell’s equations. An inclined rotating magnetic field is a Hertz Dipole and emits in vacuum the power [52]: B 2 r 6 Ω4 (3.3) Lvacuum (θ) = 0 3 sin2 θ 6c with the magnetic field B0 on the pulsar’s surface and the pulsar’s radius r. The pulsar spins with frequency Ω. The spin down luminosity LSD is a function of the observed periodicity P and its derivative P˙ = dP dt : LSD =

dErot d = dt dt



1 2 IΩ 2



˙ = 4π 2 I = IΩΩ

P˙ , P3

(3.4)

where I is the moment of inertia. We assume the pulsar’s mass to be mPulsar = 1.4M⊙ , its diameter d = 20 km and thus its moment of inertia I = 4.5 · 1038 m2 kg. For the Crab Pulsar with P = 0.033 s and P˙ = 4.23 · 10−13 , a spin down luminosity of LSD = 4 · 1031 W is expected6 . Only a small fraction of the total spin down luminosity goes in the pulsed radiation. Most of it is carried away in an ultrarelativistic particles and Poynting flux outflow. Combining equation 3.3 and 3.4, the surface magnetic field is found as a function of P and P˙ : p (3.5) B = 3.2 · 1019 P P˙ G The characteristic age of a pulsar is then computed using the following expression: τ=

Zτ 0

Zτ Zτ ZP 1 1 P P˙ P dP dt = P P˙ dt = dt = P P˙ P P˙ P P˙ 0

0

(3.6)

P0

with P dP = P P˙ dt and assuming that P P˙ (∝ B 2 ) be constant over time. Integrating above equation and assuming that P02 ≪ P 2 , it follows: τC =

P , 2P˙

(3.7)

For Crab, the surface magnetic field is then7 : B = 1.8 · 1012 G and its characteristic age: τC = 1240 years. As Crab is the remnant of the supernova in 1054, its real age is around 950 years. τC is an upper limit onto the lifetime of a pulsar, as the initial rotation frequency P0 is not known. It should be remarked that the derivation of equation 3.5 only makes sense if the spin down is dominated by magnetic dipole radiation. Nevertheless, there are other mechanisms for the pulsar to spin down: 1. Gravitational wave radiation. 2. Quadrupole magnetic radiation. 6

The whole mankind needs around 1021 J per year, so 1 second spin down luminosity could deliver power for about 20 billion years. 7 More sophisticated methods (e.g. taking into account that the magnetosphere is filled with a plasma and not a vacuum) lead to a surface magnetic field of around B ≈ 8. · 1012 G [53]

3. Pulsars

37

3. Accretion of surrounding material, e.g. in binary systems. These mechanisms predict different braking indeces n defined as: ˙ = −KΩn Ω

(3.8)

˙ its derivative. K is an arbitrary constant, depending with Ω the rotation frequency and Ω on the surface magnetic field. By differentiating equation 3.8 once and replacing Ωn−1 by ˙ Ωn /Ω = −Ω/(KΩ), one finds: ¨ ΩΩ (3.9) ˙ Ω In case of a spin down due to magnetic dipole radiation, n = 3 is expected. In case of the Crab Pulsar, a breaking index of nCrab = 2.51 ± 0.01 is measured, with most young pulsars have a breaking index < 3 [54]. The consequences of this is that either n=

• the spin down is not dominated by magnetic dipole radiation and there are additional effects causing a torque onto the pulsar denoted above [54], • or the factor K in equation 3.8 is not constant due to a changing surface magnetic field B0 with time or a change of the inclination angle α over time. The precise way how the spin down luminosity is translated to the observed emitted radiation is not well known. The coherent radio emission is believed to originate from above the polar caps. Its production mechanism is thought to be generated by coherent relativistic electrons, similar to a MASER [55]. For the electromagnetic emission of photons with wavelengths between infrared (IR) and X-rays, the emission is not coherent and the emission mechanism is hardly known [56, 57]. All theories about the emission of electromagnetic radiation from a pulsar are based on charged particles being accelerated in the strong electric field induced by the rotating magnetic field. The electric field E and the electrostatic potential Φ induced by the rotating magnetic field in a pure vacuum is found by solving the Poisson equation ∆Φ = 0 for the electric potential Φ within the pulsar’s light cylinder while assuming a vacuum outside and the magnetic field being a dipole co-rotating with the pulsar (e.g. [58]):
Bs R 5 Ω 2 3 cos θ − 1 (3.10) Φ= 6cr 3 with the neutron star’s radius R and its surface magnetic field Bs . For the Crab Pulsar with Bs ≈ 1012 G, R = d/2 ≈ 10 km and P = 0.033 ms, a potential above the polar cap of Φ ≈ 100 · 1012 V is induced. For particles with unit charge, this corresponds to field energies of 100 TeV ≫ me c2 , injecting electrons from the neutron star’s surface. Electrons are accelerated by the electric field and start the particle-radiation cascade resulting in the magnetosphere being filled with a dense plasma. Plasma currents will screen the electric field parallel8 to the magnetic field and thus maintain the co-rotation condition E · B = 0. The γ radiation from pulsars is then produced by charges accelerated in the strong electric field induced by the rotating magnetic field. The main emission mechanisms are: Synchroton emission, curvature emission and inverse compton scattering. The main difference between existing pulsar models consists of the location where the acceleration of the 8

Acceleration of electrons and positrons is only possible parallel to the magnetic field lines. If the particles are accelerated perpendicular to the magnetic field lines, they lose quickly their energy due to synchrotron radiation.

38

3. Pulsars

charged particles takes place. As mentioned above, the electric field parallel to the magnetic field is screened in regions with a plasma density above the Goldreich-Julian (GJ) limit ρGJ (r) [42]: Ω·B 1 1 (3.11) ∇·E=− ρGJ (r) = 4π 2πc 1 − |Ω×r| c

The screening does not allow particles to be accelerated. Nevertheless, if the plasma density is below the Goldreich-Julian limit, the screening is not effective enough and there is a component of E parallel to B. These regions with low plasma density – called gaps – are believed to exist either above the polar caps, extending to a few pulsar radii, or at large distances from the pulsar, close to the light cylinder at the radius: rLC =

c . Ω

(3.12)

The region called outer gap extends between the null charge surface9 and the last open field line out to the light cylinder. An extension of the polar cap to the outer magnetosphere along the null charge surface – the slot gap model – was introduced later to overcome some of the shortcomings of the polar cap model (e.g. the small pulse width of the predicted pulse profile, which contradicts the observations). Figure 3.2 shows schematically these regions and the null charge surface in the co-rotating frame. Equation 3.13 (from [59]), which assumes a mostly model-independent relation derived from simulations of γ-ray absorption as a result of magnetic pair production in a magnetic dipole field, gives a relation between the pair creation cutoff energy Emax and the height of the emission above the pulsar’s surface r/R0 (with R0 the pulsar’s radius) for a pulsar with surface magnetic field B0 and period P : (  3 ) r r r 0.1Bcrit GeV (3.13) max 1, Emax ≈ 0.4 P R0 B0 R0 with the critical magnetic field Bcrit = 4.4 · 1013 G. If one measures pulsating γ-rays with an energy above Emax , one can thus derive a lower limit for the height r above the pulsar’s surface where the emission takes place. In chapter 7 this relation will be used to constrain the emission height. In the following sections, the main pulsar models will be presented.

3.3.1 Polar Cap and Slot Gap Model The polar cap and slot gap model are related, thus both will be discussed in the same section: Polar Cap Model The region above the polar cap was first proposed in [43, 60, 61, 62] as a location of particle acceleration. A gap is formed, where the density is below ρGJ . Depending on the surface temperature of the neutron star, there are different mechanisms explaining the formation of the gap. Electrons are released from the surface of the neutron star if the surface temperature exceeds the electron thermal emission temperature and are accelerated by the electric field induced by the rotating inclined magnetic dipole. As derived above, in case of the Crab 9

The null charge surface is the area for which is Ω · B = 0. It separates the positive and negative co-rotating charges.

3. Pulsars

39

Pulsar, potential drops of up to 100·1012 V are induced in the polar cap region, accelerating electrons to Γ-factors up to Γ = E/mc2 ≈ 108 . In these models (called space-charge limited flow models (SCLF)) (e.g. [63]), the charge density at the surface of the pulsar equals the Goldreich-Julian one. The charge density ρ(r) decreases as r −3 . Thus, the accelerating electric field Ek is zero at the surface, but increases with larger radii and decreasing plasma density ρ(r) < ρGJ (r) for r larger than the neutron star’s radius [63]. If the surface temperature is below the thermal emission temperature, a vacuum gap will be also formed. Thus, the parallel component of the electric field Ek is not short out in the plasma. Accelerated charged particles will start an electromagnetic pair cascade. In these so-called vacuum gap models as well as in the SCLF models, the accelerating electric field will be screened because of the charged particles produced in the onset of the electromagnetic pair cascade [63]. Absorption mechanisms like photon splitting and magnetic pair creation inhibit γ-rays with high energies to leave the polar cap area. This efficient absorption leads to a steep cut-off of the emitted γ-ray energy of super exponential shape. For the Crab Pulsar, the model predicted cut-off lies between 1-10 GeV [32]. In general, polar cap models predict rather narrow peaks in the γ-ray pulse profile [64], due to the small emission angles.

Slot Gap Model: An extension of the polar cap model model, the slot gap model, predicts wider emission beams, needed to match the wide pulse peaks observed in γ-rays. The polar cap model was extended to higher altitudes along the last closed field lines [65, 66]. Their potential as high energy γ-ray emission site, however, has only been explored recently [67]. The most recent slot gap models [32, 64, 68] account for caustic emission, which takes into account relativistic effects like aberration and consider time of flight and phase delays, leading to a staggering of γ-rays from different altitudes within the pulsar’s magnetosphere and thus forming the characteristic pulse profile. Also the initial frame dragging (Lense Thirring effect) is considered, leading to an increase of the size of the polar cap region and thus to wider peaks in the predicted pulse profiles. For both models, the polar cap and the slot gap model, the observed emission is believed to be produced at one single pole. The main pulse P1 and the inter pulse P2 are observed when the co-rotating hollow cone beam enters the observer’s field of view and when it leaves again. The emission in between in the so called bridge comes from the insight of the cone. Due to the rotation, the absorption efficiency of the magnetic field is different for the main and the inter pulse [69], stronger for the leading peak forming the main pulse and weaker for the trailing peak forming the inter pulse. Therefore, the cut-off of the main pulse is expected to be at lower energies than the one of the inter pulse. Hints for this behavior was observed with EGRET in case of the Crab Pulsar [2]. In figure 3.8, recent flux predictions from K. Hirotani [53] and A. Harding [32] for the slot gap model are compared. The flux predicted in [53] has a flux peak value which is a factor 33 lower than the one predicted in [32] and is the same factor below the observed EGRET flux measurements [2]. The main differences between the two predictions lies in the derivation of the electric field Ek which accelerates the electrons. Harding finds a stronger electric field component parallel to the magnetic field than Hirotani. The thus derived Γe factor10 of the primary electrons is about factor 33 larger and the derived spectrum fits the measured data. 10

Γe =

E me c2

40

3. Pulsars

Figure 3.8: Slot gap model total pulsed spectrum predictions from Harding [32] (left) and Hirotani [53] (right). The dashed and solid lines in the later denote the flux predictions for different viewing angles ζ.

3.3.2 Outer Gap Model In the outer gap model, the acceleration of the charged particles takes place far out in the pulsar’s magnetosphere, between the null charge surface and the last open field line, out to the light cylinder, displayed in figure 3.2. It was first proposed in [70, 71] as a region within the pulsar’s magnetosphere where acceleration takes place. This gap is formed because of particles escaping through the light cylinder are not replaced fast enough by particles from below and the local plasma density falls below the GJ density ρGJ . For the Crab Pulsar, the electric field parallel to the magnetic field can reach around Ek ≈ 2.55 · 108 V m−1 in this region, accelerating charged particles [53]. The γ-rays in the outer gap are either produced by the upscattering of softer photons (infrared/X-rays) or by curvature radiation. As in case of the polar cap model, the γ-rays above an energy threshold of some GeV will be absorbed, but not due to magnetic pair creation, because the magnetic field strength is not sufficient enough. In the outer gap, the absorption is caused by pair production with either infrared or X-ray photons: γHE + γLE → e+ e− .

(3.14)

The produced e+ e− -pairs are accelerated in the field again, and will upscatter softer photons. The absorption results in a γ-ray spectral cut-off of exponential shape, with the cut-off energy at a few GeV or tens of GeV in case of the Crab Pulsar. The large emission angle of the outer gap has as a natural consequence the broad γ-ray pulses of the pulse profile. The photons in the main and the inter pulse originate from the same magnetic hemisphere. The main peak is produced when the hollow cone enters the line of sight, the inter pulse when it leaves the line of sight. Furthermore, the bridge emission comes from the inside of the hollow cone. Since the emission in P1 and P2 can originate from different regions within the outer gap, different spectra might be observed for P1 and P2. A tentative outer gap model spectrum of a young pulsar is shown in figure 3.9 [33]. It is composed of several components. At the lowest energies, the flux is dominated by the synchroton radiation of the electrons, only surmounted in a small energy region by the thermal emission from the pulsar’s surface. The main part of the γ-ray emission is due to curvature radiation. Another component is visible in the VHE γ-ray region: The inverse compton scattered γ-rays. It is debated whether these γ-rays escape through the light cylinder or if they are completely absorbed [53, 72]. From the experimental point of view

3. Pulsars

41

Figure 3.9: The theoretical phase-averaged spectrum for a typical young γ-ray pulsar deduced for the outer gap model. Solid lines show the curvature spectrum (CR), the synchrotron emission (Sy) and the thermal radiation from the surface (kT). The dashed curve gives the TeV pulsed spectrum from inverse Compton Scattering (CS) of the synchroton photons on the primary e± . Figure and caption adapted from [33]. it should be remarked that no TeV observation has revealed a pulsating signal from pulsars in the last 20 years, indicating a strong absorption of the inverse Compton component. In case of the Crab Pulsar, flux predictions agree acceptably with measurements done with EGRET, even in phase resolved spectra. See figure 3.10.

3.4 Conclusion Pulsars are laboratories where relativistic magnetohydrodynamics and huge magnetic fields can be investigated. Pulsating electromagnetic emission has been measured from a few thousands of pulsars in the last 40 years. Seven pulsars are even known to emit γ-radiation above 1 GeV. The mechanism leading to this emission is poorly known. Current theories assume the radiation being related to relativistic electrons, accelerated by the induced electric field in so-called gaps within the pulsar’s magnetosphere. The two main model, the polar cap (extended to the slot gap) and the outer gap model assume different regions where the particle acceleration takes place. The measurement of gamma rays at the highest energies would allow to probe these pulsar models. Depending on the emission region, different fluxes, phase diagram morphologies and spectral cut-offs are predicted. Due to their particle nature, gamma rays are ideal messengers of the mechanisms within the magnetosphere leading to their emission: They do not coherently interfere and thus any reaction on their way from the pulsar to Earth would lead to their absorption. If one measures pulsating γ-rays from Crab, it is thus certain that they were produced within the pulsar’s light cylinder.

42

3. Pulsars

Figure 3.10: Phase-averaged (top left) and phase-resolved (others) spectrum of the pulsed emission from the Crab pulsar. The thin solid, dashed, and dotted curves denote the un-absorbed photon fluxes of primary, secondary, and tertiary components, respectively, while the thick red solid one presents the spectrum to be observed, including absorption along the line of sight. Interstellar absorption is not considered. Data points are from [2] http://www.sron.nl/divisions/hea/kuiper/data.html. In the top left panel, upper limits by ground-based observations are also plotted above 50 GeV. Figure and caption adapted from [53].

3. Pulsars

43

The observation of pulsed γ-rays above a certain energy, however, turned out to be more difficult than expected. Satellite experiments were too small to measure the low expected fluxes and ground based experiments like MAGIC had their energy threshold on a too high level. In chapters 5 and 6, I will describe the observations of the Crab Pulsar with a new trigger that lowered the energy threshold of the MAGIC telescope by a factor of 2, resulting in the discovery of pulsed γ-rays above 25 GeV.

4 Observation Technique and Data Analysis

Figure 4.1: The MAGIC telescope during observations. The first successful observation of cosmic γ-rays dates back to the 1960s when the Vela satellite1 measured γ-rays from gamma ray bursts. Direct measurements of γ-rays have to be performed from space because the Earth’s atmosphere will efficiently absorb γ-rays before they reach the ground. If the γ-ray energy is high enough, however, the incident particle will start an air shower in the atmosphere, with secondary particles and radiation even reaching the surface of the Earth. By measuring these secondary particles, γ-rays can be detected indirectly, using the atmosphere as a giant calorimeter. The following chapter is dedicated to the observational techniques used in the measurements of cosmic γ-rays. The first section includes the description of so-called direct measurements which are performed on board of satellite experiments like EGRET [73], AGILE (Astrorivelatore Gamma ad Immagini ultra LEggero) [74] or FGST (Fermi Gamma-ray Space Telescope, renamed from “Glast”) [75]. Then I will discuss the ground based indirect detection of γ-rays by measuring the Cherenkov light produced in the electromagnetic air shower induced by a γ-ray, in particular the imaging air Cherenkov technique, where the 1

This military satellite was launched to monitor the “Treaty Banning Nuclear Weapon Tests In The Atmosphere, In Outer Space And Under Water”.

46

4. Observation Technique and Data Analysis

Cherenkov light is projected by a mirror on the focal plane, where the photons are registered by sensitive photodetectors (section 4.2). In the following section 4.3, the MAGIC telescope (shown during observation in figure 4.1) on the Canary Island La Palma is described, its set-up, its data acquisition, standard trigger2 and calibration system. Furthermore, I will describe the data analysis software in section 4.4. Special emphasis will be placed on a new image algorithms I developed and implemented as a part of this thesis, which allows to analyse shower images with the lowest yet measured energies down to 25 GeV. The new algorithm is more robust against fluctuations of the shower images induced by fluctuations of the night sky background (NSB) light. In the last section 4.5, the methods and algorithms used in the analysis of pulsed signals will be described.

4.1 Direct Measurement of γ-Rays with Satellites Experiments measuring the primary γ-rays outside the Earth’s atmosphere generally use techniques known from particle physics. The experiments in space have to measure the energy of the γ-ray, its direction and have to distinguish the photon from other particles3 . One technique used e.g. in the LAT detector [76] onboard of FGST is pair creation, illustrated in figure 4.2. The incoming γ-ray creates a pair in the thick pair creation foil. The traces of the electron and the positron are recorded by the particle tracking detectors, thus allowing a reconstruction of the incoming γ-ray direction. A calorimeter at the bottom of the detector then measures the energy of the e+ e− pair and thus the γ-ray energy. An anticoincidence shielding prevents triggers from charged particles and thus suppresses the background denoted above. Their energy range starts from a ≈ MeV with a higher threshold at a few tens of GeV, overlapping with the lowest energy threshold of the MAGIC telescope (see chapter 5). In case of FGST4 , the angular resolution is < 3.5◦ at 100 MeV and < 0.15◦ for energies above 10 GeV.

4.2 Imaging Air Cherenkov Telescopes (IACTs) Cosmic ray particles like protons or γ-rays above a certain threshold of a few GeV hitting the Earth’s atmosphere start an air shower, similar to the cascade formed in a calorimeter in particle physics experiments. An electromagnetic shower (initiated by a VHE γ-ray or by a relativistic electron) shows a characteristic development: an exponential increase of the number of shower particles at the beginning, and an even faster decrease after reaching its maximum. At the beginning, secondary particles are created by pair production and bremsstrahlung, which will create further particles. If the energy of those particles falls below a critical limit, where ionization losses dominate over the creation of new particles, this multiplication process dies out and the shower development stops. Energy losses for charged particles include ionization, bremsstrahlung and to a small amount Cherenkov radiation. Photons on the other hand loose energy by the photoelectric effect, Compton scattering and pair creation at energies above 1 MeV. The particle-radiation cascade forms until the shower particles reach the critical energies. For charged particles in air, this energy is at ≈ 83 MeV/Z, with Z the charge of the particle, when the energy loss by 2

The new sum trigger, which was developed as a part of this thesis, will be described in the separate chapter 5. 3 For example charged cosmic rays, which are 104 − 105 times more abundant than γ-rays. 4 FGST is able to measure γ-rays up to 300 GeV. Even though the detector is much more sensitive than EGRET, its detection area is around 1.3 m2 . All known cosmic VHE γ-ray sources show a decreasing flux above 100 GeV. Thus, satellite experiments are flux-limited to measure γ-rays in the VHE range.

4. Observation Technique and Data Analysis

47

γ

Anticoincidence shield

conversion foil

particle tracking detector

e+

e-

calorimeter

Figure 4.2: The e+ e− pair creation technique used onboard the space borne FGST experiment to detect γ-rays. The incoming γ-ray creates in the conversion foil a pair of e+ e− , which are tracked by the particle tracking detector. In the calorimeter, their energy is measured. With this setup, one can reconstruct the γ-ray energy and its incoming direction.

48

4. Observation Technique and Data Analysis

ionization dominates the bremsstrahlung and the exponential growth dies out. In case of photons, the critical energy occurs at energies of ≈ 5 MeV when the Compton radiation starts to dominate over pair creation. In case a charged cosmic ray particle hits molecules or atoms in the atmosphere, a shower will develop as well. Most of the particles produced in this cascade are pions. Neutral pions will most probably decay into two γ-rays and thus initiate the above described electromagnetic cascade. Further shower particles are kaons and antiprotons as well as muons, which can reach the ground without decaying before. A hadronic air shower therefore consists of a shower core, built up from nucleons and mesons, of electromagnetic subshowers induced by neutral pion decays and of non-interacting muons (and neutrinos). A distinct feature of these hadronic shower is the large lateral distribution of the shower particles (due to the significant transverse momentum transfer in hadronic interactions) compared to the more compact electromagnetic showers. In figure 4.3, electromagnetic and hadronic air showers are schematically displayed. In figure 4.4, the traces of the simulated shower particles for γ-rays and protons at different energies are shown. Gamma ray initiated showers show a smaller lateral distribution than hadronic showers. At the low energies simulated in these figures, hardly a shower particle reaches the ground level. Thus, the measurement of shower particles will only work for γ-rays with energies above a few hundreds of GeV. The particles in the shower, however, move with velocities faster than the speed of light in the atmosphere vparticle > cn = nc , with the refraction index n and thus emit Cherenkov radiation [77]. The light emission results from the re-orientation of the instantaneous electric dipoles induced by the particle in the medium. Because the particle passes "superluminally" through the medium, the electromagnetic impulses originated in the different trajectory elements are in phase, and a net field can be produced at distant locations. The geometric picture for this follows from Huygens’ principle (see figure 4.5). A shock-wave is created behind the particle. The wavefront of the radiation propagates at a fixed angle θ, with 1 . (4.1) cos θ = βn v

the relative velocity of the particle. Due to the change of the density of with β = particle c air with altitude the refraction index decreases with height h as: n(h) = 1 + n0 e−h/h0

(4.2)

with h0 ≈ 7.1 km and n0 ≈ 0.00029. According to equation 4.1, the Cherenkov angle changes as well. The higher the radiating particle, the smaller the angle. The Cherenkov light builds up a wave front in front of the particle with a thickness of only a few nanoseconds (10−9 s). The first measurement of Cherenkov radiation from an air shower succeeded in 1953 by using a mirror in a ton reflecting the Cherenkov light onto a photomultiplier tube (PMT) [78]. Today, in ground based VHE γ-ray astronomy, three main techniques dominate: • Sampling of the shower particles, e.g. in the TIBET [81] or MILAGRO [82] experiments. To do so, water-filled pools or scintillators are used. If a relativistic shower particle hits the water, it irradiates Cherenkov photons which are measured by photomultipliers on the ground of the pools. These experiments have large observation solid angles of nearly 2π and can be run day and night resulting in a duty cycle of

4. Observation Technique and Data Analysis

49

Primary γ Cosmic Ray (p, He, Fe,...)

γ Atmospheric Nucleus

e−

π

o

γ

π+

EM Shower

e−

π+

e+

e+

γ

γ

+ −

e− e+

γ

π

o

e−

e−

νµ

Nucleons, K, etc. Atmospheric Nucleus

π−

e+

+

e

γ e−

EM Shower

e+ γ

e+

γ

e−

e−

µ−

γ

νµ

µ+

e+ EM Shower

e−

e+

EM Shower

(a) Gamma-ray initiated shower

(b) charged cosmic ray initiated shower

Figure 4.3: Schematics of an electromagnetic (a) and a hadronic air shower (b). See text for details. Images provided by S. Commichau [79] almost 100%. Nevertheless, their energy threshold is at several hundreds of GeV and the discrimination of γ-rays against the hadronic background is not very efficient. • Measuring the photons in the Air Cherenkov light cone, as it was done for instance in the Celeste [49] or STACEE [83] experiments. These kind of experiments used modified solar towers at night to efficiently measure the Cherenkov photons sampled over a large area. They had a very low energy threshold of around 50 GeV. Nevertheless, the discrimination of γ-ray showers against the hadronic background was rather poor. • Imaging Technique, e.g. Whipple, MAGIC, H.E.S.S., VERITAS, CANGAROO. The Cherenkov photons from the air shower are focussed by a mirror onto a camera consisting of light detectors, for instance PMTs or Geiger Mode Avalanche Photodiodes (G-APD). Thus, one achieves an image of the extended air shower. Particles hitting the atmosphere parallel to the telescope’s axis will point towards the center of the camera. Due to the changing Cherenkov angle θ with height, Cherenkov photons produced at higher altitudes will be imaged closer to the center than those produced at lower altitudes. As mentioned above, the duration of the air Cherenkov flash is 2 − 3 ns. Using fast PMTs and trigger electronics therefore allows to measure the Cherenkov photons over the stable background light coming from stars and other NSB light sources. A schematics of this technique is displayed in figure 4.6. This technique allows better γ - background separation than the two above mentioned techniques, because with the imaging air Cherenkov method, one gains information about the development of the whole shower. The image of an air shower initiated by γ-rays differs in shape, extension and other parameters from the ones initiated by hadrons. The MAGIC telescope being the IACT with the lowest energy threshold, has a trigger threshold of ≈ 55 GeV, using its standard trigger system. This threshold is comparable with the ones from the non-imaging air Cherenkov telescopes.

50

4. Observation Technique and Data Analysis

(a) 50 GeV proton

(b) 50 GeV γ-ray

25 km

100 GeV gamma-ray

0 -5 km

(c) 200 GeV proton

0

+5 km

(d) 100 GeV γ-ray

Figure 4.4: The traces of the particles in air showers triggered by γ-rays and protons. The electromagnetic shower initiated by γ-rays is more compact than the hadronic shower and its lateral distribution is smaller. The red traces are from leptons and photons, the violet traces from hadrons and the green ones from muons. [80]

4. Observation Technique and Data Analysis

51

c n ∆t

wavefront

θ electron

β c ∆t

Figure 4.5: Propagation of Cherenkov light, derived from Huygens’ principle. The electron moving horizontally polarises the crossed medium. If the electron moves with a velocity larger than the local speed of light, a coherent wave is emitted by the relaxating dipoles, called Cherenkov light. In the following sections, the IACT MAGIC will be described in detail.

4.3 The MAGIC Telescope The MAGIC telescope with its 17 m tessellated mirror is currently the largest single dish IACT in the world [3]. It was designed in 1998 to measure γ-rays in the previously unexplored energy region between 30 GeV and 300 GeV [84]. Engineering work started in 2001, commissioning in 2003 and regular data taking in fall 2004. A major upgrade of the readout system was performed in 2006/2007 when the sampling frequency was increased from 300 MHz to 2 GHz, see section 4.3.4. In autumn 2007, a new additional trigger system was installed, reducing the trigger threshold from ≈ 55 GeV to ≈ 25 GeV. This trigger will be described in chapter 5. The telescope is situated on the Roque de los Muchachos on the Canary Island La Palma at 2245 m above sea level (28◦ 45′ N, 17◦ 53′ W), which is considered as one of the best sites for night time observations. Because of the high altitude, the telescope is closer to the shower maximum and thus less Cherenkov photons are absorbed by the atmosphere. These specifications – large mirror, highly efficient photodetectors, fast electronics and high altitude – lead to the lowest energy threshold of all current IACTs.

4.3.1 Mirrors, Frame and Drive An important design goal of the MAGIC telescope is its light-weight structure needed for fast repositioning in case of Gamma Ray Burst (GRB) alerts, triggered by space experiments like Swift [17]. The low weight was achieved by building the frame from carbon fiber epoxy composite tubes. The resulting total weight is as low as 60 tons, while the frame weighs even only 25 tons. The telescope stands on horizontal rails and uses an altitude-azimuth (alt-az) mount. The three servo motors of the drive system [85] provide a repositioning of 180◦ in azimuth of less than 25 seconds [84], with a tracking accuracy of < 0.022◦ . The maximal repositioning time between any two alt-az positions is < 40 seconds. To account for the bending of the light weight structure, a pointing correction model is used and regularly calibrated. To measure a possible mispointing, a so-called star guider camera is used, consisting of a CCD camera in the center of the mirror, measuring the

52

4. Observation Technique and Data Analysis

Primary particle (1 TeV) Top of atmosphere First interaction with nuclei of atmosphere at about 20 km height

Cherenkov light emission under characteristic angleθ

o θ = 0.29

~8 km a.s.l.

~5

o

θ = 0.86

km

Camera (cleaned event)

~ 20 m o

θ = 1.03 Camera

Cherenkov Telescope

Reflector

Figure 4.6: Schematics of the imaging technique. An image of the extended air shower is recorded by measuring the emitted Cherenkov light with an IACT. Image credits: S. Commichau.

4. Observation Technique and Data Analysis

53

direct light from a star and its reflection onto the camera plane, relative to 7 reference LEDs which are also located in the camera plane. This setup can also be used later in the analysis to correct a mispointing by software. The mirror of the telescope has a parabolic shape and consists of 964 49.5 cm × 49.5 cm mirror elements, with a total area of 236 m2 . Dew and ice formation is prevented by equipping each element with an internal heating system. The averaged reflectivity in the wavelength range of Cherenkov light (300 nm – 650 nm) is about 85%. The parabolic shape (instead of a spherical shape used in smaller IACTs) of the mirror prevents the timing information of photons arriving isochronously at the mirror plane. Thus, timing information of the Cherenkov light wave front is preserved, allowing for a short trigger coincidence window, which is important to increase the sensitivity to low energy γ-ray showers over the NSB, see chapter 5. Furthermore, the time development of the Cherenkov light can later be used in the analysis to distinguish between γ-ray and hadronic showers. To prevent the telescope to defocus due to bending, a so-called Active Mirror Control (AMC) [86] is used. The focussing is measured in terms of the Point Spread Function (PSF) which describes the defocussing of a point-like light source (e.g. a star or a distant flashlight). Using the AMC, a PSF of < 0.1◦ is achieved which is smaller than the size of a single PMT. To ensure a fast refocussing at different azimuth and altitude angles, look-up tables (LUT) are used, which are regularly calibrated using stars at different positions. Furthermore, the AMC allows different focal lengths. Because the shower maximum of γ-ray showers usually are at a height of 10 km above sea level, the telescope is in normal operation mode focussed to this distance.

4.3.2 Camera As light detectors, the PMTs ET9116 are used. These devices have a quantum efficiency of nearly 25%, further enhanced by 20% (relative) by coating the surface with a wavelength shifting lacquer. The camera is located at the focus point of the telescope. It has a diameter of 1.5 m and consists of 577 pixels. The inner camera is built out of 397 smaller 25 mm diameter PMTs, the outer camera out of 180 larger 39 mm PMTs, respectively5 . With the mirror’s focal length of ≈ 17 m, the camera covers a field of view (FOV) of 3.5◦ , allowing the observation of modestly extended objects. The inner pixels have a FOV of 0.1◦ , the outer ones 0.2◦ [87]. The finer pixelization of the inner camera was chosen because the images of low energy γ-rays are compact and close to the camera center (80% of the Cherenkov photons of a 25-35 GeV γ-ray shower lie on a ring between 0.25 − 0.8◦ , see figure 5.3 in chapter 5). The shower of higher energy γ-ray showers, however, develop further down into the atmosphere. Thus, the Cherenkov photons from the showers tail are imaged at a larger distance from the camera center. To contain the full shower in the camera, the large diameter is needed. Among the inner pixels, the central one is not used for the measurements of Cherenkov light, but to record the optical signal of the source being observed. For the observation of pulsars, this has the advantage to co-measure their optical pulsation. The roundish PMTs are arranged in a hexagonal structure, see figure 4.7. The PMT tubes only cover 50% of the camera plane. To capture the light hitting in between pixels, so-called Winston cones (non-imaging light collectors) are installed. Next to their nearly perfect light collection efficiency for photons reflected by the mirror, they also prevent NSB photons with an incident angle > 35◦ to reach the PMTs. The signal leaving the PMTs are AC6 coupled preamplified. In the transmitter board 5 6

In this work, “PMT” and “pixel” is used as synonyms. Capacitive coupling of the PMT signals, filtering the DC component out.

54

4. Observation Technique and Data Analysis

69cm 114cm Figure 4.7: The layout of the MAGIC camera. The smaller PMT of the inner camera are used to get a better sampling of lower energy γ-ray shower images, which rather lie in the inner camera. within the camera, the signals are translated by means of a Vertical Cavity Surface Emitting Laser (VCSEL) into optical signal, then transmitted by optical fibers to the counting house. Using optical instead of analog transmission shields the signal against noise and prevents its dispersion. In the counting house, the signals are back-transformed into electrical signals by a photodiode and then split in two parts. One part is sent to the trigger system (section 4.3.3 and chapter 5), the other to the FADC (Flash Analog Digital Converter) system (section 4.3.4).

4.3.3 The MAGIC Standard Trigger Because it is not possible to register the complete data stream from the 577 pixels sampled at 2 GHz, a Cherenkov telescope needs a trigger. A trigger pre-selects images of extended air showers. Because the Cherenkov photons of γ-ray showers arrive in a short time window of 2-3 ns, an ideal trigger uses a coincidence window of this order of magnitude. The discrimination level should be higher than the possible fluctuations of the NSB to prevent accidental triggers. Furthermore, certain preselections on the shape of a shower can be applied to distinguish γ-ray showers from hadronic ones. Due to their greater abundance of 104 − 105 , the great bulk of recorded events come from hadronic shower, however.

4. Observation Technique and Data Analysis

55

The MAGIC standard trigger [88] consists of three Levels: • Level 0 (LT0): Uses the signals from the innermost 325 pixels (displayed in figure 4.8). If the signal in a PMT reaches the discriminator threshold, a logical signal is produced. • Level 1 (LT1): In this step, a next neighbor (NN) logic is applied to the digital signals from LT0. This is done because images of γ-ray showers are rather compact. In the standard operation mode, 4 NN pixels above the LT0 discriminator threshold are required to fire LT1. Possible settings can be between 2NN and 5NN configurations. • Level 2 (LT2): Based on Field Programmable Gate Arrays (FPGA), the level 2 trigger allows complicated pattern recognition algorithms. During standard data taking, LT2 is not being used, i.e. it consists of a simple OR between the different macrocells displayed in figure 4.8. In this thesis, the term Standard Trigger refers to LT0 and LT1 with a 4NN configuration7 . The discriminator thresholds are set at 6 photo electrons (phe) and the time window to 6 ns. The calibration of the trigger system is performed several times per year. In this procedure, the calibration pulses (see section 4.3.5) are used to calibrate the discrimination level in phe and to adjust the relative arrival times of neighboring pixels. Depending on the NSB, especially during observation in moonlight, twilight or in a FOV with bright stars, the LT0 discrimination thresholds are automatically adjusted by the Individual Pixel Rate Control (IPRC). Typical individual pixel rates (LT0 rates) are between 100 MHz and 500 MHz. The total standard trigger (LT1) rate in moonless clear nights is between 200-300 Hz.

4.3.4 The Data Acquisition System If the trigger fires, it will start the Data Acquisition Queue (DAQ). Since January 2007, the signals are digitized with a Multiplexed (MUX) FADC system, increasing the sampling frequency from 300 MHz to 2 GHz. Multiplexing is used to reduce the number of readout channels and thus provides a cost effective FADC system. It is possible because even at the highest possible trigger frequencies of 2 kHz, the duty cycle of the digitizer is still very low8 [90]. As described later in section 4.4, this allows to exploit the timing resolution provided by the PMTs of < 1 ns and increases the sensitivity of the telescope by about a factor of 2. In figure 4.9, the DAQ is schematically displayed.

4.3.5 The Calibration System The number of Cherenkov photons emitted by an air shower initiated by a γ-ray is proportional to the γ-ray energy. Thus, it is important to know the conversion from the number of Cherenkov photons hitting the mirror to the ADC counts recorded by the data acquisition. In the MAGIC calibration system [91], three different methods are implemented to determine this conversion. All three make use of three very fast (3-4 ns) and high light output (109 photons per pulse) LEDs. For calibrating in different wavelengths, colored LEDs are used: 460 nm (blue), 370 nm (UV), 520 nm (green). They are situated in the 7

In [89], 3NN was used as well. As we realized then, noise triggers coming from PMT after pulses (see section 5.3) will dominate the low energy events. Furthermore, with a 3NN configuration, roundish shower images are preselected, making γ-hadron separation less efficient. 8 One event is 25 ns long, corresponding to 50 FADC slices

56

4. Observation Technique and Data Analysis

Figure 4.8: The trigger macro cells in the inner camera. The xNN groups used in LT1 lie each within the same macro cell. A trigger is raised if one xNN group within one cell reaches the trigger condition. middle of the reflector. Per pulse, each pixel receives between 50-400 phe, depending on the color of the LED. Because the light flux from the LEDs depends on the temperature and thus can vary, there are different calibration methods: Blind Pixel Method: One PMT is being ’blinded’ with a foil which has a well known attenuation factor of around 1/1000 in order to attenuate the bright light from the LED pulser. Thus, one can determine the absolute LED light flux at the camera by recording the attenuated phe spectrum, identifying the peak coming from the Gaussian electronic noise and the first phe peak and its variance and performing a Poissonian fit to the phe Winston Cones

PMT

Receiver Board Receiver

FADC

MUX controller

Optical Optical splitter delay Optical Link

Receiver

Camera Receiver

Trigger Receiver

Figure 4.9: The readout scheme of the MAGIC telescope. The signals from Cherenkov photons measured with the PMTs are transmitted with optical fibers over 60 m to the counting house. The optical signals are split. One part of the signal goes to the trigger, the other to the FADC.

4. Observation Technique and Data Analysis

57

Figure 4.10: Calibration system of the MAGIC telescope. The pulsar box located in the center of the mirror contains fast and bright LEDs. The light current from the LED can be monitored by the PIN Diode or the “Blind Pixel”, where the latter is located in the camera plane. See text for details. Image credit: S. Commichau spectrum. PIN Diode Method: With a calibrated PIN Diode at a distance of 1.5 m from the LED, the absolute light flux is measured. Excess Noise Factor: In this method, the absolute light flux is not determined, but should be stable within the time the conversion is determined. The Excess Noise Factor (F-Factor) is defined as: σ 2 − σ02 , (4.3) F = 1 + 12 µ1 − µ20

where σi2 denotes the (Poissonian) variance of the ith phe peak, and µ2i its mean value in FADC counts. For i = 0, these are the values for the (Gaussian) electronic noise (pedestal) 1 ¯ = m distribution. Using the relation Var(Q) ¯ ¯ phe · F , as deduced in [92], with Q the PMT’s Q anode output charge, Var(Q) its variance and m ¯ phe the number of photoelectrons hitting ¯ number of photons hitting the photocathode), one the first dynode (all values are for N can derive the following formula: m ¯ phe = F

σ12

¯2 Q − σ02

(4.4)

Thus, one can determine the number of phe produced at the photocathode. Furthermore, if the photon detection efficiency p of the PMT is known, the number of incident photons ¯ phe ¯ = m is N p . The first two methods provide a conversion from the light hitting the PMT to FADC counts, the third a conversion from phe produced at the first dynode of the PMT to FADC counts. Assuming that the light absorption of the winston cones and the plexiglass and the quantum efficiency (QE) of the PMT is constant, these procedures are equivalent. The set-up of the calibration system is schematically displayed in figure 4.10. In the current setup, all three methods are working. Nevertheless, in the further analysis (see section 4.4.1), only the F-factor method is used.

58

4. Observation Technique and Data Analysis

4.3.6 The Central Pixel The PMT in the center of the MAGIC camera is a standard MAGIC PMT, as they are used for the rest of the inner camera. However, its preamplifier is modified [93] to integrate visible light within 500 ns. In a further modification [94], the readout of the Central Pixel was changed such that the MAGIC readout system digitizes the signal of the Central Pixel each time the telescope is triggered. This allows to observe the optical emission from point-like objects (diameter < 0.1◦ ) in parallel to its γ-ray emission. If a pulsar with a known ephemeris is observed, this setup can be used to check the accuracy of the event time stamps [94]. During the observation of the Crab Pulsar, we recorded the signal from the central pixel as well. See chapter 6.6.

4.3.7 Pyroscope To measure the percentage of clouds covering the night sky, a pyroscope is used. This instrument measures the infrared radiation reflected by clouds between 8 − 14 µm. If the infrared spectrum follows the Planck Spectrum, the temperature of the radiating body can be computed from the integral value measured by the pyroscope. It has a FOV of 2◦ and points to the same direction as MAGIC. If a cloud enters its FOV, the sky temperature changes by a small value (typically > 1 K). Thus, the measured sky temperature by the pyroscope can be used to determine the cloud coverage in the FOV. Unfortunately, the sky temperature also depends on zenith angle, air temperature and humidity. Using a calibration depending on these values, the so-called cloudiness is determined from the sky temperature, the zenith angle, air temperature and humidity during the observation. While a low value of cloudiness means good observation conditions, a high value does not necessarily mean bad conditions. However, a fluctuating cloudiness is a strong hint for clouds passing through. One can thus exclude data taken under suboptimal conditions.

4.3.8 The Rubidium and the GPS Clock To measure pulsating γ-rays with periodicities around 30 Hz and the substructure of the pulse profile, the time stamp of each event has to be as precise as < 50 µs. This accuracy is easily achieved with a rubidium clock having a precision of 200 ns. Due to the unavoidable drift of the rubidium clock, it is synchronized with a GPS clock pulse. The time stamps (time and date) of the event is then recorded in UTC (universal time).

4.4 Data Analysis The data being recorded by the DAQ at the telescope site in La Palma contains mainly the 50 recorded FADC slices of each PMT for each triggered event. Additionally, calibration and pedestal events are stored. These kind of events are each recorded at the beginning of a new datataking run and in addition interleaved with the data with 25 Hz. These interleaved events allow to adjust the calibration constants in the offline analysis, which can change due to temperature changes and other influences. Furthermore, technical information from the subsystems (AMC, drive system, camera, weather station, cloud monitor, star guider etc.) are stored in separate files. The recorded data is written to a local RAID system and to tapes to be shipped to the external datacenter, where the standard analysis chain described below is applied. This chain displayed in figure 4.11 consists of the following steps:

4. Observation Technique and Data Analysis

charge: 4phe, arrival time: 3.5ns

PMT 1

charge: 1phe, arrival time: 5.5ns

PMT 4 PMT 5 PMT 6

signal extraction

PMT 2 PMT 3

59

image cleaning and determination of image parameters

charge: 2phe, arrival time: 2.5ns charge: 2phe, arrival time: 6.5ns charge: 3phe, arrival time: 4.0ns charge: 3phe, arrival time: 3.0ns

Figure 4.11: Schematics of the analysis chain. The charge and the arrival time of PMT signals are extracted. In the second step, the image cleaning is performed, where PMTs containing only NSB noise are excluded from the further analysis. Then, the image parameters are determined. The red line denotes the shower image’s main axis. • Signal extraction: The PMT signals measured in FADC units are calibrated. This includes the absolute calibration of the charge content measured in phe and the relative pulse arrival time calibration. In this step, for each pixel, the arrival time and the charge of the signal are stored for further analysis (section 4.4.1). • Image cleaning: Due to the NSB light, the recorded shower images are noisy. The main task of the image cleaning step is to remove pixel signals containing only noise. To distinguish pixels with NSB signals from the ones containing also signals from Cherenkov photons, an image cleaning algorithm is applied, described in section 4.4.2. The cleaning algorithm takes advantage of the fact that the images of extended air showers are compact in extension and time distribution. The resulting images are parameterized. The most important parameters are the ones introduced by Hillas, see section 4.4.3, but also additional parameters are computed, for example those related to the time development of the shower image. • Gamma-Hadron separation and energy estimation: Nearly all events surviving the image cleaning are those containing images of air showers. Nevertheless, a large fraction of these events are of hadronic origin. In the standard analysis, a Random Forest (RF) algorithm based on the image parameters is used to classify each event and a parameter, called HADRONNESS, is computed. This parameter is used to distinguish γ like events from hadronic ones, see section 4.4.4. In a further step, the energy of each event and its incoming direction will be estimated. • Flux determination: If a significant excess of signal events remains after the classification, the γ-ray flux is derived. To do so, the dead time reduced effective observation time is deduced and the energy dependent collection area (see equation 4.12) is computed, using γ-ray MC simulations. The detector and analysis (cut) efficiencies are included in the calculation of the collection area. The detector efficiency below 300 GeV and thus the collection area are strongly energy dependent. The source γ-ray spectrum9 depends on the energy, often following a power law: dN = f0 · f (E) = dA dE dt 9



E Enorm

−α

(4.5)

The spectrum corresponds to the number of γ-ray photons that arrive within a certain energy interval dE, area dA and time dt on Earth

60

4. Observation Technique and Data Analysis

with the spectral index α, f0 the flux normalization and Enorm an arbitrary normalization energy. Due to the energy resolution of below 40%, the determined spectra have to be unfolded. See section 4.4.5. Additional analysis steps, software and algorithms I developed and used for the analysis of pulsars will be described in section 4.5. In the following sections, the standard analysis chain, the modifications done by me and those that were used in the presented analysis will be explained in greater detail.

4.4.1 Extraction of the FADC Signals In the standard analysis, the calibration program callisto first computes the conversion factor from FADC slices to the phe produced at the first dynode of the PMT photocathode, using the F-factor method. To do so, the recorded calibration pulses described in section 4.3.5 are used. They are stored as special events (“calibration events”). The electronic offset of the signal is computed from pedestal events, which are random triggers and contain only NSB. From the mean arrival time of the calibration pulses, the relative timing between different pixels is calibrated. To extract the digitized signal, several methods can be used [91]. In the standard analysis, signals are extracted using the digital filter method [90]. In this method, a pulse shape is assumed for cosmic rays, MC and calibration events. The extracted charge is then the sum of the FADC slices weighed with the pulse shape value. If one assumes that the shape stays constant for different signal amplitudes, one achieves a high signal/noise ratio and a low bias of the estimated charge content. Further details can be found in [91].

4.4.2 Image Cleaning In the signal extraction step, the following values are computed in each pixel: charge (in phe), the arrival time (in ns) and the root mean square (rms) of the pedestal distribution to estimate the noise. The task of the image cleaning is now to identify the images of the Cherenkov showers and thus distinguish them from signals coming from the NSB light. Several studies were made to improve the image cleaning [94, 95, 96]. An efficient algorithm is crucial if one is interested in γ-rays below 100 GeV, where the Cherenkov light output becomes small, or for the data taken in bright nights, where the NSB can become orders of magnitudes larger. 4.4.2.1 Standard Image Cleaning In the image cleaning, one distinguishes between shower core pixels and the ones from the shower boundary. Core pixel have to fulfill stronger criteria than the boundary ones, meaning for example that their charge content has to be above a threshold of 10 phe, while the charge in boundary pixels only has to be e.g. 5 phe. However, as their name says, boundary pixels have to be close to a core pixel or another boundary pixel to be considered as such. A group of core and boundary pixel belonging to the same connected group is called island. In general, the images of γ-ray showers are expected to consist of one single island, while the hadronic ones can also consist of several islands. This is explained by the formation of sub-showers in the development of the extended air shower. However, for low energies, the Cherenkov distribution on ground of a γ-ray shower will be inhomogeneous because of Poissonian fluctuations of the Cherenkov light output and due to the fact that in the shower there are only a few particles producing Cherenkov radiation. Thus, it is highly probable that those showers consist of several islands as well.

4. Observation Technique and Data Analysis

61

Before the installation of the MUX10 FADC system, the standard procedure was to use a threshold of 10 phe for core pixels, demanding that each core pixel has a neighbor fulfilling the same condition. Pixels close to a core or a boundary pixel with charge > 5 phe become boundary pixels. Furthermore, the number of cleaning rings can be chosen. In the standard analysis, the number of rings is set to 1, meaning that all pixels at a distance further than one pixel from a core pixel will be removed from the image. This setting provides a robust and stable analysis, well suited for the analysis of data taken under different conditions. However, the analysis threshold will be rather high at above 80 GeV, significantly higher than the trigger threshold of around 55 GeV. A major improvement of the image cleaning was achieved by [94] and [97] using the time information of the signals as well. Here, a pixel becomes a core pixel when its charge is above a certain threshold and its arrival time within a small time window defined by the arrival times of the signals in the neighboring pixels. Then, the mean arrival time weighed by the charge of all core pixels is computed. A boundary pixel will be accepted as such, if its charge is above the boundary pixel threshold and its arrival time within a small time window of the core pixel mean arrival time. Since the installation of the MUX FADC with its higher time resolution of 0.5 ns, the time information is used in one of the standard settings: • cleaning with the thresholds (10/5) phe (core/boundary) and no timing information is used, or • cleaning with (6/3) phe and (4.5/1.5) ns time windows for core and boundary pixels, respectively. With these last settings, the analysis threshold approaches the standard trigger threshold. In case of data taken in sum trigger mode (see chapter 5), the analysis threshold is then around 45 GeV (for a trigger threshold of 25 GeV). In this case, a lower cleaning threshold is required. Unfortunately, the standard cleaning only allows cleaning levels of around (4.5/2) phe. Reducing the cleaning levels even further will increase the analysis threshold. This is due to the fact that for a too small cleaning level, the computation of the mean arrival time of the core pixels gets dominated by fluctuations of the NSB with an arbitrary time stamp. The tight time windows need to be increased, resulting in noisier images. The algorithm can be slightly improved by computing the mean arrival time only of the core pixel lying within the same connected cluster. 4.4.2.2 Sum Image Cleaning The sum image cleaning algorithm also takes advantage of the fact that images of extended air showers are compact in time and extension. A pixel will be marked as core, if • it lies within a group of 2, 3 or 4NN with the clipped sum11 of the group’s pixel signals is above the corresponding group’s threshold. • the arrival time of the pixel signal is within a small time window of the mean arrival time of the corresponding group. 10 11

multiplexed “Clipped” means here that the amplitude of the signal is limited at a certain threshold. This prevents large accidental signals (e.g. from afterpulses) to dominate a NN group. The clipping and cleaning levels were intensively tested.

62

4. Observation Technique and Data Analysis

Figure 4.12: A few of the NN groups used in the sum image cleaning. blue: 2NN, green: 3NN, red: 4NN.

4. Observation Technique and Data Analysis

63

The idea to form NN groups and to perform a correlated cleaning was already proposed and implemented earlier [97], see section 4.4.2.3. In figure 4.12, a few NN groups are displayed. In the new algorithm presented here, the clipped sum of the NN groups is formed. The following parameters (for x ∈ {2, 3, 4}) can be set: sumthreshxNN: clipxNN: fWindowxNN:

summed xNN cleaning threshold [phe] clipping level xNN [phe] xNN time window [ns]

For convenience, the parameters set by the user were changed to fSumThreshxNNPerPixel, where x is a value between 2 and 4. Together with the core pixel cleaning level (fCleanLvl1), the sum trigger threshold for each group is computed with: sumThreshxNN = fSumThreshxNNPerPixel * fCleanLvl1. This allows to set only the core and boundary cleaning levels and the sum cleaning levels will then be deduced with “reasonable” (default) settings for fSumThreshxNNPerPixel, without having to care for the details of the algorithm. An expert seeking to improve the xNN group cleaning levels, however, has all freedom to do so by varying fSumThreshxNNPerPixel. In figures 4.13 - 4.16, a comparison between the standard and the sum image cleaning is shown. Generally, the sum image cleaning recovers more pixels and thus the content of the shower, measured in phe and called SIZE (see also section 4.4.3), is larger. Figure 4.13 shows an example of an event, which the standard image cleaning fails to recover. Figure 4.17 shows that the energy threshold after cleaning12 for the sum image cleaning lies at 25 GeV, while the standard cleaning has a threshold of 40 GeV. The new sum image cleaning was successfully applied in the analysis of the Crab Pulsar, see chapter 6. A further analysis was applied to the Crab Nebula signal [98]. It was found that the sum image cleaning shows a better performance at low energies and the same performance as the standard cleaning at high energies. Furthermore, the after cleaning rate in nights with a large amount of NSB is more stable and more signal events can be recovered. 4.4.2.3 Calibration Image Cleaning A similar concept as the sum image cleaning was already described in [97]. In this algorithm, the extraction and the image cleaning is performed within the same step. Before a signal is extracted from the FADC signal, the algorithm searches for signals in the surrounding pixels. Only if a signal in x neighboring pixels above a certain threshold depending on the number of pixels in the group is found, it will be extracted. This algorithm also takes full advantage of the fact that the images are compact and that the Cherenkov shower arrives within a small time window. The xNN thresholds were found by considering the probability to find an accidental trigger from NSB above the tested threshold for each group. As in the case of the sum image cleaning, the cleaning level per pixel can be drastically reduced. The sum image cleaning and the calibration image cleaning show a similar performance, considering random triggers from NSB and an analysis threshold (both 25 GeV for observations with the sum trigger). A detailed description of the algorithm can be found in [97]. A direct comparison between the sum image cleaning and the calibration image cleaning 12

The energy threshold is defined as the peak of the MC γ-ray energy distribution after image cleaning, for a γ-ray source following a power law spectrum with spectral index 2.6.

64

4. Observation Technique and Data Analysis

MHCamera Mean RMS

MHCamera

MHCamera Mean 0 RMS 0

1.529 1.286

Mean RMS

9.4

8.8

8.8

8.3

8.3

7.7

7.7

7.1

7.1

6.5

6.5

5.9

5.9

5.3

5.3

4.7

4.7

4.1

4.1

3.5

3.5

2.9

2.9

2.4

2.4

1.8

1.8

1.2

1.2

0.6

0.60° 189mm

4.945 2.002

9.4

0.60 189mm

(a) MC, uncleaned image, energy: 33 GeV

0.6

0.60° 189mm

0.0

(b) standard image cleaning, no surviving pixel

0.0

(c) sum image cleaning, SIZE = 34 phe

Figure 4.13: An image of a MC γ-ray event with energy 33 GeV. While the standard image cleaning with cleaning levels (4.5/3) phe fails, the sum image cleaning recovers the image. Note that the red line in the cleaned image denotes the main axis of the shower image. MHCamera Mean RMS

Mean RMS

MHCamera

8.876 5.544

Mean RMS

7.882 5.68

21

21

21

20

20

20

19

19

19

17

17

17

16

16

16

15

15

15

13

13

13

12

12

12

11

11

11

9

9

9

8

8

8

7

7

7

5

5

5

4

4

3

3

1

0.60° 189mm

MHCamera

1.627 1.807

0

(a) MC, uncleaned image, energy: 35 GeV

1

0.60° 189mm

0

(b) standard image cleaning, SIZE = 106 phe

4 3 1

0.60° 189mm

0

(c) sum image cleaning, SIZE = 110 phe

Figure 4.14: MC γ-ray, the energy is 35 GeV, about the same energy as the one displayed in figure 4.13. However, the measured phe content (called SIZE) is 110 phe. Thus, the standard and the sum image cleaning recover the image. MHCamera Mean RMS

0.60° 189mm

MHCamera

1.054 1.268

Mean RMS

MHCamera

9.633 3.538

Mean RMS

8.052 4.473

15.2

15.2

15.2

14.2

14.2

14.2

13.3

13.3

13.3

12.3

12.3

12.3

11.4

11.4

11.4

10.4

10.4

10.4

9.5

9.5

9.5

8.5

8.5

8.5

7.6

7.6

7.6

6.6

6.6

6.6

5.7

5.7

5.7

4.7

4.7

4.7

3.8

3.8

3.8

2.8

2.8

1.9

1.9

0.9

0.9

0.0

(a) data, uncleaned image

0.60° 189mm

0.0

(b) standard image cleaning, SIZE = 38 phe

2.8 1.9

0.60° 189mm

0.9 0.0

(c) sum image cleaning, SIZE = 48 phe

Figure 4.15: An image of a recorded event. Both the standard and the sum image cleaning recover the image.

4. Observation Technique and Data Analysis

MHCamera Mean RMS

MHCamera

1.096 1.132

Mean RMS

MHCamera

5.543 1.521

Mean RMS

4.65 1.757

7.0

7.0

7.0

6.6

6.5

6.5

6.2

6.1

6.1

5.7

5.7

5.7

5.3

5.2

5.2

4.8

4.8

4.8

4.4

4.4

4.4

4.0

3.9

3.9

3.5

3.5

3.5

3.1

3.1

3.1

2.6

2.6

2.6

2.2

2.2

2.2

1.8

1.7

1.7

1.3

1.3

0.9

0.9

0.4

0.60° 189mm

65

0.0

(a) data, uncleaned image

1.3 0.9

0.4

0.60° 189mm

0.4

0.60° 189mm

0.0

(b) standard image cleaning, SIZE = 22 phe

0.0

(c) sum image cleaning, SIZE = 32 phe

Figure 4.16: A recorded event. Both the standard and the sum image cleaning recover the image. Due to the lower cleaning threshold of the sum image cleaning, it recovers more pixels and thus the phe content (SIZE) of the image is higher.

7000

standard image cleaning sum image cleaning

6000

counts

5000 4000 3000 2000 1000 10

20

30

40

50 60 energy [GeV]

70

80

90

Figure 4.17: The energy distribution of the sum and the standard image cleaning for γray MC and sum trigger simulations, see section 5.4. The energy threshold defined as the peak of the distribution of the sum image cleaning is at the trigger threshold of 25 GeV. The cleaning thresholds were chosen such that the number of pedestal event images surviving the image cleaning is below 1%.

66

4. Observation Technique and Data Analysis

in the analysis of the Crab Pulsar shows that both provide the same sensitivity towards low energy γ-rays. For the detection of pulsed γ-rays, both were used and compared. Similar results were achieved.

4.4.3 Image Parameters After the cleaning step, the image is parameterized, illustrated in figure 4.18. The so-called Hillas parameters [99] turn out to be a robust parameterization. They correspond to 1st, 2nd and 3rd moments along the main and minor axis of the image shower: • SIZE: The 1st moment of the shower image. It is the sum of the signals of the surviving pixels. It can be computed from the main island of the image only, or also include other islands. • WIDTH: The 2nd moment of the image along its minor axis. It is a measure of the transversal development of the shower. • LENGTH: The 2nd moment of the image along its major axis. It is a measure of the longitudinal development of the shower. • M3LONG: The 3rd moment of the image along its major axis. Measures the direction of the shower. This parameter is important to distinguish the shower head from its tail. Further parameters were introduced. Those include: • MEANX, MEANY: The shower’s center of gravity. • CONC: The fraction of the charge contained in the two brightest pixels. • LEAKAGE: The fraction of the total charge that is contained within the outermost ring of the camera. It is a measure for the amount of the shower not contained within the camera. • DIST: The distance between the center of gravity and the source position. It is highly correlated with the impact parameter13 of the air shower. • ALPHA: The angle between the line connecting the center of gravity of the shower with the source position and the major axis of the ellipse. The image of γ-ray showers coming from the source points towards the source position in the camera, thus having a small value of ALPHA. For the very low energy γ-ray showers analyzed in this thesis, this was the most important parameter to distinguish γ-rays from background. Further parameters derived from the time information of the shower were introduced [100]. These are: • TIME GRADIENT: The change of the arrival time along the major image axis. • TIME RMS: The RMS of arrival times of the pixels belonging to the image after image cleaning. 13

The impact parameter is the distance between the telescope’s position and the penetration point of the incident γ-ray direction with the horizontal telescope plane.

4. Observation Technique and Data Analysis

67

32 30 28 26

ALPHA

24 22 20

center of field of view

LENGTH WIDTH

DIST COG

18 16 14 12 10 8 6 4 2 0

Figure 4.18: Definition of image parameters, defined for a moment analysis of a shower image in the camera plane.

4.4.4 Event Classification with Random Forest and Estimation of the Energy The great bulk of events triggering the telescope comes from cosmic hadrons. Efficient cuts are therefore needed to separate γ-rays and hadrons and thus to reduce the background. Usually, these cuts are based on the image (and time) parameters, because these parameters are rather robust under changing observation conditions (e.g. different NSB levels, weather changes, clouds) In order to obtain optimal results, these cuts have to be determined in a high dimensional parameter space. In the standard analysis, this task is accomplished by the RF algorithm, based on multidimensional decision trees. The RF is trained with MC γ-rays as signal events. As background, a measured data sample from a position pointing to a region in the sky where no γ-rays are expected (so-called OFF data sample) is used. Each event (background and signal event) is described by a set of (image) parameters forming a vector. A RF consists of a number of independent binary decision trees based on those parameters. The algorithm performs the following steps: 1. A sample of both, signal and background data, is drawn. 2. A random (image) parameter is chosen. The cut in this parameter that minimizes the Gini index14 is applied and thus the data sample is split into two sets, forming two nodes. 3. For each of the two subsets, another random parameter is chosen and the cut applied that minimizes the Gini index. 4. The above step is repeated (the tree is “grown”) until the remaining sample after the cut either only contains events of one kind (either background or signal), or the number of remaining events is below 10. 5. Steps 1-4 are repeated N times. Each repetition forms a random tree. All trees together build the random forest. 14

The Gini index measures the node purity and is 0 if it contains only one kind of events and 1 if it contains from both the same amount.

68

4. Observation Technique and Data Analysis

This RF is then used to classify events from a sample containing γ-rays and background events. The classification is based on the HADRONNESS parameter, which is 0 if the event is γ-ray like and 1 for background like events. The HADRONNESS value for a given event is determined according to: • The parameters of the event are fed into the random tree i. h , • The event ends up in a terminal node. The value hi is computed as: hi = NhN+N γ with Nh and Nγ the number of hadrons and γs in the terminal node, respectively.

• Step 1 and 2 are repeated for all N trees. The resulting HADRONNESS h of the event is then: h=

1 X hi N

(4.6)

i

In the following analysis, a cut in HADRONNESS h depending on SIZE15 is used to separate signal from background. If one is interested to discover a new source, the cut yielding the highest sensitivity should be chosen. A cut that maximises the Q-factor fullfills this condition. This factor is defined as:

with ǫγ =

Nγ (h) Nγ (total)

Q(h) = p

ǫγ (h) ǫbackground (h)

(4.7)

the signal cut efficiency where Nγ (h) the number of signal events N

(h)

background surviving the cut in h and ǫbackground (h) = Nbackground (total) the background cut efficiency with Nbackground (h) the number of background events surviving the cut in h. The higher the Q-factor, the better the separation between signal and background. In the standard analysis, the parameter ALPHA is not used in the RF, because the different ALPHA distribution of γ-rays (which are expected to peak at ALPHA= 0) and hadrons (flat distribution) are used to extract the excess events from the data. The energy of an event is also estimated from the image parameters. Here, the parameter SIZE and DIST are the most important ones. Furthermore, the light output of a Cherenkov shower depends on the zenith angle of the shower. Also the shape of the shower and its extension, and therefore WIDTH and LENGTH, change for different energies. To estimate the energy in the MAGIC standard software, a modified version of the RF is used. Instead of separating signal from background, this RF is used to determine the energy an event belongs to. With the optimal set of parameters, an energy resolution of 25% for high energies can be achieved. For the lowest energies at around 25-45 GeV, the energy resolution, however, gets worse and reaches values around 40-50%. Further details on the RF method can be found in [101, 102, 103].

4.4.5 Determination of the Differential Flux After applying the γ-hadron separation, the excess events are determined. To do so, the distribution of an image parameter (typically ALPHA) for data taken from the source (called ON region and containing signal and background events) and data taken from an 15

The dependence on SIZE is chosen because for low SIZEs, the images from γ-ray showers and background nearly look the same, thus asking for softer cuts in HADRONNESS than for higher SIZEs where the events can be classified better.

4. Observation Technique and Data Analysis

69

ON Data

21000

OFF Data 20500

entries

Normalisation region

20000

Signal region

19500 19000 18500 18000 0

10

20

30

40

50 ALPHA

60

70

80

90

Figure 4.19: Determination of the number of excess events from the ALPHA distribution. The ON data sample contains γ-rays and background, the OFF sample only background. For ALPHA> 30◦ , the two distributions are normalized. The excess events are determined from the difference of the distribution in the signal region, ALPHA< 20◦ . OFF region (containing only background events) are drawn. From MC simulations, a parameter interval is chosen, where no signal is expected. For ALPHA, this is typically |ALPHA| > 30◦ . The distributions of ON and OFF events are normalized in this region (displayed in figure 4.19), and the excess events in the signal region defined before extracted according to: Nexcess = NON − Fnorm · NOFF (4.8) with Fnorm =

bg region NON

bg region NOFF

the normalization factor in the background interval of the distri-

bution. The significance of the detection can be computed using the assumption that ON and OFF data are independently measured data samples. Thus, the variance of the excess is: Nexcess Nexcess Nσ = 2 (4.9) =p 2 Vexcess NON + Fnorm · NOFF

Since this formula overestimates the statistical significance of the number of excess events, the following formula is used (formula (17) in [104]): s     √ NOFF NON 1 + Fnorm + NOFF ln (1 + Fnorm ) Nσ = 2 NON ln Fnorm NON + NOFF NON + NOFF (4.10) Once the excess events are determined, a differential energy spectrum is computed: f (E) =

dN dE dA dt

(4.11)

which describes the number of photons within the energy interval [E, E + dE] per time t and area A.

70

4. Observation Technique and Data Analysis

The spectrum cannot be determined directly. Thus, the data is binned in estimated energy and the number of excess events in each bin determined. This number is divided by the effective observation time teff (which is the observation time corrected for the dead time) and by the mean collection area in this bin. The collection area Aeff (E) is computed as: Z2π rZmax Aeff (E, θ) = ǫγ (E, θ, r) r dr dφ

(4.12)

φ=0 r=0

with θ the zenith angle, r the impact distance of the incoming photon and φ its azimuth angle. The collection area is determined from MC γ-rays in bins of energy Ei and zenith angle θj : N X
Nγafter cuts (Ei , θj , rk ) 2 2 rk,max − rk,min (4.13) Aeff (Ei , θj ) = π total Nγ (Ei , θj , rk ) k=1

The resulting γ-ray flux in each energy bin Ei is then computed with: f (Ei ) =

X

f (Ei , θj ) =

j

X j

Nexcess (Ei , θj ) ∆Ej teff (θj ) Aeff (Ei , θj )

(4.14)

where the sum is performed over all observed θj bins. 4.4.5.1 Unfolding Methods Due to the limited and biased energy resolution of the detector, the derived energy spectrum is generally different from the source spectrum. To correct for these effects, the spectrum needs to be unfolded. These effects are described by the migration matrix M which projects the original energy bins to reconstructed energy bins and is deduced from MC γ-rays: Yi =

X

Mij Sj

(4.15)

j

where Yi is the number of measured events with an estimated mean energy in bin i and Sj the number of events in a true energy bin j. A simple deconvolution, meaning solving the linear equations given by equation 4.15 is not an ideal solution, because the energy bins of Y and S are strongly correlated and its values fluctuate, Y = Y + ∆Y . Furthermore, in general cases, the number of estimated and true energy bins are not the same, and the matrix M therefore not square. The equations 4.15 are then solved by minimizing the summed squared differences χ0 (S; Y ) between Y and M S, weighed by the covariance matrix of Y and S [105]. More sophisticated methods use regularization, meaning that certain constraints on the true distribution Sj are posed, e.g. that Sj or its derivative to be smoothly distributed. Minimizing w (4.16) χ2 (S; Y ) = χ20 (S; Y ) + Reg(S) 2 with Reg(S) measuring the smoothness of S instead of χ20 alone, takes care of the regularization. The weights w are used to allow different relative weighs between the two terms. In the MARS package, several methods are implemented choosing different weights w and

4. Observation Technique and Data Analysis

71

regularizations Reg: Tikhonov [106], Bertero [107] and Schmelling [108]. In this work, the method of Tikhonov is used. The regularization in this method is given by: Reg =

X  d2 S 2 j

dE 2

(4.17)

j

2

d S The derivative dE 2 is either computed numerically or by using a spline interpolation between the points Sj = S(Ej ). The χ2 (S; Y ) is minimized as a function of S, given the measurements Y .

Forward Unfolding Another approach is to assume a parameterized prior distribution Sj (θ1 , ...θn ) depending on a few parameters θ1 , ...θn describing the assumed model and then minimize the squared difference between M × S and Y . This method is called forward unfolding and is used in this thesis to determine the cut-off value of the pulsed spectrum, described in section 4.5.1 and applied in section 6.8.

4.5 Analysis of Periodic Signals from Pulsars In this work, signals from the Crab source which emits a periodic signal with a frequency of ≈ 30 Hz is investigated. The analysis of these periodic sources allow for additional analysis steps that are not included in the standard chain. An important new parameter is called PHASE φ and is assigned according to the arrival time of the γ-ray at the telescope. It describes the rotational phase of the pulsar, an observer at the position of the sun would observe at a certain time t. During one rotation of the pulsar, φ goes from 0 to 1. The position of the Sun and not the Earth is chosen as a inertial reference. Otherwise, the phase φ would not only depend on the rotation of the pulsar but also on the rotation of the Earth. The drawback of this definition is the need to transform the arrival times of the γ-ray candidates to the solar barycenter (the center of mass of the solar system). The phase of the rotating and decelerating pulsar is computed as: 1 φ(t) = ν0 · (t − t0 ) + ν˙ 0 · (t − t0 )2 + ... 2

(4.18)

where t0 denotes an arbitrary time with defined phase φ0 = φ(t = t0 ) = 0. In case of Crab, φ0 is defined as the arrival phase of the radio main pulse. ν0 is the rotational frequency of d ν(t)|t=t0 its derivative, needed due to the deceleration the pulsar at the time t0 and ν˙ 0 = dt of the pulsar’s rotation. Additional Taylor series terms can be included. For the Crab pulsar, where a monthly radio ephemeris listing the values for t0 , ν0 and ν˙ 0 is available [109] and ν¨ and the following Taylor terms are small, this is not necessary. In this analysis, the program TEMPO [110] is used to perform the transformation from the telescope’s time frame tobs (measured in UTC times) to the solar barycentric arrival time tb [111]16 : r · n (r · n)2 − |r|2 + + ∆tE + ∆tS (4.19) tb = tobs + c 2cd with n the unit vector in the direction of the pulsar and r the vector connecting the solar barycenter with the telescope’s position. The first two terms perform the geometrical (Newtonian) transformation from the telescope system to the solar barycenter. ∆tE is the general relativistic (“Einstein”) delay caused by the gravitational redshift and time dilation 16

The delay caused by dispersion used in [111] is for γ-rays naturally 0.

72

4. Observation Technique and Data Analysis

Component Leading Wing Peak 1 Trailing Wing Bridge Leading Wing Peak 2 Trailing Wing Offpulse

1 1 2 2

Name LW1 P1 TW1 Bridge LW2 P2 TW2 OP

Phase Interval φ 0.88 − 0.94 0.94 − 0.04 0.04 − 0.14 0.14 − 0.25 0.25 − 0.32 0.32 − 0.43 0.43 − 0.52 0.52 − 0.88

Table 4.1: Crab Pulsar phase interval definition [112]. In figure 4.20, the phase regions are drawn.

mainly because of the motion of the Earth, and ∆tS the Shapiro delay, caused by the curved spacetime of the solar system. For Crab, the main and inter pulse P1 and P2 are for all energies, from Radio wavelengths to γ-rays at the same phase position. Thus, a sensitive test to decide if the measured phases φ are equally distributed or not, takes advantage of this fact by assuming that the two peaks are in VHE γ-rays at the same phase position as well. In [112], the phase regions denoted in table 4.1 and displayed in figure 4.20 were defined. Studying the light curve measured by EGRET at the highest energies shows that most of the emission lies within the phase regions defined as P1 and P2. If one assumes that this remains the same for the energy range of MAGIC, computing the significance of the excess found in P1+P2 together using formula 4.10 has a high sensitivity for pulsed emission. As ON phase region, (4.20)

φON ∈ [0, 0.04] ∪ [0.94, 1] ∪ [0.32, 0.43]

is used, corresponding to the P1 and P2 phase regions in table 4.1 and in figure 4.20. As OFF phase region, it is reasonable to use the offpulse phase region (OP in table 4.1): (4.21)

φOFF ∈ [0.52, 0.88]

The normalization factor is then computed from the ratio between these regions. This test is appropriate if one assumes that the positions of P1 and P2 do not change for different energies. As it can be seen in figure 3.5, this is the case for Crab for the whole electromagnetic energy range between Radio wavelengths up to the γ-rays measured by EGRET. This assumption was already made in [3, 94], where the upper limit on the pulsed VHE γ-ray emission was determined. Further periodicity tests are the H-Test [113], the Kuiper test [114], or a χ2 -test of the pulse profile.

4.5.1 Determination of the Cut-Off Energy of the Pulsed γ-Ray Emission The total γ-ray flux in the phase region P1 and P2 as well as the phase resolved γ-ray flux in each phase region defined for the Crab Pulsar measurements by EGRET, follows a power law spectrum above 100 MeV - 6 GeV: dN = f0 fγ,EGRET (E; f0 , α) = dE dt dA



E Enorm

−α

(4.22)

4. Observation Technique and Data Analysis

350

LW2

TW1 Bridge

TW2

LW1

P2

OP

TW1

P1

300

73

LW2 Bridge

TW2

LW1

OP

P2

P1

counts [a.u.]

250 200 150 100 50 -1

-0.8

-0.6

-0.4

-0.2

-0 Phase φ

0.2

0.4

0.6

0.8

1

Figure 4.20: The phase diagram of the Crab Pulsar measured by EGRET for energies above 100 MeV. The phase regions defined in table 4.1 are depicted. with spectral index α = 2.022 ± 0.014. This spectrum is expected to turn over at an energy above a few or a few tens of GeV and to show a cut-off of exponential or super exponential shape, described by: −α  β dN E fγ (E; f0 , α, Ecut , β) = · e−(E/Ecut ) (4.23) = f0 dE dt dA Enorm is used. Ecut denotes the cut-off energy and the parameter β defines the shape of the cut-off. For β = 1, the cut-off is purely exponential. The measurements of EGRET do not allow to determine the cut-off energy and the shape of the turn over of the pulsed emission from Crab. However, these two parameters can be determined from the MAGIC data described in chapter 6. To achieve a good sampling of the cut-off shape, the MAGIC data were subdivided into bins of estimated energy. In each bin i, the number of excess events Niexc with statistical error ∆Niexc is derived from the data. Using MC γ-rays, the effective area Ai (E) (after all cuts) for each estimated energy bin i is determined. The number of excess events ˆ exc (f0 , α, Ecut , β), is estimated by folding the spectrum fγ (E; f0 , α, Ecut , β) (equation N i 4.23) with the effective area in each bin i: Z ˆiexc (f0 , α, Ecut , β) = teff · fγ (E; f0 , α, Ecut , β) · Ai (E) dE N (4.24) To determine the parameters f0 , α, Ecut and β, the least squared difference between ˆ exc (f0 , α, Ecut , β) and N exc is minimized. In the fit, predicted number of excess events N i i the spectrum points measured with EGRET and COMPTEL above 100 MeV are combined with the MAGIC measurements above 25 GeV:  2 ˆ exc (f0 , α, Ecut , β) X Niexc − N X (fj − fγ (Ej ; f0 , α, Ecut , β)2 i + χ2 (f0 , α, Ecut , β) = 2 (∆fj ) (∆Niexc )2 j,HE

i,VHE

(4.25) with fj the measured flux at the energy Ej . The first sum goes over the HE flux points measured with COMPTEL and EGRET, the second sum over the VHE excess events in SIZE bins measured with MAGIC. The χ2 is minimized using ROOT’s17 TMinuit package. In chapter 6, f0 , α, Ecut and β are determined. As one expects, the power law parameters 17

ROOT is an objected oriented analysis framework, often used in high energy physics experiments. See http://root.cern.ch

74

4. Observation Technique and Data Analysis

f0 and α are mostly determined by the measurements of EGRET, while the cut-off parameters Ecut and β by the measurements of MAGIC. The method described here to compute the cut-off energy of the pulsed spectrum has several advantages over the unfolding method described above in section 4.4.5.1. It is a robust method in terms of stability against statistical fluctuations. However, an important limitation is that for this method, the shape of the cut-off needs to be parameterized, here by Ecut and β. From the χ2 value alone, it might be difficult to judge if the cut-off really has the assumed shape. Furthermore, this method is prone to possible systematic errors of the flux measurements by EGRET and to a mismatch of the energy scale calibration of MAGIC and EGRET. To overcome this disadvantages and for consistency checks, the spectrum of the pulsed emission is determined by the measurements of MAGIC alone, using the forward unfolding method with a power law and the unfolding method by Tikhonov. These methods are described in section 6.9.

5 The New Sum Trigger

Figure 5.1: The sum trigger as it is installed on place in the MAGIC counting house on the Roque de los Muchachos on La Palma. A large amount of cables was needed to adjust the relative time delay before summing the signals up. Image by Takayuki Saito. Gamma rays from various cosmic sources have been measured with the MAGIC telescope above an (analysis) threshold of 80 GeV. One type of source – Pulsars – were not discovered before 2007, however. Even though the set-up described in chapter 4 provides the lowest trigger threshold of all current IACTs, the sharp predicted cut-off of the pulsed emission prevented to measure a significant amount of γ-rays from any pulsar. In this chapter, I will describe the idea of the sum trigger which is a new trigger system designed and installed by a small group within the MAGIC collaboration in 2006/07. As a main feature, the new trigger provides an energy threshold which is about a factor of 2

76

5. The New Sum Trigger

lower than the standard trigger. In section 5.1, I will depict several ideas of new trigger systems that could lead to a lower energy threshold. The basic idea of the sum trigger is explained in section 5.2. In section 5.3, PMT noise called afterpulses, which made the implementation of the sum trigger more difficult, are depicted. The simulations will be described in 5.4. Section 5.5 then covers the electronics built in Munich and installed on La Palma, and section 5.6 the calibration of the trigger threshold. In section 5.7, the performance of the standard trigger and the sum trigger will be compared. Figure 5.1 shows a part of the sum trigger installed in the counting house at the telescope site on La Palma. The sum trigger was successfully used in the observations of the Crab Pulsar that led to its detection in VHE γ-rays above 25 GeV. This will be described in chapter 6.

5.1 Reaching the Design Goal of the MAGIC Telescope One of the main design goals of the MAGIC Telescope was to achieve an energy threshold as low as a few tens of GeV. Early plans even aimed at an energy threshold of 10 GeV [115], and thus have a sensitivity overlapping in energy with satellite experiments like EGRET. The components of the telescope were designed in order to reach an energy trigger threshold as low as possible (e.g. large mirrors, high QE PMTs, fast signal sampling). With the setup described in chapter 4, however, the resulting trigger threshold is at a rather high energy of 50− 60 GeV. One main point is that the MAGIC standard trigger as described in section 4.3.3 is not sensitive enough for the low energy showers below 50 GeV. Reasons for this are the following: • The standard trigger requests the signals from a group of compact pixels to be above the trigger threshold. The Cherenkov light of low energy showers, however, is often distributed over more pixels which are not necessarily connected. • The standard trigger uses rather large effective coincidence times between pixels of 4-8 ns. This prevents the LT0 discriminator threshold to be lowered below 6 phe as otherwise triggers coming from Poissonian fluctuations in the NSB will dominate, resulting in huge trigger rates. The DAQ of the telescope, however, can only handle data rates below 1.4 − 2 kHz. • At low energies, roundish shower images have a higher probability to trigger the standard trigger. These images have a random orientation image parameter ALPHA due to the compact next neighbor condition. Since this parameter is important to distinguish γ–rays from hadronic background events for low energies, the background discrimination is rather inefficient. It has already been tried at the beginning of the telescope operation to lower the energy threshold using special trigger configurations. The following ideas were proposed: • Accept only triggers from events lying on a donut shaped ring around the center of the camera between 0.2◦ and 0.8◦ [116]. This idea is based on the fact that more than 80% of the γ-ray Cherenkov shower images with energies below 50 GeV lie in this ring. The ring structure was implemented later in the topology of the sum trigger, see section 5.4.2. • Use a 5NN trigger setting. This allows to lower the LT0 discrimination threshold. MC simulations showed that the sum trigger finally implemented and described in this chapter, however, is more sensitive to low energy γ-rays.

5. The New Sum Trigger

77

• Based on groups of several pixels (for example 8-16 pixels), a trigger is raised if at least m out of the N pixel reach the discriminator threshold. This trigger (used for example in the H.E.S.S. experiment [117]) regards the variety of image shapes found for γ-ray showers. • Sum trigger. This idea was first proposed in [118], but not consequently followed. This trigger will be described in the next sections. • Improvement of the analysis software, in particular the image cleaning. Add image parameters that were deduced from the time information of the shower image. In the following sections, the basic idea of the sum trigger is explained, its features described and the problem of the intrinsic PMT noise (so-called afterpulses) discussed. Then, the performance of the sum trigger is compared with the MAGIC standard trigger.

5.2 The Sum Trigger: The Basic Idea Studying images of low energy Monte Carlo γ–ray events (displayed in figure 5.2) reveals the following properties: • the Cherenkov light arrives within a small time window of around 2 ns at the Camera plane, • the spatial distribution of the Cherenkov light has large fluctuations and is often not distributed within only a few PMT, as it was assumed for the MAGIC standard trigger configuration and • 85% of the Cherenkov light is distributed within a ring between 0.15◦ and 0.85◦ , relative to the source position (figure 5.3). Based on these simulations, we decided to implement an analog sum trigger. The main component of this trigger is to sum up the analog PMT signals within one of a set of overlapping patches of N connected PMTs. The DAQ trigger is raised when the summed signal in one patch reaches the sum discriminator threshold. Compared to the standard trigger, this concept has the following features: • There is a free choice on the shape or the distribution of the Cherenkov photons within one patch. • A small “natural” coincidence window of the order of a few ns. This can be changed by varying the PMT signal’s full width half maximum (FWHM) before summing up the signals. • Also small signals of a PMT below the standard trigger’s discriminator threshold contribute to the trigger signal. • Because the shower signal is distributed over several loosly connected PMTs, the summed signal will show a better signal to noise ratio than a single PMT signal. • A reduced sensitivity to Poissonian fluctuations of the NSB light. As these are only qualitative arguments, I implemented the trigger simulation for different patch shapes and overlap regions into the existing detector simulation program. The different parameters (sum trigger discrimination level, patch shape, patch size and topology, signal shape and timing) were optimized with MC trigger simulations (see section 5.4).

78

5. The New Sum Trigger

49

22

46

21

43

19

40

18

37

17

34

15

31

14

28

12

25

11

22

10

19

8

15

7

12

6

9

4

6

3

3

0.60° 189mm

0

2

0.60° 189mm

0

(a) 21 GeV γ-ray

(b) 22 GeV γ-ray

28 26 25 23 21 19 18 16 14 13 11 9 7 6 4 2

0.60° 189mm

0

(c) 24 GeV γ-ray

Figure 5.2: The distribution of Cherenkov photons (black dots) of simulated γ-ray in the Camera of the MAGIC telescope. Note that the Cherenkov photons cover a large area containing several PMTs within the camera plane.

5. The New Sum Trigger

79

Figure 5.3: The radial averaged distribution of Cherenkov photons for simulated γ-ray showers with initial energy between 20 and 30 GeV. 85% of all Photons lie on a ring between 0.15◦ - 0.85◦ . The trigger area therefore should cover this ring.

5.3 Photomultiplier Tubes Afterpulses Compared to other photodetectors like Hybrid Photomultiplier Tubes (HPD) or G-APDs, PMTs show large intrinsic noise, called afterpulses. The origin of these pulses is not completely understood. They might originate in electrons accumulating at the first dynode. From time to time the electrons are released and thus generate a large signal by hitting the anode. Another reason for afterpulsing might be due to the imperfect vacuum sealing of the tubes. While producing the PMTs, water molecules may remain inside the vacuum tube. By the acceleration of the phe between photocathode and first dynode, they are able to split the molecules and produce H+ and OH− ions. These ions are accelerated between first dynode and the photocathode. By hitting the photocathode, they can produce a large amount of phe amplified by the normal multiplication process within the PMT and thus produce large signals. The time delay of the afterpulses depends on the dimension of the tubes and on the voltage between the photocathode and first dynode, and is of the order of µs [119]. In figure 5.5, the measured and simulated rate as a function of the pulse height for a PMT exposed to NSB is displayed. The rate at pulse heights < 6 phe is dominated by NSB. Afterpulses become important above 6-7 phe, and have a rate of 100 Hz for a pulse height above 30 phe. The trigger rate of a “naive” sum trigger, which sums up the signals of a few PMTs without taking care for afterpulses, would be completely dominated by afterpulses. This could only be compensated by a rather high trigger threshold and thus reduce the sensitivity towards low energy γ-ray showers. Several ideas came up to prevent afterpulses from triggering. First we considered a veto for high pulses or demanding for at least 2 large pulses in one trigger patch. This turned out to be rather complicated and unpractical. In the end, a simple and effective solution to solve the afterpulse problem was found: Before

80

5. The New Sum Trigger

Figure 5.4: Strong signals are clipped before summed up. the summation, the analog PMT signals are clipped at a certain pulse height1 , displayed in figure 5.4. This significantly reduces the contributions to the trigger signal of afterpulses. The sensitivity for γ–rays, however, will be reduced as well, in particular at the lowest energies.

5.4 Implementation of the Sum Trigger and Afterpulses into the Detector Simulation The MC simulation software used in the MAGIC collaboration consists of three different programs. The extended atmospheric air showers are simulated with CORSIKA [121]. In the original CORSIKA code, detailed simulations of an air shower initiated by photons, protons or nucleii have been implemented. The program tracks all particles created in the avalanche, taking into account the reactions of the shower particles with the atmospheric molecules and atoms, the Earth’s magnetic field and also allows to apply different atmospheric models. In later versions of CORSIKA, the process of Cherenkov radiation emitted by relativistic charged particles (electrons, positrons, muons and hadrons) was implemented, but not the atmospheric scattering and absorption of the produced Cherenkov light. The detector simulation of MAGIC consists of the reflector and the camera program. The reflector program uses the simulated Cherenkov photons from CORSIKA and simulates the imaging of the produced Cherenkov photons by the mirror onto the camera plane, including aberration, reflectivity losses and scattering. Furthermore, it tries to adapt for atmospheric extinction and scattering of the Cherenkov light, which is not included in CORSIKA yet. The photons hitting the camera plane are then processed by the camera program. Based on the output of the reflector program, it simulates the PMT response, the electronic shaping of the signals, the standard trigger and the digitization by the FADC. Furthermore, the NSB is added to the Cherenkov signals. For the Crab observation regions, a value of 0.18 phe/ns in each pixel is used. The standard camera program did not include afterpulse simulations nor the sum trigger was set up.

5.4.1 Afterpulse Simulations The camera program uses a database in which typical NSB signals for different levels of background lights are stored. For each photon hitting the PMT, there is a small chance to 1

The idea to limit the amplitude of one pixel was originally developed for trigger concepts for G-APDs, where similar large fake pulses can occur due to optical cross-talk. [120].

5. The New Sum Trigger

layout number 3 5 6 7 8 9 10 11 12 13 14

pixels per patch 32 (4×8) 16 (2×8) 14 (2×7) 14 (2×7) 14 (2×7) 14 (2×7) 15 (3×5) 18 (3×6) 18 (3×6) 21 (3×7) 21 (3×7)

trigger thresh. Vthreshold [mV] 35 30 26 26 27 27 27 29 29 31 31

81

@10 GeV 151 85.8 264 58 102 142 150 207 316 237 263

collection area [m2 ] @20 GeV @30 GeV 2549 8394 1783 6828 3296 9876 2929 10181 1899 6938 2324 8057 2479 8558 3054 10112 5404 17134 3230 10591 2645 10184

@40 GeV 18967 17287 20874 23626 16970 18979 20080 23032 32497 23991 22512

Table 5.1: The collection areas at 10, 20, 30 and 40 GeV for different sum trigger layouts. produce an afterpulse of a certain height. The probability that an afterpulse is produced after a NSB photon hit the photocathode is given by: f (A > A0 ) = 0.0164 · e−0.25·A0 phe−1

(5.1)

with f (A > A0 ) the probability that an afterpulse is generated with amplitude A larger than A0 phe. In figure 5.5, the rate of signals with amplitudes larger than the discriminator threshold on the x-axis is displayed, when a PMT is exposed to NSB. The afterpulse rate depends on the photon intensity hitting the photocathode. The release of afterpulses will happen in the order of a few µs after the light reaches the cathode. Thus, afterpulses are mainly uncorrelated in time with the incident photon, but only dependent on the intensity of the NSB. As can be seen in figure 5.5, the rate is dominated by NSB below the kink at around 6 phe. Above this kink, afterpulses dominate. I simulated several clipping levels and found that a clipping level at around 5-7 phe discriminates the afterpulses efficiently enough. There is a trade off between the sensitivity to γ-rays, in particular at the lowest energies, reduced by the clipping and the maximum rate we can accept. As we aimed at a trigger threshold between 25-35 phe per shower, a trigger from afterpulses alone is highly improbable (corresponding to a 6 − 8σ effect). We extended the existing NSB database with the afterpulse simulations. To achieve a good statistics of afterpulses, the total simulated NSB time was extended from 10’000 ns to 500’000 ns. The simulated database therefore includes afterpulses with a pulse height of up to 45 phe.

5.4.2 Sum Trigger Simulations: Topology, Thresholds and Shape The sum trigger simulation was added to the existing trigger simulations. Before the signals are summed up, the AC coupling is simulated by first integrating the analog trigger signal and then subtracting the mean value. Several trigger configurations were simulated, with patch sizes between 8 - 32 PMTs and different overlapping regions. As most Cherenkov photons are projected on a ring between 0.15◦ − 0.85◦ in the camera, the whole trigger area covers only this ring. This structure will already provide some hardware suppression of hadronic background events. To keep a rotational symmetric configuration, the patches were designed for one sector of the camera

82

5. The New Sum Trigger

Figure 5.5: The count rate above a discriminator threshold denoted on the x-axis for one MAGIC PMT exposed to NSB. The red crosses denote the simulation with a NSB level of 0.18 phe/ns, the black ones the measured spectrum. Note that above 6 phe, afterpulses dominate the count rate.

5. The New Sum Trigger

83

and then rotated 5 times each by 60◦ . The trigger thresholds for each configuration was chosen such that random triggers from NSB and afterpulses were below 50 Hz. The clipping threshold was varied between 5 and 10 phe. It was found that for a clipping threshold between 5 and 7 phe, the resulting collection area was constant within errors. Thus, we decided to fix the clipping level at 6 phe. In the real setup, the discriminator and clipping levels fluctuate. To take this into account, I added Gaussian fluctuations (σ = 1.5 phe) to the clipping and the discriminator thresholds. Since this did not affect the resulting collection trigger area, this feature was removed from the simulations. We also simulated different bandwidths of the sum trigger electronics and found that decreasing the analog PMT pulse width efficiently suppresses the NSB. Nevertheless, if the pulse width is too short (below ≈ 2 ns), the trigger efficiency decreases as well. This has a physical reason: The duration of a small energy γ–ray event is about 2-3 ns. By decreasing the analog pulse widths, one loses showers with slightly larger time spread. Furthermore, the electronic signals are expected to jitter in time at the order of < 1 ns. To assure that the whole γ-ray signal is contained within the trigger signal, we fixed the pulse width at 3.5 ns. We also examined the influence of the optical PSF on the trigger efficiency of the telescope. The idea was to defocus the mirror slightly. Thus, the Cherenkov light would be distributed over a larger area, resulting in less peaked single PMT signals. Thereby, the Cherenkov light would not be clipped, but it would still trigger the telescope. Using simulations, we realized that this method does not increase the sensitivity, and we dropped the idea. The final sum trigger layout is displayed in figure 5.6. The trigger area is contained within 0.26◦ − 0.83◦ and contains 216 PMTs. In table 5.1, the collection areas for a few tested layouts are noted. It was found that an overlap factor of 2 (i.e. each PMT in the trigger area lies in two different trigger patches) increases the collection area significantly. An overlap factor of 3 would have increased the sensitivity even more, but hardware constraints did not allow such a configuration. These hardware constraints also led to the choice of the somewhat strangely shaped sum trigger patches.

5.5 Electronics and Installation on Site Parallel to the Monte Carlo simulations of the sum trigger, the necessary electronics was planned and designed at the electronics workshop of the MPI in Munich. To keep the design simple but flexible, the following boards, schematically displayed in figure 5.7, were chosen: 1. 35 type I boards with 8 channels each, where the clipping is done. 2. 40 type II boards. These boards can sum the signals of up to 8 pixels and form the so-called subpatches. They have 3 equal output lines, allowing each subpatch to lie within 1-3 patches. 3. 24 type III boards. These boards sum up the signals of up to 3 output signals from the type II boards Additionally, an “OR board” was needed in order to run the sum trigger and the standard MAGIC trigger in parallel. The number and characteristics of these boards led to the following constraints:

84

5. The New Sum Trigger

8.0 7.5 7.0 6.5 6.0 5.5 5.0 4.5 4.0 3.5 3.0 2.5 2.0 1.5 1.0 0.5

0.60° 189mm

0.0

(a) The sum trigger patches.

(b) The subpatches of the final sum trigger layout, each consisting of 6 pixels. Each subpatch belongs to two different patches, resulting in a overlap factor 2.

Figure 5.6: The final configuration of the sum trigger. Displayed in (a) is the total trigger area (salmon) and the base sum trigger patches (reddish colors). This configuration contains 18 pixels per sum trigger patch and occupies in total 216 pixel. Each pixel within the trigger area is contained in exactly 2 different sum trigger patches. To get the total trigger layout, each base sum trigger patch has to be rotated 5 times by each 60◦ around the center of the camera. In (b), the subpatches are displayed. Each patch consists of 3 subpatches with 6 PMTs each.

5. The New Sum Trigger

85

input signals

8 output channels

8 input channels

clipped output signals

Clipping Board I

to board II

3 output channel

8 input channels

clipped input signals

Summation Board III

1 output channel

3 input channels

Summation Board II

summed, clipped output signal

to board III

summed, clipped output signal from up to 3 subpatches

to discriminator

Figure 5.7: Schematics of the three sum trigger boards. Board I performs the clipping, board II sums the signals from up to 8 PMTs. This board defines the subpatches, which can contain maximal 8 PMTs. The output signal from 3 board II types are then summed in board III, forming the full patches. A trigger is raised if the signal from one patch exceeds the discriminator threshold. In the final sum trigger layout, one patch consists of 3 subpatches, which sum the signal of 6 PMTs. 1. The number of pixels in the trigger area is smaller than 280 (35× 8 clipping channels) 2. The total number of subpatches is smaller than 40 (i.e. ≤ 6 subpatches in each sector). 3. A subpatch can contain up to 8 pixels. 4. To provide an overlap between patches, each subpatch can be contained in up to 3 patches. 5. The number of patches is smaller than 24 (e.g. ≤ 4 patches per sector). All components were installed and tested at the MAGIC site at La Palma in October 2007. After the commissioning phase, regular data taking on Crab started by the end of October 2007.

86

5. The New Sum Trigger

5.6 The Calibration of the Sum Trigger The trigger thresholds were set in mV. To calibrate the trigger threshold and the clipping level in phe, we used 3 different methods: • Using the calibration- and afterpulses. This was done after the installation of the sum trigger. Afterpulses have a very short time spread and can be seen as δ-functions, compared to the time resolution of the MAGIC electronics [122]. Thus they show a ’perfect’ pile up of phe in time. The waveform of the afterpulses and the calibration pulses was measured with an oscilloscope. Using the F-factor method (see section 4.3.5), one can compute the average phe content of the calibration pulses. To compute the signal ratio between height (in mV) and phe content for afterpulses, one has to take into account the different time spread of both pulses. (The phe content is proportional to the area below the pulse). Thus one achieves the calibration: C = (Vcalib [mV]/N phecalib ) · F W HMcalib /F W HMafterpulse ≈ 0.037[mV/phe] • Matching the measured NSB spectrum (i.e. rate vs. discriminator threshold) with a simulated one. See figure 5.8. As the NSB can vary by a large factor depending on the sky region, this is not a very precise measurement. • Comparing the measured distribution of phe for cosmic ray events with the simulated one. See figure 5.9. The first and the second method resulted in a trigger threshold of 27 phe and a clipping level of 6 phe, while the third method showed a trigger threshold of 24 phe and a clipping level of 5.5 phe. As the first and the third method have an estimated systematic error of 10 − 15%, these measurements agree within the error. In the later analysis, a systematic error for the threshold calibration of 10% was used.

5.7 Comparison with the MAGIC Standard Trigger During each observation, the sum trigger was run in parallel to the MAGIC standard trigger. This allows to directly compare the performance of the sum trigger in comparison with the standard trigger, e.g.: the homogeneity of the COG distribution of the shower image, displayed in section 5.7.1, the energy threshold for γ-rays (section 5.7.2), the collection area in section 5.7.3, or the sensitivity towards a DC signal of the Crab Nebula for different SIZE cuts, illustrated in section 6.11.

5.7.1 Homogeneity By plotting the distribution of the COG for each shower image within the MAGIC camera, one achieves a visualization of the trigger homogeneity, see figure 5.10. While the standard trigger is expected to have a smooth distribution, the sum trigger shows a six fold structure due to the choice of the sum patches. For low SIZEs, the COG distribution is for the sum trigger more homogeneous than for the standard trigger.

5.7.2 Energy Threshold As defined above, the trigger threshold is the peak of the energy distribution for a γ-ray source with a spectral index of 2.6. Such a source was simulated with Monte Carlo γ-rays.

5. The New Sum Trigger

87

rate [Hz]

clipped NSB spectra Measured NSB rates for 1 patch. Normalisation: 632mV = 27phe

Simulated NSB rates, clipping at 6.2 phe

5

Simulated NSB rates, clipping at 6.5 phe

10

Simulated NSB rates, clipping at 7.0 phe

Simulated NSB rates, clipping at 7.25 phe

104

3

10

102

10

20

30

40

50

60 70 threshold [phe]

Figure 5.8: The measured NSB spectrum (black markers) and the simulated spectrum (colored markers) for different clipping levels. Remark: The conversion factor computed here is 0.042 mV/phe, compared to 0.037 mV/phe using the afterpulse method. The flat tail to the right in the data is due to triggers from cosmic rays (mostly hadrons). This is not implemented in this simulations.

88

5. The New Sum Trigger

proton Monte Carlo Crab OFF 200

Entries

150

100

50

0 0

50

100 150 SIZE [phe]

200

250

Figure 5.9: The distribution of SIZE in phe for data taken from a OFF region (where no γ-rays are expected from) and for a proton Monte Carlo dataset. The distributions of data and simulations agree well for SIZE> 50 phe. For lower SIZEs, the distribution fit not ideally. This is due to the rather high level of NSB in the OFF region. The energy distribution after trigger is plotted in figure 5.11. The energy threshold was lowered by about a factor of 2 from 50-60 GeV to 25 GeV. As the energy threshold of protons was reduced as well, this resulted in an estimated trigger rate of at least a factor 2.5 higher than by using the standard trigger. To reduce the energy threshold, we set the sum trigger discriminator to a rather low level, risking that Poissonian fluctuations of the NSB trigger the telescope. We estimated that thus about 50% of all trigger come from the NSB and not from cosmic rays. The total rate was increased from 200 Hz to 800-1’300 Hz, resulting in a total amount of 3.5 TB (raw) data for the 42 hours of Crab observation in sum trigger mode.

5.7.3 Collection Area and Sensitivity for Low Energies The energy dependent collection area (see section 4.4.5) characterizes the performance of an air Cherenkov telescope. Comparing the standard trigger with the sum trigger in figure 5.12 reveals an improvement of around one order of magnitude at 20-30 GeV.

5.8 Future Improvements Even though the sum trigger shows a good performance, there is a lot of room for improvement. There is no automatic setting of the trigger thresholds, resulting in a need to manually readjust the thresholds during data taking by one of the sum trigger experts using screwdrivers, a voltmeter to read the thresholds and NIM scalers to check the resulting rate. Due to changing atmospheric conditions, frequent threshold readjustments during

5. The New Sum Trigger

89

COG, Lvl 1 400 1600 300 1400 200 1200 COGY [mm]

100 1000 0

800

-100

600

-200

400

-300

200

-400 -400 -300 -200 -100 0 100 200 300 400 COGX [mm]

0

(a) standard trigger - MC

(b) standard trigger - data

(c) sum trigger - MC γ-rays

(d) sum trigger - data

Figure 5.10: The distribution of the COG of the shower images within the MAGIC camera for the standard trigger and the sum trigger. Figure a) and c) shows the expected COG distribution gained from Monte Carlo simulations for the standard trigger (also dubbed Lvl1 trigger) and the sum trigger, respectively. Remark the six fold distribution which can also be seen in the data. The inhomogeneity in the lower right corner is due to a malfunctioning sum trigger board.

90

5. The New Sum Trigger

differential rate [a.u]

Monte Carlo energy distribution, after trigger for a power law with spectral index -2.6.

18000 std MAGIC trigger

16000 sum trigger

14000 12000 10000 8000 6000 4000 2000 0 0

20

40

60

80

100

120

140

160 180 200 MC energy [GeV]

Figure 5.11: The differential rate for MC γ-rays, following a power law with spectral index of 2.6 for the MAGIC standard trigger and the sum trigger. The peak of the distribution is a measure for the energy threshold. The standard trigger has a threshold of 55 GeV, while the sum trigger has a threshold energy of 25 GeV. The distributions were normalized for MC energies above 120 GeV.

Figure 5.12: The collection area for the MAGIC standard trigger and the sum trigger. The gain in sensitivity of the sum trigger at energies below 30 GeV is about 1 order of magnitude. Remark: For energies above 200 GeV, the standard trigger is more efficient than the sum trigger. The reason for this is the smaller trigger area within the camera. (Sum trigger: ring 0.26◦ − 0.83◦ , standard trigger: full area 0◦ − 1.05◦ )

5. The New Sum Trigger

91

data taking were needed. To improve the handling, a computer controlled threshold steering was installed in summer 2008 [123]. With the program SUMO, the single patch rates are checked and the thresholds adjusted, if the rate of one patch goes below or above a certain threshold. This window is 700-800 Hz for pulsars and 400-500 Hz for extragalactic sources. To increase the sensitivity, the trigger area can be enhanced and a larger overlap between neighboring pixels is provided. An enhanced trigger area will increase the effective area for large energies and thus yield a better sensitivity in this domain. The larger overlap will increase the sensitivity for lower energy showers.

5.9 Conclusion As shown in this chapter, the sum trigger [124] is a simple, cost-effective and efficient trigger implementation for Cherenkov telescopes. Its installation at the MAGIC telescope site effectively reduced the trigger threshold by a factor of 2 from 55 GeV to 25 GeV and increased its collection area at 35 GeV by nearly an order of magnitude. At energies above 150 GeV, it is slightly less efficient than the standard trigger. This is due to the reduced trigger area. The sum trigger played a crucial role in the detection of the Crab Pulsar, described in the next chapter.

6 Observation of the Crab Pulsar using the Analog Sum Trigger In this chapter the first detection of pulsed γ-rays from Crab will be described in detail. I will start with the description of the data sample (section 6.2) and the selection criteria that were applied to ensure that only data are analyzed that was taken under optimal observation conditions (section 6.3). The following section 6.4 covers the γ-ray simulations which are adjusted for the observation of Crab with the sum trigger. The applied cuts adjusted for this data sample and the reconstruction of the pulsar phase for each event is described in section 6.5. The optical data collected in the central pixel of the camera was analyzed as well and the optical pulsation of Crab determined. This is described in section 6.6. The optical data is recorded with the same DAQ system as the γ-ray data, and the conversion of the timestamps to phases is also the same. The analysis of the optical pulsation is therefore an important cross-check to ensure that the timestamps are treated correctly in the data taking and analysis chain. In section 6.7, the first detection of pulsed γ-rays above 25 GeV from Crab is described. The spectrum in the cut-off region is modelled with a power law spectrum with exponential and super exponential cut-off, taking into account the joint measurements of COMPTEL, EGRET and MAGIC (section 6.8). This modeling is done for the total pulsed flux of the main and inter pulse together, as well as separately for each pulse. In addition, the spectrum in the cut-off region is determined based only on the measurements of MAGIC in section 6.9, again for the total flux and for each pulse separately. The pulse profile itself contains physical information on the geometry of the emission region. In section 6.10, a closer look at the morphology of the pulse profile at different energies is therefore given. The data taken with the sum trigger can also be analyzed to search for an unpulsed signal. This is demonstrated in section 6.11, where the measurement of unpulsed γ-rays from Crab at energies around 50 GeV will be described. Note that this chapter contains a listing of all obtained results. An interpretation of these results, including the implications for current pulsar models, is given in chapter 7. The detection of the pulsed emission and the determination of the cut-off energy has been published [125].

6.1 Introduction The γ-ray source candidate “Crab” contains two (potential) γ-ray sources: The Crab Pulsar and its surrounding Crab Nebula. The later is a bright and constant source of VHE γ-rays (e.g. [1, 3, 49, 126, 127, 128, 129]) and therefore taken as the VHE standard candle. The performance of an IACT is given in the percentage of the Crab flux the telescope can measure with a significance of 5σ within 50 hours. Furthermore, lacking a calibrated γ-ray test beam in the vicinity of the Earth, Cherenkov telescope have to rely on their simulations of air showers, the atmosphere, the local magnetic field and the detector. The energy spectrum from the Crab Nebula is therefore compared with measurements from

94

6. Observation of the Crab Pulsar using the Analog Sum Trigger

Figure 6.1: Measurement of pulsed VHE γ-rays from Crab with MAGIC using data taken in 2005. The excess in the marked phase region, which corresponds to the position of the two peaks in the pulse profile measured with the EGRET experiment, has a significance of ≈ 2.9σ. The (analysis) energy threshold was ≈ 60 GeV [3]. Figure provided by N. Otte [94] other instruments as a consistency check. Since the beginning of ground based VHE γ-ray observations more than 20 years ago, the Crab Pulsar in the center of the Crab Nebula has been a main target [46, 47, 48, 49, 50, 51]. Nevertheless, even the lowest energy threshold of above 50 − 60 GeV [3, 83, 130] or higher [48, 131], achieved with Cherenkov telescopes, did not allow to detect pulsed emission from any pulsar. The analysis of data taken with MAGIC in 2005 using the standard setup, however, revealed a hint for a pulsed signal for γ-ray energies above 60 GeV [3], displayed in figure 6.1. The upper limits of the pulsed VHE γ-ray flux derived from different ground based observations are displayed in figure 6.2, together with the spectrum measurements from EGRET. It was found that the 2σ upper limit for the cut-off energy was 30 GeV, assuming an exponential cut-off shape. The sharp spectral cut-off in the energy range between a few GeV - tens of GeV requires thus a very low energy threshold at about the energy of the cut-off. Due to the steepness of the spectrum in the cut-off region, lowering a detector’s energy threshold results in a large increase of the number of measured γ-rays. After the layout design of the sum trigger and the construction of the electronic parts, the sum trigger was installed and commissioned in September 2007. Data taking on Crab Pulsar started at the end of September. By February 2008, we had observed Crab for about 40 hours.

6.2 Data Sample Since the sum trigger operation was not automatized at that time, each observation was accompanied by one of the sum trigger experts to ensure a good performance under chang-

6. Observation of the Crab Pulsar using the Analog Sum Trigger

95

Figure 6.2: This figure shows on the left side (black circles) the EGRET spectrum measurements. On the right side, the flux upper limits (arrows) found with ground based Cherenkov telescopes are denoted. The dashed line denotes the fitted power law spectrum, adjusted with an exponential cut-off (see equation 3.2), where the power law part was derived for the EGRET points. For the cut-off energy, we found an upper limit of 30 GeV [3]. The red arrow denotes the energy threshold of the telescope using the standard setup. Figure provided by N. Otte [94]

96

6. Observation of the Crab Pulsar using the Analog Sum Trigger

2500

2000

accepted

counts

1500

1000

500

0 0

10

20

30

40 50 60 cloudiness [%]

70

80

90

100

Figure 6.3: The distribution of cloudiness for the whole observation period.

ing observation conditions. The single patch rates were counted by NIM scalers, and the thresholds of the 24 patches set by a screw driver. Because other sources were also observed with the sum trigger, the patch discriminator thresholds had to be changed within a few minutes when the telescope was moved from one source to the next. All data were taken at zenith angles below 20◦ . At larger zenith angles, the shower develops further away and less Cherenkov light arrives at the telescope, resulting in an increased trigger threshold. During all observations, the standard trigger was run in parallel to the sum trigger, connected by an OR-module. Every event was flagged with a bit pattern describing the type of event (calibration, pedestal, data) and which trigger fired (LT1, LT2, Sum). This allows to compare the performance of the sum and the standard trigger directly. Great care was taken to ensure that only data under optimal conditions was taken. Still, about 50% of all data did not fullfill the strong criteria depicted in the next section. To push the trigger threshold as low as possible, we lowered the discriminator thresholds to values such that the total trigger rate was at 800 − 1000 Hz, close to the maximal rate of 1.4-2 kHz allowed by the DAQ system. The data taking in La Palma is subdivided into intervals of roughly 3 weeks, corresponding to the period when observation without moon is possible. We took sum trigger data between period 57 in October 2007 and period 63 in February 2008.

100

100

80

80 cloudiness [%]

cloudiness [%]

6. Observation of the Crab Pulsar using the Analog Sum Trigger

60

40

20

60

40

20

0

54394.2 54394.4 54394.6 54394.8 MJD [days]

54395

0

54395.2 54395.4

54412

54414

100

100

80

80

60

40

20

54420

54422

60

40

20

54438

54440

54442 54444 MJD [days]

54446

0

54448

54466

(c) Period 61

100

100

80

80

60

40

20

0 54490

54468

54470 54472 54474 MJD [days]

54476

54478

(d) Period 62

cloudiness [%]

cloudiness [%]

54416 54418 MJD [days]

(b) Period 60

cloudiness [%]

cloudiness [%]

(a) Period 59

0

97

60

40

20

54492

54494 54496 MJD [days]

(e) Period 63

54498

0

54444.02

54444.04

54444.06 54444.08 MJD [days]

54444.1

(f) Night 20071210

Figure 6.4: The figures (a) - (e) show the cloudiness curve as a function of MJD. As an example, figure (f) shows the detailed cloudiness for the night 20071210. The large fluctuations during this night indicate bypassing clouds or other changing weather conditions, leading to the decision to remove this night from the analysis.

98

6. Observation of the Crab Pulsar using the Analog Sum Trigger

140 accepted

120

counts

100 80 60 40 20 0 150

200

250

300 350 400 cleaning rate [Hz]

450

500

550

Figure 6.5: The distribution of the event rates after cleaning and removal of nights not passing the cloudiness cut.

25

SIZE>200 100 0.9 · 1016 GeV was achieved, while the measurements with AGN revealed EQG1 > 0.21 · 1018 GeV [153] on a confidence level of 95%. Compared to these sources, the Crab Pulsar is rather close (L ≈ (2.0 ± 0.5) kpc), but the maximum γ-ray energy is higher than for γ-ray bursts at tens of GeV and the emission is on extremely fast time scales. The width of the main pulse P1 is approximately 3 ms. As summarized in table 6.6., the position of the peak can be determined with a precision of ∆ΦP1 = 0.017, corresponding to approximately 560 µs. Using the measurements of EGRET, Kaaret [155] determined the phase delay of P1 for different energy bins and found a lower limit on the linear mass scale of quantum gravity of EQG1 > 1.8 · 1015 GeV.1 Based on combined measurements from EGRET and MAGIC, a new lower limit on the quantum gravity effects can be derived, using the arrival phases of P1 and P2 in the 1

Note that it is assumed that there is not intrinsic energy dependent phase delay of the emission of pulsed photons at different energies.

7. Summary, Conclusion and Outlook

pulse . P1 P2

∆φ mu.l. = ∆E [1/GeV] 0.0019 0.0004

∆t mu.l. = ∆E [µs/GeV] 63 14

145

El.b. (Q.G.) [GeV] 0.4 · 1016 1.6 · 1016

Table 7.1: The first and the second column denote the 95% u.l. on the slope of the fitted linear function in figure 6.32, in terms of [phase]/GeV and µs/GeV, respectively, given the pulsar rotation period. The third column denotes the 95% lower boundary on the quantum gravity mass scale, assuming a linear vacuum dispersion: EQ.G. = Lc ∆E ∆t , where c denotes the speed of light, L = 2.2 kpc the distance ∆t to the pulsar and m = ∆E the upper limit on the slope, given in the second column. dφ

(E)

different energy bins. The slope mφ = Pj does not significantly differ from 0. Thus, dE its 2σ upper limit can be derived. With this upper limit and the pulsar rotation period dT (E) = T0 mφ , where T0 = 0.033 s, the time delay of each pulse can be computed. mT = Pj dE TPj (j = 1, 2) denotes the arrival time of Pj relative to the radio main pulse arrival time. From equation 7.5, the following (linear) QG energy scale can be derived:   L dTPj (E) −1 L EQG1 = = (7.6) c dE mT c In table 7.1., the corresponding upper limits on mφ and mT and the lower boundary on the QG mass scale are displayed. The highest upper limit is derived for P2 and is: EQG1 > 1.6 · 1016 GeV

(7.7)

In a second approach, I used only the optical and γ-ray pulse profiles measured with MAGIC. This prevents the introduction of a possible systematic error due to the different radio ephemeris used in the EGRET and MAGIC measurements. To counteract a possible intrinsic time delay between the emission of γ-rays and optical photons, the upper limit P2 ) pulse arrival phase was computed. The energy on the delay of the averaged ( TP1 +T 2 difference between the photons measured in optical and in γ-ray energies by MAGIC is ≈ 33 GeV.2 The averaged arrival phase in visible wavelengths is φoptical = 0.193 ± 0.003 and at 33 GeV φγ = 0.20 ± 0.01, yielding a phase difference of ∆φ = 0.010 ± 0.011, not significantly different from 0. The lower limit on EQG1 becomes: EQG1 > 0.8 · 1016 GeV

(7.8)

The lower limits deduced here are up to one order of magnitude better than the best limit so far deduced from pulsars [155] and at the same order of magnitude as the limits from γ-ray bursts. By increasing the amount of data in a future measurement of the Crab pulsar and thus increasing the phase resolution of the VHE γ-ray pulse profile, the lower boundary might be improved by another factor 2-5. 7.1.2.4 The P2/P1 Ratio: Comparison with the Predictions of a Multicomponent Model In figure 7.4, our results on the energy dependent P2/P1 ratio are compared to predictions from a multicomponent model fit to the phase resolved optical- γ-ray spectrum [156]. The 2

This is the mean energy of the γ-ray photons, assuming that their spectrum follows a power law with spectral index −3.1. The optical photons have an energy of 10−9 GeV.

146

7. Summary, Conclusion and Outlook

Figure 7.4: The multicomponent model prediction for the P2/P1 ratio. The data points are from [2] (black) and from MAGIC (red). The colored lines denote different model fits. multicomponent model assumes that the optical - X-ray Crab spectrum follows two different log parabolic shapes. Motivated by the Compton scattering mechanism, it is assumed in [156] that the γ-ray spectrum follows the same spectral curvature as the optical - X-ray spectrum but is shifted by a constant energy. The peaks of the γ-ray log parabolic spectra are at 300 MeV and 2 GeV. By assuming a cut-off energy at 15 GeV, a prediction for the HE γ-ray emission is obtained. Our results do not allow to constrain the multicomponent model predictions.

7.2 Outlook The discovery of the Crab Pulsar in VHE γ-radiation is one of the “Holy Grails of ground based γ-Astronomy”, as it was called by T. Weekes in the opening session of the GAMMA conference 2008 in Heidelberg. It led to the conclusion that for the Crab Pulsar the polar cap region above the magnetic polar region is excluded as an emission region of VHE γ-rays. Instead, the emission region is constrained to the outer magnetosphere – an important achievement in the research of pulsar physics. The FGST experiment, launched in summer 2008, aims to measure γ-rays between 20 MeV to 300 GeV,3 thus providing for the first time an energetic overlap of VHE γray observations performed by satellite and ground based experiments. While the energy resolution of FGST was tested and calibrated by MC and a test beam [157], flux measurements of γ-ray sources can be used to cross-calibrate the energy scale of Cherenkov 3

Note that the FGST instrument is due to its small detector area of ≈ 1 m2 limited to detect the small pulsed flux of Crab above an energy of a few tens of GeV. Using our flux measurements, and assuming that FGST has a duty cycle of 100%, it will measure in one year around 30 photons with an energy above 50 GeV. In order to measure the same amount of photons with an energy above 100 GeV, it has to measure at least 4 years.

7. Summary, Conclusion and Outlook

147

Telescopes like MAGIC [158], where the energy calibration relies solely on MC simulations. This will reduce MAGIC’s systematic error on the energy scale. Due to its steep spectral turnover, the flux measurement from a source like the Crab Pulsar is an ideal opportunity to perform this task. The results obtained for the cut-off energies of the total pulsed flux and for the phase resolved flux are statistically very solid, but might suffer from a systematic mismatch between the energy scale of the EGRET experiment and MAGIC. A spectrum measured by FGST overlapping in energy with MAGIC could be used to tune the energy scale of the MAGIC telescope and thus reduce this mismatch. As stated in chapter 6 and here, I found a strong hint on the 3.4σ level that the Crab Pulsar emits pulsed radiation above an energy of 60 GeV. With an observation duration of around 50 hours, this signal should reach the detection limit of 5σ, thus allowing to deduce highly significant spectrum between 25 GeV and 100 GeV. The measurements of non-pulsed γ-rays at 50 GeV from the Crab Nebula encourages the use of the sum trigger also for non-pulsating sources. In particular, the observation of AGN with steep source spectra could be ideal targets. More data will also help to determine the P2/P1 ratio as a function of energy and time. This ratio may vary on a cycle of 14 years [159]. Furthermore, a small hint at 1.4σ for an additional peak in the off pulse region was found. More data will show whether this a statistical fluctuation or a real signal. The phase resolution will also be improved with more data, allowing more precise studies on the pulse profile morphology. In the first half year after its launch, FGST has already found pulsating signals from a dozen of pulsars [160]. Among those there might be pulsars having an even higher cut off energy and a harder spectrum than the one we measured from the Crab Pulsar. Measurements of MAGIC and FGST together could then determine the cut-off spectrum on a very precise level.

A Gamma/Hadron Separation Using Shower Autocorrelations During my thesis, I tried out several methods to improve the Gamma/Hadron separation, which in the end was not used for the analysis of the Crab Pulsar. One of these methods, making use of the autocorrelation function of the shower image, is presented here. The images of γ-ray and hadronic background showers differ in spatial and time distribution. Hadronic showers are usually more spread out in space and time than γ-rays and they often have islands, originating from subshowers initiated by e.g. pions. Identifying the images from these subshowers over the NSB light is a challenging task.

A.1 General Idea Assuming that the NSB causes Poissonian signal fluctuations in a PMT, one wants to distinguish these signals from the ones caused by the shower itself. A robust method to do so is the image cleaning algorithm, described above in section 4.4.2. This algorithm, however, removes not only NSB fluctuations, but also the small signals coming from subshowers. Assuming that these subshowers are correlated in time with the main shower, I used the autocorrelation function Af (x, y, t) to identify those. For a given function f (x, y, t), its autocorrelation function Af (x, y, t) is defined as: Z Z Z Af (x, y, t) = f (x + x′ , y + y ′ , t + t′ ) · f (x′ , y ′ , t′ ) dx′ dy ′ dt′ (A.1)

The function f is thus used as a filter for itself. An important property of the autocorrelation function is that it follows a δ-distribution if f (x, y, t) is Poissonian or Gaussian noise: Af, noise (x, y, t) = δ(x, y, t) (A.2) displayed in 1D in figure A.1. If f (τ ) consists of noise and two signal peaks, the autocorrelation function can be used as an efficient noise filter. This is displayed in figure A.2. The function f (x, y, t) corresponds to the FADC value in each slice (time t) and pixel (spatial coordinate (x, y)). With the main axis (LENGTH) of the shower as abscissa and the minor axis (WIDTH) as the ordinate, I introduced cartesian coordinates along the shower main axis. The √ P ′ ′ fFADC (x′ , y ′ , t) · e− (x−x )+(y−y )/σ fcartesian (x, y, t) =

NN(x′ ,y ′ )

P

e−

√

(x−x′ )+(y−y ′ )/σ

(A.3)

NN(x′ ,y ′ )

fFADC (x′ , y ′ , t) is the time and charge calibrated FADC value of the pixel at the position (x′ , y ′ ). The autocorrelation function can be computed much faster (O(2n log n))1 if one 1

n denotes the number of sampling points of f (x, y, t). Using the the direct computation in equation A.1 needs O(n2 /2) computational steps

150

A. Gamma/Hadron Separation Using Shower Autocorrelations

A f ( τ)

f(t) 4

1

3 0.8

2 1

0.6 0 -1

0.4

-2 0.2

-3 -4

0

-5 0

2000

4000

6000

8000

10000

-4000

-2000

0

2000

4000 τ

t

(a) f (τ ): Gaussian Noise

(b) Af (τ )

Figure A.1: Figure (a) shows an example of Gaussian noise f (τ ) and figure (b) its autocorrelation function Af (τ ). A f ( τ)

f(t)

1.18

5 4

1.16 3 2

1.14

1 0

1.12

-1 1.1

-2 -3

1.08

0

2000

4000

6000

8000

10000

0

1000

2000

3000

4000

5000 τ

t

(a) f (τ ): Noise and two signals

(b) Af (τ )

Figure A.2: The function f (τ ) consists of Gaussian noise and two hardly visible peaks at the positions τ = 5000 and τ = 7500. Its autocorrelation function shows a clear peak at τ = 0 at τ = 2500, coming from the hardly visible two peaks in f (τ ). Because the shape of the two signal peaks is similar, the autocorrelation function filters the noise. Entries Mean

Entries 0 Mean 0.1612 RMS 1.194

0 -0.0004181

RMS

0.7906

11.1 10.2 9.3 8.4 7.5 6.6 5.7 4.8 3.9 3.0 2.1 1.2 0.4 -0.5 -1.4 -2.3 -3.2

0.60° 189mm

0.60° 189mm

(a) t = 0 slices

0.2

ShowerRect Entries 200 Mean x 0.08309 Mean y 0.003736 RMS x 0.3264 8 0.1673 RMS y

0.1

6

0

4

0.3

-0.1

2

-0.2

0

-0.3 -0.6

-0.4

-0.2

0

0.2

0.4

0.6

(g) t = 0 slices

Entries Mean

200 0.06594 -0.003022 0.2829 8 0.1364

0.1

6

0

4

-0.1

2

-0.2

0

-0.3 -0.6

-0.4

-0.2

0

0.2

0.4

0.6

(h) t = 1 slices

RMS

0.2

ShowerRect Entries 200 Mean x 0.0779 Mean y 0.01089 RMS x 0.2911 8 0.1381 RMS y

0.1

6

0

4

0.3

-0.1

2

-0.2

0

-0.3 -0.6

-0.4

-0.2

0

0.2

0.4

0.6

(i) t = 2 slices

0

Entries

-0.01284

Mean

0.7487

0.2

ShowerRect Entries 200 Mean x 0.07052 Mean y 0.02799 RMS x 0.3352 8 0.152 RMS y

0.1

6

0

4

-0.1

2

-0.2

0

-0.3 -0.6

-0.4

-0.2

0

0.2

0.4

0.6

(j) t = 3 slices

Entries Mean

0.7625

0.2

ShowerRect Entries 200 Mean x 0.07296 Mean y 0.02859 RMS x 0.3371 8 0.1616 RMS y

0.1

6

0

4

-0.1

2

-0.2

0

-0.3 -0.6

-0.4

-0.2

0

0.2

0.4

0.6

(k) t = 4 slices

0.7731

11.1 10.2 9.3 8.4 7.5 6.6 5.7 4.8 3.9 3.0 2.1 1.2 0.4 -0.5 -1.4 -2.3 -3.2

0.60° 189mm

(e) t = 4 slices 0.3

0 0.002419

RMS

11.1 10.2 9.3 8.4 7.5 6.6 5.7 4.8 3.9 3.0 2.1 1.2 0.4 -0.5 -1.4 -2.3 -3.2

0.60° 189mm

(d) t = 3 slices 0.3

0 -0.05848

RMS

11.1 10.2 9.3 8.4 7.5 6.6 5.7 4.8 3.9 3.0 2.1 1.2 0.4 -0.5 -1.4 -2.3 -3.2

0.60° 189mm

(c) t = 2 slices

ShowerRect Entries Mean x Mean y RMS x RMS y

0.2

Entries Mean

0.9166

11.1 10.2 9.3 8.4 7.5 6.6 5.7 4.8 3.9 3.0 2.1 1.2 0.4 -0.5 -1.4 -2.3 -3.2

0.60° 189mm

(b) t = 1 slices 0.3

0 0.07037

RMS

11.1 10.2 9.3 8.4 7.5 6.6 5.7 4.8 3.9 3.0 2.1 1.2 0.4 -0.5 -1.4 -2.3 -3.2

(f) t = 5 slices ShowerRect Entries Mean x Mean y RMS x RMS y

0.3

0.2

200 0.02835 -0.007839 0.3491 8 0.1673

0.1

6

0

4

-0.1

2

-0.2

0

-0.3 -0.6

-0.4

-0.2

0

0.2

0.4

0.6

(l) t = 5 slices

Figure A.3: A shower in time slices in Camera and Rectangular coordinates.

A. Gamma/Hadron Separation Using Shower Autocorrelations

CorrRect Entries 200 Mean x 0.001574 Mean y 0.0007087 RMS x 0.3409 RMS y 0.1709

0.03 0.028 0.026 0.024 0.022 0.02 0.018 0.016 0.014 0.012 0.3

CorrRect Entries 200 Mean x -0.0006977 Mean y 0.0003805 RMS x 0.3442 RMS y 0.1719

0.03 0.028 0.026 0.024 0.022 0.02 0.018 0.016 0.014 0.012 0.3

0.2

CorrRect

0 -0.1 -0.2 -0.3 -0.6

-0.4

-0.2

0

0.2

0.4

0 -0.1 -0.2 -0.3 -0.6

(a) t = 0 slices

-0.4

-0.2

0

0.2

0.4

-0.1 -0.3 -0.6

0 -0.1 -0.2 -0.3 -0.6

-0.4

-0.2

0

0.2

0.4

-0.4

-0.2

0

0.2

0.6

0 -0.1 -0.2 -0.3 -0.6

(e) t = 4 slices

-0.4

-0.2

0

0.2

0.4

-0.3 -0.6

0.6

-0.4

-0.2

0

0.2

0.4

0.6

(d) t = 3 slices

CorrRect Entries 200 Mean x -0.0009884 Mean y 0.0003638 RMS x 0.3453 RMS y 0.1723

0.03 0.028 0.026 0.024 0.022 0.02 0.018 0.016 0.014 0.012 0.3

CorrRect Entries Mean x Mean y RMS x RMS y

0.03 0.028 0.026 0.024 0.022 0.02 0.018 0.016 0.014 0.012 0.3

200 -0.002419 0.000912 0.3458 0.1723

0.2 0.1 0 -0.1 -0.2 -0.3 -0.6

(f) t = 5 slices

0 -0.1

0.2 0.1

0.1

0.6

-0.2

CorrRect Entries 200 Mean x -8.287e-05 Mean y 0.0001425 RMS x 0.3461 RMS y 0.1723

0.03 0.028 0.026 0.024 0.022 0.02 0.018 0.016 0.014 0.012 0.3

0.4

(c) t = 2 slices

0.2 0.1

0 -0.2

CorrRect

0.2

CorrRect Entries 200 Mean x -0.001297 Mean y 0.0001663 RMS x 0.3458 RMS y 0.1725

0.03 0.028 0.026 0.024 0.022 0.02 0.018 0.016 0.014 0.012 0.3 0.2

0.1

0.6

(b) t = 1 slices

Entries 200 Mean x -0.001375 Mean y 0.0001163 RMS x 0.3462 RMS y 0.1721

0.03 0.028 0.026 0.024 0.022 0.02 0.018 0.016 0.014 0.012 0.3

200 -0.001821 3.065e-05 0.3451 0.1723

0.2 0.1

0.6

Entries Mean x Mean y RMS x RMS y

0.03 0.028 0.026 0.024 0.022 0.02 0.018 0.016 0.014 0.012 0.3

0.2 0.1

151

-0.4

-0.2

0

0.2

0.4

(g) t = 6 slices

0.6

0.1 0 -0.1 -0.2 -0.3 -0.6

-0.4

-0.2

0

0.2

0.4

0.6

(h) t = 7 slices

Figure A.4: The autocorrelation function of the shower displayed in figure A.3. It peaks at the time slice 0. In this example, after t = 4 slices, the autocorrelation function at (x, py) = (0, 0) is reduced to 1/e of its value at t = 0. In the spatial direction r = x2 + y 2 , Af (r, t = 0) is reduced to 1/e at a distance of r = 0.14◦ of the value at r = 0. uses the fourier transform: Af (x, y, t) =

Z Z Z

S(kx , ky , ω) · eiωt eikx t eiky t dkx dky dω

(A.4)

where Z Z Z 2 1 −iωt −ikx t −iky t S(kx , ky , ω) = f (x, y, t) · e e e dx dy dt 2π

(A.5)

denotes the power spectrum of f (x, y, t). The fourier transform was computed with the FFTW algorithm [161]. For each event, its autocorrelation function is computed. For the event displayed in figure A.3, its autocorrelation function is displayed in figure A.4. As in the 1D examples above, the autocorrelation function Af (x, y, t) peaks at (x = 0, y = 0, t = 0) due to its self similarity.

A.2 Conclusion In an additional step, the autocorrelation function is parameterized. Its spatial extension in x and y direction is characterised by the variance σx2 and σy2 , and its temporal extension by σt2 , derived for each event. Comparing the distribution of these parameters for MC γ-rays and hadronic background reveals slight differences and allows to reduce the background by 5%, additionally to the cut applied by the random forest, see figures A.5 and A.6. Unfortunately, the computing cost is huge, even though FFTW is used, and not justified by the additional 5% gain in background subtraction.

152

A. Gamma/Hadron Separation Using Shower Autocorrelations

(a) background

(b) MC γ-rays

q Figure A.5: The spatial variance σr = σx2 + σy2 of the autocorrelation function Af (x, y, t). The unit on the x-axis is arbitrary, the y-axis denotes the number of counts per bin. The σr for background reaches higher values.

(a) background

(b) MC γ-rays

Figure A.6: The variance in time σt of the autocorrelation function, measured at (x, y) = (0, 0). The unit on the x-axis is arbitrary, the y-axis denotes the number of counts per bin. The autocorrelation of hadronic background is in average larger than for γ-rays.

A.3 Future Development: Implementation of the Algorithm using a GPU and CUDA An idea to shorten the computing time for this algorithm is to use a graphics processor unit (GPU) for the fast fourier transform and for the computation of the shower image in rectangular coordinates. First tests I performed with a NVIDIA GeForce GTX 280 GPU using the CUDA2 environment, showed an acceleration of at least a factor 100 for the computation of the rectangular coordinates. Further studies are in process.

2

Compute Unified Device Architecture, see http://www.nvidia.com/object/cuda_get.html

B List of Acronyms and Abbreviations

AC AGILE AGN AMC CCD CGRS COG COMPTEL CUDA DAQ DC

EGRET eV FADC FGST FPGA FOV G-APD GeV GJ GPU H.E.S.S. HE HPD IACT IPRC LE LED LT0 LT1

“Alternate Current”, is also used for capacitive coupling of a signal Astrorivelatore Gamma ad Immagini ultra LEggero Active Galactic Nuclei Active Mirror Control Charge Coupled Device Compton Gamma Ray Satellite Center Of Gravity COMpton TELescope Compute Unified Device Architecture Data Acquisition Queue “Direct Current”. Used in this thesis for a constant or steady signal, e.g. the expected DC γ-ray flux from the Crab Nebula compared to the pulsed flux from the Crab Pulsar Energetic Gamma Ray Experiment Telescope electron Volt, 1eV = 1.6 · 10−12 J Flash Analog Digital Converter Fermi Gamma-ray Space Telescope Field Programmable Gate Array Field Of View Geiger Mode Avalanche Photodiode Giga electron Volt 1 GeV = 109 eV Goldreich-Julian (density) Graphics Processor Unit High Energy Stereoscopic System High Energy Hybrid Photo Detector Imaging Air Cherenkov Telescope Individual Pixel Rate Control Low Energy Light Emitting Diode Level 0 Trigger Level 1 Trigger

154

LT2 LUT MAGIC MC MeV MJD MUX ndf NN NSB phe PMT PSF PULSAR QE QG RF rms SCLF SED STACEE TeV UV VCSEL VERITAS VHE

B. List of Acronyms and Abbreviations

Level 2 Trigger Look Up Table Major Atmospheric Gamma-ray Imaging Cherenkov Telescope Monte Carlo Mega electron Volt 1 MeV = 106 eV Mean Julian Day MUltipleXed number of degrees of freedom Next Neighbor Night Sky Background photoelectron Photomultiplier Tube ("Pixel" and "PMT" are synonyms in this work) Point Spread Function PULSating StAR Quantum Efficiency Quantum Gravity Random Forest root mean square Space Charge Limited Flow Spectral Energy Density Solar Tower Atmospheric Cerenkov Effect Experiment Terra electron Volt 1 TeV = 1012 eV UltraViolet Vertical-Cavity Surface-Emitting Laser Very Energetic Radiation Imaging Telescope Array System Very High Energy

List of Figures 0.1

The Crab Nebula as seen in radio, infrared and optical wavelengths, in Xand γ-rays . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

15

1.1 1.2

The cosmic ray spectrum . . . . . . . . . . . . . . . . . . . . . . . . . . . . Illustration of an extended air shower initiated by a cosmic ray particle . . .

18 19

2.1 2.2 2.3 2.4 2.5

A selection of scientific targets of very high energy γ-ray astronomy The location of the sources detected with EGRET . . . . . . . . . The VHE γ-ray sources for energies above 100 GeV . . . . . . . . . The Kifune Plot . . . . . . . . . . . . . . . . . . . . . . . . . . . . Different interaction mechanism between γ-rays and leptons . . . .

. . . . .

21 22 23 23 26

The discovery plot of the first pulsar in 1967 by J. Bell . . . . . . . . . . . . The locations of the vacuum gaps within the pulsar’s magnetosphere . . . . The size of a canonical pulsar . . . . . . . . . . . . . . . . . . . . . . . . . . The pulsar period P and its derivative P˙ for the known radio pulsars . . . . The pulse profiles of the seven pulsars detected with EGRET . . . . . . . . The pulse profile of the Crab Pulsar measured with the EGRET experiment for energies above 100 MeV . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.7 The phase resolved spectrum measurements of the Crab Pulsar . . . . . . . 3.8 Slot gap model flux predictions for the Crab Pulsar . . . . . . . . . . . . . . 3.9 Outer gap model flux predictions for a young pulsar . . . . . . . . . . . . . 3.10 Outer gap model flux predictions for the Crab Pulsar . . . . . . . . . . . . .

30 31 32 33 34 34 35 40 41 42

4.1 4.2 4.3 4.4 4.5 4.6 4.7 4.8 4.9 4.10 4.11 4.12 4.13 4.14 4.15 4.16 4.17 4.18 4.19

45 47 49 50 51 52 54 56 56 57 59 62 64 64 64 65 65 67 69

. . . . .

. . . . .

. . . . .

. . . . .

3.1 3.2 3.3 3.4 3.5 3.6

The MAGIC telescope during observations. . . . . . . . . . . . . . . . . . . Pair creation technique used in satellite γ-experiments . . . . . . . . . . . . Schematics of an electromagnetic and a hadronic air shower. . . . . . . . . . The traces of the particles in air showers triggered by γ-rays and protons . . Propagation of Cherenkov light . . . . . . . . . . . . . . . . . . . . . . . . . Schematics of the imaging technique . . . . . . . . . . . . . . . . . . . . . . The layout of the MAGIC camera . . . . . . . . . . . . . . . . . . . . . . . . The trigger macro cells . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . Readout scheme of the MAGIC telescope . . . . . . . . . . . . . . . . . . . . Calibration system of the MAGIC telescope . . . . . . . . . . . . . . . . . . Schematics of the analysis chain . . . . . . . . . . . . . . . . . . . . . . . . . NN groups used in the sum image cleaning . . . . . . . . . . . . . . . . . . . An image of a MC γ-ray event with energy 33 GeV . . . . . . . . . . . . . . An image of a MC γ-ray event with energy 35 GeV . . . . . . . . . . . . . . An image of a recorded event . . . . . . . . . . . . . . . . . . . . . . . . . . An image of a recorded event . . . . . . . . . . . . . . . . . . . . . . . . . . The energy distribution of the sum and the standard image cleaning . . . . Definition of the image parameters . . . . . . . . . . . . . . . . . . . . . . . Determination of the number of excess events from the ALPHA distribution

156

List of Figures

4.20 The phase diagram of the Crab Pulsar measured by EGRET for energies above 100 MeV . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 5.1 5.2

The installed sum trigger in the MAGIC counting house . . . . . . . . . . . The distribution of Cherenkov photons of simulated γ-rays in the Camera of the MAGIC telescope . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 5.3 The radial averaged distribution of Cherenkov photons . . . . . . . . . . . . 5.4 Strong signals are clipped before summed up . . . . . . . . . . . . . . . . . 5.5 The count rate above a discriminator threshold for one MAGIC PMT exposed to NSB . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 5.6 The final configuration of the sum trigger . . . . . . . . . . . . . . . . . . . 5.7 Schematics of the three sum trigger boards . . . . . . . . . . . . . . . . . . . 5.8 The measured and simulated NSB spectra . . . . . . . . . . . . . . . . . . . 5.9 The distribution of SIZE in phe for simulated and measured protons . . . . 5.10 The distributions of the shower centers for sum trigger and standard trigger 5.11 The differential rate for MC γ-rays for the standard trigger and the sum trigger . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 5.12 The collection area for the MAGIC standard trigger and the sum trigger . . 6.1 6.2 6.3 6.4 6.5 6.6 6.7 6.8 6.9 6.10 6.11 6.12 6.13 6.14 6.15 6.16 6.17 6.18 6.19 6.20 6.21 6.22 6.23 6.24 6.25

Measurement of pulsed VHE γ-rays from Crab with MAGIC using data from 2005 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . EGRET spectrum measurements from Crab and flux upper limits deduced with ground based instruments . . . . . . . . . . . . . . . . . . . . . . . . . The distribution of cloudiness for the whole observation period . . . . . . . The cloudiness as a function of time . . . . . . . . . . . . . . . . . . . . . . The distribution of the event rates after cleaning . . . . . . . . . . . . . . . The sensitivity for γ-rays from the Crab Nebula . . . . . . . . . . . . . . . . The track of Crab along the sky . . . . . . . . . . . . . . . . . . . . . . . . . The SIZE dependent cut in ALPHA . . . . . . . . . . . . . . . . . . . . . . The optical lightcurve of the Crab pulsar . . . . . . . . . . . . . . . . . . . . The pulse profile for γ-rays above 25 GeV . . . . . . . . . . . . . . . . . . . Pulse profile for the inverted cut in ALPHA . . . . . . . . . . . . . . . . . . The pulse profiles in four SIZE bins . . . . . . . . . . . . . . . . . . . . . . . Comparison of measured and simulated γ-rays (LE) . . . . . . . . . . . . . . Comparison of measured and simulated γ-rays (HE) . . . . . . . . . . . . . The number of excess events in P1 and P2 normalized to the total number of excess events in each pulse . . . . . . . . . . . . . . . . . . . . . . . . . . The probability of the fit to the measured data EGRET, COMPTEL and MAGIC data . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . Measured and predicted number of excess events for three different models . The joint fit to the MAGIC, EGRET and COMPTEL data . . . . . . . . . The number of measured excess events the prediction based on the power law cut-off . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . The measured and predicted number of excess events in P1 and P2 . . . . The probability map for the spectral fit to the single pulses . . . . . . . . . The number of measured excess events in the main and inter pulse and the prediction from a power law fit to the MAGIC data . . . . . . . . . . . . . . The correlation between generated and reconstructed energy . . . . . . . . . The bias of the reconstructed energy . . . . . . . . . . . . . . . . . . . . . . The energy resolution at low energies . . . . . . . . . . . . . . . . . . . . . .

73 75 78 79 80 82 84 85 87 88 89 90 90 94 95 96 97 98 98 99 101 103 103 104 107 108 109 111 111 112 113 114 116 117 119 120 121 121

List of Figures

157

6.26 The differential total pulsed spectrum spectrum in P1+P2 before unfolding 6.27 The unfolded spectrum of the total pulsed flux . . . . . . . . . . . . . . . . 6.28 The predicted and measured number of excess events for P1 and P2 assuming a power law spectrum . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 6.29 The single pulse spectra of P1 and P2 before unfolding . . . . . . . . . . . . 6.30 The unfolded single pulse spectra . . . . . . . . . . . . . . . . . . . . . . . . 6.31 The measured pulse profiles from EGRET and MAGIC for different energies 6.32 The arrival phases of the main and inter pulse . . . . . . . . . . . . . . . . . 6.33 The arrival phase difference between the main and inter pulse as a function of energy . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 6.34 The pulse profiles above 25 GeV and 60 GeV, derived with a Gaussian Kernel Density Estimator . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 6.35 The ratio of P2 to P1 for four energy bins . . . . . . . . . . . . . . . . . . . 6.36 The ALPHA plots for estimated energies between 50 GeV and 80 GeV . . . . 6.37 The ALPHA plots for estimated energies between 30 GeV and 50 GeV . . . . total pulsed flux for energies between 1 MeV and 200 GeV pulsed spectrum in the main pulse P1 . . . . . . . . . . . . pulsed spectrum in the inter pulse P2 . . . . . . . . . . . . multicomponent model prediction for the P2/P1 ratio . . .

. . . .

. . . .

122 123 123 125 126 128 130 131 132 134 136 136

7.1 7.2 7.3 7.4

The The The The

. . . .

. . . .

. . . .

. . . .

. . . .

140 141 142 146

A.1 A.2 A.3 A.4 A.5 A.6

Autocorrelation function of Gaussian noise . . . . . . . . . . . . . . . Autocorrelation function of two peaks distorted with Gaussian noise A shower in time slices in Camera and Cartesian coordinates . . . . The autocorrelation function of the shower above . . . . . . . . . . . The spatial variance of the autocorrelation function . . . . . . . . . . The variance in time of the autocorrelation function . . . . . . . . .

. . . . . .

. . . . . .

. . . . . .

. . . . . .

150 150 150 151 152 152

Bibliography [1] T. C. Weekes et al. Observation of TeV gamma rays from the Crab nebula using the atmospheric Cerenkov imaging technique. ApJ, 342:379–395, July 1989. [2] L. Kuiper, W. Hermsen, G. Cusumano, R. Diehl, V. Schönfelder, A. Strong, K. Bennett, and M. L. McConnell. The Crab pulsar in the 0.75-30 MeV range as seen by CGRO COMPTEL. A coherent high-energy picture from soft X-rays up to highenergy gamma-rays. A&A, 378:918–935, November 2001. [3] J. Albert et al. VHE γ-Ray Observation of the Crab Nebula and its Pulsar with the MAGIC Telescope. ApJ, 674:1037–1055, February 2008. [4] CHANDRA. Chandra discovers x-ray ring around cosmic powerhouse in crab nebula. http://chandra.harvard.edu/photo/1999/0052/. [5] S. Swordy. Cosmic ray spectrum. http://astroparticle.uchicago.edu/cosmic_ray_spectrum_picture.htm. [6] J. Abraham et al. Correlation of the Highest-Energy Cosmic Rays with Nearby Extragalactic Objects. Science, 318:938, November 2007. [7] F. Aharonian et al. A New Population of Very High Energy Gamma-Ray Sources in the Milky Way. Science, 307:1938–1942, March 2005. [8] M. Pohl. Gamma ray astronomy. Proc. of the 27th ICRC, 2001. [9] E. Fermi. On the Origin of the Cosmic Radiation. Phys. Rev., 75:1169 – 1174, 1949. [10] W. Bednarek et al. Pulsar contribution to the cosmic rays between the knee and the ankle. Proc. of the 29th ICRC, Pune, India, 2005. [11] I. Büsching, C. Venter, and O. C. de Jager. Contributions from nearby pulsars to the local cosmic ray electron spectrum. Advances in Space Research, 42:497–503, August 2008. [12] D. Hooper, P. Blasi, and P. Dario Serpico. Pulsars as the sources of high energy cosmic ray positrons. Journal of Cosmology and Astro-Particle Physics, 1:25, January 2009. [13] R. C. Hartman et al. The Third EGRET Catalog of High-Energy Gamma-Ray Sources. ApJS, 123:79–202, July 1999. [14] R. Wagner. VHE γ-ray Sky Map and Source Catalog. http://www.mpp.mpg.de/ rwagner/sources/index.html. [15] J. Hinton. Gamma-ray Astronomy. astro-ph/0712.3352, December 2007. [16] A. Reimer, M. Joshi, and M. Böttcher. The Blazar 3C 66A in 2003-2004: hadronic versus leptonic model fits. American Institute of Physics Conference Series, 2008.

160

Bibliography

[17] N. Gehrels et al. The Swift Gamma-Ray Burst Mission. ApJ, 611:1005–1020, August 2004. [18] F. Aharonian et al. 3.9 day orbital modulation in the TeV γ-ray flux and spectrum from the X-ray binary LS 5039. A&A, 460:743–749, December 2006. [19] J. Albert et al. Variable Very-High-Energy Gamma-Ray Emission from the Microquasar LS I +61 303. Science, 312:1771–1773, June 2006. [20] A. F. Iyudin, R. Diehl, H. Bloemen, W. Hermsen, G. G. Lichti, D. Morris, J. Ryan, V. Schoenfelder, H. Steinle, M. Varendorff, C. de Vries, and C. Winkler. COMPTEL observations of Ti-44 gamma-ray line emission from CAS A. A&A, 284:L1–L4, April 1994. [21] B. J. Teegarden et al. INTEGRAL SPI Limits on Electron-Positron Annihilation Radiation from the Galactic Plane. ApJ, 621:296–300, March 2005. [22] F. A. Aharonian. Very high energy cosmic gamma radiation : a crucial window on the extreme Universe. World Scientific Publishing, 2004. [23] S. L. Adler, J. N. Bahcall, C. G. Callan, and M. N. Rosenbluth. Photon splitting in a strong magnetic field. Physical Review Letters, 25:1061–1065, 1970. [24] M. G. Baring and A. K. Harding. Photon Splitting and Pair Creation in Highly Magnetized Pulsars. ApJ, 547:929–948, February 2001. [25] Z. Bialynicka-Birula and I. Bialynicki-Birula. Nonlinear Effects in Quantum Electrodynamics. Photon Propagation and Photon Splitting in an External Field. Phys. Rev. D, 2:2341–2345, November 1970. [26] A. A. Kozlenkov and I. G. Mitrofanov. Two-photon generation of electron-positron pairs in a strong magnetic field. Journal of Experimental and Theoretical Physics, 91:1978–1989, December 1986. [27] A. E. Shabad and V. V. Usov. Gamma-quanta capture by magnetic field and pair creation suppression in pulsars. Nature, 295:215–217, January 1982. [28] L. Woltjer. The Stability of Force-Free Magnetic Fields. ApJ, 128:384, September 1958. [29] A. K. Harding and A. G. Muslimov. Pulsar Polar Cap Heating and Surface Thermal X-Ray Emission. I. Curvature Radiation Pair Fronts. ApJ, 556:987–1001, August 2001. [30] Q. Luo and D. B. Melrose. Curvature emission and absorption - Single particle treatment. Proceedings of the Astronomical Society of Australia, 10:45–47, 1992. [31] I. P. Ochelkov and V. V. Usov. Curvature radiation of relativistic particles in the magnetosphere of pulsars. I - Theory. Ap&SS, 69:439–460, May 1980. [32] A. K. Harding, J. V. Stern, J. Dyks, and M. Frackowiak. High-Altitude Emission from Pulsar Slot Gaps: The Crab Pulsar. ApJ, 680:1378–1393, June 2008. [33] R. W. Romani. Gamma-Ray Pulsars: Radiation Processes in the Outer Magnetosphere. ApJ, 470:469, October 1996.

Bibliography

161

[34] A. K. Harding and D. Lai. Physics of strongly magnetized neutron stars. Reports on Progress in Physics, 69:2631–2708, 2006. [35] A. Hewish, S. J. Bell, J. D. H. Pilkington, P. F. Scott, and R. A. Collins. Observation of a Rapidly Pulsating Radio Source. Nature, 224:472, November 1969. [36] J. M. Comella, H. D. Craft, Jr., R. V. E. Lovelace, J. M. Sutton, and G. L. Tyler. Crab Nebula Pulsar NP 0532. Nature, 221:453, 1969. [37] F. Pacini. Energy Emission from a Neutron Star. Nature, 216:567, November 1967. [38] F. Hoyle, J. V. Narlikar, and J. A. Wheeler. Electromagnetic Waves from Very Dense Stars. Nature, 203:914, August 1964. [39] W.-Y. Chau. Gravitational Radiation from Neutron Stars. ApJ, 147:664, February 1967. [40] D. Lorimer. Handbook of pulsar astronomy . Cambridge University Press, 2004. [41] A. G. Lyne. Pulsar astronomy. Cambridge : Cambridge University Press, 2005. [42] P. Goldreich and W. H. Julian. Pulsar Electrodynamics. ApJ, 157:869, August 1969. [43] P. A. Sturrock. A Model of Pulsars. ApJ, 164:529, March 1971. [44] E. P. J. van den Heuvel. Pulsar Magnetospheres and Pulsar Death. Science, 312:539– 540, 2006. [45] D. J. Thompson. Gamma ray pulsars. Astrophysics and Space Science Library, 2004. [46] F. Aharonian et al. Search for pulsed VHE gamma-ray emission from young pulsars with HESS. A&A, 466:543–554, May 2007. [47] D. E. Alexandreas et al. A search for UHE cosmic rays from the Crab pulsar/nebula. Proc. of the 22nd ICRC, 1991. [48] O. Celik et al. Observations of the Crab Nebula and Pulsar with VERITAS. astroph/0709.386, September 2007. [49] M. de Naurois et al. Measurement of the Crab Flux above 60 GeV with the CELESTE Cerenkov Telescope. ApJ, 566:343–357, February 2002. [50] R. W. Lessard et al. Search for Pulsed TEV Gamma-Ray Emission from the Crab Pulsar. ApJ, 531:942–948, March 2000. [51] S. Oser et al. High-Energy Gamma-Ray Observations of the Crab Nebula and Pulsar with the Solar Tower Atmospheric Cerenkov Effect Experiment. ApJ, 547:949–958, February 2001. [52] J. D. Jackson. Classical Electrodynamics. 1998. [53] K. Hirotani. High Energy Emission from Rotation-Powered Pulsars: Outer-gap vs. Slot-gap Models. astro-ph/0809.1283, September 2008. [54] M. A. Livingstone, V. M. Kaspi, F. P. Gavriil, R. N. Manchester, E. V. G. Gotthelf, and L. Kuiper. New phase-coherent measurements of pulsar braking indices. Ap&SS, 308:317–323, April 2007.

162

Bibliography

[55] Y. E. Lyubarsky. Physics of pulsars. Amsterdam: Harwood Academic Publishers, |c1995, 1995. [56] W. Becker and J. Truemper. The X-ray luminosity of rotation-powered neutron stars. A&A, 326:682–691, October 1997. [57] F. Pacini. The Secular Decrease of Optical and X-Ray Luminosity of Pulsars. ApJ, 163:L17, January 1971. [58] V. S. Beskin, A. V. Gurevich, and I. N. Istomin. The electrodynamics of a pulsar magnetosphere. Journal of Experimental and Theoretical Physics, 85:401–433, August 1983. [59] M. G. Baring. High-energy emission from pulsars: the polar cap scenario. Advances in Space Research, 33:552–560, 2004. [60] J. Arons and E. T. Scharlemann. Pair formation above pulsar polar caps - Structure of the low altitude acceleration zone. ApJ, 231:854–879, August 1979. [61] J. K. Daugherty and A. K. Harding. Electromagnetic cascades in pulsars. ApJ, 252:337–347, January 1982. [62] M. A. Ruderman and P. G. Sutherland. Theory of pulsars - Polar caps, sparks, and coherent microwave radiation. ApJ, 196:51–72, February 1975. [63] A. K. Harding. Gamma Rays From Rotation-Powered Pulsars. astro-ph/0208421, August 2002. [64] K. Hirotani. Outer-Gap versus Slot-Gap Models for Pulsar High-Energy Emissions: The Case of the Crab Pulsar. ApJ, 688:L25–L28, November 2008. [65] J. Arons. Pair creation above pulsar polar caps - Geometrical structure and energetics of slot gaps. ApJ, 266:215–241, March 1983. [66] E. T. Scharlemann, J. Arons, and W. M. Fawley. Potential drops above pulsar polar caps - Ultrarelativistic particle acceleration along the curved magnetic field. ApJ, 222:297–316, May 1978. [67] A. G. Muslimov and A. K. Harding. Extended Acceleration in Slot Gaps and Pulsar High-Energy Emission. ApJ, 588:430–440, May 2003. [68] J. Dyks and B. Rudak. Two-Pole Caustic Model for High-Energy Light Curves of Pulsars. ApJ, 598:1201–1206, December 2003. [69] B. Rudak, J. Dyks, and T. Bulik. High-energy radiation from pulsars: A challenge to polar-cap models. In W. Becker, H. Lesch, and J. Trümper, editors, Neutron Stars, Pulsars, and Supernova Remnants, page 142, 2002. [70] K. S. Cheng, C. Ho, and M. Ruderman. Energetic radiation from rapidly spinning pulsars. I - Outer magnetosphere gaps. ApJ, 300:500, January 1986. [71] K. S. Cheng, C. Ho, and M. Ruderman. Energetic Radiation from Rapidly Spinning Pulsars. II. VELA and Crab. ApJ, 300:522, January 1986. [72] K. Hirotani. private communication. 2008.

Bibliography

163

[73] D. J. Thompson. Gamma ray astrophysics: the EGRET results. Reports on Progress in Physics, 71(11):116901, November 2008. [74] M. Tavani et al. The AGILE space mission. Nuclear Instruments and Methods in Physics Research A, 588:52–62, April 2008. [75] S. M. Ritz et al. Fermi (formerly GLAST) Mission Overview. American Astronomical Society Meeting, 2009 . [76] T. Reposeur, D. Dumora, L. Guillemot, M. Lemoine-Goumard, D. Parent, and D. Smith. Early Pulsar Observations and GLAST LAT Performance. The First GLAST Symposium, 2007. [77] P. A Cerenkov. Visible Emission of Clean Liquids by Action of γ Radiation. Doklady Akad. Nauk SSSR 2, page 451, 1934. [78] J. Galbraith and G. Jelley. Light Pulses from the Night Sky Associated with Cosmic Rays. Nature, February 1953. [79] S. Commichau. Observation of Very High Energy Gamma-Rays from the Galactic Center with the MAGIC Telescope. PhD thesis, March 2007. ETH-17118. [80] D. Hrupec. Extragalactic sources of rapidly variable high energy gamma radiation. PhD thesis, 2008. [81] Tibet As Gamma Collaboration. Diffuse Gamma-Ray Search from the Galactic Plane Using the TIBET AS Data. Proc. of the ICRC, 1997. [82] C. Sinnis et al. The MILAGRO Gamma-Ray Observatory. Bulletin of the American Astronomical Society, 1996. [83] D. A. Williams et al. Astrophysics around 100 GeV with STACEE. New Astronomy Review, 48:359–366, April 2004. [84] A. Moralejo et al. The MAGIC telescope for gamma-ray astronomy above 30 GeV. Chinese Journal of Astronomy and Astrophysics Supplement, 3:531–538, December 2003. [85] T. Bretz, D. Dorner, R. M. Wagner, and P. Sawallisch. The drive system of the major atmospheric gamma-ray imaging Cherenkov telescope. Astroparticle Physics, 31:92–101, March 2009. [86] A. Biland et al. The Active Mirror Control of the MAGIC Telescope. ph/0709.1574, September 2007.

astro-

[87] J. Cortina et al. Camera Control and Central Control of the MAGIC Telescope. Proc. of the 28th ICRC, 2003. [88] A. Stamerra, D. Corti, M. Mariotti, R. Paoletti, R. Pegna, A. Stammera, and N. Turini. The Trigger System of the MAGIC Telescope: On-Line Selection Strategies for Cherenkov Telescopes. Proc. of the 28th ICRC, 2003. [89] M. Rissi. Search for Pulsed High Energy γ-Rays from the Direction of the Central Pulsar in the Crab Nebula with the 17m Ground Based Air Cherenkov Telescope MAGIC. diploma thesis, February 2005.

164

Bibliography

[90] H. Bartko, F. Goebel, R. Mirzoyan, W. Pimpl, and M. Teshima. Tests of a prototype multiplexed fiber-optic ultra-fast FADC data acquisition system for the MAGIC telescope. Nuclear Instruments and Methods in Physics Research A, 548:464–486, August 2005. [91] M. Gaug et al. Calibration of the MAGIC Telescope. Proc. of the 29th ICRC, Pune, India, 2005. [92] R. Mirzoyan. Conversion factor calibration for MAGIC based on the use of measured F-factors of PMTs. MAGIC-TDAS 00-15, 2000. [93] F. Lucarelli et al. Development and first results of the MAGIC central pixel system for optical observations. Proc. of the 29th ICRC, Pune, India, 2005. [94] A. N. Otte. Observation of VHE Gamma-Rays from the Vicinity of magnetized Neutron Stars and Development of new Photon Detectors for Future Ground based Gamma-Ray Detectors. PhD Thesis, 2008. [95] M. Gaug. Calibration of the MAGIC telescope and observation of gamma-ray bursts. PhD Thesis, 2006. [96] N. Tonelli. Study of the VHE gamma-ray emission from the Active Galactic Nucleus 1ES1959+650. PhD Thesis, 2006. [97] M. Shayduk and others. A New Image Cleaning Method for the MAGIC Telescope. Proc. of the 29th ICRC, Pune, India, 2005. [98] R. de los Reyes, I. de Oya, and M. Gaug. private communication. 2008. [99] A. M. Hillas. Cerenkov light images of EAS produced by primary gamma. Proc. of the 3rd ICRC, 1985. [100] E. Aliu et al. Improving the performance of the single-dish Cherenkov telescope MAGIC through the use of signal timing. Astroparticle Physics, 30:293–305, January 2009. [101] J. Albert et al. Implementation of the Random Forest Method for the Imaging Atmospheric Cherenkov Telescope MAGIC. astro-ph/0709.3719, September 2007. [102] L. Breiman. Random forests. Machine Learning, 45:5, 2001. [103] T. Hengstebeck. Measurement of the energy spectrum of the BL Lac object PG1553+113 with the MAGIC telescope in 2005 and 2006. PhD Thesis, 2006. [104] T.-P. Li and Y.-Q. Ma. Analysis methods for results in gamma-ray astronomy. ApJ, 272:317–324, September 1983. [105] W. Wittek. Unfolding (revised). MAGIC-TDAS 05-05, 2005. [106] A.-N. Tikhonov and V.-J. Arsenin. Methods of Solution of Ill-posed Problems. 1979. [107] M. Bertero. Linear inverse and ill-posed problems. Advances in Electronics and Electron Physics, 75:1–120, 1989.

Bibliography

165

[108] M. Schmelling et al. The method of reduced cross-entropy. A general approach to unfold probability distributions. Nuclear Instruments and Methods in Physics Research A, 340:400, 1994. [109] A. Lyne et al. Jodrell bank crab pulsar monthly ephemeris. http://www.jb.man.ac.uk/~pulsar/crab.html. [110] D. Manchester and D. Nice. Tempo. http://www.atnf.csiro.au/research/pulsar/tempo/. [111] J. H. Taylor and J. M. Weisberg. Further experimental tests of relativistic gravity using the binary pulsar PSR 1913 + 16. ApJ, 345:434–450, October 1989. [112] J. M. Fierro, P. F. Michelson, P. L. Nolan, and D. J. Thompson. Phase-resolved Studies of the High-Energy Gamma-Ray Emission from the Crab, Geminga, and VELA Pulsars. ApJ, 494:734, February 1998. [113] O. C. de Jager, B. C. Raubenheimer, and J. W. H. Swanepoel. A powerful test for weak periodic signals with unknown light curve shape in sparse data. aap, 221:180– 190, August 1989. [114] N. H. Kuiper. Tests concerning random points on a circle. Proceedings of the Koninklijke Nederlandse Akademie van Wetenschappen, Series A 63: 38. [115] D. Petry et al. The MAGIC Telescope - Bridging the Gap at 10 - 200 GeV. Astronomical Society of the Pacific Conference Series, 1999. [116] Razmik Mirzoyan. private communication. 2006. [117] S. Funk, G. Hermann, J. Hinton, D. Berge, K. Bernlöhr, W. Hofmann, P. Nayman, F. Toussenel, and P. Vincent. The trigger system of the H.E.S.S. telescope array. Astroparticle Physics, 22:285–296, November 2004. [118] E. de Oña Wilhelmi. The Optimization of the MAGIC Telescope for Pulsar Observations. PhD thesis, December 2003. [119] C. Hsu. private communication. 2007. [120] E. Lorentz. Clipping of PMT Afterpulses. private communication, 2007. [121] D. Heck, G. Schatz, T. Thouw, J. Knapp, and J. N. Capdevielle. Corsika: A monte carlo code to simulate extensive air showers. 1998. FZKA-6019. [122] Maxim Shayduk. private communication. 2007. [123] R. Paoletti. Sumo - Control Program For the Sumtrigger System. MAGIC internal document, 2008. [124] M. Rissi, N. Otte, T. Schweizer, and M Shayduk. A new trigger provides a lower energy threshold for the MAGIC Cherenkov telescope. Nuclear Science Symposium Conference Record, 2008. NSS ’08. IEEE. [125] E. Aliu et al. Detection of pulsed gamma-ays above 25 gev from the crab pulsar. Sep 2008.

166

Bibliography

[126] P. Goret, T. Palfrey, A. Tabary, G. Vacanti, and R. Bazer-Bachi. Observations of TeV gamma rays from the Crab nebula. A&A, 270:401–406, March 1993. [127] F. Aharonian et al. The Crab Nebula and Pulsar between 500 GeV and 80 TeV: Observations with the HEGRA Stereoscopic Air Cerenkov Telescopes. ApJ, 614:897– 913, October 2004. [128] M. Amenomori et al. Observation of Multi-TEV Gamma Rays from the Crab Nebula using the Tibet Air Shower Array. ApJ, 525:L93–L96, November 1999. [129] G. Vacanti et al. Gamma-ray observations of the Crab Nebula at TeV energies. ApJ, 377:467–479, August 1991. [130] H. Manseri et al. Results from the CELESTE experiment. In SF2A-2004: Semaine de l’Astrophysique Francaise, 2004. [131] F. Aharonian et al. Observations of the Crab nebula with HESS. A&A, 457:899–915, October 2006. [132] S. C. Commichau, A. Biland, D. Kranich, R. de los Reyes, A. Moralejo, and D. Sobczynska. Geomagnetic Field Effects on the Imaging Air Shower Cherenkov Technique. astro-ph/0709.1251. [133] T. Oosterbroek et al. Simultaneous absolute timing of the Crab pulsar at radio and optical wavelengths. A&A, 488:271–277, September 2008. [134] A. K. Harding. private communication, 2008. [135] F. W. Stecker, S. D. Hunter, and D. A. Kniffen. The likely cause of the EGRET GeV anomaly and its implications. Astroparticle Physics, 29:25–29, February 2008. [136] S. J. Sheather and M. C. Jones. A reliable data-based bandwidth selection method for kernel density estimation. Journal of the Royal Statistical Society. Series B, 460:743–749, December 1991. [137] A. Pellizzoni et al. AGILE observations of Crab (and other pulsars). In Polarimetry days in Rome: Crab status, theory and prospects, 2008. [138] D. A. Moffett and T. H. Hankins. Polarimetric Properties of the Crab Pulsar between 1.4 and 8.4 GHZ. ApJ, 522:1046–1052, September 1999. [139] K. Hirotani. private communication. 2009. [140] S. S. Eikenberry and G. G. Fazio. The Infrared to Gamma-Ray Pulse Shape of the Crab Nebula Pulsar. ApJ, 476:281, February 1997. [141] J. W. Percival et al. The Crab pulsar in the visible and ultraviolet with 20 microsecond effective time resolution. ApJ, 407:276–283, April 1993. [142] J. Alfaro, H. A. Morales-Técotl, and L. F. Urrutia. Loop quantum gravity and light propagation. Phys. Rev. D, 65(10):103509, May 2002. [143] G. Amelino-Camelia and D. V. Ahluwalia. Relativity in Spacetimes with ShortDistance Structure Governed by an Observer-Independent (Planckian) Length Scale. International Journal of Modern Physics D, 11:35–59, 2002.

Bibliography

167

[144] G. Amelino-Camelia, J. Ellis, N. E. Mavromatos, D. V. Nanopoulos, and S. Sarkar. Tests of quantum gravity from observations of γ-ray bursts. Nature, 393:763–765, June 1998. [145] J. Ellis, K. Farakos, N. E. Mavromatos, V. A. Mitsou, and D. V. Nanopoulos. A Search in Gamma-Ray Burst Data for Nonconstancy of the Velocity of Light. ApJ, 535:139–151, May 2000. [146] J. Ellis, N. E. Mavromatos, and D. V. Nanopoulos. String theory modifies quantum mechanics. Physics Letters B, 293:37–48, October 1992. [147] J. Ellis, N. E. Mavromatos, D. V. Nanopoulos, and A. S. Sakharov. Quantum-gravity analysis of gamma-ray bursts using wavelets. A&A, 402:409–424, May 2003. [148] L. Smolin. How far are we from the quantum theory of gravity? ArXiv High Energy Physics - Theory e-prints, March 2003. [149] A. A. Abdo et al. Fermi Observations of High-Energy Gamma-Ray Emission from GRB 080916C. Science, in print, 2009. [150] S. D. Biller et al. Limits to Quantum Gravity Effects on Energy Dependence of the Speed of Light from Observations of TeV Flares in Active Galaxies. Physical Review Letters, 83:2108–2111, September 1999. [151] O. Blanch, J. Lopez, and M. Martinez. Testing the effective scale of quantum gravity with the next generation of gamma ray telescopes. Astroparticle Physics, 19:245–252, May 2003. [152] S. Sarkar. Possible Astrophysical Probes of Quantum Gravity. Modern Physics Letters A, 17:1025–1035, 2002. [153] J. Albert et al. Probing quantum gravity using photons from a flare of the active galactic nucleus Markarian 501 observed by the MAGIC telescope. Physics Letters B, 668:253–257, October 2008. [154] F. Aharonian et al. Limits on an Energy Dependence of the Speed of Light from a Flare of the Active Galaxy PKS 2155-304. Physical Review Letters, 101(17):170402, October 2008. [155] P. Kaaret. Pulsar radiation and quantum gravity. A&A, 345:L32–L34, May 1999. [156] E. Massaro, R. Campana, G. Cusumano, and T. Mineo. The optical to γ-ray emission of the Crab pulsar: a multicomponent model. A&A, 459:859–870, December 2006. [157] L. Baldini et al. Preliminary results of the LAT Calibration Unit beam tests. In The First GLAST Symposium, volume 921 of American Institute of Physics Conference Series, 2007. [158] D. Bastieri et al. Using the photons from the Crab Nebula seen by GLAST to calibrate MAGIC and the imaging air Cherenkov telescopes. Astroparticle Physics, 23:572–576, July 2005. [159] W. F. Tompkins, B. B. Jones, P. L. Nolan, G. Kanbach, P. V. Ramanamurthy, and D. J. Thompson. EGRET Gamma-Ray Observations of the Crab P1/P2 Ratio. ApJ, 487:385, September 1997.

168

Bibliography

[160] P. M. Saz Parkinson et al. Fermi Gamma-ray Space Telescope Observations of Gamma-ray Pulsars. ArXiv e-prints, January 2009. [161] M. Frigo and S. Johnson. FFTW - the fastest fourier transform of the west. http://www.fftw.org .

Acknowledgement I am deeply grateful for the excellent working conditions I have in Felicitas Pauss’ group at the Institute for Particle Physics at the ETH in Zurich. I appreciate the freedom to study subjects I am interested in and the confidence in my success. I want to express my gratitude to Felicitas Pauss for the opportunity to write my PhD thesis as a member of the MAGIC collaboration and for her support during my dissertation. I am thankful to Christoph Grab, who kindly accepted to act as co-referent and from whom I learnt interesting insights in the world of particle physics and statistics. I want to thank Adrian Biland who was my supervisor. We had many productive discussions, and his critical objections helped a lot to achieve and improve the results presented in this thesis. I am much indebted to Daniel Kranich who carefully read this thesis and helped to improve the manuscript, and to Johanna Amberg Sayed who proof-read the entire thesis. My thanks also to Dorothee Hildebrand, Quirin Weitzel and Sabrina Stark, who also read parts of the manuscript. It was a great pleasure to work together with the MAGIC sum trigger team: Nepomuk Otte, Maxim Shayduk and Thomas Schweitzer, from whom I learnt a lot about physics, simulations, programming, electronics and determination against all odds, as well as about hiking trails and hidden beaches on La Palma. Florian Goebel, Eckart Lorenz, Razmick Mirzoyan and Masahiro Teshima provided us with vital inputs for the design and the construction of the sum trigger as well as for the difficult data analysis. Last but not least, I want to thank Marcos Lopez who joined the sum trigger analysis team in 2007. Working in the MAGIC collaboration has been a very fruitful experience. I particularly appreciated the team work during the long night shifts at the MAGIC telescope on La Palma and the discussions at the numerous meetings. I appreciated very much the collaboration with my colleagues at ETH, Sebastian Commichau and Sabrina Stark, who introduced me to the experiment and helped me to understand the programs and physical topics with patience. I enjoyed working together with Isabel Braun, Ilia Britvitch, Dorothee Hildebrand, Daniel Kienzler, Thomas Krähenbühl, Daniel Kranich, Carmelo Marchica and Quirin Weitzel. I want to thank Gert Viertel, Günther Dissertori, Hans Anderhub, Ulf Röser, Hanspeter von Gunten and Urs Horisberger. My thanks to the secretaries Johanna Amberg Sayed and Nadia Sigrist-Müller. The time I spent here at ETH would not have been the same without my friends Martin Bernet, Mirko Birbaumer, Simon Bruderer, and Michel Sauter. I want to thank them and my family for their support and the different point of view they could provide. My appreciation to the Omori family. This thesis is dedicated to my love Regina.

170

Bibliography