nexoncn.com

文档资料共享网 文档搜索专家

文档资料共享网 文档搜索专家

当前位置：首页 >> >> A Search for Gamma-Ray Bursts and Pulsars, and the Application of Kalman Filters to Gamma-R_图文

arXiv:astro-ph/0202088v1 4 Feb 2002

A SEARCH FOR GAMMA-RAY BURSTS AND PULSARS, AND THE APPLICATION OF KALMAN FILTERS TO GAMMA-RAY RECONSTRUCTION

a dissertation submitted to the department of physics and the committee on graduate studies of stanford university in partial fulfillment of the requirements for the degree of doctor of philosophy

By Brian Jones October 1998

c Copyright 1999 by Brian Jones All Rights Reserved

ii

I certify that I have read this dissertation and that in my opinion it is fully adequate, in scope and quality, as a dissertation for the degree of Doctor of Philosophy.

Peter Michelson (Principal Adviser)

I certify that I have read this dissertation and that in my opinion it is fully adequate, in scope and quality, as a dissertation for the degree of Doctor of Philosophy.

Vah? e Petrosian

I certify that I have read this dissertation and that in my opinion it is fully adequate, in scope and quality, as a dissertation for the degree of Doctor of Philosophy.

Robert Wagoner

Approved for the University Committee on Graduate Studies:

iii

Abstract

High-energy γ -ray astronomy was revolutionized in 1991 with the launch of the Energetic Gamma-Ray Experiment Telescope (EGRET ) on board the Compton GammaRay Observatory. In addition to unprecedented instrument e?ective area and a narrow point-spread function, EGRET provided photon time-tagging to an absolute accuracy of 100 ?s. The opportunity to analyze high-quality γ -ray data requires sophisticated statistical and analytic tools. Part I describes the analysis of periodic and transient signals in EGRET data. A method to search for the transient ?ux from γ -ray bursts independent of triggers from other γ -ray instruments is developed. Several known γ -ray bursts were independently detected, and there is evidence for a previously unknown γ -ray burst candidate. Statistical methods using maximum likelihood and Bayesian inference are developed and implemented to extract periodic signals from γ -ray sources in the presence of signi?cant astrophysical background radiation. The methods allow searches for periodic modulation without a priori knowledge of the period or period derivative. The analysis was performed on six pulsars and three pulsar candidates. The three brightest pulsars, Crab, Vela, and Geminga, were readily identi?ed, and would have been detected independently in the EGRET data without knowledge of the pulse period. No signi?cant pulsation was detected in the three pulsar candidates. Furthermore, the method allows the analysis of sources with periods on the same order as the time scales associated with changes in the instrumental sensitivity, such as the orbital time scale of CGRO around the Earth. Eighteen Xray binaries were examined. None showed any evidence of periodicity. In addition, methods for calculating the detection threshold of periodic ?ux modulation were developed. iv

The future hopes of γ -ray astronomy lie in the development of the Gamma-ray Large Area Space Telescope, or GLAST . Part II describes the development and results of the particle track reconstruction software for a GLAST science prototype instrument beamtest. The Kalman ?ltering method of track reconstruction is introduced and implemented. Monte Carlo simulations, very similar to those used for the full GLAST instrument, were performed to predict the instrumental response of the prototype. The prototype was tested in a γ -ray beam at SLAC. The reconstruction software was used to determine the incident γ -ray direction. It was found that the simulations did an excellent job of representing the actual instrument response.

v

Acknowledgements

The acknowledgements page of the modern physics thesis has become a stage upon which seasoned and world-weary graduate students perform a tragicomic stand-up routine to amuse their family and friends. It is altogether ?tting and proper that they should do so, since it is unlikely that those family and friends will read any of the rest of the thesis. Besides, they hope they will be mentioned. Nevertheless, it is a good medium, and a useful way to point out the people and institutions without whom, for better or worse, this work would never have been performed. The number of people who have been a part of this cause is enormous. There is no way I can mention all of the friends and colleagues whom I appreciate, so I will just touch on a few. Beginning with the most general, I must thank the taxpayers of the United States who, whether they knew it or not, spent a vast amount of money to advance the frontiers of γ -ray astrophysics. I can only hope to have done them proud. I was fortunate enough to pro?t from data taken by the EGRET instrument, on board the Compton Gamma-Ray Observatory. I have deep appreciation for the work done to get EGRET into orbit and taking data by the scientists and sta? of Goddard Space Flight Center, especially Dave Thompson, who was always particularly supportive of my work, and the Max Planck Institut f¨ ur Extraterrestriche Physik. This collaboration has done an outstanding job in designing, building, and operating an instrument which has radically changed our understanding of the γ -ray sky. I’m proud to have been a small part of that group. I’m also proud to be a part of a second international collaboration, this one responsible for GLAST . I have bene?tted from the chance to work with scientists from Goddard, the Naval Research Laboratory, UC Santa Cruz, and Columbia University vi

? in the United States, as well as the Ecole Polytechnique and CEA in France, the University of Trieste in Italy and the University of Tokyo, as well as dozens of other institutions worldwide too numerous to mention. Being on the inside of a developing Big Science project has been eye-opening, to say the least. Speci?cally, I’d like to thank Steve Ritz, Neil Johnson, and Eric Grove for listening to the brash opinions of a grad student as if he were a Real Scientist. The fact that no Stanford (or SLAC) people are included above re?ects not their absence, but their proximity. A number of people in and around Stanford have made an indelible impression on me. My advisor Peter Michelson has given me tremendous freedom to pursue my own interests, even when they were not what most astrophysicists would call mainstream. Bill Atwood has been a tremendous motivator, and taught me a great deal. His methods are generally quite harsh, by his own admission, and he is never satis?ed. For some reason which I still do not understand, this had the e?ect of spurring me on to better, more careful, and more thorough work, completed more rapidly, rather than to hatred and bitterness. Joe Fierro’s Pucklike sensibilities drew me to the group in the ?rst place, in the belief that it would be a fun place to work. He was a great help once I pestered him enough. Tom Willis graciously put up with a newcomer in his o?ce, and helped me through the steep learning curve in EGRET minuti?. Although he’s not at Stanford, Prof. Dave Wilkinson along with the rest of the Princeton faculty showed me what I am really capable of, and gave me the best undergraduate physics education available. After that, grad school was a breeze. And of course, without Marcia Keating in the Varian o?ce, I would never have ?led a study list, gotten paid, or managed to ?le the appropriate graduation forms. It doesn’t matter how many Nobel Laureates we have wandering around; without Marcia, this place would fall apart. She will be greatly missed as she moves on to new challenges. Pat Nolan, my de facto advisor, deserves a paragraph to himself. His patience, advice, and knowledge have been a prime factor in whatever success I may have had. He also deserves rich thanks for maintaining the computer systems and knowing arcane details about UNIX so that I didn’t have to. It is a great privilege to have had the opportunity to learn from him. vii

Also critical have been my fellow graduate students in the Physics Department, especially Doug Natelson and Bill Tompkins. Doug, as my erstwhile roommate, had to put up with a lot and endured it well. Bill, as roommate and o?cemate, found he could never avoid me. He therefore sought to educate me, so that I would be more bearable. If I have given half as good as I’ve gotten from him, I consider myself a success. On a more personal note, I would like to thank all of the members of TMBS for their support and prayers. I know the Lord would have given me strength without them, but they knocked at the door persistently, and it has paid o?. Primarily, however, support for an undertaking such as a Ph.D. has to come from the family. Fortunately, I have been blessed with the best family I could imagine. My parents have supported me in my quest “to be a scientist” from the very beginning, and always believed I could do it, even if they weren’t sure why I wanted to. Plus, they paid for college. My brother Bill, aside from the non-trivial accomplishment of being a fantastically supportive brother, decided to study physics just so I wouldn’t feel alone. And last, but by no means least, my beautiful wife Kim met me while I was yet a graduate student, and saw some kernel of value in me even then. Her love, support, and con?dence in me have given me the strength to ?nish this work. My great blessing is that although this thesis is now part of the past, we have many wonderful years together in our future.

viii

O LORD, our Lord, how majestic is your name in all the earth! You have set your glory above the heavens. When I consider your heavens, the work of your ?ngers, the moon and the stars, which you have set in place, what is man that you are mindful of him, the son of man that you care for him?

Psalm 8 (NIV)

The heavens declare the glory of God; the skies proclaim the work of his hands. Day after day they pour forth speech; night after night they display knowledge. There is no speech or language where their voice is not heard. Their voice goes out into all the earth, their words to the ends of the world. In the heavens he has pitched a tent for the sun, which is like a bridegroom coming forth from his pavilion, like a champion rejoicing to run his course.

Psalm 19 (NIV)

ix

Contents

Abstract Acknowledgements iv vi

I

Methods for Time-Series Analysis

1

2 3 5 5 7 8 10 11 12 13 14 14 15 16 17 19

1 Introduction 1.1 1.2 The γ -Ray Observatories . . . . . . . . . . . . . . . . . . . . . . . . . The EGRET instrument . . . . . . . . . . . . . . . . . . . . . . . . . 1.2.1 1.2.2 1.3 Instrumental Design . . . . . . . . . . . . . . . . . . . . . . . Instrumental Calibration and Performance . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

Successes with EGRET

2 Statistical Methods in γ -Ray Astronomy 2.1 The Nature of the γ -Rays . . . . . . . . . . . . . . . . . . . . . . . . 2.1.1 2.1.2 2.2 2.2.1 2.2.2 2.2.3 2.3 2.4 Di?use Background . . . . . . . . . . . . . . . . . . . . . . . . Spectral Di?erences . . . . . . . . . . . . . . . . . . . . . . . . Point-Spread Function . . . . . . . . . . . . . . . . . . . . . . Sensitive Area . . . . . . . . . . . . . . . . . . . . . . . . . . . Energy Dispersion . . . . . . . . . . . . . . . . . . . . . . . .

Instrumental E?ects . . . . . . . . . . . . . . . . . . . . . . . . . . .

Likelihood Analysis . . . . . . . . . . . . . . . . . . . . . . . . . . . . Applying Likelihood Analysis to EGRET . . . . . . . . . . . . . . . .

x

2.4.1 2.4.2 2.5 2.6

Binned Likelihood

. . . . . . . . . . . . . . . . . . . . . . . .

20 23 24 25 26 28 29 30 31 32 33 34 36 39 40 41 43 45 46 48 49 53 55 60 62 63 65 67 69 70

Maximizing the Likelihood with LIKE . . . . . . . . . . . . . .

Parameter Estimation vs. Hypothesis Testing . . . . . . . . . . . . . Bayesian Methods . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2.6.1 2.6.2 2.6.3 2.6.4 Bayes’ Theorem . . . . . . . . . . . . . . . . . . . . . . . . . . The Odds Ratio . . . . . . . . . . . . . . . . . . . . . . . . . . Marginalization and Con?dence Regions . . . . . . . . . . . . Advantages and Disadvantages . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . Upper Limits from LIKE . . . . . . . . . . . . . . . . . . . . . Calculating More Accurate Upper Limits . . . . . . . . . . . . Implications for Extant Conclusions . . . . . . . . . . . . . . .

2.7

Calculating Upper Limits 2.7.1 2.7.2 2.7.3

3 γ -Ray Bursts 3.1 3.2 Recent Observations . . . . . . . . . . . . . . . . . . . . . . . . . . . γ -Ray Burst Models . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.2.1 3.2.2 3.2.3 3.3 3.4 Energy production mechanisms . . . . . . . . . . . . . . . . . Blast wave theories . . . . . . . . . . . . . . . . . . . . . . . . Observable consequences . . . . . . . . . . . . . . . . . . . . .

EGRET observations . . . . . . . . . . . . . . . . . . . . . . . . . . . Possible EGRET -only Bursts . . . . . . . . . . . . . . . . . . . . . . 3.4.1 3.4.2 3.4.3 Statistical Methods . . . . . . . . . . . . . . . . . . . . . . . . Results . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . Discussion . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

3.5

Conclusions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

4 Periodic Time-Series Analysis 4.1 Pulsars . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4.1.1 4.1.2 4.2 4.2.1 EGRET Contributions . . . . . . . . . . . . . . . . . . . . . . Models . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . Previous Methods . . . . . . . . . . . . . . . . . . . . . . . . . xi

Statistical Methods . . . . . . . . . . . . . . . . . . . . . . . . . . . .

4.2.2 4.2.3 4.2.4 4.2.5 4.2.6 4.3 4.3.1 4.3.2 4.3.3 4.4 4.4.1 4.4.2

Maximum Likelihood . . . . . . . . . . . . . . . . . . . . . . . Application to EGRET . . . . . . . . . . . . . . . . . . . . . Bayesian Inference . . . . . . . . . . . . . . . . . . . . . . . . Upper Limits and Thresholds . . . . . . . . . . . . . . . . . . Bin-free Maximum Likelihood . . . . . . . . . . . . . . . . . . Measurement of Known Pulsars . . . . . . . . . . . . . . . . . Searches for Geminga-like Pulsars . . . . . . . . . . . . . . . .

75 78 83 85 90 91 92 95

Searching for Pulsars . . . . . . . . . . . . . . . . . . . . . . . . . . .

Results . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 100 Thresholds and Searches . . . . . . . . . . . . . . . . . . . . . 102 Results . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 103 107

X-Ray Binaries . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 101

5 Conclusions

II

The October 1997 GLAST -prototype Beam Test

111

112

6 GLAST : The Next Generation 6.1 6.2 Potential Improvements

. . . . . . . . . . . . . . . . . . . . . . . . . 113

The Baseline GLAST Instrument . . . . . . . . . . . . . . . . . . . . 115 118

7 Testing the GLAST Science Prototype 7.1 7.2

?

The SLAC e Beam . . . . . . . . . . . . . . . . . . . . . . . . . . . 118 The Beam Test Instrument . . . . . . . . . . . . . . . . . . . . . . . . 121 7.2.1 7.2.2 7.2.3 Tracker . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 121 Calorimeter . . . . . . . . . . . . . . . . . . . . . . . . . . . . 124 Anti-Coincidence Detector . . . . . . . . . . . . . . . . . . . . 125

7.3 7.4

Summary of Data . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 126 Simulations . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 127 129

8 Reconstructing Events 8.1

Pair Production . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 130 xii

8.2 8.3

Track Reconstruction . . . . . . . . . . . . . . . . . . . . . . . . . . . 132 The Annotated Kalman Filtering Formulae . . . . . . . . . . . . . . . 136 8.3.1 8.3.2 8.3.3 8.3.4 The Filtering Equations . . . . . . . . . . . . . . . . . . . . . 138 The Smoothing Equations . . . . . . . . . . . . . . . . . . . . 139 Goodness of Fit . . . . . . . . . . . . . . . . . . . . . . . . . . 141 Kalman Filter Implementation for the Beam Test . . . . . . . 142 The Exhaustive Search . . . . . . . . . . . . . . . . . . . . . . 143 Beam Test Algorithm . . . . . . . . . . . . . . . . . . . . . . . 144

8.4

Track Finding Algorithm . . . . . . . . . . . . . . . . . . . . . . . . . 143 8.4.1 8.4.2

8.5 8.6 8.7

Measurement Error and Multiple Scatter Estimates . . . . . . . . . . 146 Noise . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 148 Energy . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 150 8.7.1 8.7.2 Energy Splitting Energy Dispersion . . . . . . . . . . . . . . . . . . . . . . . . . 151 . . . . . . . . . . . . . . . . . . . . . . . . 151 . . . . . . . . . . . . . . . . . . . . . . . . . 153

8.8

Potential Improvements 8.8.1 8.8.2 8.8.3

Track Fitting . . . . . . . . . . . . . . . . . . . . . . . . . . . 153 Track Finding . . . . . . . . . . . . . . . . . . . . . . . . . . . 154 Energy Estimations . . . . . . . . . . . . . . . . . . . . . . . . 156

8.9

Calculating the Point-Spread Function . . . . . . . . . . . . . . . . . 157

8.10 Extended Kalman Filters . . . . . . . . . . . . . . . . . . . . . . . . . 159 9 Instrument Response 9.1 9.2 9.3 9.4 160

Alignment . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 160 Cuts . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 161 Expected Beam Test Point-Spread Widths . . . . . . . . . . . . . . . 163 Conclusions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 166 9.4.1 9.4.2 Comparison of Beam Test and Monte Carlo Results . . . . . . 166 Implications for GLAST . . . . . . . . . . . . . . . . . . . . . 169 171 173

A SSB Arrival Time Corrections B Summary of Kalman Filtering Equations xiii

C Track Finding Algorithm Bibliography

174 177

xiv

List of Tables

1.1 3.1 3.2 3.3 3.4 3.5 3.6 3.7 4.1 4.2 4.3 4.4 4.5 4.6 4.7 4.8 4.9 6.1 6.2 7.1 7.2 γ -Ray telescope performance . . . . . . . . . . . . . . . . . . . . . . . Three minute burst time scale results . . . . . . . . . . . . . . . . . . Ten minute burst time rcale results . . . . . . . . . . . . . . . . . . . Thirty minute time scale results . . . . . . . . . . . . . . . . . . . . . One hour burst time scale results . . . . . . . . . . . . . . . . . . . . Triggered bursts not independently detected . . . . . . . . . . . . . . Hardness ratios for burst candidates . . . . . . . . . . . . . . . . . . . Joint burst probabilities . . . . . . . . . . . . . . . . . . . . . . . . . Typical pulsar parameters . . . . . . . . . . . . . . . . . . . . . . . . Detection of known γ -ray pulsars by timevar . . . . . . . . . . . . . Thresholds for pulsed sources . . . . . . . . . . . . . . . . . . . . . . 7 54 54 55 55 56 57 59 65 93 98

GLAST pulsed signal thresholds . . . . . . . . . . . . . . . . . . . . . 100 Results of period search with timevar. . . . . . . . . . . . . . . . . . 101 Low mass X-ray binaries . . . . . . . . . . . . . . . . . . . . . . . . . 104 High mass X-ray binaries . . . . . . . . . . . . . . . . . . . . . . . . . 105 Low-mass X-ray binary results . . . . . . . . . . . . . . . . . . . . . . 106 High-mass X-ray binary results . . . . . . . . . . . . . . . . . . . . . 106 GLAST baseline characteristics . . . . . . . . . . . . . . . . . . . . . 115 Institutions in the GLAST collaboration . . . . . . . . . . . . . . . . 117 Radiation lengths of Pb available . . . . . . . . . . . . . . . . . . . . 123 Energy bands used for analysis . . . . . . . . . . . . . . . . . . . . . 124 xv

9.1 9.2

Cut e?ciencies . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 163 Calculated point-spread widths . . . . . . . . . . . . . . . . . . . . . 170

xvi

List of Figures

1.1 1.2 1.3 2.1 2.2 2.3 2.4 2.5 3.1 3.2 4.1 4.2 4.3 4.4 4.5 4.6 4.7 4.8 4.9 The electromagnetic spectrum . . . . . . . . . . . . . . . . . . . . . . The Compton Gamma-Ray Observatory . . . . . . . . . . . . . . . . The EGRET Instrument . . . . . . . . . . . . . . . . . . . . . . . . . Full sky photon counts . . . . . . . . . . . . . . . . . . . . . . . . . . Full sky galactic di?use map . . . . . . . . . . . . . . . . . . . . . . . Full sky exposure map . . . . . . . . . . . . . . . . . . . . . . . . . . Full sky intensity map . . . . . . . . . . . . . . . . . . . . . . . . . . Bayesian upper limits . . . . . . . . . . . . . . . . . . . . . . . . . . . γ -Ray bursts observed by BATSE . . . . . . . . . . . . . . . . . . . . Overall γ -ray burst signi?cance from the Monte Carlo simulations. . . Geminga likelihoods . . . . . . . . . . . . . . . . . . . . . . . . . . . Geminga light curves . . . . . . . . . . . . . . . . . . . . . . . . . . . Geminga phase o?sets . . . . . . . . . . . . . . . . . . . . . . . . . . Marginalized geminga likelihood . . . . . . . . . . . . . . . . . . . . . Pulsar periods vs. period derivatives . . . . . . . . . . . . . . . . . . Likelihood of ?ux modulation from Crab . . . . . . . . . . . . . . . . Light curves of Crab . . . . . . . . . . . . . . . . . . . . . . . . . . . Likelihood of ?ux modulation from Vela . . . . . . . . . . . . . . . . Light curves of Vela . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2 4 6 11 13 22 24 34 38 52 82 83 85 86 92 94 94 95 96 96 97

4.10 Likelihood of ?ux modulation in PSR 1706-44 . . . . . . . . . . . . . 4.11 Light curves of PSR 1706-44 . . . . . . . . . . . . . . . . . . . . . . .

xvii

2 ˙ 4.12 E/D vs. pulsar period . . . . . . . . . . . . . . . . . . . . . . . . . .

99

6.1 7.1 7.2 7.3 7.4 8.1 8.2 8.3 8.4 8.5 8.6

The GLAST satellite . . . . . . . . . . . . . . . . . . . . . . . . . . . 113 Location of End Station A at SLAC . . . . . . . . . . . . . . . . . . . 120 Beam test experimental scheme . . . . . . . . . . . . . . . . . . . . . 121 Pancake and Stretch tracker con?gurations . . . . . . . . . . . . . . . 122 CsI(Tl) Calorimeter . . . . . . . . . . . . . . . . . . . . . . . . . . . . 125 Kalman ?ltering . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 133 Kalman smoothing . . . . . . . . . . . . . . . . . . . . . . . . . . . . 135 Measurement error distributions . . . . . . . . . . . . . . . . . . . . . 147 Strip occupancies . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 149 Monte Carlo energy dispersion . . . . . . . . . . . . . . . . . . . . . . 152 Probability distributions of e? e+ energy split as a function of total incident energy . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 158

9.1 9.2 9.3 9.4

Point-spread width dependence on total energy and energy splitting fraction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 164 Pancake X-Projection PSF vs. Monte Carlo . . . . . . . . . . . . . . 166 Stretch X-Projection PSF vs. Monte Carlo . . . . . . . . . . . . . . . 168 Projected point-spread distributions . . . . . . . . . . . . . . . . . . . 169

xviii

Part I Methods for Time-Series Analysis

1

Chapter 1 Introduction

The desire to experiment with the extremes of nature is innate. While it is almost always amusing, it can be informative as well. Such is the case in the study of some of the most energetic photons produced in space: the γ -rays (Figure 1.1). With at least ten million times the energy of ordinary optical photons, γ -rays represent a unique window into the most energetic processes in astrophysics—the (electromagnetically) roaring jets from active galactic nuclei, the arcing plasmas in the intense gravitational ?elds of pulsars, and the enigmatic and inordinately powerful explosions known as γ -ray bursts. The immense value of γ -rays for astrophysics lies both in their role as telltale markers of large energy-generation processes and in their likelihood of passing unperturbed through vast reaches of intergalactic space. Measuring the γ -ray spectra of astrophysical objects sharply constrains estimates of their total energy output.

Energy (eV) -7 -9 10 10

-5 -3

10

10

-1 10 1 eV

10

1

10

3

10

5

10

7

10

9

Radio

Microwave

Infrared Optical

Ultraviolet

X-ray

Gamma-ray

Figure 1.1: The Electromagnetic Spectrum. γ -Rays occupy the highest energy extreme of the electromagnetic spectrum, from 10 MeV to over 300 GeV.

2

CHAPTER 1. INTRODUCTION

3

Since γ -rays are so energetic, many astrophysical sources emit the bulk of their total power output at these high energies. In addition, γ -rays travel relatively unimpeded through space. Since they carry no charge, they are nearly una?ected by galactic and intergalactic magnetic ?elds. Their small interaction cross section means that they are relatively una?ected by dust and gas in the intervening space between the source and the detector. A high-energy γ -ray can travel through the central plane of the Galaxy with only a 1% chance of being absorbed [43]. γ -Rays may be observed from Earth essentially unchanged since they left the distant violence in which they were created. However, every silver lining has its cloud. The Earth’s atmosphere is very good at absorbing γ -rays. Unfortunately, it means that precise astrophysical γ -ray observations must be done in space. The second di?culty in γ -ray astronomy is intrinsic to the energy production mechanisms that produce γ -rays in the ?rst place. Because γ rays are so energetic, most sources produce very few of them. Detectors must be very e?cient to collect these rare photons, and at the same time be able to discriminate against the sea of undesirable charged particles trapped in the Earth’s magnetic ?eld and albedo γ -rays which are generated in the Earth’s atmosphere. This discrimination requires background rejection on the order of one part in ?105 or better.

1.1

The γ -Ray Observatories

Experience with terrestrial accelerator-based γ -ray detectors suggested that a spark chamber might be an e?ective astrophysical γ -ray detector. In the mid-1970s, SAS 2 [39] and COS B [12] proved the concept of a γ -ray satellite telescope, while discovering several of the brightest γ -ray sources. Simultaneously, NASA envisioned the Great Observatories for Space Astrophysics program: a series of satellite telescopes designed to give unprecedented insights into electromagnetic emission from infrared to γ -rays. Under the auspices of this program, the Space InfraRed Telescope Facility (SIRTF ) infrared telescope has been designed, the Advanced X-ray Astrophysics

CHAPTER 1. INTRODUCTION

4

Figure 1.2: The Compton Gamma-Ray Observatory. EGRET is the dome on the right end of the spacecraft. It is coaligned with COMPTEL. Rounding out the instruments aboard CGRO are OSSE , which is very sensitive to e? e+ annihilation lines, and the omni-directional BATSE γ -ray burst detector. Facility (AXAF) has been designed and built, the Hubble Space Telescope has revolutionized optical astronomy, and the Compton Gamma-Ray Observatory (CGRO ) has explored the γ -ray sky from less than 0.1 MeV to more than 10 GeV (Figure 1.2). The ?ve orders of magnitude in energy of the electromagnetic spectrum observed by CGRO require four di?erent instruments on the satellite. The lowest energy γ rays interact primarily through the photo-electric e?ect. The Oriented Scintillation Spectrometer Experiment (OSSE ) covers the energy range from 0.05–10 MeV with a ?eld of view of 3? . 8×11? . 4 [78]. The Compton Telescope (COMPTEL) detects Compton scattered electrons, the most signi?cant γ -ray interaction in the energy range between 1 MeV and almost 30 MeV, to image the γ -ray sky with a ?eld of view of ?1 sr [178]. The Energetic Gamma-Ray Telescope Experiment (EGRET ) measures pairconversion events in a spark chamber, like SAS 2 and COS B . It is sensitive to energies between 20 MeV and 30 GeV, with a ?eld of view of ?1 sr [71, 92]. In addition, a fourth instrument aboard CGRO is optimized to detect γ -ray bursts. The Burst and

CHAPTER 1. INTRODUCTION

5

Transient Source Experiment (BATSE ) consists of eight uncollimated detectors, one on each corner of the CGRO spacecraft, sensitive to 25 keV–2 MeV γ -rays with nearly uniform coverage of the sky [47]. The EGRET instrument is the focus of Part I of this work. The instrument was built and operated by a collaboration of scientists at Stanford University, Goddard Space Flight Center (Greenbelt, Maryland), the Max Planck Institut f¨ ur Extraterrestrische Physik (Garching, Germany), and the Grumman Aerospace Corporation (Bethpage, New York). It was launched aboard CGRO on the Space Shuttle Atlantis (STS-37) on April 5, 1991 and was deployed two days later. It was activated on April 15, and began taking data on April 20. The instrument and its characteristics have been extensively documented [71, 92, 93, 150, 193]; we will brie?y touch on the highlights relevant to data analysis in §2.2. Future γ -ray telescopes will further extend our understanding of astrophysical γ -

ray processes. The GLAST instrument, scheduled for launch in 2005 [134, 2, 14, 15], will improve upon the successes of EGRET as it brings γ -ray astronomy to the 21st century. The test of a GLAST science prototype will be the focus of Part II of this work. Details of the current instrument baseline design are given in §6.2.

1.2

The EGRET instrument

Spark chambers detect γ -rays via pair production. Pair production refers to the process whereby a γ -ray converts to an electron–positron pair in the presence of matter. The detection process is more fully described in §8.1. The resulting electrons and positrons are easily detected because they are charged particles.

1.2.1

Instrumental Design

To optimize the detection and resolution of γ -rays, the EGRET instrument consists of a series of thin tantalum (Ta) sheets interleaved with planes of conducting wires spaced by 0.8 mm. Below this multilayer spark chamber is a NaI(Tl) calorimeter known as the TASC (Total Absorption Shower Counter). Surrounding the spark

CHAPTER 1. INTRODUCTION

6

ANTI-COINCIDENCE COUNTER

CLOSELY SPACED SPARK CHAMBERS

WIDELY SPACED SPARK CHAMBERS

TIME OF FLIGHT COINCIDENCE SYSTEM PRESSURE VESSEL

NaI(Tl) ENERGY MEASUREMENT COUNTER

Figure 1.3: The EGRET instrument. The total height of the spark chambers is approximately 1 m. chamber is a monolithic plastic scintillator to reject charged cosmic ray particles. The schematic design is shown in Figure 1.3. 99.5% of γ -rays will pass undetected through the anticoincidence scintillator into the 28 closely-spaced spark chamber modules, where it has a ?33% chance of conin the spark chamber along their trajectories. Below the closely-spaced modules is a time-of-?ight system designed to measure whether the particles are upward- or downward-moving. The system consists of two layers of 4×4 arrays of plastic scintillator tiles, spaced 60 cm apart. By combining measurements from the two layers, the time-of-?ight delay can be measured to within ?1.5 ns, and the general direction of the particles can be estimated. At any given time, a limited number of general direcorbital parameters, and are designed to exclude the Earth’s limb, as well as limiting the ?eld-of-view when desired. If the time-of-?ight system registers the passage of a downward particle from a valid region of the sky, and the anticoincidence system has not been triggered by a tions are considered valid for an instrument trigger. The valid directions depend on verting to an e? e+ pair. If it does so, the pair will ionize the (mostly) neon gas

CHAPTER 1. INTRODUCTION

7

Field of View E?ective Area >100 MeV Angular Resolutiona Energy Resolutionb Point Source Sensitivityc (photons cm?2 s?1 )

a b

SAS 2 0.25 sr 100 cm2 1? .5 ?100% 10?6

COS B 0.25 sr 70 cm2 1? .5 42% 10?6

EGRET 1.0 sr 1200 cm2 0? .6 18% 10?7

GLAST 2.6 sr ?7000 cm2 ?0? .1 10% <4 ×10?9 d

RMS at 500 MeV full-width, half maximum at 100 MeV c >100 MeV, 106 s exposure, unless noted d 1 year, high-latitude, >100 MeV, 5σ

Table 1.1: Performance of four high-energy γ -ray telescopes. Continually improving technology is re?ected in improving performance from the earliest instruments, SAS 2 and COS B , through EGRET and on to the proposed GLAST instrument. charged particle, a high voltage pulse is applied to the wires in the spark chamber modules. Ionized paths short the wires to ground, and the a?ected wires are recorded digitally in ferrite cores. The Total Absorption Spectrometer Calorimeter (TASC), located below the spark chamber, measures the total energy of the γ -ray event. In consists of 8 radiation lengths of NaI(Tl), and has an energy resolution of ?20% FWHM from a few tens of on-board clock to an absolute accuracy of 100 ?s and a relative accuracy of 8 ?s. The energy measurements made by the TASC are corrected on the ground for energy lost in the spark chamber and shower leakage. MeV to several tens of GeV. Events are tagged with an arrival time by the CGRO

1.2.2

Instrumental Calibration and Performance

Good calibration of the EGRET instrument was critical to the proper understanding of its data. The calibration was as extensive as its literature [193, 92, 93, 115, 150]; only the results will be stated here. There are three areas in which we will need to know the instrument performance: the point-spread function, or the distribution of the measured γ -ray incident angles as a function of the true incident angle; the sensitive (or e?ective) area, or the physical area for collecting γ -rays multiplied by

CHAPTER 1. INTRODUCTION

8

the e?ciency, as a function of position on the sky at any given time; and the energy dispersion, or the distribution of measured energy as a function of the true energy. These three functions were measured and recorded in tabular form as a function A reasonable approximation to the point-spread width assumes a relatively simple from a point on the sky may be taken as of aspect angle and energy. Their use in data analysis will be discussed in §2.2.

functional form. The half-angle which de?nes a cone containing ?68% of the γ -rays

θ68 = 5? . 85(E/100 MeV)?0.534 where E is the energy in MeV [193].

(1.1)

The sensitive area and energy dispersion are not easily expressible in functional form. Tables of their values were created in machine readable form, and analysis programs access them directly. The performance of the EGRET instrument compares very favorably with its predecessors, SAS 2 and COS B . The order of magnitude increase in e?ective area and improved point-spread function lead to the order of magnitude improvement in the point source sensitivity. A comparison of the telescopes, along with the proposed GLAST telescope, is shown in Table 1.1.

1.3

Successes with EGRET

The Compton Gamma-Ray Observatory has been very fortunate in successfully achieving and exceeding its design goals. Still operating some 7 years after launch, it has almost quadrupled its planned lifetime. While the entire observatory has been critical in advancing our understanding of astrophysics from 15 keV to 30 GeV (most notably the shocking revelation from BATSE that the mysterious γ -ray bursts are isotropically distributed on the sky), we will concentrate here on the contributions made by EGRET , with a view toward future advancements to be made by GLAST . EGRET has signi?cantly improved our understanding of pulsars [43, 191]. Six γ -ray pulsars have been identi?ed by EGRET , and their pulse periods have been

CHAPTER 1. INTRODUCTION

9

measured. EGRET observations of pulsars will be discussed in Chapter 4. Significant advancements have been made through observations of the Galactic [73] and extragalactic [212] di?use background emission. The largest number of identi?ed EGRET sources are the γ -ray blazars. Roughly 60 blazars have been identi?ed above 100 MeV, leading to new insights into blazar emission mechanisms [62]. While the BATSE γ -ray burst measurements were the most revolutionary discovery, EGRET has made signi?cant contributions to our understanding of γ -ray bursts at the highest energies [122]. γ -Ray bursts will be discussed in Chapter 3. CGRO has yielded a wealth of information about the γ -ray sky. Much of that information has directly increased our understanding of astrophysical systems. In particular, EGRET has given us an unprecedented view of the high-energy γ -ray sky. EGRET has identi?ed a great number of new sources; the launch of GLAST will give us a tool to understand what exactly it is that we have found.

Chapter 2 Statistical Methods in γ -Ray Astronomy

Optical astronomy has long been famous for breathtaking images of distant galaxies, star-forming clouds, and beautiful nebul?. While γ -ray astronomy can produce equally beautiful results, the nature of the photons and the instrument yield data which must be analyzed very di?erently than data from other wavelengths. The ?rst major di?erence is the sparsity of the photons; integration times of days are usually required to observe all but the brightest sources. A multitude of astrophysical conditions also a?ect the nature of the data. Cosmic ray interactions with Galactic dust and gas produce di?use background γ -rays. Di?erent energy generation mechanisms produce di?erent spectral pro?les across the energy range observed by EGRET . Several speci?c instrumental responses also shape the data analysis process. A ?nite point-spread function means that γ -ray directions will be determined to within a statistical distribution around the true source direction. Instrumental sensitivity to three orders of magnitude in energy across the ?eld of view of more than a steradian is not uniform as the platform orbits at 17,000 miles per hour. Meanwhile, our estimates of all these parameters depend on the γ -ray energy, which is only known to limited accuracy. Despite these impediments to observation, EGRET has been very successful in

10

CHAPTER 2. STATISTICAL METHODS IN γ -RAY ASTRONOMY

11

Figure 2.1: All γ -rays above 100 MeV measured by EGRET in phase I, II, and III, binned at a scale of 0? . 5. making γ -ray observations. This chapter will examine the nature of the EGRET data, and the statistical methods used to analyze it.

2.1

The Nature of the γ -Rays

Several major features of the high-energy γ -ray sky are evident from the simplest possible examination. Figure 2.1 shows the raw photon counts for the whole sky in 0? . 5 ×0? . 5 bins in Galactic coordinates. The Galactic center is evident at the center of the map, as is the Galactic disk. There are several bright spots in the plane, as well as some evident high-latitude sources. We will ?rst examine the di?use background, then see how statistical methods along with understanding of the background will help us identify point sources and estimate their locations and ?uxes.

CHAPTER 2. STATISTICAL METHODS IN γ -RAY ASTRONOMY

12

2.1.1

Di?use Background

Almost two decades before the launch of SAS 2 , Hayakawa [66] predicted that highenergy γ -rays would be produced as cosmic rays interacted with interstellar gas yielding pions, which would decay, directly or indirectly, into γ -rays. Indeed, experience with SAS 2 , COS B , and EGRET has shown strong correlations of the di?use γ -ray background at low Galactic latitudes with known Galactic structural features such as the spiral arms [39, 63, 119]. Based on indications that the di?use γ -ray ?ux observed by SAS 2 and COS B was approximately the intensity and shape expected from cosmic ray interactions with interstellar matter, a model of the di?use ?ux was made for the purposes of EGRET data analysis [10, 75]. There are three main processes by which di?use γ -rays are generated. The domiray protons interact with dust particles or gas. These pions then decay to high-energy γ -rays. Below ?70 MeV, γ -rays can be produced by either bremsstrahlung of cosmic rays in interstellar clouds or by inverse-Compton upscattering of low-energy photons. A good model of the Galactic di?use γ -ray intensity thus requires a good model of the distribution of interstellar matter in the Galaxy, and a good model of the cosmic ray ?ux throughout the Galaxy. The ?rst may be well approximated using maps of the distribution of hydrogen, which comprises most of the interstellar matter in the Galaxy. Atomic hydrogen has been carefully mapped with observations of the 21 cm hyper?ne transition line. Molecular hydrogen is more di?cult to map, but may be approximated by assuming that CO is a good tracer, easily identi?ed by its 2.6 mm emission line. A constant ratio of CO to molecular hydrogen is typically assumed throughout the galaxy. Much more di?cult is the estimation of the Galactic cosmic ray ?ux. Since the ?ux cannot be directly measured, assumptions about the distribution of cosmic rays must be made. For the purposes of EGRET analysis, Bertsch [10] and Hunter [75] assume that cosmic rays are in dynamic equilibrium with the interstellar magnetic pressure and the gravitational pressure of the galactic matter. These assumptions, convolved with the instrument point-spread function (§2.4), result in the map of diffuse Galactic γ -rays shown in Figure 2.2, which is in good agreement with the observed nant process above ?70 MeV is the decay of pions. Pions are produced when cosmic-

CHAPTER 2. STATISTICAL METHODS IN γ -RAY ASTRONOMY

13

Figure 2.2: The assumed Galactic contribution to the di?use γ -ray background, taken from measurements of atomic and molecular hydrogen, and convolved with the EGRET point-spread function. This is a map of Gij , as de?ned in equation (2.20). intensity map shown in Figure 2.4. Of course, at high galactic latitudes the di?use background is primarily extragalactic [182], although there is signi?cant galactic diffuse background at high latitude. Some or all of this extragalactic background is due to a large number of weak sources, while some of it may be due to a truly di?use background [212].

2.1.2

Spectral Di?erences

EGRET ’s broad energy range allows the spectra of sources and the di?use background to be measured. The spectrum of almost every EGRET source is well ?t by a power law. Since γ -rays are so sparse, the power laws are usually quoted as a function of the number of photons, instead of a function of energy as is sometimes used for other wavelengths. The form is then I (E ) = I0 E ?α photons cm?2 s?1 MeV?1 (2.1)

CHAPTER 2. STATISTICAL METHODS IN γ -RAY ASTRONOMY

14

where I (E ) is the di?erential photon ?ux. The spectral index α is close to 2.0 for most sources, though it can be as low as 1.42 for pulsars. The Galactic di?use background has a bit softer spectrum—about 2.1.

2.2

Instrumental E?ects

The quality and nature of the data depend equally on the photons and the instrument which observes them. Any attempt to analyze the data must consider the instrumental response as the basis for an analysis method. The three main aspects of the instrumental response that we must consider are the point-spread function, the sensitive area, and the energy dispersion. There is no reason, a priori, to expect that the instrument response function should separate cleanly into these three functions. However, it o?ers great simpli?cation to the data analysis, and in practice seems to be a good approximation.

2.2.1

Point-Spread Function

It is important to be able to quantify the ability of a γ -ray telescope to correctly reconstruct the true incident direction of a γ -ray. To precisely de?ne this, we will distinguish between the “point-spread density” and the “point-spread function.” The point-spread density, or PSD, refers to the probability density distribution of incident γ -ray directions measured by the instrument from a point source. This distribution may in general be a function the true position of the point source (the inclination and azimuth relative to the centerline of the telescope) and the energy of the γ -ray: PSD = PSD(θ, φ; θ0 , φ0 , E0 ) (2.2)

where θ0 , φ0 represents the true source position and E0 is the true γ -ray energy. The apparent γ -ray direction is (θ, φ). Often, we will require a probability, as opposed to a probability density. The di?erential probability of measuring a photon in a di?erential

CHAPTER 2. STATISTICAL METHODS IN γ -RAY ASTRONOMY

15

element dθ dφ is given by dP = PSD(θ, φ; θ0 , φ0 , E0 ) sin θ dθ dφ (2.3)

The point-spread function is the di?erential probability, often integrated over azimuth and energy to yield a function of inclination, as in equation (2.6). The EGRET PSD was measured at the Stanford Linear Accelerator Center in 1986. A beam of electrons with tunable energies between 650 MeV and 30 GeV was back-scattered o? pulsed laser photons. γ -Rays were produced between 15 MeV and 10 GeV by inverse-Compton scattering [115]. The point-spread density was measured as a function of apparent γ -ray position for 10 discrete energies, 5 inclination angles (θ0 ) and 3 azimuthal angles (φ0 ). The resulting tables yield the relative probability of detecting a photon at (θ, φ) assuming values of the other three parameters. In addition, EGRET operates in up to 87 di?erent “modes,” corresponding to di?erent triggering criteria.? These modes are designed to maximize the operating ?eld of view, even when part of the geometric ?eld of view is obscured by the Earth or its limb.

2.2.2

Sensitive Area

A second function which is clearly critical for data analysis is the sensitive area of instrument. The sensitive area (or e?ective area) is the projected area of the detector multiplied by its e?ciency. It too is a function of incident γ -ray parameters: SA = SA(θ0 , φ0 , E0 ) (2.4)

Clearly, the sensitive area of the instrument is also dependent on the instrument mode.

These di?erent modes correspond to allowing only photons from certain broad regions of the sky as de?ned by coincidence of di?erent combinations of time-of-?ight tiles.

?

CHAPTER 2. STATISTICAL METHODS IN γ -RAY ASTRONOMY

16

2.2.3

Energy Dispersion

Finally, the analysis must consider energy dispersion. The energy dispersion function gives the distribution of measured energy for a given true energy. The measured energy varies from the true energy because of noise in the photomultipliers, ?uctuations in the shower leakage from the calorimeter, and incomplete correction for energy losses elsewhere in the instrument. ED = ED(E ; E0 , θ0 , φ0 ) (2.5)

Taken together, these three functions yield the point-spread width approximation in equation (1.1) in the following way. We ?rst notice from the calibration data that the point-spread function is roughly azimuthally symmetric, and that it does not vary widely with the true inclination angle. Then we can ?nd PSF(θ) = 2π N

Emax E =Emin ∞ E ′ =0

E′

?α

PSF(θ, E ′ )ED(E ′ , E )dE ′dE

(2.6)

where E is the measured γ -ray energy, E ′ is the true γ -ray energy, and N is a normalization factor. The deviation between the apparent incident angle and the true incident angle is θ. The point-spread function PSF(θ) is the integral of the trueenergy dependent point-spread function, weighted by the spectrum, integrated over the measured energy band from Emin to Emax , and integrated over all true energies, weighted by the energy dispersion function. This re?ects the fact that there is some probability that a γ -ray of any given true energy will have a measured energy between Emin and Emax . This integral was done numerically for a number of energy bands, and a Gaussian ?t to the results led to equation (1.1) [193].

CHAPTER 2. STATISTICAL METHODS IN γ -RAY ASTRONOMY

17

2.3

Likelihood Analysis

The sparsity of astrophysical γ -rays and the complicated instrumental response of the EGRET instrument suggests statistical data analysis that functions at the photonby-photon level, taking into account backgrounds and the instantaneous instrument state to extract the most information from each photon. Early analyses of COS B data were based on a cross-correlation method [69]. However, this method could not easily handle the highly structured background that is typical of high-energy γ -ray astrophysics. Later, a maximum likelihood technique was brought to bear on COS B data with much greater success [164]. Based on this success, maximum likelihood techniques were adopted for use with EGRET [117]. The central idea of likelihood analysis is very simple. Given a set of models, we wish to ?nd the model which is most likely to be responsible for the observed data. The likelihood is de?ned as the probability of the observed data, given a choice of model. The likelihood is written as follows: L(D |M ) (2.7)

where D represents the observed data, and M the model. Quantities to the right ability. The maximum likelihood method determines the best model by maximizing this likelihood function. A special, but very common, case is that of a parameterized model. For example, consider that we wish to measure the ?ux of a γ -ray source. We will imagine an idealized detector with 100% e?ciency and an angular area of α sr from the source. The number of γ -rays emitted in a unit time is Poisson distributed, so the likelihood of measuring n′ photons from a source with intensity ? is given by α e?? ?n L (D | ? ) = 4π n! (2.8) of the | sign are taken as given and ?xed, making the likelihood a conditional prob-

where n is the number of photons emitted from the source, and n′ = (α/4π )n. To ?nd our best estimate of ?, we maximize this likelihood. In fact, we will ?nd it more convenient (and mathematically equivalent) to maximize the logarithm of the

CHAPTER 2. STATISTICAL METHODS IN γ -RAY ASTRONOMY

18

likelihood. ln L(D |?) = ?? + n ln ? + ln the derivative of equation (2.9) to zero, we ?nd 0 = ?1 + n/? =? ? = n =

α 4πn!

(2.9)

The last term is constant, and thus for maximization purposes can be ignored. Setting

4π ′ n α

(2.10)

Unsurprisingly, in this simple example, we ?nd that the number of photons measured per unit time divided by the subtended angle fraction is the best estimate of the ?ux. One other example is illustrative. Let us assume we have some Gaussian process with a constant mean and a known variance σ 2 . A series of observations x is made, and the most likely mean value ? is desired. The likelihood function is L(x|?) = Taking the logarithm, we have ln L(x|?) = ?1/2 (xi ? ?)2 σ2 (2.12) e?

i

(xi ??)2 2σ 2

(2.11)

i

We notice that ?2 ln L is formally equivalent to χ2 . In fact, this is an example of a general result: in the limit that all distributions involved are Gaussian, the maximum likelihood result is the same as the χ2 minimizing result. This is an example of a more general result known as Wilks’ Theorem [211]; it will be described more thoroughly in §2.5. Maximum Likelihood Con?dence Regions. Maximum likelihood methods also yield con?dence regions. Following Eadie, et al. [36], we ?rst consider the case of a Gaussian distribution with unit variance and unknown mean ?. The likelihood is

L(? x|?) = (2π )?N/2 exp ?

1 N (xi ? ?)2 2 i=1

CHAPTER 2. STATISTICAL METHODS IN γ -RAY ASTRONOMY

19

= (2π )?N/2 exp ?

N 1 N (xi ? x ? )2 (? x ? ?)2 exp ? 2 2 i=1

(2.13)

ln L ≥ ?1/2. This corresponds to the interval ?1 ≤ ? ? x ? ≤ +1. From the properties of the normal distribution, we know that this must contain 68.3% of the distribution. This would be merely a curiosity if not for the following. Suppose that the likeliSimilarly, the interval corresponding to ln L ≥ ?2 contains 95.5% of the distribution. hood function is a continuous function of ?, with only one maximum. In that case, we may reparameterize our observed variable in terms of a new variable g (?) such that the likelihood as a function of g is parabolic. We may now ?nd the con?dence region in g as we did above. Given that region, we may invert the reparameterization

The ln L is thus a parabola in ? of the form ? N (? ? x ?)2 . In the case that N = 1, let 2

to ?nd a con?dence interval in ?. Furthermore, we notice that the function is the same, whether it is parameterized by ? or g . That is, ln L(? x|?) = ln L(? x|?(g )) (2.14)

Thus, we can ?nd the con?dence region in ? directly by determining the point at which the ln L has decreased by 1/2 for 68% or 2 for 95.5%, without ever ?nding likelihood as a function of g is a Gaussian, but it is not necessarily central in ?. the reparameterization g . However, note that the interval is central in g since the

2.4

Applying Likelihood Analysis to EGRET

We can see now how to proceed in EGRET data analysis. Take all the data accumulated in some time period. The likelihood of that data is the product over di?erential elements of angle and energy of the Poisson probability density of detecting the photons, given the rate of photons times the probability of detecting each photon. The rate ? can be expressed as ?(?, b, E ; ?0, b0 ) =

∞ E =0

(2.15)

[I (?0 , b0 )PSF(?, b, E ′ ; ?0 , b0 ) + B (?, b)] SA(?, b, E ′ )ED(E ′ , E )dE ′

CHAPTER 2. STATISTICAL METHODS IN γ -RAY ASTRONOMY

20

The rate is a function of the measured energy and position on the sky. The source intensity I depends only on the true source position (?0 , b0 ). The background B is assumed to have already been convolved with the point-spread function. The three instrument functions depend on the instrument mode m as well; this dependence will be suppressed for clarity. The likelihood is then the product over all di?erential parameter elements of the Poisson probability of measuring (or not measuring, as was the case) a photon in that element, given the rate in that element from equation (2.15). The appearance of derivatives of products encourages us to use the logarithm of the likelihood. Denoting the integrated rate over measured energies as ? ?= likelihood becomes ln L(?0 , b0 ) = [?? ?(?, b; ?0 , b0 ) + n ln ? ?(?, b; ?0 , b0 )] d? db (2.16)

Emax Emin

?(?, b, E ) dE , the log

?

b

where n is the number of photons observed in a di?erential element d? db. Since the element is di?erential, this must be either 1 or 0. The integral thus divides into an integral and a sum:

N

ln L(?0 , b0 ) = ?

?

b

? ?(?, b; ?0 , b0 )d? db +

i=1

ln ? ?(?i , bi ; ?0 , b0 )

(2.17)

where the sum is evaluated for the parameters of each photon. While this looks fairly simple, ? ? is quite a complicated object, implicitly containing four integrals. We would like to maximize equation (2.17) over ?0 and b0 , the source position. Computationally, this is a herculean task, requiring evaluation of a six-dimensional integral at each trial point (?0 , b0 ). Clearly, some simpli?cation is necessary.

2.4.1

Binned Likelihood

The most obvious simpli?cation is to give up on individual photons, and create binned maps. While binning is always undesirable, as it loses information in the data, in this case binning is minimally undesirable. Our derivation above considered di?erential elements in the parameters; essentially, we let the bin size go to zero. It has been

CHAPTER 2. STATISTICAL METHODS IN γ -RAY ASTRONOMY

21

shown that using a ?nite but small bin size speeds computation dramatically at very little expense of accuracy [199]. A standard analysis program called LIKE was developed for the map-based likelihood analysis of EGRET data [117]. LIKE considers the total expected γ -ray rate in each pixel, typically 0? . 5 ×0? . 5 on the sky.? This rate is the sum of the expected rates from the isotropic background, the galactic background, and a point source. In order to estimate rates for binned maps, the photons from any desired observation interval are binned as in Figure 2.1. In addition, the instrument exposure to the sky must be calculated. This is a function of the amount of observing time and the sensitive area during that time. Let T (θ, φ, m) be the amount of observing “livetime,” (that is, elapsed viewing time that the instrument was active, excluding occultations and instrument dead time) for a location (θ, φ) on the sky in an instrument observing mode m between time t1 and t2 . Only the total observing time in each mode is relevant; it need not be contiguous. Then the exposure E (?E ; θ, φ) for a given measured energy range ?E to a point (θ, φ) on the sky is E (?E ; θ, φ, t1 , t2 ) = where ?(?E ; θ, φ, m) = A

Emax E =Emin ∞ E ′ =0

?(?E ; θ, φ, m) T (θ, φ, m; t1 , t2 )A

m

(2.18)

E′

?α

SA(E ′ ; θ, φ, m)ED(E, E ′ )dE ′ dE

(2.19)

This exposure explicitly depends on the spectral index of the source or background, whichever is being observed. This is the ?rst example of a continuing di?culty; the spectral index is required to calculate the exposure, but it is not known until after the ?ux is determined. In principle, the spectral index should be allowed to vary everywhere during the maximization process. An approximation would be to iterate. However, since most spectral indices are very nearly 2.0, and varying the index

? This scale is somewhat arbitrary. It was based largely on the scale to which the Galactic hydrogen has been mapped, as well as to ensure su?cient photons in each bin. It was not chosen to optimize the likelihood method. Nevertheless, for most energies the bin size is much smaller than the instrument point-spread width.

CHAPTER 2. STATISTICAL METHODS IN γ -RAY ASTRONOMY

22

Figure 2.3: Combined EGRET exposure in phase I, II, and III. Phase I was an all-sky survey, with roughly even exposure. Phases II and III consisted of pointed observations for sources of interest. Lighter shade corresponds to more exposure. requires recalculating the exposure maps (a computationally expensive prospect), a constant index is assumed. The exposure can be calculated for each pixel in the standard map; an example is shown in Figure 2.3. To allow for the degrees of freedom in the background maps described above, the background is assumed to be a linear combination of isotropic and galactic background whose coe?cients will be optimized. The isotropic di?use count estimate for cm?2 s?1 sr?1 , and Eij is the instrument exposure to pixel i, j in units of cm2 s sr. a given bin is gb Eij , where gb (“g-bias”) is the coe?cient to be ?t in units of photons The Galactic component of the rate will depend on both the di?use map and the

instrument point-spread function, as well as the exposure. Denoting the binned radio di?use gas map as fij , the rate due to the Galactic di?use background is Gij =

kl

fkl Ekl PSF(θij,kl ) kl PSF(θij,kl )

(2.20)

where θij,kl is the angular distance between pixel i, j and pixel k, l. The denominator is a normalization factor, necessary since the sums over k and l may not be over

CHAPTER 2. STATISTICAL METHODS IN γ -RAY ASTRONOMY

23

the entire sky. Pixels far from the point of interest will contribute negligibly to the estimation of source ?ux and location, but it is critical to good ?ux measurements to have a normalized point-spread function. Therefore, we allow analysis of only the pixels within some radius of analysis Ranal , but renormalize the point-spread function. In addition, to allow for the unknown cosmic ray ?ux and the CO/H2 ratio, we will allow a constant multiplier gm (“g-mult”) to be optimized as well. So, then, given k sources in the ?eld of view, the expected number of counts in bin i, j is ?ij = gm Gij + gb Eij + ck PSF(θij,k )

k

(2.21)

where θij,k is the angular distance between the position of source k and pixel i, j . This count estimate is a function of 3k + 2 parameters: the k source strengths, latitudes, and longitudes; gm , and gb . The ?t values of gm are consistently in the range 0.92–1.08. The ?t level of gb is usually around 2 ×10?5 photons cm?2 s?1 sr?1 [43].

2.4.2

Maximizing the Likelihood with LIKE

Just as in the exact case, the distribution of photon counts in each bin is Poisson. The likelihood of a given map, then, is: L(D |ck , ?k , bk , gm , gb ) = ?ijij e??ij nij !

n

(2.22)

ij

with ?ij given by equation (2.21). The maximum likelihood estimates of ck , ?k , bk , gm , and gb are simultaneously solving the the set of equations ? = 0 ln L(?) ?ck ?=? ? ? = 0 ln L(?) ??k ?=? ? ? = 0 ln L(?) ?bk ?=? ? ? = 0 ln L(?) ?gm ?=? ? (2.23) (2.24) (2.25) (2.26)

CHAPTER 2. STATISTICAL METHODS IN γ -RAY ASTRONOMY

24

Figure 2.4: γ -Ray intensity as measured by EGRET . The counts map in Figure 2.1 is divided by the exposure map in Figure 2.3 to yield the γ -ray intensity. ? = 0 ln L(?) ?gb ?=? ? (2.27)

where ? ? is optimized to satisfy all these conditions. This is just the multi-dimensional maximization of the likelihood function—there are 3k + 2 equations to solve simultaneously. For this reason, it is usually computationally much faster (although less accurate) to ?nd the brightest source ?rst, then ?x its parameters and ?t the next source. Once the locations of the sources are established, equation (2.21) divided by the exposure map yields an intensity map of the sky in Figure 2.4; equivalent to the optical image observed simply by peering through the eyepiece of an optical telescope.

2.5

Parameter Estimation vs. Hypothesis Testing

The foregoing discussion concerned the estimation of the position and ?ux of a source which was assumed to exist. The maximum likelihood technique not only allows for parameter estimation as described above but also for hypothesis testing. Hypothesis testing estimates the signi?cance of the source detection. Highly signi?cant sources are reliable; low signi?cance sources may be statistical ?uctuations.

CHAPTER 2. STATISTICAL METHODS IN γ -RAY ASTRONOMY

25

Consider two hypotheses, M0 and M1 . For concreteness, let us assume that M0 , the null hypothesis, states that there is no source in the ?eld of view. Further, we assume that M1 states that there is a source with non-zero ?ux in the ?eld of view. If there are n degrees of freedom in M0 , and m degrees of freedom in M1 , then asymptotically distributed as χ2 (m?n) : ? 2 ln Wilks’ Theorem [211] states that ?2 times the logarithm of the likelihood ratio is L(D |M0 ) ? χ2 (m?n) L(D |M1 )

(2.28)

In the case of a source whose position is known from other experiments, there are 3 degrees of freedom in M1 (ck , gm , gb ) and 2 degrees of freedom in M0 (gm , gb ). Thus, the likelihood is distributed as χ2 1 . We de?ne the EGRET test statistic TS ≡ ? 2(ln L0 ? ln L1 ). From the de?nition of the χ2 distribution, this implies that the √ signi?cance of a detection is given by TSσ in the standard nomenclature. Thus, a

TS = 16 source is detected at the 4σ level. Of course, that is for a single trial, and in fact many sources observed by other telescopes have been searched for with EGRET . The case of previously unknown sources is a little more complicated. M1 now has 5 degrees of freedom: ck , ?k , bk , gm , gb . But the null hypothesis is degenerate in the position—the likelihood is independent of position when the ?ux is zero. In that case, Wilks’ Theorem does not hold, so the distribution is not known. Some theoretical work has been done to determine this distribution [198]. Based on Monte Carlo simulations, EGRET sources are accepted as detections if, for sources at Galactic √ √ latitudes |b| > 10? , TS > 4, and for sources at Galactic latitudes |b| < 10? , TS > 5. one spurious excess with TS > 16 will be detected in an analysis of the entire sky [117, 43].

Monte Carlo simulations have shown that for a perfect background model, roughly

2.6

Bayesian Methods

Despite their formal similarity, conceptual di?erences between maximum likelihood analysis and Bayesian methods have kept the latter relegated, for the most part, to the

CHAPTER 2. STATISTICAL METHODS IN γ -RAY ASTRONOMY

26

statistical backwaters of astrophysics. Some recent work using Bayesian methods has yielded useful results [174, 110, 111]; since Bayesian methods will be an appropriate alternative framework for later chapters, a brief overview will be presented here. Traditional, or frequentist, statistics are predicated on the notion that extreme values of some function of the data given a null hypothesis indicate that the null hypothesis is probably wrong, and therefore the test hypothesis is probably true. Much of the confusion of the general population about “Statistics” is due to the two twists involved in frequentist analysis: ?rst of all, the goal is to disprove the thing we suspect false, rather than ?nd evidence for what we suspect true; and second, this is done by considering all possible outcomes from an ensemble of data sets. In to the distribution expected from the ensemble of data sets that might be generated if the null hypothesis were true. If the measured value of the likelihood is extreme, according to this distribution, we claim that the null hypothesis has been excluded to some con?dence level. the case of likelihood statistics, we calculate the statistic (?2 ln L) and compare it

2.6.1

Bayes’ Theorem

In contrast, the Bayesian method demands that we stay at all times in the realm of probability: speci?cally, the probability that the test hypothesis is true. To develop the mathematics, let us begin with what Scargle [174] calls the “obviously true” form of Bayes’ Theorem: P (M | D )P (D ) = P (D | M )P (M ) (2.29)

Formally, this follows from the de?nition of conditional probabilities. The probability of some statement M being true, given that another statement D is true is the joint probability of M and D . Any joint probability may be expressed as the probability of one statement times the conditional probability of the other. Or stated another way, the probability of A and B equals the probability of B and A. Thus Bayes’ Theorem is proved. Clearly, the notation in equation (2.29) is suggestive. As we have identi?ed above, the probability of the data, given a model, is known as the likelihood of the data.

CHAPTER 2. STATISTICAL METHODS IN γ -RAY ASTRONOMY

27

Equation (2.29) points out that while we have been calculating likelihoods, what we that a given model is true? This is the question which Bayesian methods set out to answer. To that end, we rearrange equation (2.29) into a more useful form. P (M | D ) = L (D | M )P (M ) P (D ) (2.30) really desire is P (M |D ); that is, given the data that we have, what is the probability

Equation (2.30) gives us exactly what we want: the probability of a model being true, given the data that we have observed. We do not resort to any hypothetical ensembles of data, and more importantly, we make direct statements about the model in question. one of a discrete series, or it may be parameterized. In the case that the model is a function of a continuous parameter, then the left side of equation (2.30), known as the posterior probability, becomes a function P (M (θ)|D ). It is interpreted as the probability that the true value of the parameter θ0 is given by θ. Given the importance of the posterior probability, it behooves us to understand the right side of equation (2.30). The likelihood has been fully discussed. P (M ) is known as the prior probability of the model. Except (apparently) in the case of quantum mechanics, probabilities are used to represent ignorance. We do not know the true value of a parameter, or which model is correct, so we assign probabilities to represent the knowledge that we do have. P (M ) represents the knowledge we have about the system before we receive the data D . A common example is that of a parameterized model in which we know that the true value of the parameter θ lies somewhere between θmin and θmax ? . The prior is then ?at over the interval [θmin , θmax ] and zero elsewhere. In other situations, it may be more appropriate to take a scale-invariant prior, which would have a logarithmic form. Such a prior is generally referred to as a “least informative prior.” It re?ects only information about the structure of the experiment, and does not favor any speci?c outcome. Of course,

?

Analogously to our likelihood calculations in §2.3, the model in question may be

It will often be possible to take the limit of our results as θmin → ?∞ and θmax → ∞.

CHAPTER 2. STATISTICAL METHODS IN γ -RAY ASTRONOMY

28

if there is speci?c information about the true parameter value, the prior should re?ect that information. Finally, we must address the denominator of equation (2.30). This term, the probability of the data, serves as a normalization. It expresses the probability of measuring the data regardless of which model is actually true. If this can be rigorously calculated, then the left side of equation (2.30) will be a well-normalized probability distribution. In practice, it is often more practical to sidestep the issue by forming the odds ratio.

2.6.2

The Odds Ratio

Usually it is the case that we compare two discrete models, or that we compare a discrete model with a class of models characterized by a ?nite number of continuous parameters. In that case, the odds ratio conveniently handles our normalization issues. The odds ratio is formed by comparing the posterior probabilities of the di?erent models. Consider two discrete models, M0 and M1 . For concreteness and comparison with §2.5, let us assume that M0 states that there is no source, and that ratio is P (M0 |D ) L(D |M0 )P (M0 ) = P (M1 |D ) L(D |M1 )P (M1 ) M1 states that there is a source of a given ?ux at a given position. Then the odds (2.31)

This gives the probability that there is no source relative to the probability that there is a source. Note that it only allows those two possibilities. A value of 1/3, for example, would mean it was three times more likely that there was a source than that there was no source. We know nothing about the probability that there were two sources. M1 (?) was a class of models, parameterized by the source strength. The odds ratio in that case is also a function of the parameter: O (? ) = P (M1 (?)|D )P (M1(?)) P (M0 |D )P (M0) (2.32) Of course, the example in §2.5 was more complicated than this. The model M1 =

CHAPTER 2. STATISTICAL METHODS IN γ -RAY ASTRONOMY

29

The function O (?) gives the odds that there is a source of strength ? at a given position versus that there is no source. In most cases, we will be interested in the total odds ratio; that is, the odds that there is a source regardless of its strength. ratio using the procedure known as marginalization. This is akin to the detection signi?cance calculated in §2.5. We may calculate this

2.6.3

Marginalization and Con?dence Regions

We may eliminate uninteresting parameters by marginalization. This process acquired its odd name through a historical accident, when “integrations” were carried out numerically by adding columns of numbers into the margins of the page. In essence, the method is mathematically simple. Given any conditional probability, and the probability of the condition, we may integrate to eliminate the condition: P (A) =

Bmax B =Bmin

P (A|B )P (B )dB

(2.33)

The application to the odds ratio is immediately clear. We simply integrate the numerator over all possible ? to ?nd the total odds ratio. A very similar process may be used to ?nd con?dence intervals. Consider the situation when a source is known to exist. We then wish to ?nd the best estimate of its ?ux, and a con?dence interval for that estimate. Since the ?ux ? must take on some positive value, we may evaluate the denominator of equation (2.30): P (D ) =

∞ 0

L(D |M (?) )P (?)d?

(2.34)

Then the probability that the true value of ? is between ?? and ?+ is P (? ? ≤ ? 0 ≤ ? + ) =

?+ ?? L(D |M (?) )P (?)d? ∞ 0 L(D |M (?) )P (?)d?

(2.35)

For a given con?dence level, there are an in?nite number of choices of con?dence intervals, as there are with frequentist statistics. There is no requirement that con?dence intervals be contiguous. However, useful intervals can often be found by requiring

CHAPTER 2. STATISTICAL METHODS IN γ -RAY ASTRONOMY

30

resulting interval is central by de?nition, and in the limit that the probability density it is inside the interval, it is minimal.

?? = ? ? ? δ? and ?+ = ? ? + δ?, where ? ? is the most likely value of ? [36]. The

is symmetric about the maximum and is smaller everywhere outside the interval than

2.6.4

Advantages and Disadvantages

Advocates of frequentist and Bayesian methods have unfortunately been polarized into two extreme camps, with only the most acrimonious communication between them. In fact, both methods are rather like a powerful hunting ri?e. Used properly, they are e?cient and successful at doing their job. Improper use may result in permanent catastrophic injury and/or death. We will brie?y examine the objections to both methods, and in doing so ?nd that the disagreements are actually objections to using the methods improperly. The major objection to the Bayesian method is the use of a prior. It is said that the prior subjectivizes what should be an objective procedure, and therefore reduces the results to a sort of “modi?ed best guess” of the experimenter. The Bayesian responds that that is exactly true; indeed, if we use probabilities to express our ignorance, then we should hope that the results of an experiment reduce our ignorance. Furthermore, the Bayesian claims that all assumptions and prior knowledge are made explicit under the Bayesian formulation. It is certainly clear that the prior is the Achilles’ heel of the Bayesian method. There is no objective, prescribed method for obtaining a prior. However, there is also no objective, prescribed method for choosing a traditional statistic. The use of χ2 , for example, implicitly assumes that the errors in the measurements are Gaussian, which may or may not be the case. The primary objection to traditional statistics is the ad hoc procedure of selecting a statistic. Statistics are chosen based on their power and appropriateness, to be sure, but the choice is also largely guided by experience. The justi?cation for the choice of statistics is generally its past success, rather than any a priori reasoning. In contrast, the Bayesian method prescribes the calculation for any experiment: form the likelihood and weight by the prior.

CHAPTER 2. STATISTICAL METHODS IN γ -RAY ASTRONOMY

31

A secondary objection is philosophical in nature. Bayesians prefer to treat the true value of the parameter as a random variable, with the experimental data as the ?xed and unchanging measure of reality. The parameter then has some probability of falling within the con?dence region. The traditional approach treats the true value of the parameter as ?xed and unchanging, and the data as only one possible outcome in an imagined ensemble, a shadow or projection of the true reality. The con?dence region that we calculate from this instantiation of the data then has some probability of covering the true value of the parameter.

§

Despite these objections, the two methods are actually compatible, and under the right circumstances, equivalent. In most circumstances, the maximum likelihood method is equivalent (for parameter estimation) to the Bayesian result with a ?at prior. Whenever the implicit assumptions of a traditional analysis are matched by the explicit assumptions of a Bayesian analysis, the results will be the same.

2.7

Calculating Upper Limits

A signi?cant positive detection of a source is the goal of all telescopes. However, null results can also be useful [133, 200]. It is often valuable to set an upper limit on the γ ray ?ux of a source known from X-ray, optical, or radio observations. Determining the value of the upper limit is a statistical endeavor that requires a very careful de?nition of the goal. It has been noted that “the question of how to calculate an upper limit in the vicinity of a physical boundary is one of the most divisive in high-energy physics” [5]. Unfortunately, this is precisely the limit in which we ?nd ourselves. A negative source ?ux is unphysical; nevertheless, a measurement of a weak (or non-existent) source in the presence of a large background may easily result in a data set best ?t with a negative source ?ux. In order to combine such a result with other results, the ?ux must be reported as negative, with the con?dence region found as in §2.3

The philosopher will note a certain correspondence to historically important epistemological viewpoints. Sartre may debate Plato on the true nature of reality, but both are crushed by the rock of Sisyphus.

§

CHAPTER 2. STATISTICAL METHODS IN γ -RAY ASTRONOMY

32

under the maximum likelihood method [5]. Otherwise, any combination of results from di?erent observations or di?erent experiments would be biased. Once we have agreed upon a point estimate of the parameter, we must consider the con?dence region that we wish to comprise the upper limit. In analogy to the con?dence regions around point estimates, we will de?ne the “1σ upper limit” (in frequentist terms) to be the top of an interval which will, when constructed from an ensemble of data sets, contain the true value of the ?ux 68.3% of the time. In Bayesian parlance, this means that the integral of the posterior probability distribution from zero to the 1σ upper limit will be 68.3%.

2.7.1

Upper Limits from LIKE

The upper limits generated by LIKE do not ful?ll this de?nition. There are three situations for which LIKE must generate an upper limit. The ?rst is when the ?ux measurement is positive, but the con?dence region extends to negative ?ux; e.g., but the con?dence region extends into positive territory. The third is when the ?ux measurement and the entire con?dence interval are negative. LIKE handles these situations in the following way. An upper limit is always quoted, and the ?ux measurement is quoted only if certain conditions are met. If TS > 1 and the ?ux is positive, the ?ux and con?dence regions are quoted. The upper limit is the top end of the con?dence interval. If TS < 1 or the ?ux is negative, only an upper limit is quoted. (Note that this immediately introduces the bias discussed above.) If the ?ux is positive, the upper limit is the top end of the con?dence interval. If the ?ux is negative, LIKE ?nds the width of the con?dence region, then shifts it so that it is centered on zero. The upper limit is then the top end of the shifted con?dence interval. The upper limit calculated for strong sources with positive ?ux clearly does not ful?ll our de?nition. The con?dence interval on the point estimate will contain the true value 68.3% of the time. The interval between zero and the upper limit is larger a measurement of 5 ± 10. The second is when the ?ux measurement is negative,

CHAPTER 2. STATISTICAL METHODS IN γ -RAY ASTRONOMY

33

than the con?dence interval, and so must contain the true value at least 68.3% of the time. The upper limit for sources with positive ?ux whose con?dence intervals extend to negative ?uxes also does not ful?ll our de?nition. The 68.3% upper limit is found from the one-sided con?dence interval; that is, the value σu.l. such that

σu.l. ?∞

f (? ?|?0 )d? ? = 68.3%

(2.36)

a Gaussian function f , the con?dence of the integral evaluated to σu.l. = ? ? + σ may ? ? + σ depends strongly on the shape of f . The most di?cult situation is when the maximum likelihood source strength is

where f (? ?|?0 ) is the probability of measuring ? ? given a true ?ux of ?0 . In the case of

be evaluated. In general, the con?dence of an interval extending to an upper limit of

negative. Barrett, et al. [5] suggest “lifting up” the measured ?ux to zero, evaluating the likelihood function and taking the upper limit. Instead, LIKE ?nds the width of the con?dence interval about the measured (negative) ?ux, then centers that con?dence interval about zero. It is very di?cult to estimate the probability that a con?dence interval obtained in this way would cover the true source strength; one must consider, for each possible value of the true source strength, all possible data sets. Any one of those data sets could fall into any of the categories we have outlined.

2.7.2

Calculating More Accurate Upper Limits

The upper limits calculated in this way are not only confusing and non-intuitive, they also do not (in general) ful?ll our requirement: a con?dence interval from zero to the upper limit should cover the true value 68.3% (or some other speci?ed fraction) of the time. A Bayesian approach is in this case more intuitive and ful?lls our requirements. Instead of considering the ensemble of data sets that might produce a given ?ux measurement, we start directly with the data, and consider the range of true ?ux values that might produce the data. The requirement that all ?uxes must be positive Our posterior probability distribution is identical to the likelihood, except that it is is easily ful?lled; we form a prior that is ?at for ? ≥ 0 and zero for negative ?.

CHAPTER 2. STATISTICAL METHODS IN γ -RAY ASTRONOMY

34

Figure 2.5: Upper limits calculated by Bayesian methods. This hypothetical likelihood function represents a measurement that gives a most likely value of some parameter to be less than zero. If the parameter is constrained by the physics to be positive, the Bayesian formalism suggests forming a prior which is zero for negative parameters values, and constant for positive parameter values. The posterior distribution would then look like the shaded portion of the curve. zero for all negative values of ?, and renormalized so that the integral over all values is unity. (Figure 2.5). The interpretation is then straightforward. The most likely value of the ?ux is zero; the 68.3% upper limit is found by integrating the posterior distribution from zero to some value σu.l. where the integral equals 68.3%.

2.7.3

Implications for Extant Conclusions

There are two sorts of conclusions typically drawn from EGRET upper limits. The ?rst is generally qualitative. A known X-ray or radio source is not detected in the EGRET data. A “low” upper limit is taken as evidence that γ -ray emission is small or non-existent. A “high” upper limit is taken as weak evidence that there is little γ -ray emission, but that there was not su?cient data to draw a conclusion. For these sorts of conclusions, the LIKE upper limits are adequate (e.g., [109, 77]). The second sort of conclusion is generally quantitative and statistical in nature. Often, possible source variability is examined in a number of observations. Some EGRET observations of a source (usually an AGN) yield signi?cant detections, and

CHAPTER 2. STATISTICAL METHODS IN γ -RAY ASTRONOMY

35

some yield only upper limits. Given ?ux measurements and con?dence regions, it is a simple matter to formulate a χ2 or likelihood test to determine if a constant ?ux model is compatible with the data. It is critical for such a test that upper limits have well-de?ned statistical properties. The upper limits generated by LIKE and quoted in the EGRET catalogs do not have these properties [195, 196, 65]. Unfortunately, catalogs of variability have been compiled based precisely on these upper limits [121]. All variability conclusions which involve upper limits are suspect, and should be treated as qualitative suggestions rather than quantitative results. Work by W. F. Tompkins is in progress to compile variability catalogs which have been calculated with statistically meaningful upper limits [199]; these should be used as soon as they are available.

Chapter 3 γ -Ray Bursts

In the midst of the Cold War, in the late 1960s, a number of satellites were launched carrying γ -ray detectors. Sensitive to γ -rays between ?200 keV and ?1.4 MeV, the Vela satellites were designed to detect the testing or use of nuclear weapons. Between July of 1969 and July of 1972, four Vela satellites, equally spaced in the same circular orbit, detected 16 bursts of γ -ray energy. Comparisons of the γ -ray arrival times in di?erent satellites determined that the origin of the bursts was more than 10 orbital diameters away. Thus, the ?rst theory of γ -ray bursts, Soviet nuclear testing, was ruled out. Although national security concerns delayed the publication of these intriguing results by Klebesadel, Strong, and Olson [97] for a number of years, they would mark the beginning of the longest-standing mystery in astrophysics since the Shapley-Curtis debates. A number of other bursts were observed over the next two decades [129, 208]. A consensus emerged fairly early that the source of the mysterious γ -ray bursts was Galactic in origin [53].? Preliminary detections of ?ux over 1 MeV strengthened this conclusion [149]; Schmidt “showed” that detection of emission over 1 MeV required a Galactic origin, since the source luminosity required for bursts at more than ?100 pc would imply an energy density that would result in γ –γ pair production [175]— the optical depth to Earth for 1 MeV γ -rays would be >1. While there was never

This work will deal only with the so-called “classical γ -ray bursts.” Another class of transient γ -ray sources, the soft gamma repeaters, are clearly a di?erent type of object.

?

36

CHAPTER 3. γ -RAY BURSTS

37

signi?cant evidence as expected that the bursts were preferentially located in the Galactic disk, it was presumed that the BATSE experiment on board CGRO would clear up any remaining ignorance about γ -ray bursts. Despite the general consensus on the location of γ -ray bursts, the ?eld attracted a wide variety of speci?c theories. Nemiro? identi?es 99 distinct γ -ray burst theories put forth between 1968 and the end of 1991 [148]. Early theories placed bursts both locally and cosmologically distant, with energy generation mechanisms ranging from relativistic dust grains in the solar neighborhood to cometary collisions with neutron stars to white hole emission to vibrations of cosmic strings at z = 1000 [4]. Nevertheless, by 1981, most models involved neutron stars in the Galactic plane. (Note, however, that Paczy? nski [157] put forth a theory of cosmological bursts with an optically thick e? e+ plasma out?ow in 1986. However, Paczy? nski also put forth a number of other unrelated theories of γ -ray bursts; thus it is unclear whether his apparent success is due to prescience or judicious covering of theory space.) In September 1991, Gerald Fishman, representing the BATSE team at the ?rst Compton Symposium, summarized the results of the γ -ray burst search since the launch of CGRO in April of that year. Everything that had been believed about the origin of bursts was shaken. The data, published ?rst by Charles Meegan, Fishman, and others in Nature [123], showed an isotropic distribution of γ -ray bursts across the sky. Furthermore, there was already evidence of the so-called “edge” of the γ -ray burst distribution. That is, there were fewer low-?ux bursts than would be expected if bursts were standard candles uniformly distributed in Euclidean space. The distribution was incompatible with a galactic disk population. Old theories die hard, and much e?ort went into searches for the expected spatial and temporal correlations between bursts in the BATSE catalogs. While it appeared in the ?rst BATSE catalog that some correlations might exist, additional bursts in the second catalog made isotropy much more likely [37]. Such statistical tests must be performed very carefully, due to the biases acquired from variable thresholds and exposures [160, 162, 107]. Similarly, only very marginal evidence could be found for repeating γ -ray bursts [161]. The distribution of all bursts observed to date (Figure 3.1) is compatible with a uniform distribution across the sky [125].

CHAPTER 3. γ -RAY BURSTS

38

Figure 3.1: Locations of γ -ray bursts observed by BATSE as of July 3, 1998, in Galactic coordinates [125]. The BATSE data quickly winnowed away many γ -ray burst theories. But like opportunistic species after an ecological disaster, new theories quickly sprung up to ?ll the void (e.g., [186]). The new theories were classi?ed primarily by the location of the bursts: locally (in the Oort cloud), in the hypothesized large Galactic halo, or at cosmological distances. The ?rst of these was essentially abandoned when no suitable energy generation mechanisms could be envisioned.

?

Halo models generally involved

interactions with old neutron stars that had been ejected from the plane into the halo. While the presence of a dark Galactic halo has been inferred from the rotation curves of other galaxies, its extent has not been conclusively measured [80]. The isotropy of the γ -ray burst distribution suggested that the radius of the halo would need to be many times the distance between the Earth and the Galactic center. Similarly, the halo of M31 (Andromeda) would be evident in the distribution once BATSE detected a large enough sample of bursts. The isotropy of the ?rst 1122 γ -ray bursts detected by BATSE [124] suggests if the bursts were in the Galactic halo, the halo would have a radius greater than 100 kpc. To be consistent with the lack of any excess toward Andromeda, the scale would have to be less than 400 kpc; a rather tight constraint. If the bursts were at cosmological distances, evolution e?ects conveniently explain the

C. D. Dermer proposed a mechanism at a conference in Santa Barbara, California in 1995 involving antimatter comets annihilating with normal comets in the Oort cloud. This theory was apparently never published.

?

CHAPTER 3. γ -RAY BURSTS

39

observed “edge” of the distribution. The energy required to produce the observed ?ux of the binding energy of a neutron star. Such an energy scale leads naturally to catastrophic theories involving neutron stars. It had been realized quite early [175] that the large energy release required for distant γ -ray bursts would produce an environment that was optically thick to γ –γ pair production. An assumption of cosmological origin of the γ -ray bursts led to a bifurcation of theories into those describing the energy generation mechanism and those describing the propagation of a large amount of compact energy, regardless of its source. The relevance of the latter theories for EGRET observation will be discussed in §3.2. from such distances (?1051 –1054 ergs, modulo any beaming factor) is a few percent

3.1

Recent Observations

In 1981, Fichtel & Trombka speculated that “identi?cation of the [γ -ray burst] objects with observations at other wavelengths will probably be required before signi?cant progress can be made in determining their origin” [41]. The day for which they were waiting arrived February 28, 1997, when the Italian satellite BeppoSAX discovered an X-ray afterglow to a γ -ray burst [30, 31]. The BeppoSAX satellite has a collection of instruments, including a γ -ray burst monitor, a wide-?eld imaging X-ray camera with a positional resolution of about 3 arcminutes, and a narrow-?eld X-ray camera with a position resolution of about 1 arcminute. Within 8 hours of the detection of the burst and its imaging with the wide-?eld camera, the narrow-?eld camera observed a source, consistent with the burst error circle, at α2000 = 05h 01m 57s , δ2000 = 11? 46′ 42′′ [50]. Optical and radio transients at the same position were quickly found [52, 159]. Some time later the optical transient was also seen by the Hubble Space Telescope [173]. The optical, radio, and X-ray sources were observed to decay as a power law consistent with the ?reball blast wave theories of M? esz? aros and Rees [127] and Wijers et al. [210]. In the following months several other bursts would be found to be consistent with transient sources in other wavelengths. Another revelation came from GRB 970508,

CHAPTER 3. γ -RAY BURSTS

40

observed May 8, 1997 by the wide-?eld camera on BeppoSAX and by its narrow-?eld camera 5.7 hours later. An optical counterpart was quickly found by Bond [17] and a radio counterpart found by Frail et al. [49] consistent with the narrow-?eld BeppoSAX position. Metzger and his colleagues at Caltech identi?ed a set of absorption features well, at the same redshift [130]. For the ?rst time since their discovery 30 years before, there was direct experimental evidence that γ -ray bursts were cosmological [131]. Another notable burst, GRB 971214 [67], was detected at a number of di?erent wavelengths: optical [59], infrared [55], and ultraviolet [16], as well as being detected by a number of γ -ray and X-ray instruments including BATSE , RXTE , and Ulysses [95]. Absorption and emission lines have been detected from this burst as well, yielding a redshift measurement of 3.43—?rmly establishing the cosmological nature of the γ -ray bursts [101]. Most recently, the redshift of the host galaxy of GRB 980703 has been measured at the Keck Observatory. Absorption and emission lines yield a redshift of z = 0.966 [35]. A good review of the current state of observational a?airs is given by Wijers [209]. Now that the observational mechanism for establishing γ -ray burst counterparts has been established, the ?eld is changing rapidly. As of this writing, 11 γ -ray bursts have been observed in radio and optical wavelengths [108]. These observations and others which will certainly be made in the near future have strongly a?ected, and will continue to a?ect, the leading theories of γ -ray bursts. implying a redshift of ≥ 0.835 [132]. By early June, they detected emission lines as

3.2

γ -Ray Burst Models

Until the launch of BATSE , γ -ray burst theories outnumbered the bursts themselves. In light of the data provided by BATSE , BeppoSAX , and other satellites, it is worthwhile to examine the leading theories with an eye to understanding the consequences observable by EGRET . We will ?rst look at the theories of the energy source powering γ -ray bursts, and then examine some of the theories of the γ -ray generating shock waves created by the bursting source.

CHAPTER 3. γ -RAY BURSTS

41

3.2.1

Energy production mechanisms

Gigantic explosions hold fascination for the theoretical physicist as much as for the schoolchild. Cosmological γ -ray bursts require the largest explosions known, and thus attract theorists in droves (e.g., [60]). As long ago as 1975, Malvin Ruderman noted [172], “For theorists who may wish to enter this broad and growing ?eld, I should point out that there are a considerable number of combinations, for example, comets of antimatter falling onto white holes, not yet claimed.” Mercifully, the considerable experimental data amassed since 1975 has largely narrowed the ?eld to two general mechanisms: neutron star–neutron star mergers and massive stellar collapses known as “collapsars” or “hypernov?.” Nevertheless, at this stage in our understanding of γ -ray bursts it is unreasonable to think that all theories could be placed in one of only two categories; these represent only the most prominent of the theories. Neutron star mergers. The most popular theory at the time of this writing is the neutron star merger model [145]. The basic premise is very simple. Neutron star binary systems slowly lose energy as a result of gravitational radiation. Eventually, the neutron stars will spiral into each other, presumably resulting in a large release of energy. A variant of the model replaces one of the neutron stars with a black hole. Narayan [145] cites several advantages of such a model, as well as a number of issues to be resolved. First, the source population is known to exist. Neutron star binaries have been observed, and their energy loss corresponds with that expected from gravitational radiation to better than 1% [188]. Second, the energetics are of the right order. The energy release from such a merger would exceed 1053 ergs in ?1 ms within 100 km. Finally, the frequency of such mergers may be approximately right. Estimates of the merger rate include the range from 10?6 to 10?4 yr?1 per standard galaxy, which matches the observed burst rate under a cosmological scenario and isotropic emission from the energy source. An early objection to the neutron star merger model was the so-called “no host galaxy” problem. Optical observations of γ -ray burst error circles before the launch of BeppoSAX failed to ?nd the galaxies expected as the hosts of the colliding neutron stars. However, it was realized that binaries can often acquire substantial recoil

CHAPTER 3. γ -RAY BURSTS

42

velocities as a result of their two supernova explosions. For reasonable parameters, binaries may travel 1–1000 kpc before merging [201]. Therefore, by the time a binary neutron star system becomes a γ -ray burst, it may have left its progenitor galaxy. The diverse character of observed γ -ray bursts has been pointed out as a problem for the merger model. Neutron stars are expected to have a very narrow mass range. Merging binaries should almost always consist of about 3 M⊙ collapsing in a very clean gravitational system. Such a homogeneous energy generation mechanism should result in a fairly homogenous population of γ -ray bursts, though beaming e?ects, magnetic ?elds, and an inhomogenous environment may explain the observed di?erences. Hypernov?. A much newer idea is the hypernova model due to Paczy? nski [158]. The idea is similar to Woosley’s “failed supernova” model [214]. A very massive star then the angular momentum of the star requires the formation of a rotating dense undergoes core collapse to form a ?10 M⊙ black hole. If the star is rapidly rotating,

torus of material around the Kerr black hole. Any previously existing magnetic ?eld will be signi?cantly a?ected by the collapse; most possibilities strengthen the local magnetic ?eld. Paczy? nski estimates that the rotational energy of the black hole should be of order 1054 ergs. He also ?nds that the maximum rate of energy extraction is Lmax ≈ 1051 ergs s?1 B 1015 G

2

MBH 10 M⊙

2

(3.1)

Achieving the required energy release requires ?elds on the order of 1015 G. No mechanism for generating such ?elds is o?ered, although many theorists are currently working on the details of the e?ects of the core collapse on the ambient magnetic ?elds. The observed spectrum from either model of γ -ray bursts depends more on the details of the ?reball model discussed in the next section than on the energy generation mechanism. The primary di?erence between these mechanisms is their location. While neutron star mergers are often expected far from galaxies, hypernov? are found in star-forming regions, since their progenitors are rapidly evolving massive stars. The lack of optical observation of GRB 970828, despite fairly deep searches, combined with

CHAPTER 3. γ -RAY BURSTS

43

the reported large hydrogen column density [142] has led Paczy? nski to conclude that the optical emission was not observed due to extinction by dust. A substantial number of detections of γ -ray bursts at multiple wavelengths should illuminate this question in the very near future.

3.2.2

Blast wave theories

While very early experiments seemed to detect thermal radiation from γ -ray bursts, better instruments soon observed a characteristic power law emission in the X-ray and γ -ray regime. Theoretical explanations have centered on the blast-wave model, where the injection of a large amount of energy into a very small volume leads to an expanding ?reball, optically thick to pair production, expanding into some surrounding medium. The simplest model—also due to Paczy? nski [157]—assumes the spherical expansion of an optically thick e? e+ plasma with no surrounding medium. This model results in a blackbody spectrum. A series of more and more complicated models were based on this simplistic one, ?nally culminating in the model of M? esz? aros, Rees, and Papathanassiou [128], which calculated spectra expected for a variety of magnetic ?eld con?gurations and particle acceleration e?ciencies, including the shock fronts and reverse shocks which arise from the expansion of the burst ejecta into the surrounding medium. Shock fronts in blast waves have become the baseline from which various energy generation mechanisms diverge. Their popularity stems from their success: they can naturally explain how energy can be reconverted into γ -rays from the particles that must emerge from the ?reball, they can naturally explain the observed burst time scales, and, assuming the medium into which the ?reball expands varies slightly from burst to burst, they can explain the wide variety of burst pro?les and time scales observed [165, 204]. The observed afterglows may be explained as emission from the slowing, cooling shock [205, 203]. However, several theoretical issues remain. The mechanism of particle acceleration in relativistic shocks is not well understood, though it is widely assumed that such acceleration results in a power law energy

CHAPTER 3. γ -RAY BURSTS

44

spectrum [6]. The magnitude and nature of the magnetic ?eld is unknown. Preexisting magnetic ?elds, if su?ciently strong, will provide signi?cant synchrotron radiation in the particle blast wave. In addition, turbulent mixing in the shock can generate signi?cant ?elds due to charge separation. Finally, the coupling between electrons, baryons, and the magnetic ?eld is not well understood. Energy radiation from elections is very e?cient, but most of the energy is in the baryons, or perhaps the magnetic ?eld, which do not radiate their energy so readily. The nature of their coupling a?ects the total radiated power as well as the time scale of the radiation [28, 168]. While a detailed investigation of the various blast-wave models is beyond the scope of this work, a general overview is given, based on that found in Chiang & Dermer [28]. The burst begins with the deposition of 1051 –1055 ergs of energy in a radius of about 100 km. The nature of this energy is not known. It is assumed that most of this energy will be transformed into kinetic energy of baryons, which expand adiabatically. The baryons soon become cold; that is, the baryon speeds in the comoving frame of the bulk ?ow become sub-relativistic. The bulk Lorentz factor is then given by Γ0 ? E0 /M0 c2 where M0 is the rest mass of baryons. As this sphere freely expands into the surrounding medium, it sweeps up material.

A shock front begins to form. At some “deceleration radius” rd the shell can no longer be approximated as freely expanding, and the bulk kinetic energy of expansion begins to be reconverted into internal energy of the baryons. This radius is the point at which the integrated momentum impulse of the swept-up matter equals the original baryonic

3 rest mass: rd ≈ (3/4πρΓ0 )M0 where ρ is the density of the surrounding material.

Predictions of what happens at this deceleration radius form the core of most

blast-wave models. A function Γ(r ) is derived which expresses essentially the rate at which baryonic kinetic energy is converted to radiation. Since cosmological scenarios require high initial Γ0 of order 102 –103, much of the energy release in all scenarios is compressed into the ?rst few tens of seconds in the observer’s frame. Nevertheless, models of the blast wave at the deceleration radius and beyond make predictions about the observed spectrum at various times as well as the observed burst time scales.

CHAPTER 3. γ -RAY BURSTS

45

Given the recent evidence that γ -ray bursts really are cosmological in origin, the energy required for the ?uences observed on Earth make some sort of ?reball almost inevitable. Blast-wave models have qualitatively reproduced some aspects of γ -ray burst observations, and appear to be a promising theoretical road to pursue. However, it should be noted that blast waves as they have been described cannot be the whole story. In either the neutron star merger model or the hypernova model, there is a natural symmetry axis associated with the source. A symmetric blast wave then seems rather unlikely. Asymmetry in the blast wave will lead to beaming on some scale. Beaming will of course reduce the energy required per burst, but increase the rate at which the bursts must occur. It is not clear how beaming would a?ect the radiation produced by shock waves. The radiation reaching the Earth has already been beamed into a cone of order 1/Γ2 0 . This tight beaming means that all of the γ -rays observed at the Earth are emitted from a very small piece of the blast wave; therefore, the amount of large-scale symmetry in the wave may be irrelevant to the spectrum observed.

3.2.3

Observable consequences

Without the compass of observations, we are doomed to drown in an ever-deepening sea of theories. Fortunately, the tools of the experimental astrophysicist are becoming more and more powerful. A number of X-ray satellites are currently operating, with more planned, which produce error boxes small enough to allow radio astronomers to bring their substantial instruments to bear. GLAST will revolutionize high-energy γ -ray astronomy in the next decade just as EGRET did in this decade. And new types of astronomy, previously relegated to science ?ction or crackpot dreams, are ˇ beginning to provide useful data. Air Cerenkov detectors can localize high-energy γ rays from a few hundred GeV up to many TeV to much less than a degree. B¨ ottcher & Dermer calculate the high-energy emission from proton synchrotron radiation and ˇ photopion-induced pair cascades, and ?nd that future high-sensitivity Cerenkov telescopes with low energy cuto?s (or, it should be noted, good very-high energy response

CHAPTER 3. γ -RAY BURSTS

46

from GLAST ; see Chapter 6) could measure the level of the infrared background radiation [18], since high-energy γ -rays will interact with the infrared background to produce electron-positron pairs. While no current neutrino detectors have su?cient sensitivity, a measurement of the neutrino ?ux from a γ -ray burst could shed light on the energy generation mechanism [206]. Additionally, an intense gravitational collapse will produce gravitational waves, which could be detectable with the detectors currently under construction [99].

3.3

EGRET observations

While BATSE has enjoyed the spotlight for most of the contributions of CGRO to the γ -ray burst problem, EGRET has made some important observations as well. At least 16 bursts occurred outside the ?eld of view of the spark chamber, but were nevertheless detected in the calorimeter [22, 103, 177]. Five γ -ray bursts have been detected in the EGRET spark chamber [34]. Each has characteristics worthy of some examination. The ?rst burst detected in the EGRET spark chamber occurred on May 3, 1991. Six photons were detected in the spark chamber in two seconds, coincident with the signal received from BATSE [176]. Measurements of the background before the burst suggested that 0.18 photons per two seconds were expected in the spark chamber. Additionally, the burst was evident in the TASC calorimeter data, as well as the anticoincidence trigger rate. The anticoincidence dome is sensitive to photons above about 20 keV via Compton scattering, while the TASC has four triggering levels, corresponding approximately to energies above 1.0 MeV, 2.5 MeV, 7.0 MeV, and 20.0 MeV. Firm detection of this event with all three systems, coincident with the BATSE trigger not only demonstrated γ -ray burst detection with EGRET but also measured the ?rst γ -ray burst emission above 100 MeV: one of the spark chamber photons had an energy of approximately 230 MeV. The next burst was detected a month later, on June 1, 1991 [103]. This burst was also seen in the TASC and anticoincidence dome, and four photons were detected in the spark chamber when 1.5 background photons were expected. The detection of

CHAPTER 3. γ -RAY BURSTS

47

high-energy photons from these bursts began to shift the evidence toward power-law emission instead of thermal emission. Furthermore, the last photon to arrive in the spark chamber was nearly 70 seconds after the initial BATSE trigger. Bu?alo Bills quarterback Jim Kelly remembers January 31, 1993 as the day the Bills lost their third consecutive Super Bowl, this time to the Dallas Cowboys by a score of 52–17. High-energy astrophysicists worldwide remember that date for the detection of what still stands as the most intense BATSE burst ever seen [100]. BATSE count rates exceeded 2×106 counts s?1 . Signi?cant structure to the burst was observed at scales below 10 ms, and the ?rst intense spike lasted 200 ms. Fortunately, the burst was in the ?eld of view of the EGRET spark chamber, which detected 16 γ -rays (compared to 0.04 expected by chance) in 25 seconds [181]. Two of the γ -rays had energies of nearly 1 GeV. EGRET measured a power-law photon spectrum with index ??2. The hard photon spectrum lent further credence to shock acceleration closer than ?50 pc. A total ?uence measurement could not be made by EGRET , since the dead time per spark chamber event is approximately 100 ms. It is likely that a number of high-energy γ -rays passed though the spark chamber in the ?rst 100 ms. It would be more than a year before EGRET would detect another burst, but that burst would provide plenty of fuel for the raging γ -ray burst debate. The burst, which arrived on February 17, 1994 and whose initial pulse lasted for 180 s, was detected by COMPTEL [96], Ulysses , BATSE , and EGRET . Of primary theoretical importance was a single EGRET photon, detected some 4500 s after the initial burst, with an energy of 18 GeV [76]. The probability of this photon originating in the 940217. Ten photons were detected during the ?rst 180 s, concurrent with the BATSE detection, including a 4 GeV photon and a 3.5 GeV photon. The remaining 18 photons were detected over the next 1.5 hours. Hurley et al. [76] point out a number of conclusions which can be immediately drawn. While the universe is optically thin at γ -ray energies, the cosmic microwave background becomes an e?cient medium for γ –γ pair production for su?ciently high γ -ray energies. For a 25 GeV photon, this background was 5×10?6 . In fact, EGRET detected a total of 28 photons from GRB models, although the high-energy photons were taken as evidence that the burst was

CHAPTER 3. γ -RAY BURSTS

48

consideration limits the source to z<55. To avoid attenuation from the intergalactic infrared background, the source distance is further limited to z<2.5, ruling out early universe theories [4]. Finally, they note that if the spectrum of the ?rst 180 s of GRB 940217 is extrapolated to 1 TeV, the predicted ?uxes would be detected by air ˇ Cerenkov telescopes and the Milagro air-shower array. Unfortunately, bursts with such high-energy emission are evidently somewhat rare, and the ?elds of view of air ˇ Cerenkov detectors are quite small. Nevertheless, the current network of fast burst position detections may allow observation of delayed emission like that of the February 17 burst. The very long delay of emission has also been the source of much theoretical speculation. The existence of such delayed emission is a natural consequence of blastwave models. While the initial γ -ray burst is caused by the formation of a shock with the surrounding medium, the blast front will continue to radiate as it decelerates. The most recent burst detected in the EGRET spark chamber arrived only a few weeks later, on March 1, 1994. This burst was similar to the ?rst two which had been detected, with 7 photons of maximum energy 160 MeV, arriving within 20 s [177]. The spectrum of this burst was softer (??2.5) than the February 17, 1994 burst. It is unlikely that any more bursts will be observed in the EGRET spark chamber,

as instrumental lifetime concerns have reduced EGRET to Target of Opportunity observation only. Nevertheless, the observations made by EGRET create signi?cant constraints on the models used to explain γ -ray bursts. The GLAST instrument, with its very wide ?eld of view (§6.2) and large sensitive area up to very high energies (300 GeV), should further constrain the high energy behavior of γ -ray bursts.

3.4

Possible EGRET -only Bursts

Observations of γ -ray bursts in the last decade have shown conclusively that the high-energy emission has a power-law spectrum. If, for some bursts, the burst has a harder spectrum than the background, the signal-to-noise ratio should increase at higher energies. Some blast wave theories predict that a spectral break may occur at or just below EGRET energies [33]. In such models, the total peak spectral power early in the burst is emitted above 100 MeV.

CHAPTER 3. γ -RAY BURSTS

49

It is thus possible that some γ -ray bursts will be detectable by EGRET , but not by lower-energy instruments such as BATSE . The evidence of such bursts would therefore be lurking in the EGRET photon database. As was discussed in §2.3, the best way to deal with the photon-by-photon nature of EGRET data is through the use data for γ -ray bursts. While they did ?nd the known bursts in the data they searched, they apparently did not perform a comprehensive search of the EGRET database. The di?culty with such a search is that the exposure for each photon has to be calculated individually. This process is very time-consuming, due to the the extended nature of the EGRET calibration ?les. Therefore, a comprehensive search required some adaptation of statistical methods [87]. Those methods, and their results, will be described below. of likelihoods. Buccheri et al. [20] developed a likelihood method to search EGRET

3.4.1

Statistical Methods

The γ -ray burst search algorithm was designed to be fast and e?cient at ?nding bursts. From a statistical standpoint, a burst may be de?ned as a time interval with a measured rate which is, to some speci?ed con?dence, incompatible with Poisson ?uctuations. The method, then, is fairly straightforward: ?rst, ?nd the background rate; second, ?nd any time intervals with su?ciently high rates to rule out Poisson ?uctuations at the given con?dence level; and ?nally, verify that the spatial distribution of the photons is consistent with a point source. The sequential nature of this method makes it fast compared with a full likelihood analysis; however, it also involves some binning, and as we will see, complicates our estimates of the signi?cance of detections. To further speed the algorithm, we search for signi?cant intervals in two steps. The ?rst compares photon arrival rates only. The rate is de?ned as the number of photons from some area on the sky per second. It requires only the amount of instrument live time and a photon count. The second step will measure the photon ?uxes (photons cm?2 s?1 ). This is certainly the more physically relevant quantity, but it is also computationally more expensive to determine the state of the EGRET

CHAPTER 3. γ -RAY BURSTS

50

instrument for the entire interval. Therefore, promising candidate intervals will be found by comparing rates, and then checked for signi?cance by comparing ?uxes. All photons in the viewing period with energies over 100 MeV and zenith angles less than 100? . 5 are sorted into these bins. The total instrument live time to the center of each square is calculated. This yields an array of approximate background count rates. The critical number of photons required for a signi?cant interval, Ncrit , is given implicitly by the Poisson formula α=

Ncrit ?1 n=0

To establish a background rate, the ?eld of view is binned into 5? × 5? squares.

e??t (?t)n n!

(3.2)

where ? is the average rate, t is the time interval being searched, and α is the con?dence level. The probability of observing Nobs >Ncrit is 1 ? α. To acquire a set of candidate events while avoiding speci?c time binning, each

photon is considered to be the start of an interval. All photons arriving within the standard interval length (either 1 hour, 30 minutes, 10 minutes, or 3 minutes) and within 5? of the initial photon are considered part of the same event and counted. If the number of observed photons exceeds Ncrit for that area of the sky, the candidate event is accepted for further evaluation. Because the number of candidate events which pass the ?rst cut is relatively small, it is practical to calculate instrument exposure and compare the candidate ?ux to the expected background ?ux. The single trial probability of such an event, PS , is rate ? replaced with the ?ux. calculated as in equation (3.2), with the time t replaced by the exposure E and the Unfortunately, PS cannot be interpreted as the probability that the ?ux in a given interval is not a Poisson ?uctuation. In fact, we have searched many di?erent intervals, and the appropriate con?dence level must take into account the number of trials. If the intervals had been ?xed and preselected so as not to overlap, then each interval would be statistically independent. The total probability PN that a given ?ux could not be attributed to Poisson ?uctuations would be given by the complement of

CHAPTER 3. γ -RAY BURSTS

51

the product of the chances that it would not be seen in each of any one interval:

N PN = 1 ? PS

(3.3)

where N is the number of independent trials. However, preselecting non-overlapping time bins would not be a good way to search for γ -ray bursts. If a burst did not happen to fall entirely within one bin, but instead split its ?ux between two bins, it would probably not be signi?cant in either bin. The sliding interval described above was designed to avoid this problem. Nevertheless, the sliding interval introduces its own problems: namely, that the intervals are no longer independent. In order to evaluate the actual signi?cance of an interval PT as a function of PN , a Monte Carlo simulation was performed. It was expected (and veri?ed) that the signi?cance of an interval should be monotonic in PN (Figure 3.2). Approximately data set searched contained 1.2×106 photons. Each Monte Carlo photon was assigned 5.7×108 simulated photons were generated for the Monte Carlo data set; the actual

an arrival time, drawn from a Poisson distribution with a given background rate. Each photon was also assigned a uniformly distributed x and y coordinate between 0? and 40? , simulating EGRET ?eld of view. The data were simulated using actual EGRET exposure information calculated for a typical point on the sky over many viewing periods, yielding an exposure set representative of the windowing and noncontinuous exposure actually obtained with EGRET . This exposure set was used as many times as necessary to generate the entire Monte Carlo data set. To determine appropriate background rates, 104 actual ?uxes were measured from random points in di?ering viewing periods. These ?uxes were sorted, and the highest and lowest 20% were discounted in order to exclude nonrepresentative outliers. The range of this tightened distribution of ?uxes was used as the range for the uniformly distributed ?uxes in the Monte Carlo simulation. This was done to ensure the use of typical, though not rigorously representative, background ?uxes. Since all probabilities are found from the di?erence between measured and expected ?uxes, they depend only

CHAPTER 3. γ -RAY BURSTS

52

Figure 3.2: Overall signi?cance PT versus raw probability PN from the Monte Carlo simulations. PN is calculated for each Monte Carlo event from equation (3.2) and equation (3.3). A cumulative count of the Monte Carlo events yielding such a PN yields the cumulative probability of that events, PT , given by the vertical axis. weakly on the precise values of the background ?uxes. The range of backgrounds found was 7.0×10?7 –4.8×10?6 photons cm?2 s?1 for a 5? circle on the sky. In some cases, the search algorithm may have detected the same event more than

once; that is, the interval following the ?rst photon in the event yielded a signi?cant rate increase, and the interval following the second photon also yielded a signi?cant rate increase. In these cases, the two probabilities must be independently combined, because the Monte Carlo distribution was found by counting all signi?cant events, regardless of whether or not they were part of the “same” event. In addition to exhibiting a ?ux increase, γ -ray burst photons should be statistically consistent with a point source origin. Candidate bursts found in the time series analysis were examined spatially using a likelihood analysis technique [36]. Likelihoods were calculated for all photons within approximately 20? of the ?rst photon of the interval. A null model of smooth background was compared to a model with a variable position point source plus a background rate taken to be the average background found above. To simplify computation, errors in photon position were taken from the width of the energy dependent best ?t Gaussian containing 68% of the EGRET

CHAPTER 3. γ -RAY BURSTS

53

point-spread function given by equation (1.1). The usual EGRET test statistic is likelihoods of the source and null models, respectively. The spatial analysis required then de?ned [117] as in §2.5: TS ≡ ?2(ln Ls ? ln Ln ), where the Ls and Ln are the

likelihood calculation for each individual photon. The standard EGRET likelihood software, LIKE [117], is designed to evaluate likelihoods based on maps of photon counts. Separate likelihood software was thus developed for this study. Simpli?cations in this implementation make direct comparison of these values of TS with those from LIKE inexact. For the purposes of source detection in EGRET analyses, a TS of 16 is considered to be a signi?cant (4σ ) detection. This is based on the measured distribution of TS. However, the likelihood statistic calculated here su?ers from a very low count rate; often the likelihoods are computed from less than 10 photons. The statistics of TS with very few counts are not well characterized. Furthermore, the background rate across the whole 20? is taken to be the same as that found for the central 5? circle. In regions where the background is spatially varying, this may distort the measured TS. Nevertheless, the TS measurement adds valuable information to the evaluation of burst candidates. Candidates with a very low TS, corresponding to little or no point-like structure, should be discarded. A sharp variation in photon arrival rate across a large area of sky is probably not due to an astrophysical process. It should be remembered that there is already a selection e?ect of candidate events due to the fact that the time series analysis considers photons con?ned to a ?ve degree circle. This selection e?ect as well as the low count statistics make translation of TS into a con?dence level problematic. Note also that the spatial and time series probabilities are not completely independent.

3.4.2

Results

We analyzed 182 viewing periods using these methods, corresponding to observations between 1991 April 22 and 1996 March 21. Viewing periods after this time used new instrument modes to compensate for various instrument malfunctions and narrow?eld viewing. This data was not searched, since the ?eld of view was signi?cantly

CHAPTER 3. γ -RAY BURSTS

54

Date 1993 Jan 31 1994 Feb 17 1994 Apr 27 1993 Mar 17 1992 Oct 8 1991 May 4

Time 18:57:12 23:03:05 01:31:01 06:40:42 04:35:04 05:16:33

? 287 152 121 65 201 204

b 51 -55 -0.7 21 31 13

TS 33.1 18.8 26.7 17.1 15.7 18.7

Number Expected 0.048 0.074 0.120 0.009 0.047 0.011

Number Observed 6 5 6 3 4 3

Max Energy (MeV) 1240 3382 680 361 508 960

PT >99.994% 98.8% 95.4% 52.2% 32.9% 29.0%

Table 3.1: Three minute time scale results. Burst candidates are listed with their times (UT), galactic coordinates ? and b of the maximum likelihood position, spatial TS, expected and observed numbers of photons, highest photon energy, and signi?cance, as derived from the Monte Carlo simulations.

Date 1993 Jan 31 1993 Jul 20 1993 Feb 22 1992 Mar 22 Time 18:57:12 07:07:55 06:46:44 08:25:50 ? 287 121 20 333 b 51 42 -52 0.4 TS 33.1 23.2 21.8 29.6 Number Expected 0.160 0.097 0.099 0.623 Number Observed 6 5 5 8 Max Energy (MeV) 1240 878 615 11100 PT 99.85% 74.0% 72.2% 28.1%

Table 3.2: Ten minute time scale results, as in Table 3.1. smaller, and the systematics of the new modes are not as well understood. Each viewing period was searched on four time scales: one hour, 30 minutes, 10 minutes, and 3 minutes. Results are presented in Tables 3.1—3.4. In each table, PT is the probability that no such events would be detected in the entire EGRET data set for that time scale in the null hypothesis. All ?uxes quoted below are found by dividing the number of photons by the total exposure in the scale time, even if the photons evidently arrived in a much shorter time. For events with very few photons, the ?ux is not very well de?ned, since an arbitrary change in cuto? time may reduce the exposure without changing the number of photons. Thus, the ?ux for the purposes below is always found by considering all of the exposure in the period. Two of the bursts triggered by BATSE, the Super Bowl Burst of 1993 January 31 and the 1994 February 17 burst, were also independently detected with this algorithm. The other three EGRET -detected BATSE-triggered bursts were not strong enough to be independently detected. For comparison, these bursts and their characteristics are listed in Table 3.5. Several independent detections were made, the most signi?cant occurring on 1994 April 27.

CHAPTER 3. γ -RAY BURSTS

55

Date 1993 Jan 31 1993 Jul 20

Time 18:57:12 06:50:45

? 287 120

b 51 42

TS 33.1 28.8

Number Expected 0.387 0.205

Number Observed 8 6

Max Energy (MeV) 1240 880

PT 97% 77%

Table 3.3: Thirty minute time scale results, as in Table 3.1.

Date 1993 Jan 31 1994 Apr 27 Time 18:57:12 21:34:18 ? 287 169 b 51 3 TS 33.1 26.0 Number Expected 0.387 1.146 Number Observed 8 10 Max Energy (MeV) 1240 522 PT 97% 70%

Table 3.4: One hour time scale results, as in Table 3.1.

3.4.3

Discussion

Each PT is the chance that measured photon arrival times are not due to random Poisson noise in the entire EGRET data set for that particular time scale. Thus, an event PT = 90% in the one hour time scale would be a spurious detection in 10% of the ensemble of possible EGRET data sets searched on a one hour time scale. However, four time scales have been searched. If each time scale were independent, this would be counted as four trials. But the same data were searched in each time scale, so the trials are not independent. Rigorous evaluation of the overall number of independent trials is problematic, and is not attempted for fear of producing meaningless results. The calculated probabilities also do not take into account the photon energies. However, high energy photons are rarer than lower energy photons, so the arrival of a similar number of high energy photons would constitute a more signi?cant event. If the detected bursts exhibited a plethora or paucity of high energy photons, the event energy spectrum would be expected to di?er appreciably from that of the background. Because each event has too few photons to meaningfully determine a spectrum, a hardness ratio was calculated. All photons from all bursts in each time scale were collected, and a hardness ratio for the set was calculated. Even such a collection su?ers from poor counting statistics, as evidenced by the large errors (Table 3.6). The hardness ratio is de?ned as the number of photons with energies greater than 300 MeV divided by the number with energies between 100 and 300 MeV. A typical background spectral index of 2.0 corresponds to a hardness ratio of 0.5. The measured

CHAPTER 3. γ -RAY BURSTS

56

Scale 3 min 10 min 30 min 1 hour 1991 June 1 3 min 10 min 30 min 1 hour 1994 March 1 3 min 10 min 30 min 60 min

Date 1991 May 3

Number Number Expected Observed PS PT 0.202 3 99.88% < 10?300 0.672 5 99.93% < 10?300 1.288 5 98.97% < 10?300 1.288 5 98.97% < 10?300 0.396 4 99.93% < 10?300 0.633 5 99.95% < 10?300 0.633 5 99.95% < 10?300 2.064 9 99.97% < 10?300 0.102 3 99.98% < 10?300 0.333 4 99.96% < 10?300 0.569 4 99.72% < 10?300 0.569 4 99.72% < 10?300

Table 3.5: Known BATSE -triggered γ -ray bursts not detected independently by EGRET in each time scale. The expected and observed numbers of photons are given, along with the single trial probability PS and the signi?cance PT as determined by the Monte Carlo. hardness ratios are consistent with the background hardness ratios, implying that the PT would not be signi?cantly modi?ed if energy information were taken into account. All probabilities are dependent on the statistical distribution found from the Monte Carlo simulation. Thus, any possible factors which might make the Monte Carlo an imperfect simulation of the real data must be examined. Edge e?ects may slightly a?ect the ?nal distribution, since the 5? circles drawn around photons near the edge of the region will reach past the edge. However, a short simulation was done in a 20? by 20? square, and the resultant distribution did not di?er signi?cantly from that calculated with a larger area, indicating that edge e?ects play a small part. The Monte Carlo simulation also assumed spatially uniform exposure and background rate in order to make computation time reasonable. These assumptions will be poor only when the exposure or rate varies signi?cantly over the 5? circle. Exposure generally varies smoothly, so this is probably a good approximation. The background rate can vary extensively if a strong point source is nearby, but the di?use background

CHAPTER 3. γ -RAY BURSTS

57

Time Scale Hardness Ratio 3 min 0.723 ± 0.338 10 min 0.385 ± 0.202 30 min 0.500 ± 0.433 1 hour 0.500 ± 0.306 Table 3.6: Hardness ratios for each time scale. The hardness ratio is de?ned as the number of photons above 300 MeV divided by the number of photons between 100–300 MeV. does not vary rapidly on this scale. Thus, we may need to worry about the distribution if events in the real data are detected near steady point sources; otherwise, the constant background approximation should be relatively good. Spatial correlations must also be considered in the evaluation of burst candidates. Unfortunately, rigorous treatment of spatial correlations is fraught with di?culty. The maximum likelihood distribution is not well characterized for low counts. Furthermore, the 5? search radius used will select for spatially correlated events. It was hoped that spatial analysis would add signi?cantly to our understanding of these events; unfortunately, because of these biases and correlations the spatial analysis has yielded little additional insight. The most interesting candidate event is clearly the 1994 April 27 01:31 event in the 3 minute time scale. The signi?cance of this event, while much lower than the BATSE -independent detection of the Super Bowl Burst, is almost as high as the BATSE -independent detection of the 1994 February 17 burst. This detection occurred while EGRET was in its most common, largest e?ective area mode. The observed 27? o? the instrument axis. It thus seems unlikely that it is a spurious to the six photons within three minutes of elapsed time classi?ed as good events, there were seven additional events which were rejected by the standard EGRET data analysis for a variety of reasons. One of these generated too few sparks in the spark chamber, one had its vertex in the wall of the chamber, and the rest failed to produce an acceptable time of ?ight measurement. Estimated trajectories of these events Earth zenith angle to the event was ? 40? , well away from the horizon. It was

detection caused by pushing the operating envelope of the instrument. In addition

CHAPTER 3. γ -RAY BURSTS

58

suggest that they originated from the same area on the sky as the good events. Although a precise calculation of the signi?cance of these events is impossible, there appear to be more such events from that direction on the sky than would otherwise be expected. The time interval corresponding to the 1994 April 27 candidate was examined in the 30 - 100 MeV range as well. No photons were detected during that period. However, EGRET sensitive area to this energy range is comparatively small, so that only 0.046 background photons were expected (compared to 0.120 in the 100 MeV and above range.) The TASC and anti-coincidence dome rates did not show a signi?cant change in rate. With such small numbers, it is di?cult to assess the signi?cance of the lack of low energy photons. Furthermore, a check was done of BATSE detections nearby. No trigger occurred at that time, although BATSE triggering was enabled. The last previous trigger occurred some 20 hours earlier. No other bursts were detected so signi?cantly. However, there were several less signi?cant detections, making it worthwhile to consider the probability that at least one γ -ray burst has been detected. In the three minute time scale, the April 27, 1994 burst is fairly signi?cant by itself. However, it is interesting to calculate the joint probability that all the events in each time scale are due to random Poisson noise. We will exclude the Super Bowl Burst and the February 17 Burst, since they are already known to be a γ -ray bursts. Table 3.7 shows the joint probabilities in each time scale from the arrival time data. There is apparently some evidence for burst occurrence in the short time scales, with no signi?cant evidence for occurrence in the half-hour and hour long time scales. The detected candidate events were compared to the Second BATSE Burst Catalog to see if any coincident detections were made. With the exception of the bursts actually triggered by BATSE, that is, the Super Bowl and February 17 bursts, no coincidences were found. While the Poisson error analysis may be encouraging, it is important to consider the contribution of systematic errors before coming to any conclusions. The search for γ -ray bursts by the method used is most sensitive to systematic errors in the

CHAPTER 3. γ -RAY BURSTS

59

Time Scale Probability 3 min 99.0% 10 min 94.8% 30 min 77% 1 hour 70% Table 3.7: Joint probabilities of detecting all the marginal candidates in each time scale. While systematic errors could be responsible, there is apparently evidence for burst activity at or below the three minute time scale. determination of background rates. If the background is consistently underestimated relative to the burst photons, many apparently signi?cant events are not actually signi?cant. The most likely source of systematic errors in the background rate is the error in the estimate of the instrument exposure. A subtle point should be explored here. EGRET switches observing modes every few minutes. Each mode has a di?erent amount of sensitive area to a given position on the sky, as well as a di?erent ?eld of view. Approximately 80% of the time EGRET is in the single mode with the largest ?eld of view and sensitive area. In the limit that the instrument is always in the same observing mode, any exposure errors will cancel each other, since the error induced in the background ?ux will be exactly the same as the error induced in the burst ?ux. However, if the exposure error is di?erent for each instrument mode, then a burst candidate detected in a rare mode might exhibit considerable systematic error. The sensitivity of EGRET diminishes as the gas in the spark chamber ages. This e?ect has been partially characterized and compensated for in the EGRET exposure data. The time scale for this drift is typically several months, much longer than the burst search scales and also at least twice as long as the longest background averaging time. The nondetection of the 1994 April 27 candidate burst by BATSE must be considered in light of current models for γ -ray burst production. The relativistic ?reball model of M? esz? aros, Rees, and Papathanassiou [128] suggests the possibility of a high energy ?at spectrum tail. If the spectrum is ?at down to the BATSE energy range, the signal would be lost in the soft-gamma background. However, that model does

CHAPTER 3. γ -RAY BURSTS

60

not predict a ?at spectrum at BATSE wavelengths. Signi?cant redshift e?ects might move the break of the ?reball model spectrum to below BATSE energies, but then x-ray detection would be expected. Nevertheless, other models [33] suggest that at the earliest times after the burst, the spectral break could be at EGRET energies.

3.5

Conclusions

The mystery of the γ -ray bursts has recently experienced its third major revelation. The ?rst occurred in the late 1960s and early 1970s, when the energetic bursts were discovered. The next occurred in 1992, when BATSE made clear the isotropic, nonhomogeneous distribution of the bursts—pointing to the almost inconceivable conclusion that the bursts were of cosmological origin. The third was the association, ?rst by BeppoSAX , and then by many others of a few bursts with counterpart afterglows in many wavelengths, resulting in the measurement of cosmological redshifts. One of the outstanding questions about γ -ray bursts has thus been answered. Nevertheless, we are far from a complete understanding of the burst mechanism. The central engine is probably either a merger of neutron stars or a hypernova. Energetic concerns suggest fairly strongly that the initial energy release goes into baryonic kinetic energy, and that the observed burst is radiation from a shock front created when the out?owing baryons sweep up surrounding material. The existence of shock fronts suggests high-energy power law emission, but the exact mechanism is still unknown. Given the uncertainty of the ?eld, it was judged prudent to examine the EGRET data for the presence of γ -ray bursts not detected in other wavelengths. One possible detection of a high-energy γ -ray burst event has been found in EGRET data independent of a trigger from BATSE . Two previously detected γ -ray bursts were independently detected, verifying the search method. The new event, occurring on 1994 April 27, was detected with a statistical signi?cance of 95.4% on a 3 minute time scale although this does not include systematic errors. Pointlike structure and the presence of several additional photons which converted in the walls of the EGRET instrument lend qualitative support to the detection.

CHAPTER 3. γ -RAY BURSTS

61

Taken together, several less signi?cant events suggest burst activity on the three minute time scale with probability 99.0%. Activity is suggested at a somewhat lower probability on the ten minute time scale, and at signi?cantly lower probability on thirty minute and one hour time scales.

Chapter 4 Periodic Time-Series Analysis

The full-sky map of the γ -ray sky that EGRET made in its ?rst year of operation provided unprecedented discovery and localization of γ -ray sources (Table 1.1). However, EGRET also provides excellent arrival time information about individual photons. Combined with good directional information, this timing data can be used to look for coherent periodic variation. Six pulsars have been positively identi?ed in EGRET data, and the precise rotation of pulsars makes them an obvious candidate for periodic analysis. X-ray binary systems also exhibit regular periodic behavior in the X-ray bands; it may be hoped that some of these could be detected by EGRET as well. Since there have been extensive EGRET observations of pulsars, we will begin by looking at pulsar models for which EGRET (or GLAST ) may have discriminating power, as well as the contributions made to pulsar understanding by EGRET so far. We will develop the necessary statistical apparatus to search for periodicity in a γ ray signal. We will then see how to apply likelihood statistics to pulsar analysis, and the results of that analysis. Finally, we will see how the statistical tools for periodic time-series analysis may be applied to X-ray binary searches.

62

CHAPTER 4. PERIODIC TIME-SERIES ANALYSIS

63

4.1

Pulsars

Two years before the Vela satellites would begin to secretly detect γ -ray bursts, Jocelyn Bell and Anthony Hewish detected a periodic radio signal that pulsed every 1.3 s [70]. The signal remained coherent for days at a time. However, the mysterious source of the signal would be divined before the γ -ray bursts would be declassi?ed. The discovery of the Vela pulsar, with a period of 89 ms [106], and the Crab pulsar, with a period of 33 ms [183], would narrow the ?eld of prospective sources to one: a rapidly rotating neutron star. Such an object is the only astrophysical object capable of radiating the power observed at the pulse frequencies observed without ?ying apart under its own rotation. The basic description of the pulsar was worked out over the next few years. It had long been realized [24] that a su?ciently massive body would overcome its electron degeneracy pressure by self-gravitation. This situation can come about when a massive star (>M⊙ ) exhausts its nuclear fuel. Without the energy released by the fusion reactions in the core, the star collapses, the gravitational potential condenses the electrons and protons, and a central core of neutrons develops. The remaining outer material is blown o? in a catastrophic explosion: the supernova. Conservation of some fraction of the original stellar angular momentum results in the central neutron star rotating with a frequency of 0.1–100 Hz. Conservation of the magnetic ?ux in the star (assuming it remains a perfect conductor) results in a magnetic ?eld at the neutron star surface of order 1012 G. The observed ?ux from pulsars, then, is the result of a rapidly rotating magnetic dipole (e.g., [43]). Some basic pulsar parameters may be inferred from elementary mechanics. If we make the assumption that all radiated pulsar power comes from the rotational energy,

1 then the power is given by the time derivative of the rotational energy 2 I ?2 , where

I is the moment of inertia and ? is the rotation frequency. Converting this to the pulsar period and period derivative, we have ˙ = (2π )2 I p/p E ˙ 3 (4.1)

Making some basic assumptions about a pulsar as a rotating dipole, Ostriker &

CHAPTER 4. PERIODIC TIME-SERIES ANALYSIS

64

Gunn [153] ?nd that the radiated power can be related to the magnetic ?eld:

2 6 4 ˙ ?2 B a? E 3 c3

(4.2)

where B is the surface magnetic ?eld and a is the radius of the pulsar (typically 10-15 km). Combining these two, we arrive at an estimate of the ?eld strength in terms of the period and period derivative 2 c3 B 2 = (2π )2 6 Ipp ˙ 3 a (4.3)

Typical values of pulsar parameters are given in Table 4.1. For any given pulsar, we assume that its radius and moment are constant. If we assume also that its magnetic ?eld remains constant, we arrive at the di?erential equation pp ˙ = k governing the pulsar’s age. This equation is integrable, yielding 1 2 (p ? p 2 0 ) = kτ 2 (4.4)

where p0 was the pulsar period at time zero, and k is a constant formed from the parameters in equation (4.3). By convention, we assume that the ?nal period is much larger than the initial period, and, solving for τ above, designate the pulsar’s characteristic age by 1 τ = p/p ˙ 2 (4.5)

The characteristic age seems to be a reasonable approximation to the actual age for the only pulsar we can verify. The Crab pulsar was born in a supernova which was observed at Earth on July 4, 1054 [118, 197]. The Crab rotates with a period of ?33.39 ms, and has a period derivative of ?4.21 × 10?13 , yielding a characteristic spindown is probably least accurate at very short times after the supernova explosion. This e?ect will diminish with pulsar age, while glitches in the pulsar period and its derivative will add error to our age estimate. One possible mechanism is gravitational wave emission due to mass multipoles [153]. Ostriker & Gunn ?nd that to fully age of about 1280 yrs. This is correct to within 30%, but the assumption of linear

CHAPTER 4. PERIODIC TIME-SERIES ANALYSIS

65

Parameter Period p Period derivative p ˙ Magnetic Field Characteristic Age Moment of Inertia ˙ Radiated Power E

Typical Values 1.6 ms–6 s 10?15 –10?13 s/s 1012 –1014 G 103 –106 yr ?1045 g cm2 ?1031 –1034 erg s?1

Table 4.1: Typical pulsar parameters for γ -ray pulsars. account for the age di?erence, approximately 1/6 of the Crab’s current luminosity would be in gravitational waves. Some elementary physics gives a reasonable ?rst approximation to a pulsar model. Nevertheless, physics is a fractal subject, and the number of questions still to be answered is seemingly independent of the number of questions already answered. Fortunately for EGRET , γ -rays turn out to be an excellent window on pulsar energygeneration mechanisms.

4.1.1

EGRET Contributions

The EGRET instrument along with the other instruments on board the CGRO has revolutionized our understanding of at least one class of pulsars. Extensive reviews of those contributions [43, 46, 146, 152, 191, 194] have been published elsewhere; we will only review their highlights. Six pulsars have been unambiguously discovered in high-energy γ -rays. All six share some common features. First of all, they tend to have double-peaked light curves, with the peaks separated by somewhat less than 180? . γ -Ray pulsars tend compared to the average γ -ray source spectral index of ?2.3). Five of the six have 2 ˙ among the highest measured values of E/D , the radiated power divided by the pulsar may be calculated by assuming that all of the kinetic energy lost from the pulsar spin-down is radiated isotropically, and that the radio dispersion measure gives a good distance estimate to the pulsar. to have harder spectra than most sources (spectral indices range from ?1.4 to ?1.8,

2 ˙ distance squared. E/D should be proportional to the observed power at Earth. It

CHAPTER 4. PERIODIC TIME-SERIES ANALYSIS

66

The Crab pulsar was detected by both SAS 2 [189] and COS B [8, 29] in γ -rays some ten years after its initial radio discovery [183]. It was strongly detected in many di?erent EGRET viewing periods, yielding high quality light curves and spectra [151, 202, 43, 44]. All γ -ray pulsars are observed to have very steady luminosities; in fact, they have been used as constant calibration sources for some variability studies [121]. Early observers suggested, however, that the ratio of the luminosity in the two γ -ray light curve peaks of the Crab pulsar varied with a 14-year cycle [213, 155]. However, further observations ruled out signi?cant variation on timescales of order 14 years [200]. When it was ?rst discovered in radio observations [106], no one expected the Vela pulsar to be the brightest object in the γ -ray sky. Nevertheless, it was the strongest source detected by SAS 2 [190], and remains the brightest steady source in the EGRET catalog [65]. Although there is signi?cant timing noise (that is, the pulsar period ?uctuates on short time scales), Vela is easily detected in a number of EGRET viewing periods [94]. The most fascinating γ -ray pulsar is Geminga. The third brightest object in the γ -ray sky after Vela and the Crab, it was unique when discovered by SAS 2 in 1975 in that it had no obvious counterpart in other wavelengths [98, 40]. Most notably, it had no radio counterpart—a prominent feature of all other known pulsars. It would be 17 years before pulsation was discovered in the X-ray emission [58]. Using the X-ray ephemeris, pulsation was quickly detected in EGRET [11], SAS 2 [114], and COS B [13] data. Indeed, the γ -ray signal is so strong that various period searching techniques would have independently discovered the pulsation in Geminga without the aid of an X-ray ephemeris (§4.2.3, [116]). Recently, there has been an uncon?rmed report of pulsation detected in radio [102]; whether or not this detection can be con?rmed, the radio emission is clearly very weak. The remaining three γ -ray pulsars are less well-known, but no less important to our understanding (or lack thereof) of γ -ray pulsar physics. PSR 1706-44 was only discovered as a radio pulsar in 1992 [79] and soon thereafter in the EGRET data [192]. Unlike the other ?ve γ -ray pulsars, which have two peaks separated by 140? –180? in phase, PSR 1706-44 exhibits a single peak in its γ -ray light curve.

CHAPTER 4. PERIODIC TIME-SERIES ANALYSIS

67

Unexpected in its own way is PSR 1055-52. It was discovered as a γ -ray pulsar in 2 ˙ EGRET data in 1993 [45], although it ranks approximately 29th in known E/D [43]. Why PSR 1055-52 is observed when other, apparently more luminous pulsars are not 2 ˙ remains a mystery. As a ?nal confusion, PSR B1951+32—fourth in E/D —was not observed as a pulsed source in EGRET Phase I. Further pointed observations in Phase III resulted in a weak pulsed signal [167]. EGRET observations of γ -ray pulsars pose a unique challenge for theorists. The small number of pulsars and the wide variety of observed behavior makes it di?cult for any one model to match all the observations while remaining simple enough to make concrete predictions. Nevertheless, two main classes of models dominate the thinking of most pulsar theorists today.

4.1.2

Models

The theory and rami?cations of γ -ray pulsar models are extensively discussed elsewhere [112, 43, 191]; we will review the highlights of the two major models. Most current γ -ray pulsar models may be categorized as either “polar-cap” or “outer-gap” models. These two classes di?er mainly in the region where the energy is emitted. Polar-cap models The polar-cap models build on the early pulsar models of Goldreich & Julian [54], which de?ned the standard pulsar framework. Although many researchers have extended their work, Daugherty and Harding are largely responsible for developing the model for γ -ray emission [61, 32]. The “polar cap” is a region at the magnetic pole of the pulsar, which in general is not aligned with the rotation axis. Some magnetic ?eld lines emerge from one pole and re-enter the pulsar at the other pole. However, since the pulsar is rapidly rotating, there is some radius (given by r = ?/c) at which the ?eld lines can no longer corotate with the pulsar without exceeding the speed of light. This radius de?nes the light cylinder. The closer the point from which the ?eld line emanates is to the magnetic pole, the closer the ?eld lines will approach the light cylinder. Field lines which emanate from within some critical radius will not return

CHAPTER 4. PERIODIC TIME-SERIES ANALYSIS

68

to the pulsar; instead, they will close at in?nity. The polar cap is that region around the magnetic pole inside the critical radius. The polar cap is important because free charges placed on open ?eld lines will escape to in?nity. Potential di?erences across the polar cap can rip charge o? of the pulsar surface. This charge then ?ows along the curved ?eld lines, emitting curvature radiation. At the lowest altitudes above the pulsar surface, this curvature radiation is so energetic that it pair-produces in the magnetic ?eld. Eventually, the charged particles have lost enough energy that the photons (γ -rays) no longer pair produce. These γ -rays are observed by EGRET . The details of the ?eld shapes and the emission regions determine the γ -ray and radio light curves. Currently, it seems as if many γ -ray pulsar characteristics can be explained with the polar cap model. However, it has been suggested that general-relativistic frame-dragging e?ects may more strongly in?uence the γ -ray emission than the precise shape of the magnetic ?eld [144, 143]. Outer-gap models A substantially di?erent model of energy emission is the “outer-gap” model, ?rst proposed by Cheng, Ho, & Ruderman [26, 27]. In these models, vacuum gaps may arise along the last closed ?eld lines in the pulsar magnetosphere. The gaps separate charge of di?erent sign on opposite sides of the last closed ?eld line, causing signi?cant electric potentials. Charged particles may be accelerated across these potentials, following the magnetic ?eld and emitting curvature radiation. The details of the shapes of the emitting regions and the resultant light curves are not immediately obvious, but extensive work has been done to calculate the impact of the magnetosphere geometry on the γ -ray observations [171, 215]. The outer-gap models also seem to describe many of the qualitative features observed in γ -ray pulsars. Remaining Questions It is still not clear which model of pulsar emission more accurately describes γ -ray pulsars. In general, γ -rays emitted from outer-gap regions will be emitted in a more fan-like pattern, while γ -rays emitted from polar-caps will be more tightly beams.

CHAPTER 4. PERIODIC TIME-SERIES ANALYSIS

69

However, both models can reproduce a subset of the pulsar characteristics observed in EGRET data [152, 194, 44, 43]. Unfortunately, the sample of observed γ -ray pulsars is still small, and it is di?cult to determine which observed characteristics are generally representative and which are particular to the observed pulsars. It is therefore critical to increase the number of observed γ -ray pulsars. The GLAST mission will bring a superior γ -ray telescope to bear on the problem (Chapter 6). However, it is possible that there are additional Geminga-like radio-quiet pulsars which have already been observed by EGRET . Our task is to ?nd likely pulsar candidates, and determine their pulsation parameters.

4.2

Statistical Methods

Searching for pulsation in EGRET sources poses a unique set of problems to the pulsar researcher. Dispersion, the scourge of radio pulsar searches, is non-existent in γ -rays. However, since the ?ux is quantized into γ -rays, EGRET might detect a photon from a pulsar only once per thousand rotational periods. Long integrations of coherent pulsar emission are required to resolve the periodic behavior. Furthermore, the particular structure of EGRET data, organized by individual photon information, is signi?cantly di?erent than the ?ux information found in other wavelengths. Low count rates mean that the instrument state can vary signi?cantly between the arrival of individual photons. The low rates also require that as much information as possible be extracted from each photon. For these reasons, special statistical care must be taken in analyzing EGRET data for pulsar signals. Extensive coherent pulsation analysis has been already been done, and it will be instructive to examine the challenges speci?c to EGRET previously overcome (§4.2.1). First, however, we must approach some di?culties basic to all pulsation analysis. In order to be sensitive to pulsed emission from the fastest pulsing sources known (a few milliseconds), we must have timing accuracy to a fraction of the pulse period. EGRET records photon arrival times to an absolute accuracy of 100 ?s, two and a half orders of magnitude smaller than the period of the Crab. However, the radius

CHAPTER 4. PERIODIC TIME-SERIES ANALYSIS

70

of the Earth’s orbit is 8 light-minutes. Photons from source observations separated by six months time will have light travel times di?ering by hundreds of thousands of pulse periods. In order to compensate for the motion of the Earth, we will translate all measured arrival times into Solar System Barycenter (SSB) times. The SSB time is the time when a photon would have arrived at the solar system barycenter if none of the mass of the solar system were present. The conversion between measured time (tUTC ) and barycenter time tb has been worked out in detail elsewhere [43, 188]. The expression will be motivated and examined in Appendix A; the schematic result is: tb = tUTC + ?convention + ?location + ?Einstein + ?Shapiro (4.6)

where adjustments to Universal Time are made to correct for bookkeeping conventions such as leap seconds, for the current location of the spacecraft with respect to the barycenter, for the gravitational redshift (“Einstein delay”) e?ects, and for the delay due to the gravitational potential (“Shapiro delay”). The full result is given by equation (A.5).

4.2.1

Previous Methods

Two methods have been used to examine coherent periodicity in EGRET data. The ?rst one, used for analysis of pulsars with known ephemerides, was implemented by Joe Fierro in a program called PULSAR [42]. PULSAR selects photons within an energydependent cone about the source position. The acceptance radius, chosen to maximize the signal-to-noise ratio, is given by equation (1.1). The pulsar period, along with the ?rst and second period derivatives, are used to assign to each photon a phase, corresponding to the fraction of a rotation through which the pulsar has turned since a reference time t0 . Recall that a γ -ray will arrive, on average, only once per hundred or thousand rotations of the pulsar; by retaining only the fractional rotation, the emission from many rotations is added together coherently. This process is known as epoch folding. Plotting the resulting photon rate as a function of phase yields a light curve—the average pro?le of emission through a rotation. Photons can then be selected as a function of phase for further analysis.

CHAPTER 4. PERIODIC TIME-SERIES ANALYSIS

71

The phase of a photon is determined from the pulsation period and its derivatives, as measured from some reference epoch: φ(t) = φ(t0 ) + ν (t ? t0 ) + 1/2ν ˙ (t ? t0 )2 + 1/6¨ ν (t ? t0 )3 + . . . (4.7)

where ν = 1/p is the pulsation frequency and p the period. In practice, it is only rarely that the second-derivative term of the Taylor expansion is important; higher orders are negligible. The constant term in the expansion is kept only to allow γ -ray light curves to be compared with radio light curves. Of primary importance in such an analysis is the determination of the presence of intensity modulation at the known frequency. Ephemerides for radio pulsars have been collected into a database stored at Princeton University [187]. γ -Ray photon phases based on these ephemerides may be tested for periodicity. The χ2 test The simplest test is the χ2 test. Photons are sorted into m bins, based on their phase. Since photon count statistics are distributed as a Poisson, the variance in the number of photons in each bin is Nm , the number of photons in the bin. The usual χ2 statistic is calculated, and the probability of the observation under the null model is computed. A su?ciently low probability of observation under the null model is taken as evidence for variability. While this method is simple, it has a number of undesirable aspects. First, the signi?cance of the detection depends strongly on accidents of bin boundaries—a problem inherent to binning schemes. Second, the χ2 is invariant under permutations of bins. We expect pulsar light curves to smoothly vary, with one or two main peaks, and would like our statistical method to have some power to identify such light curves.

2 The Zm test 2 The Zm test was developed to try to address some of these di?culties [9, 21, 43].

The basis of this test is the realization that the light curve is well represented by the phase density f (φ), which is the fraction of the integrated intensity in each dφ. In

CHAPTER 4. PERIODIC TIME-SERIES ANALYSIS

72

its simplest form, corresponding to m = 1 and known as the Rayleigh statistic [113], this phase density will be compared to a sinusoid of unknown amplitude and phase. The general form uses more terms of the Fourier expansion. The best estimate of the phase density function is a series of δ -functions: f (φ ) = 1 N

N

δ (φ i )

i=1

(4.8)

We will de?ne fm (φ) as the function composed of the ?rst m Fourier components of

2 f (φ). To calculate the Zm statistic, we determine a comparison (or model) phase ? ? density function f m (φ). We then ?nd the squared di?erence from fm (φ) and fm (φ) 2π 0 2 ? fm (φ) ? f m (φ) dφ

2 Zm ≡ 2πN

(4.9)

If we take the comparison density function to be constant, this is simply computed; it can be shown to be equal to 2N odd Fourier coe?cients from fm (φ). degrees of freedom [7]. This approach clearly mitigates some of the problems with the χ2 test. Primarily, photon binning has been eliminated. In addition, the Fourier decomposition expects a periodic signal—the degeneracy to permutation of ?ux peak location is removed.

2 Unfortunately, the Zm statistic has shortcomings of its own. The signi?cance of a m 2 2 k =1 (αk + βk ), where αk and βk are the even and 2 Zm is asymptotically distributed as χ2 with 2m

detection is a strong function of the number of harmonics chosen. It is not di?cult to see why this should be the case. Retaining only m harmonics of a sum of δ -functions amounts to smoothing the observed phase density function. Too much smoothing (a very small m) washes out any signal that may be present. Too little smoothing (a very large m) results in the signal being swamped by noise from the higher harmonics. To ameliorate this di?culty, another statistic was invented. The H -test In the quest to obtain more and more signi?cant results, it naturally occurred to astrophysicists to try a variety of values of m, and see which one results in the most

CHAPTER 4. PERIODIC TIME-SERIES ANALYSIS

73

2 2 signi?cant value of Zm . Fortunately, de Jager et al. realized that the set of Zm

obtained in this way would no longer be distributed as χ2 2m . At this point, there is no further guiding insight; an ad hoc limit of 20 is set on m, and the distribution of the term involving m favors models with fewer Fourier components, and the constant search for a pulsed signal in EGRET data [43]. Limitations of these statistics Statistics for analyzing pulsar data have evolved, with each generation of statistics alleviating some of the problems of the previous generation, but creating problems of its own. The discussion of likelihood statistics in §2.3 sheds some light on why these problem will be given in §4.2.2; here we will quickly examine the above statistical measures. Since the H -test is the most sophisticated measure, we will focus on it. If we permit ourselves to use some Bayesian language, we see that the goal of the H -test is to eliminate the nuisance parameter m. The method is really to take H as a

2 function of m, and maximize this H (m) with respect to m. Each Zm has the form 2 H ≡ max1≤m≤20 Zm ? 4m +4 is found from a Monte Carlo simulation. The inclusion of

ensures that the statistic will be positive. This test has become the primary tool to

di?culties arise. A full discussion of the application of likelihood statistics to this

of the squared di?erence between the null model and a representation of the data comprised of m harmonics? . This is just the logarithm of the likelihood of the “data representation,” assuming a null model and Gaussian errors of variance Nm . The constant term only insures positivity; clearly it does not a?ect the distribution of the statistic. The only remaining term contains an m. To see its signi?cance, let us write the exponentiation of H (m)—that is, instead of the logarithm of the likelihood, the

Note that each representation, labeled by m, is derived from the data, which is not permitted in Bayesian analysis. Therefore, the likelihood formed is not “the likelihood of the data, given the model,” but “the likelihood of the representation, given the model.” Since the representation is a (lossy) smoothing of the data, we have no reason to expect a one-to-one relationship between the likelihoods. Nevertheless, the qualitative argument is still valid—this only means that multiplying 2 the Zm by a prior and integrating over m would still give incorrect numerical results.

?

CHAPTER 4. PERIODIC TIME-SERIES ANALYSIS

74

likelihood itself.

2 exp H (m) = exp(Zm ? 4m + 4) = (e4 e?4m ) exp 2 fm (φ) ? f? m (φ) dφ

(4.10)

The integral is the likelihood; it is the exponential of a sum of squared di?erences, which of course is a Gaussian. The ?rst term, e?4m , is readily identi?ed as the prior! The Bayesian would now integrate? H (m) over m to eliminate the unnecessary parameter m as described is §2.6.3. However, if H (m) is at all peaked as a function of m, then eH (m) is even more peaked, and the integral of eH (m) over m is very nearly m of the likelihood times a prior of e?4m . This example makes the comparison of the frequentist and Bayesian schools quite explicit. The Bayesian makes arbitrary, subjective choices of priors. The frequentist makes arbitrary, subjective choices in the method of smoothing—for example, it is not clear that the basis of sines and cosines is the optimal smoothing for the H -test. The advantage to the Bayesian method is that it makes explicit when the subjective choices are to be made, and completely prescribes the rest of the analysis process. As we have seen, the further disadvantage is that the scientist creating new statistics cannot always consider all the rami?cations of his choices. Unfortunately, the H

2 statistic cannot be easily adapted to a Bayesian analysis; for while Zm may appear

equal to the maximum. The H -test, then, is an approximation to the integral over

to be a true likelihood, it is not. The likelihood must represent the probability of the

2 data. Zm represents the probability of fm (φ), which in general, does not have the

same probability density as the actual data. Furthermore, it is not at all clear that e?4m is the most suitable prior. If de Jager et al. had claimed to use such a prior, they would have had to spend some time justifying their choice. However, under the guise of creating an ad hoc statistic, they have surreptitiously inserted that prior into the analysis—probably fooling even themselves. Is the H -test then worthless, or worse yet, misleading? No. Within the constraints of the simulations done to quantify its signi?cance, and within the assumptions that it implicitly makes, it is completely valid. However, the likelihood approach o?ers a

? 2 if only the Zm were a well-formed likelihood.

CHAPTER 4. PERIODIC TIME-SERIES ANALYSIS

75

method which is both simpler to understand, and at least as powerful. The correct generalization of the H -test will be discussed in §4.2.6. Fourier Transforms The statistical tests described above have been used extensively in analyzing EGRET data, but have only had success when a pulsar ephemeris from radio or X-ray observations was used to ?nd the photon phases φi . The Geminga pulsar is the ?rst known “radio-quiet” pulsar, with little or no radio emission despite a strong pulsed signal in X-ray and γ -ray observations [11, 120]. It is natural to search for other such radio-quiet pulsars in the EGRET data. A likelihood method of pulsar searching will be discussed in §4.3; here we will brie?y touch on a Fourier transform method. The traditional method of searching for periodic behavior at an unknown fre-

quency is the fast Fourier transform (FFT). The main advantage of the FFT is its speed. Indeed, such an algorithm has been implemented to search for periodic signals from pulsar candidates [116]. It has been shown to be successful in detecting Geminga—although as we will later see, the signal from Geminga is so strong that this is not a de?nitive test. The drawbacks to the FFT are also signi?cant. First of all, the FFT cannot detect a period derivative. Trial period derivatives must be removed from the photon phases, and a new FFT performed for each period derivative. The FFT also cannot naturally deal with noncontiguous variable exposure and backgrounds. In order to detect the weakest pulsar signals, we will need to identify, to the extent possible, the photons most likely to be source photons, and account for the instrument exposure on a photon-by-photon basis.

4.2.2

Maximum Likelihood

The ideal statistical method for pulsed-signal analysis would be sensitive to variations in pulsed ?ux from the candidate source, rather than pulsed count rates in a region around the source. Given the success of the maximum likelihood method for spatial analysis. Again, we will ?nd that the mathematics may be interpreted in a maximum analysis (§2.4), it seems natural to apply likelihood to the problem of periodic signal

CHAPTER 4. PERIODIC TIME-SERIES ANALYSIS

76

likelihood framework or a Bayesian framework. Following a suggestion by Gregory & Loredo [56], we will generalize the binned χ2 method to a maximum likelihood method. In §4.2.6 we will examine the possibility of generalizing an unbinned method. The ?rst thing we must do is to construct the likelihood—the probability of ob-

serving our data, given our model and the state of the instrument. Since the state of the instrument and the background are di?erent for every photon, we ?rst divide our time interval and spatial intervals into su?ciently small subintervals that there is either zero or one photon in each subinterval. We then have an array in parameter space of volume elements ?t ?? ?b. The probability that a given element has no photons in it will be the Poisson probability of detecting zero photons during a time ?t from a direction element ?? ?b with rate r (t, ?, b), which depends in general on time and observed position. The probability of detecting one photon follows similarly: p0 (t, ?, b) = e?r(t,?,b)?t???b p1 (t, ?, b) = r (t, ?, b) ?t???b e?r(t,?,b)?t???b (4.11)

We identify each element in t, ?, and b with a sequential, unique label. Furthermore, let us suppose that N of these parameter-space volume elements contain one photon, and Q of them contain no photons. N + Q is then the total number of parameter-space volume elements, and N is the total number of photons. We may ?nd the likelihood of this con?guration by multiplying the probabilities of each element: L= p1 (ti , ?i , bi )

i∈α1 j ∈α0

p0 (tj , ?j , bj )

(4.12)

where α0 is the set of Q volume elements with zero photons, and α1 is the set of N volume elements with one photon. Plugging in our expressions from equation (4.11), we may factor out the exponential and arrive at L = (?t???b)N r (ti , ?i , bi ) exp ??

?

N +Q j =1

i∈α1

r (tj , ?j , bj ) ?t???b?

?

(4.13)

The sum in the exponential is over all volume elements. It is simpler to separate this

CHAPTER 4. PERIODIC TIME-SERIES ANALYSIS

77

sum into three independent summations, one for each parameter. In the limit that ?t, ??, and ?b are very small, these sums become integrals. Now, the i are just arbitrary labels. We renumber the volume elements so that elements 0 through N contain one photon, and we have L = (?t???b)N

N i=1

r (ti , ?i, bi ) exp ?

t

?

b

r (t, ?, b) dt d? db

(4.14)

Note that nowhere have we assumed that the ?t are contiguous. Thus, this likelihood expression is valid even if there are gaps in the data. The rate r (ti , ?i, bi ) is a function of both the source strength and background. We assume that the background is constant in time, and varies in space according to the standard EGRET gas map [10, 75]. The positional dependence of the source rate will be taken from the instrument point-spread function. We will thus take the rate to be ?(ti ) + B (ti , ?i , bi ) r (ti , ?i , bi ) = A(ti , ?i , bi )f (4.15)

? is the time-dependent model source ?ux as before, A(ti , ?i , bi ) represents the where f instrument response function, and B (ti , ?i, bi ) is the background ?ux. We will worry about the detailed forms of A and B below; however, we note here that the time dependence of A and B is a result of the time dependence of the instrument sensitive area. To simplify the products and exponents, we take the logarithm of the likelihood:

N

ln L = N ln(?t???b) +

i=1

?(ti ) + Bi ) ? ln(Ai f

?

b

t

?(ti ) + Bi ) d? db dt (4.16) (Ai f

where Ai and Bi are shorthand for A(ti , ?i , bi ) and B (ti , ?i , bi ). Now, eventually we will either maximize the likelihood, or ?nd a ratio of likelihoods. Thus, constant additive terms in ln L may be dropped. We drop the ?rst term and separate the integral to get

CHAPTER 4. PERIODIC TIME-SERIES ANALYSIS

78

N

ln L =

i=1 ?

?(ti ) + Bi ) ? ln(Ai f

b t

?(t) d? db dt ? A(t, ?, b)f

t

?

b

B (t, ?, b) dt d? db

(4.17)

Although the last term is (implicitly) a function of source position, it will be constant for our purposes, as we take source position to be given. Since it is constant, we drop it. In the second term, we may split A into PSD(?, b)SA(t), the point-spread density and the sensitive area. The point-spread function is normalized, by construction, for all times, and SA(t) depends only on the source position, not the photon positions. Therefore we may split the integral into ? A′ (?, b) d? db ?(t)SA(t)dt f (4.18)

?

b

t

The spatial integrals are normalized, and the resulting likelihood function is given by

N

ln L =

i=1

?(ti ) + Bi ) ? ln(Ai f

t

?(t)SA(t)dt f

(4.19)

? is the model source ?ux. It is this model that we will vary to Remember that f ?nd the model most likely to have produced our data. The data itself appears only implicitly in equation (4.19), in the Ai and Bi , and most importantly in the ti at which the model ?ux is sampled.

4.2.3

Application to EGRET

?. It is tempting to follow the lead We have not yet selected a functional form for f 2 ? be a function composed of some number of harmonics, statistic, and let f of the Zm perhaps with parameters to be determined. We will delay such an attempt to §4.2.6, and ?rst develop a method based on χ2 which we will ?nd to be less computationally demanding. ? that corresponds to χ2 binning is stepwise constant, with The functional form of f m di?erent steps. The immediate bene?t of such a choice comes in the simpli?cation

CHAPTER 4. PERIODIC TIME-SERIES ANALYSIS

79

of the second term of equation (4.19): ?(t)SA(t)dt ?→ f

m j =1

t

? f j τj

(4.20)

where τj is the sensitive area integrated over the livetime in bin j ; that is, the total exposure in bin j . Now equation (4.19) may be stated as a sum over bins

m

ln L =

j =1

? ?

i∈bin j

?j + Bi ) ? f ?j τj ? ln(Ai f

?

(4.21)

? Thus the maximum of ln L may be found by independently maximizing the f j for ?—an enormous each bin, instead of simultaneously maximizing all the parameters of f computational savings. The analysis is thereafter quite simple, and quite slow. A candidate pulsar is selected for analysis. All photons within some radius, usually 15? , are collected. Each photon with su?cient energy (>100 MeV) which arrives from a position su?ciently far from the Earth’s limb (<95? from the instrument axis) during a valid instrument mode is accepted for analysis. For each of these photons, the arrival time is corrected to the Solar System barycenter, and the phase-independent likelihood factors Ai and Bi are calculated. For periods less than a minute, it may reasonably be assumed that the instrument exposure will be evenly distributed across the bins. For longer periods, such as X-ray binary orbital periods, all instrument exposure during the observation is calculated as a function of SSB time. To ameliorate the primary objection to binning, namely, that the arbitrary placement of the bin boundaries may obscure a signal, we may allow an o?set to the bin boundary. In practice, this means rebinning the data for a number of di?erent boundary o?sets. Of course, such an operation increases the number of trials. Sampling Density A subset of parameter space is identi?ed to be searched for pulsed signals. Some boundaries of the range of period and period derivative are set by the physics of the

CHAPTER 4. PERIODIC TIME-SERIES ANALYSIS

80

system being examined; those will be motivated in §4.3 and §4.4. The sampling size required is determined by the statistics. To determine the number of steps required to sample the period space, we consider the di?erence in phase of the last photon in the observation as a function of the change in trial periods. If two trial periods di?er by so little that the change in phase of the last photon as computed by equation (4.7) is less than the width of a bin, then the most likely fj ’s computed by maximizing equation (4.21) will be unchanged. So, we choose the minimum spacing so that the phase of the last photon will change by an amount of order the bin width, and calculate the number of steps to take in period across the desired range. Np = (νmax ? νmin ) m δt (4.22)

where δt is the length of the observation. Note that a longer observation means more photons will be observed, increasing the signal-to-noise ratio, but the period space must be searched more densely. Similarly, we can ?nd the minimum period derivative we expect could make a di?erence. In this case, our condition is that φ(δt) ? νδt is an appreciable fraction frequency derivative ν ˙ = p/p ˙ 2 , we have

1 of a bin. Neglecting higher derivatives, this happens when 2 νδt ˙ ? p/m. Since the

p ˙min =

2 p3 m(δt)2

(4.23)

The maximum period derivative must be determined by the physics. With the range of period derivatives for a given period in hand, we may compute the number of samplings in period derivative required. Again, we ?nd the step size by equating the di?erence in computed phase with the bin size: 1 1 ν ˙ (δt)2 ? (ν ˙ + δν ˙ )(δt)2 = 1/m 2 2 and thus the number of steps required in the period derivative: 1 ˙ max ? ν ˙ min ) m (δt)2 Np ˙ = (ν 2 (4.25) (4.24)

CHAPTER 4. PERIODIC TIME-SERIES ANALYSIS

81

Signi?cance Estimates For each point in parameter space, the ?ux in each bin is estimated by maximizing the likelihood. Since the sampling of parameter space was chosen to make each point roughly independent, the number of points sampled in parameter space should be a good estimate of the number of independent trials. The signi?cance is then the probability of measuring no likelihoods greater than or equal to that observed in the number of trials. Automated analysis with timevar These methods have been automated into a tool for searching for periodicity in EGRET data called timevar [81]. Capable of running either interactively or in batch mode, timevar searches a given region of parameter space, for a given number of bins. It maximizes the likelihood for each period and period derivative, and reports the signi?cance of the most likely period and period derivative. In addition, it yields a light curve, both as the most likely ?ux in each bin, and for comparison with PULSAR, the number of photons within the 68% containment radius in each bin. To demonstrate the capabilities of timevar, we examine the Geminga pulsar. We will check viewing period 413.0, corresponding to observations between 7 March 1995 16:44:00 and 21 March 95 14:08:39. A small range around the known period was searched, and the results are presented in Figure 4.1. The light curves are shown in Figure 4.2. The probability of the observation in the null model is less than 10?25 . The great advantage of timevar over previous analysis methods is that it combines period and period derivative searching ability with a proper likelihood treatment of the source and background as perceived through the instrument point-spread function and sensitive area. PULSAR has been useful to examine known pulsars with radio ephemerides, although it does not optimally use the point-spread and background information about each photon. The disadvantage to timevar is its computational expense. Examining a single pulsar candidate over a reasonable range of period and period derivatives currently takes on the order of months of continuous running on a Sun Microsystems Sparc 10. Mattox et al. [116] have sidestepped this problem ?rst by

CHAPTER 4. PERIODIC TIME-SERIES ANALYSIS

82

Figure 4.1: The likelihood of periodic modulation from Geminga, as a function of period, for a 30-bin light curve. Plotted on the y -axis is 2 ln L, which is distributed as χ2 29 , since the average intensity in the null model is a free parameter. The probability of the peak in the null model is less than 10?25 . The maximum likelihood value of the period is 0.237097785 s, assuming a period derivative of 1.09744 × 10?14 s/s and a truncated Julian day epoch 8750.0. turning to a Fourier algorithm—which is much faster, but cannot properly account for backgrounds and point-spread functions—and second by using a massively parallel computer. Further optimization may be possible to speed timevar, and certainly additional computing power will reduce the required search time. A similar program designed for GLAST will be slowed by the additional order of magnitude in the number of photons to process, but since the sensitive area will be larger, the elapsed time required to observe a su?cient number of photons will be shorter, allowing period space to be sampled less densely. In addition, faster computers will be available to do the calculations. The net e?ect will be that period searching with the GLAST equivalent to timevar will be quite feasible (Table 4.4).

CHAPTER 4. PERIODIC TIME-SERIES ANALYSIS

83

Figure 4.2: Light curves of the most likely period as found in Figure 4.1. The light curve on the left is the most likely ?ux in each bin, after accounting for the instrument point-spread function and sensitive area and the background. The light curve on the right is the raw photon count inside the 68% containment radius, as given by equation (1.1). This is the same as the light curve used by PULSAR.

4.2.4

Bayesian Inference

A further advantage of the likelihood method employed by timevar is the ease by which results can be evaluated under a Bayesian framework. Strictly speaking, the plot shown in Figure 4.1 is meaningless in a maximum likelihood framework—only the maximum value, and not the entire likelihood function, is relevant. However, intuitively we expect the rest of the function to have some meaning. Local maxima should represent frequencies more likely to be present in the data. We would expect formulation includes all this information into the posterior probability distribution. We simply multiply the likelihood function given in Figure 4.1 by the prior probability distribution. Then we may integrate the posterior distribution over interesting intervals, or simply present it as it is. Appropriate Priors The least-informative prior probability distribution for the pulsation period is scaleinvariant. That is, it should not depend on the units used for measuring period, and it should be the same whether the search is done in period or frequency space. The appropriate solution is ?at in the logarithm of the period. Therefore, we have harmonics of the main pulse frequency to be present. As we saw in §2.6, the Bayesian

CHAPTER 4. PERIODIC TIME-SERIES ANALYSIS

84

Pprior (p) ∝ 1/p. The upper and lower bounds of the period search de?ne the cuto?, so that the prior can be normalized: pmin pmax 1 p

Pprior (p) = ln

(4.26)

The least informative prior for the period derivative is uniform. The period derivative is a dimensionless number, and thus is the most natural unit in which to work. The last parameter of interest is the phase o?set of the bin boundaries. This is also dimensionless, and thus the least informative prior is ?at. Gregory & Loredo [56] point out that marginalization over m yields the best light curve independent of the number of bins. This is exactly the spirit of the H -test. However, this is computationally impossible for a period search given the present state of the art. It is interesting to note here that the proper treatment of this problem is very di?cult using frequentist statistics. Two models with di?erent m have a di?erent number of free parameters. The Bayesian framework takes care of this naturally; each free parameter (that is, each of the m ?uxes) has its own prior; thus the complete prior for all the ?uxes (assuming a ?at prior between zero and the large cuto?) will ?m ? be Pprior (f . This prior is a quantization of Ockham’s Razor: models with m) = f

max

more free parameters (that is, larger m) are discouraged. The frequentist method compares each model to the null model, but cannot directly compare the two. Figure 4.1 is actually one of three di?erent likelihood functions measured for different phase o?sets. It is the slice of the likelihood corresponding to the most likely o?set. The perpendicular slice through the likelihood as a function of o?set ?xed at the best period is shown in Figure 4.3. The phase o?set is a classic nuisance parameter; it arose by the arbitrary choice of t0 in equation (4.7). We therefore marginalize over the o?set and arrive at the likelihood function shown in Figure 4.4. To ?nd the overall signi?cance of a detection? in the Bayesian framework, we multiply the likelihood function shown in Figure 4.4 by the period prior given in equation (4.26). The ratio of this probability with the probability of the null model

In our case, this will be the overall signi?cance of a detection for all models with m bins. In principle, however, the parameter m could be marginalized over as well, yielding the true overall detection signi?cance.

?

CHAPTER 4. PERIODIC TIME-SERIES ANALYSIS

85

Figure 4.3: The likelihood of pulsation in Geminga as a function of the phase o?set of the bin boundaries, for ?xed period. The period chosen is that of the maximum likelihood for any o?set, and corresponds to the peak seen in Figure 4.1. The width of the plot corresponds to one bin; therefore, the leftmost point and the rightmost point are the same. (m=1) gives the odds ratio. In this case, the probability of periodic modulation is completely dominated by the peak, and the signi?cance of the detection is 1?7×10?46 . The Geminga pulsar is clearly detected by either method. The detections are

so strong that even if a very large region of phase space had been searched, the detections would still be signi?cant. Finding the signi?cance that would have resulted from such a search may be calculated by assuming that the peak we found would be the maximum likelihood over the entire search, and that the number of trials would be the number of steps found in §4.2.3. A thorough search of viewing period 413.0 from 50 ?s to 500 ?s, over the appropriate period derivative ranges would require approximately 7 × 108 trials. The signi?cance of the detection in Figure 4.1 would

harsh attack. It is therefore also useful to develop some formalism for upper limits and detection thresholds.

still be 1 ? 10?20 . However, it is clear that weaker detections will not survive such a

4.2.5

Upper Limits and Thresholds

desirable quality of an upper limit is that it be statistically well-behaved; that is, that

The problem of upper limits was brie?y touched upon in §2.7. There we noted that a

CHAPTER 4. PERIODIC TIME-SERIES ANALYSIS

86

Figure 4.4: The likelihood of pulsation in Geminga. This likelihood function is similar to that shown in Figure 4.1, except that the likelihood has been marginalized over the phase o?set. a “1σ upper limit” be greater than the true value in 68% of the ensemble of possible data sets. The value of this upper limit depends on the details of the data. In contrast, in some cases it may be more useful to calculate a threshold. We de?ne a threshold to be the intensity (of a source, or of pulsation) that would be detected by our criteria in some fraction of the possible data that might be measured from that source. This may be quanti?ed: an x% detection threshold is that true parameter value for which x% of possible data sets would result in detection. We make a statement about the outcome of a hypothesis test (is there pulsation?) instead of a statement about a point-estimate (is the pulsed fraction less than some value?). There are three advantages to this. First of all, it more closely re?ects the intuitive “upper limit” concept—that is, it is now reasonable to say, “The threshold for detection is f . Since we do not detect, then the actual value is probably less than f .” We will quantify this below. Second, we may calculate thresholds from the underlying distributions of the data; upper limits can only be calculated with the speci?c data

CHAPTER 4. PERIODIC TIME-SERIES ANALYSIS

87

observed. This means that thresholds can be precalculated. This is a great advantage, especially for instruments which are still being designed or built. Finally, since we can use distributions, we can calculated thresholds (at least partially) analytically. For the purposes of pulsation analysis, we will de?ne the “pulsed fraction threshold” as the fraction of the ?ux (background plus source) that must pulse such that a fraction ξ of data sets observed with true sinusoidal variation with threshold amplitude will result in detection. We will assume that any data set with signi?cance greater than α will be a detection and the variation will be sinusoidal. α-point of the integral χ2 distribution. That is, we seek C such that For simplicity, we de?ne C ≡ 2 ln L. We will take ξ = 0.5, so that we may ?nd the

C 0

χ2 (d, t) dx =

α, where d is the number of degrees of freedom, and t is the number of trials.

更多相关文档:

Prospects *of* observing pulsed radiation from *gamma-ray* *pulsars* with H.E.S.S_专业资料。Observations *and* theoretical studies have demonstrated that *the* pulsed...

Constraints on *the* Galactic Corona Models *of* *Gamma-Ray Bursts* From *the* 3B ...have revolutionized our understanding *of* *the* birth velocities *of* radio *pulsars*...

Cordes 2000) *and* EGRET detections, upper limits *and* diffuse background measurements, we *find* *a* best-fit luminosity law *for* *the* *gamma-ray* *pulsar* population...

This brief review describes *the* motivation *for* *gamma-ray* *pulsar* studies*,* *the*...*gamma-ray bursts,* it will have *gamma-ray* capabilities useful *for* *pulsar* ...

ApJ, in press NSF-ITP/95-128 ARE *GAMMA-RAY BURSTS* DUE TO ROTATION POWERED HIGH VELOCITY *PULSARS* IN *THE* HALO ? Dieter H. Hartmann1;3 *and* Ramesh ...

including *gamma-ray bursts*, *pulsars,* blazars, interstellar clouds, *and* *a* ...2. X-Rays: *Searches* *for* Counterparts from *the* Highest Energies Down As ...

Long *and* short *gamma-ray bursts,* *and* *the* *pulsar* kicks_专业资料。One *of* *the* mysteries that surround *gamma-ray bursts* (GRB) is *the* origin *of* two ...

PulsarSpectrum simulating *gamma-ray* *pulsars* *for* *the* GLAST mission_专业资料。...*the* three CGRO experiments*,* *the* *Burst* *and* Transient Source Experiment (BATSE...

Dynamos, Super-*pulsars* *and* *Gamma-ray bursts* *The* remnant *of* *a* neutron star binary coalescence is expected to be temporarily stabilised against gravitational ...

Very High Energy *Gamma* *Rays* from *the* Vela *Pulsar*Nebula_专业资料。We have observed *the* Vela *pulsar* region at TeV energies using *the* 3.8 m imaging ...

更多相关标签:

相关文档

- A search for x-ray counterparts of gamma-ray bursts with the ROSAT PSPC
- The search for possible mesolensing of cosmic gamma-ray bursts. Double and triple bursts in
- Dynamos, Super-pulsars and Gamma-ray bursts
- Attributes of flares in Gamma Ray Bursts sample I
- Scattered Emission from A Relativistic Outflow and Its Application to Gamma-Ray Bursts
- OG.2.3.17 Status and Results of the Search for Gamma-Ray Bursts above 1 TeV with the HEGRA
- STATUS AND RESULTS OF THE SEARCH FOR GAMMA-RAY BURSTS ABOVE 1 TEV WITH THE HEGRA EXPERIMENT
- Supernovae, Pulsars and Gamma-Ray Bursts A Unified Picture
- An Off-line Scan of the BATSE Daily Records and a Large Uniform Sample of Gamma-Ray Bursts
- Prompt Ultraviolet-to-Soft-X-Ray Emission of Gamma-Ray Bursts Application to GRB 031203

文档资料共享网 nexoncn.com
copyright ©right 2010-2020。

文档资料共享网内容来自网络，如有侵犯请联系客服。email:zhit325@126.com