# The Flux Variability of Markarian 501 in Very High Energy Gamma Rays

## 1 Introduction

Three active galactic nuclei (AGN) have been discovered to be very high energy (VHE, E$`>`$300 GeV) $`\gamma `$-ray sources by the Whipple Observatory $`\gamma `$-ray collaboration: Markarian 421 (Mrk 421) (Punch et al. (1992)), Markarian 501 (Mrk 501)(Quinn et al. (1996)) and 1ES 2344+514 (Catanese et al. (1998)). These are the three closest BL Lacertae objects (BL Lacs), with redshifts in the range 0.0308 - 0.044, and are among the brightest at X-ray energies. They are all classified as X-ray selected BL Lacs (XBLs) as their synchrotron spectra extend into the X-ray range. A fourth BL Lac, PKS2155-304, has been detected in VHE $`\gamma `$-rays by the University of Durham group (Chadwick et al. (1999)).

The Energetic Gamma-Ray Experiment Telescope (EGRET) on board the Compton Gamma-Ray Observatory has detected at least 51 AGN at energies $`>`$100 MeV (Thompson et al. (1995); Mukherjee et al. (1997)). They are all members of the blazar class of AGN, which include flat spectrum radio quasars and BL Lacs. Of the EGRET-detected blazars, 14 are BL Lacs with 12 Radio Selected BL Lacs (RBLs) and only two XBLs. Mrk 421 is the only VHE source in this catalogue and it is among the weakest. However, Mrk 501 has recently been detected at the 4$`\sigma `$ level with EGRET (Kataoka et al. (1999)). One of the most striking characteristics of the EGRET-detected blazars is variability; 42 of the 51 AGN exhibit variability (Mukherjee et al. (1997)). Variability time-scales as short as 4 hours have been observed (Mattox et al. (1997)).

VHE observations of Mrk 421 and Mrk 501 have also revealed extreme variability (e.g. Gaidos et al. (1996); Quinn et al. (1996)). The VHE flux has been measured to vary by nearly a factor of 100 in Mrk 421 (McEnery et al. (1999)) and, as we show in the following sections, Mrk 501 has been measured with fluxes ranging from 0.1 to 5 times the Crab Nebula flux with the Whipple Observatory atmospheric Čerenkov telescope. The large collection area ($`3.5\times 10^5`$ m<sup>2</sup>) of the Whipple Observatory telescope permits sensitive studies of variability on time-scales inaccessible to space-based telescopes. Indeed, the shortest observed variability of any blazar at any $`\gamma `$-ray energy, a 30 minute duration flare observed from Mrk 421 (Gaidos et al. (1996)), was measured with this telescope. Both Mrk 421 and Mrk 501 have been closely monitored by the Whipple Collaboration since their discovery, with an $``$0.5 hr exposure/night being sufficient for detection of flaring activity. Prior to 1997, the $`\gamma `$-ray emission from Mrk 421 was generally observed to have a higher mean flux (Schubnell et al. (1996)) and to have been more frequently variable (Buckley et al. (1996)) than that of Mrk 501.

In Spring 1997, the Whipple Collaboration observed Mrk 501 to be in an unprecedented high emission state at VHE energies, as subsequently confirmed by several independent Čerenkov imaging groups (Protheroe et al. (1998)). A public Compton Gamma-Ray Observatory Target of Opportunity was initiated in response to a request by the Whipple Collaboration in 1997. Evidence of correlated variability in data from the Whipple Observatory $`\gamma `$-ray telescope, the Oriented Scintillation Spectrometer Experiment (OSSE) and the All-Sky Monitor (ASM) of the Rossi X-ray Timing Explorer is presented in Catanese et al. (1997). During this observing campaign the energy output of Mrk 501 in VHE $`\gamma `$-rays was comparable to that in the 2-100 keV range but the variability amplitude was larger. The correlations seen may imply some relativistic beaming of the emission, given that the spectrum extends to $`>`$7 TeV (Samuelson et al. (1998)). There was also some indication that the optical U-band flux was higher on average in the month of peak $`\gamma `$-ray activity.

Here, based on 4 years of data, we present a study of the variability of the $`\gamma `$-ray flux above $``$350 GeV from Mrk 501. Details of the observations are given in §2 and the analysis methodology, including $`\gamma `$-ray selection criteria and the test for variability, is described in §3. The results of the analysis are presented in §4 and their implications briefly discussed in §5.

## 2 Source Observations

Observations were made with the Whipple Observatory 10 m atmospheric Čerenkov imaging telescope (Cawley et al. (1991)), located on Mt. Hopkins in southern Arizona. The 10 m reflector images the Čerenkov radiation from cosmic-ray and $`\gamma `$-ray initiated air-showers onto a high resolution camera mounted in the focal plane. Subsequent off-line analysis of the images (described below) facilitates the selection of candidate $`\gamma `$-ray events.

The high resolution camera utilises fast photomultiplier tubes (PMTs) arranged in a hexagonal array, with inter-tube spacing of $`0\stackrel{}{\mathrm{.}}25`$. During 1995 and 1996 the camera consisted of 109 PMTs with a resulting field of view (FOV) of $`2\stackrel{}{\mathrm{.}}8`$. An event was recorded when any two of the inner 91 PMTs registered a signal $`>`$40 photoelectrons within an effective resolving time of 15 ns. For the 1997 observations a further 42 PMTs were added, resulting in a FOV of $`3\stackrel{}{\mathrm{.}}4`$. The trigger condition remained the same as for the 109 PMT camera. In Summer 1997 a new camera, containing 331 pixels, was installed. This camera has a FOV of $`4\stackrel{}{\mathrm{.}}8`$. The trigger condition for this enlarged camera required that any two of the 331 pixels produce a signal $`>`$40 photo-electrons within an effective resolving time of 8 ns. The telescope was also triggered artificially, after every 24 events for the 109/151 pixel cameras and once every second for the 331 pixel camera, to determine the background sky-brightness level in each PMT. Light-cones, which minimise the dead-space between PMTs and reduce the albedo effect, were used on the 109 and 151 pixel cameras, but were not yet installed on the 331 pixel camera.

In general, two modes of observation are used: on/off and tracking. With the on/off mode the source is tracked continuously for 28 minutes and then, to estimate the background, a region offset in right ascension (RA) by 30 minutes (allowing 2 minutes slew time) is tracked. This has the disadvantage that an equivalent amount of observation time is spent looking away from the source. Alternatively, the tracking mode, where the background is estimated from the on source run itself (see §3.2), can be used. In this case an off source run is not required for each on source run, allowing continuous monitoring of an object. However, off runs are still needed to determine the response of the telescope to background events.

Observations are typically made when the source zenith angle is less than 35 and are referred to as small zenith angle (SZA) observations. Observations at large zenith angles (LZA, typically 55 to 70) may also be made. Increasing the zenith angle has the effect of increasing the energy threshold, but has the benefit of increasing the collection area. Thus, it is an excellent method for increasing photon statistics to facilitate the determination of the energy spectrum at higher energies. For a detailed description of the LZA technique see Krennrich et al. (1997).

Since its discovery as a $`\gamma `$-ray source in 1995 (Quinn et al. (1996)), the VHE $`\gamma `$-ray emission from Mrk 501 has been monitored intensively with the Whipple Observatory 10 m telescope. Only SZA observations taken under good sky conditions have been considered in the analysis for variability presented here. Our selection includes data from 56, 50, 55 and 49 nights of observation in the Spring – Summer periods of 1995, 1996, 1997 and 1998 respectively. Table 1 summarises the resulting database.

## 3 Data Analysis

### 3.1 $`\gamma `$-ray Selection

The vast majority of events detected by Čerenkov telescopes are cosmic rays. Candidate $`\gamma `$-ray events are selected on the basis of the shape and orientation of the Čerenkov images: $`\gamma `$-ray images are typically more compact and elliptical than background hadronic images and tend to have their major axes aligned with the source location in the FOV. Background cosmic ray events, on the other hand, have random orientations.

Each image is first subjected to a cleaning procedure (Fegan (1997)) which suppresses pixels which are dominated by light from fluctuations of the night-sky background. A moment-fitting routine is then used to calculate various image parameters. The shape of each image is characterised by the parameters *length* and *width* and the orientation by $`\alpha `$, the angle between the major axis of the image and the line joining the source location in the FOV to the centroid of the image. In addition, for the data taken in 1998, an *asymmetry* cut has been included. $`\gamma `$-ray images have a cometary shape with a tail which points away from the source location in the FOV (Buckley et al. (1998)) and thus their intensity profiles have positive asymmetry. The larger FOV of the camera used in 1998 allows this parameter to be accurately determined, something which was not possible with the smaller FOV cameras.

Before the application of shape and orientation cuts a *software trigger* cut is also applied to eliminate events close to threshold, some of which are induced by noise fluctuations. The software trigger involves cuts on the image *size* (i.e. the total number of photoelectrons recorded), the counts in each of the brightest two tubes (*max1, max2*), and a requirement that at least three tubes above a low noise threshold (2.25$`\sigma `$, where $`\sigma `$ is the RMS sky-noise in a PMT, as determined from the artificially triggered events) be neighbours (*NBR3*). A *distance* cut is applied to eliminate images which are too close to the center of the camera and will have poor $`\alpha `$ reconstruction and also those events which have occurred too close to the edge of the FOV and may be truncated.

Due to the continuous evolution of the high resolution camera, and changes such as deterioration of mirror reflectivity through weathering and the presence or absence of light-cones, the optimum data analysis cuts differ for each year. The cuts used for a given telescope configuration are optimised on an independent data set, usually data taken on the Crab Nebula or, in the case of the 1998 data, Mrk 421. Table 2 lists the cuts used for each years’ analyses.

### 3.2 $`\gamma `$-ray Rate and Flux Calculation

There are slightly different analysis methods for the on/off and tracking observation modes. For the on/off observations the background is estimated from the off source run, which is assumed to be on a sky region which does not include a $`\gamma `$-ray source. This analysis mode has been discussed at length elsewhere (e.g. Kerrick et al. (1995); Catanese et al. (1998)).

For tracking observations the background is estimated from the on source run itself. All of the $`\gamma `$-ray selection criteria apart from orientation ($`\alpha `$) are applied to the data. The background is then estimated from events which are not oriented towards the source. In this analysis, background events with values of $`\alpha `$ between 20 and 65 are used. Images having values of $`\alpha `$ between 65 and 90 are discarded because of possible systematic effects due to truncation at the camera’s edge. Once the number of events with orientations in the 20 to 65 range is known then the expected number of background events in the *signal* domain ($`\alpha `$ of 0 to 15, or to 10 for data taken in 1998) can be estimated. Off source data recorded for this source and others can be combined to calculate a ratio, $`r\pm \mathrm{\Delta }r`$, of the number of events in the signal region to those in the 20 to 65 region in the absence of a source. In this case, the significance of a $`\gamma `$-ray excess, $`S`$, is given by:

$$S=\frac{N_{on}rN_{off}}{\sqrt{N_{on}+r^2N_{off}+(\mathrm{\Delta }r)^2N_{off}^2}}$$
(1)   
and the $`\gamma `$-ray rate ($`R\pm \mathrm{\Delta }R`$) is calculated from:

$$R\pm \mathrm{\Delta }R=\frac{N_{on}rN_{off}}{t}\pm \frac{\sqrt{N_{on}+r^2N_{off}+(\mathrm{\Delta }r)^2N_{off}^2}}{t}$$
(2)   
where $`N_{on}`$ is the number of counts in the $`\gamma `$-ray domain ($`\alpha <10^{}`$ or $`15^{}`$), $`N_{off}`$ is the number of counts in the 20 to 65 $`\alpha `$ range and $`t`$ is the duration of the observation.

The inclusion of the statistical error on the tracking ratio effectively limits the amount (duration) of tracking data which can be usefully analysed for an excess. Once the duration of the tracking data exceeds that of the off source data used in the calculation of the tracking ratio the error on the tracking ratio starts to dominate and to limit the significance of a detection.

For the purpose of investigating possible differences in results produced by the on/off and tracking analyses, the Crab Nebula data were analysed with both methods (Quinn (1997)). The results demonstrated that the $`\gamma `$-ray rates derived using both methods were consistent and stable; the tracking analysis did not introduce any apparent variability in the rate. For the analysis presented in this paper, all of the on source data were combined with the tracking data and the resulting database analysed using the tracking analysis.

The data presented here were taken with different telescope configurations having different sensitivities. It is therefore necessary to normalize when comparing the different data sets. To a first approximation, this can be achieved by converting the rates to fractions of the $`\gamma `$-ray rate from the Crab Nebula taken with the same telescope configuration and analysed with the same cuts. The Crab Nebula is believed to be a steady source of VHE $`\gamma `$-rays, as has been observed by the Whipple Collaboration over the past decade (Cawley et al. (1999)). To convert a given $`\gamma `$-ray rate to an integral flux, the rate as a fraction of the Crab Nebula flux was multiplied by $`(1.05\pm 0.24)\times 10^{10}`$cm<sup>-2</sup>s<sup>-1</sup>, which represents the integral Crab Nebula flux above 350 GeV (Hillas et al. (1998)), the threshold of the analysis presented here.

### 3.3 Test for Variability

To search for temporal variability in the $`\gamma `$-ray flux we apply a $`\chi ^2`$ test for a constant rate. The $`\chi ^2`$ sum is converted into a probability ($`P_{\chi ^2}`$) that the emission is constant about the mean using the incomplete gamma function *gammq*(a,x) (Press et al. (1988))

$$P_{\chi ^2}=gammq(\frac{N1}{2},\frac{\chi ^2}{2})$$
(3)   
where N-1 is the number of degrees of freedom. The number of trials is taken into account by calculating the probability ($`P_{trials}`$) of $`P_{\chi ^2}`$ occurring in N trials from:

$$P_{trials}=1(1P_{\chi ^2})^N.$$
(4)   
This method is used to test whether the distribution of measured $`\gamma `$-ray rates is consistent with statistical fluctuations about the mean for a range of time-scales.

## 4 Results

The data have been analysed using the tracking analysis described above. This differs from that applied to the 1995 data by Quinn et al. (1996), in that a statistical error on the tracking ratio is now included and that a 10% systematic error is no longer added to tracking results as a careful study showed that the results of the on/off and tracking analysis methods are in close agreement (Quinn (1997)). In fact, these changes tend to cancel each other.

For the 1995 data set as a whole we obtain a $`\gamma `$-ray rate of approximately 10% of that of the Crab Nebula ($`0.18\pm \mathrm{\hspace{0.17em}0.02}`$min<sup>-1</sup>, giving a 9.1$`\sigma `$ excess). This rose to approximately 20% of the Crab Nebula rate for the following season (analysis of 1996 data reveals an excess of 11.1$`\sigma `$ and a $`\gamma `$-ray rate of $`0.26\pm \mathrm{\hspace{0.17em}0.02}`$min<sup>-1</sup>), indicating that the average emission level had doubled since the previous year.

The monthly and nightly average rates, in fractions of the Crab Nebula rate, over all four years of Mrk 501 observation are shown in Figure 1. The rate in 1995 appears to have been constant with the exception of one night, MJD 49920, when the it was approximately 4.6$`\sigma `$ above the average. A $`\chi ^2`$ test gives a chance probability of $`1.2\times 10^3`$ (after accounting for trials) that the daily averages were constant during that month. The $`\chi ^2`$ probability that the daily averages are constant over the entire 5 months of observation is 0.17 (after trials). The probability that the emission is constant when averaged on monthly time-scales is 0.06 (after trials). Conversely, in the 1996 data set there were no obvious flaring episodes, but when averaged on time-scales of a month there is significant variability. The $`\chi ^2`$ probability that the emission is constant for the monthly averages is $`3.8\times \mathrm{\hspace{0.17em}10}^6`$, after accounting for trials. When each month is examined for variability with the rates averaged on time-scales of a day no significant variability is found. However, the $`\chi ^2`$ probabilities are smaller than for a similar analysis of the 1995 data, suggestive of increased day-scale flickering in the $`\gamma `$-ray emission. The results of the $`\chi ^2`$ test for variability of the monthly and daily averages are shown in Tables 3 and 4 respectively.

Observations of Mrk 501 in 1997 began in January using the LZA technique (Krennrich et al. (1998)). The initial observations suggested that the flux was much higher than in previous years. Conventional SZA observations commenced as soon as possible in February. Other ground-based $`\gamma `$-ray experiments were notified and high emission levels were verified. In March a joint IAU circular by the CAT, HEGRA and Whipple groups (Breslin et al. (1997)) announced preliminary results. The Whipple observations continued through June and included observations made with the LZA technique. The results of the LZA data are presented elsewhere (Krennrich et al. (1998)).

Our SZA observations revealed that the VHE $`\gamma `$-ray emission was very strong throughout Spring – Summer 1997. Significant variability was observed in the monthly averages. The emission appeared to increase steadily from February through May and then to level off in June. The daily rates also exhibit dramatic variability. Day-to-day changes in the flux by factors $`>`$4 were observed and on eight occasions the flux more than doubled between consecutive nights. On four occasions there were equally rapid decays in the rate. The average flux level for the season is 1.4 times that of the Crab Nebula - an increase by a factor of 14 from the level in 1995. The peak rate, observed on MJD 50554, is 3.7 times the Crab rate. The average rates for the 1997 data were calculated from only one run (the first with elevation above 55) per night. This was done to remove a bias in the observing strategy whereby observations of Mrk 501 on a given night continued only if the source was in a very active state.

The strong signal-to-background level in the 1997 data allowed a search to be made for variability on time-scales shorter than one day. For this test, data from each night on which there were three or more runs (approx. 1.5 hours, see §2) were analysed to test for run-to-run variability. A total of 24 nights satisfied this criterion. The $`\chi ^2`$ probability that the emission was constant was calculated for each of the nights; the distribution of these probabilities is shown in Figure 2a. For a statistically-variable source this distribution should be flat but there is an excess of nights with small probabilities. Assuming that the variations are purely statistical, the probability of getting seven nights in this first bin (width = 0.025) out of 24 trials is 1.45$`\times `$10<sup>-6</sup>. Of these seven nights, two show statistically significant variations within themselves. The probability for constant emission on MJD 50577 is $`5.2\times \mathrm{\hspace{0.17em}10}^6`$ while for MJD 50607 the probability is $`5.8\times \mathrm{\hspace{0.17em}10}^8`$ (after accounting for trials). The flux on MJD 50607 has a doubling time of $``$2 hours. We thus identify these two nights as exhibiting significant hour-scale variability while the five other nights exhibit marginal variability. Figure 3 shows the $`\gamma `$-ray rates for the two nights with significant variability.

A search of the 1997 data for variability on time-scales of less than half an hour has also been performed. For this test each of the 28 minute runs, 144 in total, was divided into three equal length intervals. Each triplet was then analysed for variability. No significant variations were found and the distribution of $`\chi ^2`$ probabilities (Figure 2b) does not indicate any excess of low probabilities. Hence, we see no evidence for significant sub-hour scale variability.

The average $`\gamma `$-ray rate for the 1998 data set was 0.42$`\pm `$0.04 min<sup>-1</sup>, approximately 20% of the rate obtained from the Crab Nebula i.e. on average the emission was much lower than in 1997. In fact, the average rates for March, April and May are comparable to the level of the initial detection in 1995. There were however two significant flaring events. The first occurred in early March where an apparent rise and decay were observed. Unfortunately this flare is poorly sampled due to bad weather. The flux was observed to be relatively high ($``$1.3 times the Crab Nebula flux) on MJD 50876 and on the following night a flux of $``$5.0 times that of the Crab Nebula was recorded. This is the largest flux detected to date from Mrk 501 by the Whipple Observatory 10 m telescope. For the next observation on MJD 50880 the measured flux was still relatively high ($``$1.3 times the level of the Crab Nebula flux). A second flare occurred in June. The flux increased on two consecutive nights, peaking at $``$1.1 times the Crab Nebula flux on MJD 50991, then decayed on a similar time-scale. The average flux level immediately to either side of this flare was below the sensitivity of the telescope.

## 5 Discussion

We have demonstrated that rapid variability, a common characteristic of blazars at all observed energies, is also present in the VHE $`\gamma `$-ray emission from Mrk 501. Our 4 year data set spans a remarkable change in the flux level from Mrk 501: the average yearly emission level exhibited a fourteen-fold increase between 1995 and 1997 and the average daily flux varied by a factor of $``$50 (see Figure 1). In 1997 large amplitude day-scale flares occurred frequently and were usually followed by equally rapid decays. Day-scale changes in flux by factors as large as 4.7 were observed. Episodes of significant hour-scale variability were detected, with one having a doubling time of 2 hours. In addition, there is evidence of consistent hour-scale variability which is not resolved in individual episodes. A variable flux from Mrk 501 was reported by at least four other atmospheric Čerenkov observatories in 1997 (see, e.g. Protheroe et al. (1998)).

The data presented here suggest that Mrk 501 was more variable when the flux level was higher. However, this effect could also be due to the sensitivity of the telescope. At low flux levels it takes longer to accumulate a significant signal, so the search for short-term variability may be limited by poor statistics. To address this issue we performed a test to see if the day-scale variability observed in 1997 would have been detected in 1996 and/or 1995 and if the month-scale variability seen in 1996 and 1997 would have been detected in 1995. We calculated the percentage deviations about the mean level from a period where significant variability was observed and then, using the mean signal and background level from another period, calculated the signal (and statistical error) which would have been observed given the same percentage deviations about that mean. The results are that neither the 1997 month-scale nor day-scale variability would have been detected at a significant level if present in 1995. However, the 1996 month-scale variability would have been significant if present in 1995 (chance probability of $`10^6`$), while the 1997 degree of day-scale variability, if present in 1996, would have been detectable (chance probability $`10^7`$). Thus we conclude that there was a change in the flaring characteristics, in addition to the change in the mean flux level between the different observing seasons.

The VHE $`\gamma `$-ray emission from Mrk 501 exhibits rapid variability similar to that of the emission from Mrk 421. A major difference between the two objects is that the emission from Mrk 501 seems to have a base-level, which changes on monthly and yearly time-scales, whearas the VHE $`\gamma `$-ray emission from Mrk 421 has been described as consisting of a series of rapid flares with no underlying baseline (Buckley et al. (1996)). The variability of the VHE $`\gamma `$-ray emission of Mrk 501 in 1995 and 1996 was similar to that at other wavelengths, with small amplitude, slow variations being more common than fast, large amplitude flares.

The increase in the VHE $`\gamma `$-ray power and variability in 1997 was accompanied by an increase in the hard X-ray power and an extension of the synchrotron spectrum to at least 100 keV (Catanese et al. (1997)). This is consistent with an inverse-Compton mechanism for producing VHE $`\gamma `$-rays. However, our observations of Mrk 501 to date cannot discriminate between the inverse-Compton (electron) models and those where the dominant particles producing $`\gamma `$-rays in the jet are protons. More densely sampled light curves, covering as broad a wave-band as possible, are needed to provide more insight into the mechanisms responsible for the VHE radiation from BL Lac objects. Future, more sensitive $`\gamma `$-ray observations of AGN with proposed detectors such as GLAST, HESS, MAGIC and VERITAS will allow the structure of flares on shorter time-scales to be determined.

We acknowledge the technical assistance of K. Harris and E. Roache. This research is supported by grants from the US Department of Energy and by NASA, by PPARC in the UK and by Forbairt in Ireland.