A&A 457, 899-915 (2006)
DOI: 10.1051/0004-6361:20065351

Observations of the Crab nebula with HESS

F. Aharonian1 - A. G. Akhperjanian2 - A. R. Bazer-Bachi3 - M. Beilicke4 - W. Benbow1 - D. Berge1 - K. Bernlöhr1,5 - C. Boisson6 - O. Bolz1 - V. Borrel3 - I. Braun1 - F. Breitling5 - A. M. Brown7 - R. Bühler1 - I. Büsching8 - S. Carrigan1 - P. M. Chadwick7 - L.-M. Chounet9 - R. Cornils4 - L. Costamante1,21 - B. Degrange9 - H. J. Dickinson7 - A. Djannati-Ataï10 - L. O'C. Drury11 - G. Dubus9 - K. Egberts1 - D. Emmanoulopoulos12 - P. Espigat10 - F. Feinstein13 - E. Ferrero12 - A. Fiasson13 - G. Fontaine9 - Seb. Funk5 - S. Funk1 - Y. A. Gallant13 - B. Giebels9 - J. F. Glicenstein14 - P. Goret14 - C. Hadjichristidis7 - D. Hauser1 - M. Hauser12 - G. Heinzelmann4 - G. Henri15 - G. Hermann1 - J. A. Hinton1,12 - W. Hofmann1 - M. Holleran8 - D. Horns16 - A. Jacholkowska13 - O. C. de Jager8 - B. Khélifi9,1 - Nu. Komin13 - A. Konopelko5 - K. Kosack1 - I. J. Latham7 - R. Le Gallou7 - A. Lemière10 - M. Lemoine-Goumard9 - T. Lohse5 - J. M. Martin6 - O. Martineau-Huynh17 - A. Marcowith3 - C. Masterson1,21 - T. J. L. McComb7 - M. de Naurois17 - D. Nedbal18 - S. J. Nolan7 - A. Noutsos7 - K. J. Orford7 - J. L. Osborne7 - M. Ouchrif17,21 - M. Panter1 - G. Pelletier15 - S. Pita10 - G. Pühlhofer12 - M. Punch10 - B. C. Raubenheimer8 - M. Raue4 - S. M. Rayner7 - A. Reimer19 - O. Reimer19 - J. Ripken4 - L. Rob18 - L. Rolland14 - G. Rowell1 - V. Sahakian2 - L. Saugé15 - S. Schlenker5 - R. Schlickeiser19 - U. Schwanke5 - H. Sol6 - D. Spangler7 - F. Spanier19 - R. Steenkamp20 - C. Stegmann5 - G. Superina9 - J.-P. Tavernet17 - R. Terrier10 - C. G. Théoret10 - M. Tluczykont9,21 - C. van Eldik1 - G. Vasileiadis13 - C. Venter8 - P. Vincent17 - H. J. Völk1 - S. J. Wagner12 - M. Ward7


1 - Max-Planck-Institut für Kernphysik, PO Box 103980, 69029 Heidelberg, Germany
2 - Yerevan Physics Institute, 2 Alikhanian Brothers St., 375036 Yerevan, Armenia
3 - Centre d'Étude Spatiale des Rayonnements, CNRS/UPS, 9 Av. du Colonel Roche, BP 4346, 31029 Toulouse Cedex 4, France
4 - Universität Hamburg, Institut für Experimentalphysik, Luruper Chaussee 149, 22761 Hamburg, Germany
5 - Institut für Physik, Humboldt-Universität zu Berlin, Newtonstr. 15, 12489 Berlin, Germany
6 - LUTH, UMR 8102 du CNRS, Observatoire de Paris, Section de Meudon, 92195 Meudon Cedex, France
7 - University of Durham, Department of Physics, South Road, Durham DH1 3LE, UK
8 - Unit for Space Physics, North-West University, Potchefstroom 2520, South Africa
9 - Laboratoire Leprince-Ringuet, IN2P3/CNRS, École Polytechnique, 91128 Palaiseau, France
10 - APC, 11 place Marcelin Berthelot, 75231 Paris Cedex 05, France [*]
11 - Dublin Institute for Advanced Studies, 5 Merrion Square, Dublin 2, Ireland
12 - Landessternwarte, Universität Heidelberg, Königstuhl, 69117 Heidelberg, Germany
13 - Laboratoire de Physique Théorique et Astroparticules, IN2P3/CNRS, Université Montpellier II, CC 70, Place Eugène Bataillon, 34095 Montpellier Cedex 5, France
14 - DAPNIA/DSM/CEA, CE Saclay, 91191 Gif-sur-Yvette, Cedex, France
15 - Laboratoire d'Astrophysique de Grenoble, INSU/CNRS, Université Joseph Fourier, BP 53, 38041 Grenoble Cedex 9, France
16 - Institut für Astronomie und Astrophysik, Universität Tübingen, Sand 1, 72076 Tübingen, Germany
17 - Laboratoire de Physique Nucléaire et de Hautes Énergies, IN2P3/CNRS, Universités Paris VI & VII, 4 place Jussieu, 75252 Paris Cedex 5, France
18 - Institute of Particle and Nuclear Physics, Charles University, V Holesovickach 2, 180 00 Prague 8, Czech Republic
19 - Institut für Theoretische Physik, Lehrstuhl IV: Weltraum und Astrophysik, Ruhr-Universität Bochum, 44780 Bochum, Germany
20 - University of Namibia, Private Bag 13301, Windhoek, Namibia
21 - European Associated Laboratory for Gamma-Ray Astronomy, jointly supported by CNRS and MPG

Received 3 April 2006 / Accepted 5 July 2006

Abstract
Context. The Crab nebula was observed with the HESS stereoscopic Cherenkov-telescope array between October 2003 and January 2005 for a total of 22.9 h (after data quality selection). This period of time partly overlapped with the commissioning phase of the experiment; observations were made with three operational telescopes in late 2003 and with the complete 4 telescope array in January-February 2004 and October 2004-January 2005.
Aims. Observations of the Crab nebula are discussed and used as an example to detail the flux and spectral analysis procedures of HESS. The results are used to evaluate the systematic uncertainties in HESS flux measurements.
Methods. The Crab nebula data are analysed using standard HESS analysis procedures, which are described in detail. The flux and spectrum of \ensuremath {\gamma }-rays from the source are calculated on run-by-run and monthly time-scales, and a correction is applied for long-term variations in the detector sensitivity. Comparisons of the measured flux and spectrum over the observation period, along with the results from a number of different analysis procedures are used to estimate systematic uncertainties in the measurements.
Results. The data, taken at a range of zenith angles between 45 \ensuremath {^{\circ }} and 65 \ensuremath {^{\circ }}, show a clear signal with over 7500 excess events. The energy spectrum is found to follow a power law with an exponential cutoff, with photon index $\Gamma = 2.39$ $\pm$ $0.03\ensuremath{_{{\rm\tiny {stat}}}} $ and cutoff energy $E_{\rm c} = (14.3 \pm 2.1\ensuremath{_{{\rm\tiny {stat}}}} )~{\rm TeV}$ between 440 GeV and 40 TeV. The observed integral flux above 1 TeV is $(2.26 \pm 0.08\ensuremath{_{{\rm\tiny {stat}}}} ) \times 10^{-11}~\ensuremath{{\rm cm}^{-2}~{\rm s}^{-1}} $. The estimated systematic error on the flux measurement is estimated to be 20%, while the estimated systematic error on the spectral slope is 0.1.

Key words: gamma rays: observations - ISM: individual objects: Crab nebula - ISM: supernova remnants

1 Introduction

The Crab supernova remnant (SNR) is an exceptionally well studied object, with extensive observations of the system existing across the entire accessible spectrum. At a distance of 2000 parsecs, with an age of 950 years, it is a prototypical centre-filled SNR, or plerion, as defined by Weiler & Panagia (1980). Within the supernova remnant lies the Crab pulsar, with a rotational period of 33 ms and a spin-down luminosity of L = 5 $\times $ $10^{38}~{\rm erg}~{\rm s}^{-1}$. This energy source powers a surrounding synchrotron nebula, and polarization measurements exist from radio to hard X-ray wavelengths (Wilson 1972), indicating the non-thermal origin of the radiation detected. The total energy available from the pulsar to power the system is of the order of 1049 erg. This is believed to be the power source for production of very high energy (VHE)  \ensuremath {\gamma }-rays.

The rotational energy of the pulsar is thought to be mostly carried away by a relativistic wind of electrons and positrons. Interaction of this wind with the surrounding medium causes a standing termination shock wave (Kennel & Coroniti 1984; Rees & Gunn 1974). Electron acceleration may be due to a Fermi-type process (Achterberg et al. 2001) or to driven reconnection of the alternating magnetic field at this termination shock (Michel 1994; Coroniti 1990). The interaction of accelerated electrons with ambient photon fields (in this case mostly synchrotron photons) can produce VHE  \ensuremath {\gamma }-rays via the inverse Compton process.

The Crab nebula was discovered at VHE energies in 1989 (Weekes et al. 1989) and emission has been confirmed by a number of other experiments, (Masterson et al. 1999; Goret et al. 1993; Baillon et al. 1993; Atkins et al. 2003; Aharonian et al. 2000; Smith et al. 2000). Due the high flux from the source relative to other known TeV sources, and its expected flux stability, it is conventionally used as a standard reference source for VHE astronomy. At the latitude of the HESS experiment the Crab nebula culminates at 45 \ensuremath {^{\circ }}, so observations of this source must always be made at large zenith angles; the greater effective optical depth of the atmosphere increases the energy threshold and affects the sensitivity of the detector. As the size of the Crab nebula is small compared to the HESS point-spread function, it may be treated as a point source for this analysis.

A detailed description of the HESS detector is given here, along with a discussion of the principal sources of systematic error in the atmospheric-Cherenkov technique. The HESS observations of the Crab nebula are then detailed, and used as an example in a discussion of the data calibration and analysis methods used in source reconstruction and flux and spectral measurements. The stability of the \ensuremath {\gamma }-ray flux and energy spectrum is measured using a number of analysis methods, and a correction for variations in the long-term optical efficiency of the detector is described. Using these results the systematic uncertainties on flux measurements with HESS are quantified. The sensitivity of the detector for source analysis is also discussed. It should be noted that while the analysis presented here is generally used to analyse targets observed by HESS, other techniques are also used, e.g. Lemoine-Goumard et al. (2006).

2 The HESS experiment

HESS is situated in the Khomas highlands of Namibia (23 \ensuremath {^{\circ }}16'18'' South, 16 \ensuremath {^{\circ }}30'00'' East), at an elevation of 1800 m above sea level. The four HESS telescopes are placed in a square formation with a side length of 120 metres. This distance was optimised for maximum sensitivity at the planned energy threshold of 100 GeV.

2.1 The detectors

The HESS telescopes are of steel construction, with altitude/azimuth mounts capable of precisely tracking any source from 0.0 \ensuremath {^{\circ }} to 89.9 \ensuremath {^{\circ }} in elevation, with a slew rate of 100 \ensuremath {^{\circ }} per minute (Bolz 2004). The dishes have a Davies-Cotton style hexagonal arrangement (Davies & Cotton 1957) with a flat-to-flat diameter of 13 m, composed of 382 round mirrors, each 60 cm in diameter. Thus the effective mirror surface area is $107~{\rm m}^2$. Further details of the optical structure are given by Bernlöhr et al. (2003). The mirrors are remotely adjustable under computer control, and a fully automated procedure is used, in conjunction with a CCD camera mounted in each dish, for optimal alignment onto the focal plane of each telescope camera, which is 15 m distant. Due to the rigidity of the dishes, this alignment is stable over time scales in excess of one year. The stability has been verified by regular monitoring of the optical point spread function. Details of the mirror alignment system and the optical point spread function are discussed by Cornils et al. (2003).

The HESS cameras each consist of a hexagonal array of 960 Photonis XP2960 photo-multiplier tubes (PMTs). Each tube corresponds to an area of 0.16 \ensuremath {^{\circ }} in diameter on the sky, and is equipped with Winston cones to capture the light which would fall in between the PMTs, and also to limit the field of view of each PMT in order to minimise background light. The camera is of modular design, with the PMTs grouped in 60 drawers of 16 tubes each (Vincent et al. 2003), which contain the trigger and readout electronics for the tubes, as well as the high voltage (HV) supply, control and monitoring electronics. The total field of view of the detector is 5 \ensuremath {^{\circ }} in diameter.

The trigger system of the HESS array consists of three levels. First, a single pixel trigger threshold is required, equivalent to 4 photo-electrons (p.e.) at the PMT cathode within an interval of 1.5 nanoseconds. Second, a coincidence of 3 triggered pixels is required within a sector - a square group of 64 pixels - in order to trigger a camera. Each camera has 38 overlapping sectors. Third, when the detector is operating in stereo mode, a coincidence of two telescopes triggering within a window of (normally) 80 nanoseconds is required. Only cameras which have individually triggered are read out in a stereo event. The stereo trigger system and the trigger behavior of the HESS array is described by Funk et al. (2004).

Table 1: Details of the observations of the Crab nebula with HESS between October 2003 and January 2005. For the purposes of this study, the data has been divided into 3 subsets.

During the first and second level trigger formation, the individual signals from each pixel are stored in two analogue ring sampler (ARS) circuits. A high gain and a low gain circuit are used to give optimal signal dynamic range. The signals captured by each tube are digitised in the drawer, before being collected by a central CPU in the camera and sent to the central data acquisition system (DAQ) by optical ethernet connection (Borgmeier et al. 2003).


  \begin{figure} \par\subfigure[]{\includegraphics[width=7.4cm,clip]{5351f1a.eps} ... ...8mm} \subfigure[]{\includegraphics[width=7.4cm,clip]{5351f1b.eps} } \end{figure} Figure 1: The system rate vs. time within a run for a) a 4 telescope run passing the run selection, b) a 3 telescope run failing the run selection. The run is rejected due to the instability in the rate caused by clouds passing through the field of view. The difference in absolute rate between the two runs is due to the differing zenith angles of the observations, as well as the number of telescopes active in each run.
Open with DEXTER

The HESS experiment commenced observations in June 2002 with the first telescope. The second telescope was installed in February 2003, the hardware level stereo coincidence trigger was added in July 2003. In the interim the two telescopes observed in parallel and an off-line coincidence was used for stereo analysis. The third telescope was installed in September. The full array was completed in December of 2003 and has been operational since then.

   
2.2 Systematic uncertainties

The imaging atmospheric-Cherenkov technique depends on a form of electromagnetic calorimetry to estimate the energy of observed particles. In order to accurately measure the energy of the primary particle which gives rise to an air shower, it is necessary to understand the relationship between the particle energy and the signal recorded in the cameras. Monte Carlo simulations of air showers in the atmosphere are used to predict the light yield in the detector, and thus the recorded signal, as a function of energy and shower position relative to the observer.

There are three main contributions to uncertainties in the measured air shower information, and thus the absolute flux calibration of the detector:

1.
The camera response. The single photo-electron response of each PMT varies strongly with the detector voltage and is measured using an LED system mounted in front of the camera. A second LED system mounted on each dish provides a uniform illumination across the camera and is used to correct for relative quantum efficiency variations of the PMTs, as well as in the reflectivity of the Winston cones in front of each PMT. The complete calibration of the HESS telescope is described by Aharonian et al. (2004a).

2.
The optical response of the instrument, including the mirrors, Winston cones, shadowing by structural components of the system and the quantum efficiency of the photo-cathodes of the PMTs. This response can be measured by studying the Cherenkov light from single muons passing close to the telescope, assuming the camera response is well measured. The use of muons in monitoring the telescope efficiencies is detailed by Bolz (2004) and Leroy et al. (2004). The optical response of the instrument degrades over a timescale of years as, for example, the mirror reflectivity decreases. This decrease in optical response, relative to that used in the Monte Carlo simulations, is taken into account in estimating the flux from a  \ensuremath {\gamma }-ray source. This is described in detail in Sect. 3.4.

3.
The interactions of particles and light in the atmosphere. The atmosphere is the largest and least well understood component of a Cherenkov detector, being subject to variations in pressure, temperature and humidity. Two important effects of variability in the atmosphere are density profile variations, which affect directly the height of the shower maximum in the atmosphere and thus the intensity of the light seen at the telescope, and absorption of Cherenkov light in the atmosphere by clouds and dust, which leads directly to a reduction in the telescope trigger rates, and incorrect \ensuremath {\gamma }-ray energy reconstruction. These effects are discussed further by Bernlöhr (2000). Atmospheric monitoring devices are used to understand the local conditions under which the data has been recorded; such measurements are discussed by Aye et al. (2005). Variations in the atmosphere can lead to rapid variability in the detector response, on a timescale of hours, so runs taken under variable conditions are rejected, as discussed in Sect. 3.2.

   
3 Observations

3.1 The Crab data set

The observations of the Crab nebula discussed in this paper are summarised in Table 1. Data set I was taken between October 2003 and January 2004 with three telescopes, and 5.3 h of data are included in this sample. These data range in zenith angles from 45 \ensuremath {^{\circ }} to 55 \ensuremath {^{\circ }}. A further 6.3 h (data set II) are included from data taken in January 2004 with the complete array of four telescopes, at a range of zenith angles from 45 \ensuremath {^{\circ }} to 65 \ensuremath {^{\circ }}. Data set III includes 11.3 h, taken between October 2004 and January of 2005, also with 4 telescopes. All observations were taken using a method (wobble mode) whereby the source is alternately offset by a small distance within the field of view, alternating between 28 min runs in the positive and negative declination (or right ascension) directions (Fomin et al. 1994). This observation mode allows the other side of the field of view, which does not contain the source, to be used as a control region for estimation of the background level. The wobble offsets were varied from 0.5 \ensuremath {^{\circ }} to 1.5 \ensuremath {^{\circ }} for the data discussed here. The range of time over which the observations have been taken allows us to study the long-term stability of the HESS system.

Observations of the Crab were also taken in 2002 with a single telescope, as discussed by Masterson et al. (2003). These are not included in this analysis, as the sensitivity is much lower than that in stereo mode, and the systematic uncertainties of single-telescope observations are greater than for stereo mode.

   
3.2 Data quality selection

Systematic effects on the measured flux and energy spectrum (as discussed in Sect. 2.2) may be ameliorated by rejecting data recorded in non-optimal conditions. Although the Crab was observed for a total of nearly 45 h, only 22.9 h are included in this analysis. The remaining observations have been rejected as not meeting the run quality criteria. Some were made when the sky conditions were less than optimal; the presence of clouds or excessive dust in the atmosphere can lead to the absorption of Cherenkov light and thus fluctuations in the system trigger efficiency, causing systematic uncertainties in the measured \ensuremath {\gamma }-ray flux. Figure 1 shows the trigger rate as a function of time for two runs, one with a stable system trigger rate close to the predicted level for this zenith angle, the other exhibiting variability due to the presence of clouds.

Runs for which the mean trigger rate is less than 70% of the predicted value (as discussed by Funk et al. (2004)), or for which the rms variation in the trigger rate is above 10%, are rejected. As the Crab is mainly visible during the rainy season in Namibia when weather conditions are not optimal, the rejection rate for these data is not representative of all HESS observations. The mean system rate is 240 Hz for the four telescope data and 180 Hz for the 3 telescope data.

Quality checks are also routinely carried out in order to reject runs in which the array tracking system may not be functioning correctly, leading to errors in the reconstructed position of the source, and thus the flux. The nominal performance of the HESS tracking system is discussed by Hofmann et al. (2003). The tracking errors reported from the DAQ are monitored and runs with rms deviations of more than 10 arcsec in altitude or azimuth are rejected, in order to exclude runs in which the tracking system malfunctioned. However, no observations on the Crab in this study (and in general very few runs) have failed this test. As an independent check of the tracking quality the DC PMT currents (Aharonian et al. 2004a) are used to estimate the amount of light impingent on each pixel during every run as a function of time. A map is then made of the sky brightness in the field of view of each telescope. The positions of known stars are then correlated with this map, giving a measure of the pointing position of each telescope, independent of the tracking system and standard pointing corrections. A sample image for the region surrounding the Crab nebula is shown in Fig. 2. Runs are rejected if the pointing deviation is greater than 0.1 \ensuremath {^{\circ }}. Again no runs are found to have failed this test, which acts as an auxiliary check of the tracking system.


  \begin{figure} \par\includegraphics[width=7.8cm,clip]{5351f2.eps} \end{figure} Figure 2: Sky map of DC pixel current (in units of $\mu $A, measured in a window of $16~\mu$s) for the region surrounding the Crab nebula. The positions of stars (read from the Hipparcos catalogue - Perryman & ESA 1997) match well with the peaks in the map, which correspond to peaks in sky brightness. The B and V magnitudes of each star are also given for comparison. The brightest star, $\zeta $ Tauri (B magnitude 2.8), has caused the camera to disable PMTs, thus there is no current measure for the region close to this star.
Open with DEXTER

The presence of bright stars in the field of view may trigger the over-current protection of individual PMTs, causing them to be turned off. Over-current due to bright star images is predicted and the relevant PMTs are turned off in advance, for the duration of the star transit through the PMT. Occasionally other bright, transient light normally meteorites, however lightning, airplanes and satellites may also cause problems. The events cannot be predicted and thus trigger the over-current protection, causing the PMT to be turned off for the remainder of the run. The mean number of PMTs inactive during each run for this and other hardware related reasons is monitored, and telescopes with more than 10% of the PMTs missing at any one time are rejected from the analysis. The PMT quality cut only rejects occasional runs from each data set, and the mean number of pixels per telescope disabled during the Crab nebula observations is 66, with an rms of 7.

The observation time after run selection is corrected by taking into account the dead-time of the system, when the trigger is not sensitive to air showers. Estimation of the live-time (the time during an observation in which the telescope system is sensitive to triggers from the sky) is discussed by Funk et al. (2004). The total live-time for each data set is listed in Table 1. The systematic error on this estimation is less than 1%.

3.3 Photoelectron to pixel amplitude calibration


  \begin{figure} \par\includegraphics[width=7.45cm,clip]{5351f3.eps} \end{figure} Figure 3: Distribution of the intensity of the 3rd highest pixel for 4 telescopes (points), compared with Monte Carlo simulations (solid lines). The overall flux of the simulated data has been adjusted to fit the real data.
Open with DEXTER

All of the sources of uncertainty discussed in Sect. 2.2 affect the pixel intensities, and thus the energy reconstruction of each shower. The trigger efficiency and energy threshold are also affected, contributing to uncertainty in absolute flux measurements. Part of this effect, namely the camera response, is continuously calibrated using the LED system; the accuracy of this calibration system can be evaluated by comparing the reconstructed pixel intensities from Monte Carlo simulations of air showers with similar measurements on real data. Figure 3 compares the intensities in the pixel with the third highest signal amplitude, as seen in real data for all triggering events, with those from Monte Carlo simulations of protons (which comprise a large majority of events detected) at the same zenith angle. As the trigger requirement is 3 pixels above a certain amplitude, the amplitude distribution of the 3rd pixel shows the effective pixel trigger level of each telescope. The peak in the distribution is determined by the trigger threshold of the telescope and the energy spectrum of the simulated and real events. It can be seen that the Monte Carlo simulations reproduce well the data ( \ensuremath {\chi ^{2}}/d.o.f. = 15.5/23).

   
3.4 Optical efficiency correction

The optical efficiencies of the HESS telescopes change over time, probable causes include degradation in the mirror reflectivity. This effect happens on a timescale of years and can be monitored using images of local muons in the field of view, for which the light yield can be predicted, as discussed in Sect. 2.2.

The effect manifests itself as a reduction in the image intensity for each event, compared to the intensity expected from Monte Carlo simulations. This causes a shift in the absolute energy scale of the detector, as events are reconstructed with energies which are too low. This effect is corrected by incorporating a scaling factor into the energy estimation for each event. The image intensity used in the energy estimation is scaled by the ratio of the mean optical efficiency (over the telescopes) for the run ( ${\rm Eff}_{\rm run}$) to the mean optical efficiency as derived from the Monte Carlo simulations ( ${\rm Eff}_{\rm mc}$). The corrected energy is then derived from this scaled image intensity in the standard manner. The distribution of the relative optical efficiency ( $\frac{{\rm Eff}_{\rm run}}{{\rm Eff}_{\rm mc}}$) for the runs included in this analysis is shown in Fig. 4. The application of the optical efficiency correction in flux estimation is discussed further in Sect. 6.3.


  \begin{figure} \par\includegraphics[width=6cm,clip]{5351f4.eps} \end{figure} Figure 4: Distribution of the relative optical efficiency per run for data sets I (solid line), II (dashed line) and III (dotted line). It can be seen that the relative optical efficiency is significantly decreased in data set III relative to the other data sets.
Open with DEXTER

4 Analysis

After a set of images of an air shower has been recorded, they are processed to measure Hillas parameters based on the second moments of the image (Hillas 1985). These parameters are then used for event selection and reconstruction. A diagram illustrating the parameter definitions is shown in Fig. 5.

4.1 Image cleaning and moment analysis

The first step in the moment analysis procedure is image cleaning. This is required in order to select only the pixels containing Cherenkov light in an image. Other pixels, which contain mainly night sky background (NSB) light are not used in the analysis. Images are cleaned using a two-level filter, requiring pixels in the image to be above a lower threshold of 5 p.e. and to have a neighbour above 10 p.e., and vice versa. Cleaning thresholds of 4 p.e. and 7 p.e. have also been shown to work satisfactorily, but may be more sensitive to uncertainties due to NSB light variations. This method selects spatially correlated features in the image, which correspond to air shower Cherenkov light. This method tends to smooth out shower fluctuations in a simple and repeatable manner.

After image cleaning, an image of a  \ensuremath {\gamma }-ray shower approximates a narrow elliptical shape, while images of background hadronic events are wider and more uneven. The Hillas parameters are then calculated for each cleaned image; these parameters are the basis for event selection. The total amplitude of the image after cleaning is also calculated, along with the mean position of the image in the camera, which corresponds to the centroid of the ellipse.


  \begin{figure} \par\includegraphics[width=7.9cm,clip]{5351f5.eps} \end{figure} Figure 5: Definition of simple Hillas parameters, calculated for a  \ensuremath {\gamma }-ray image, which may be approximated as an ellipse. Important parameters for this analysis are the width, length, distance. An image from a second telescope is superimposed to demonstrate the geometrical technique for source position reconstruction. The parameter $\theta $, which is the magnitude of the angular offset in shower direction reconstruction, is also shown.
Open with DEXTER

4.2 Stereo reconstruction

The arrival direction of each event is reconstructed by tracing the projected direction of the shower in the field of view (which corresponds to the major axis of the image) to the point of origin of the particle. For stereo observations it is possible to intersect the major axes of the shower images in multiple cameras, providing a simple geometric method of accurately measuring the shower direction; more details, including methods to further improve the reconstruction accuracy are given by Hofmann et al. (1999), method I from that paper is used here. Images are only used in the stereo reconstruction if they pass the selection cuts on distance (to avoid camera-edge effects) and image intensity. If less than two telescope images pass these cuts the event is rejected.


  \begin{figure} \par\includegraphics[width=6.5cm,clip]{5351f6.eps} \end{figure} Figure 6: Distribution of excess events in \ensuremath {\theta ^{2}} for the complete Crab data set, after event selection and background subtraction. The Monte Carlo derived point-spread function described in Eq. (1) is also shown, normalised to the excess distribution. The vertical lines denote the \ensuremath {\theta ^{2}} selection cuts listed in Table 2.
Open with DEXTER

Table 2: Selection cuts optimised for various purposes, as described in the text. Cuts are applied on MRSW and MRSL, as well as on the distance ($\theta $) from the reconstructed shower position to the source. A minimum of two telescopes passing the per-telescope cuts, on image amplitude and distance from the centre of the field of view, are also required. Standard cuts, as well as hard, loose and extended cuts, as described in the text, are listed.

Figure 6 shows the excess distribution of \ensuremath {\theta ^{2}} for data sets I-III, including events with two, three and four telescopes; $\theta $ is defined in Fig. 5, it is the angular offset between the reconstructed shower direction and the true direction of the Crab nebula. The distribution of reconstructed shower directions is usually expressed in units of  \ensuremath {\theta ^{2}}, as this ensures a constant solid angle on the sky per bin. The value of the cut on reconstructed shower direction is thus given in units of  ${\rm degrees}^{2}$ in Table 2 for various sets of selection cuts, and plotted in Fig. 6. A strong excess is seen close to zero, corresponding to events coming from the direction of the Crab nebula. This distribution defines the accuracy in the reconstructed arrival directions for \ensuremath {\gamma }-ray events from a point source and is described by the point spread function (PSF). This function can be approximated by the sum of two, one-dimensional Gaussian functions:

\begin{displaymath}% {\rm PSF} = {A} \left(\exp \left(\frac{-\theta^{2}}{2 \sig... ...eft( \frac{-\theta^{2}}{2 \sigma_{2}^{2}} \right) \right)\cdot \end{displaymath} (1)

This function is fitted to the \ensuremath {\theta ^{2}} distribution for simulated Monte Carlo \ensuremath {\gamma }-rays. For simulations at 50 \ensuremath {^{\circ }} zenith angle the standard deviation parameters  $\sigma_{1}$ and  $\sigma_{2}$ are 0.046 \ensuremath {^{\circ }} and 0.12 \ensuremath {^{\circ }} respectively. The relative amplitude of the second Gaussian, $A_{\rm rel}$ is 0.15, while the absolute amplitude (A) is proportional to the number of events in the fit. The 68% containment radius for 50 \ensuremath {^{\circ }} zenith angle is 0.12 \ensuremath {^{\circ }}, while that for 10 \ensuremath {^{\circ }} zenith angle is 0.10 \ensuremath {^{\circ }}. The point-spread function is shown in Fig. 6 with the amplitude parameter (A) fit to the Crab nebula data. The \ensuremath {\chi ^{2}}/d.o.f. of this fit is 53/29.

The position of the centre of the Cherenkov light pool, which corresponds to the projected impact point of the original particle track on the ground, can also be reconstructed by intersecting shower axes, projected into the plane perpendicular to the system observing direction. It is vital to reconstruct this position in order to accurately measure the amount of light originally emitted by the shower and thus the shower energy. The rms error on the reconstructed impact parameter, which is the projected distance of the extrapolated shower track to a telescope, for Monte Carlo simulations is less than 10 m for events falling within 200 m of the centre of the array.

4.3 Scaled parameter analysis

The mean scaled width method, similar to that used by the HEGRA collaboration (Daum et al. 1997), is used to classify images as either \ensuremath {\gamma }-ray like or hadron like, in order to reject non \ensuremath {\gamma }-ray background events. The main difference to the HEGRA method is in the definition of the scaled parameter itself, in the HEGRA case this is defined as ${p}_{\rm sc} = {p}/ \langle {p} \rangle$. A lookup table is used to predict the mean width and length for a  \ensuremath {\gamma }-ray as a function of the amplitude of the shower image in the camera and impact parameter. Then the value for a particular event (p) can be compared with the expected value  $\langle {p} \rangle$ according to the formula:

\begin{displaymath}% {p}_{\rm sc} = ({p} - \langle {p} \rangle)/\sigma_{{p}}. \end{displaymath} (2)

The mean value $\langle {p} \rangle$ and the scatter  $\sigma_{{p}}$ for an event vary with the image amplitude and impact distance, as well as the zenith angle. Lookup tables are generated for 13 zenith angles (Z) from 0 \ensuremath {^{\circ }} to 70 \ensuremath {^{\circ }}, based on Monte Carlo simulations. The true impact parameter of the simulated shower is used in filling the table.

When analysing real data, the reconstructed impact parameter is used along with the image amplitude for each telescope image to find $\langle {p} \rangle$ and  $\sigma_{{p}}$ in the lookup table. Linear interpolation (in  $\cos~ ({Z})$) between the two nearest simulated values is then done to find the correct value for a particular observation zenith angle. The mean reduced scaled width (MRSW) and the mean reduced scaled length (MRSL) are then calculated by averaging over the telescope images passing the image amplitude selection cut for each event: ${\rm MRSW} = \left(\sum_{{\rm tel}}{{p}_{\rm sc}}\right) / {N}_{{\rm tel}}$.

Figure 7a shows a comparison between the MRSW from Monte Carlo simulations of protons and \ensuremath {\gamma }-rays and from real data at a zenith angle of 50 \ensuremath {^{\circ }}. It can be seen that the data (before selection cuts) correspond well to Monte Carlo simulated protons, as expected, while there is good separation between the data and Monte Carlo simulated \ensuremath {\gamma }-rays, which are chosen to have a photon index ($\Gamma $) of 2.59, similar to that previously measured for the Crab nebula (Aharonian et al. 2000).


  \begin{figure} \par\mbox{\subfigure[]{\includegraphics[height=6.2cm,width=6.5cm,... ...]{\includegraphics[height=6.25cm,width=7.5cm,clip]{5351f7b.eps} }} \end{figure} Figure 7: The distributions of mean reduced scaled width (MRSW) for Monte Carlo \ensuremath {\gamma }-ray simulations ($\Gamma $ = 2.59) a) with Monte Carlo proton simulations ($\Gamma $ = 2.70) and actual off data before selection cuts, b) with real events from the direction of the Crab nebula (data set II) after selection cuts and background subtraction. All distributions are for zenith angle = 50 \ensuremath {^{\circ }}. The vertical lines denote the standard cuts described in Table 2.
Open with DEXTER

4.4 Selection cuts

Selection cuts on the mean scaled parameters, image intensity and \ensuremath {\theta ^{2}} are simultaneously optimized to maximise the detection significance ($\sigma $, as defined by Li & Ma (1983)) for sources with typical fluxes and energy spectra. The optimisation population consists of a mixture of \ensuremath {\gamma }-ray simulations (selected to give the desired flux and spectrum for optimisation) and real background data. In the presence of background, the significance achieved for a given source increases with the square root of the observation time; instrument performance is therefore characterised by \ensuremath{\sigma/\!\sqrt{t~({\rm h}^{-1})}}. The optimised cuts yield the maximum \ensuremath{\sigma/\!\sqrt{t~({\rm h}^{-1})}} for a source of that type. It should be noted that the optimum selection cuts in any analysis depend on the energy spectrum of the Monte Carlo simulations used in the optimisation procedure, and it may be necessary to optimise selection criteria separately for much harder or much softer energy spectra. As a rule this is not done however in source searches, in order to preserve the a priori nature of the analysis. The background data used in the optimisation is then not further used, in order to avoid the possibility of optimising on background fluctuations and compromising the statistical independence of the results. The selection cuts are summarized in Table 2. A number of alternative sets of cuts are presented here, which are used for analysis of HESS sources. It will be shown that the reconstructed flux and spectrum of the Crab are consistent for these various selection criteria.

The standard set of selection cuts has been optimised to give the maximum \ensuremath{\sigma/\!\sqrt{t~({\rm h}^{-1})}} for a flux 10% of the Crab (standard cuts), with a similar spectrum. The Hard cuts are optimised for a source with a flux 1% of the Crab flux, and a $\Gamma $ of 2.0. These cuts give a higher significance for weak, hard spectrum sources, at the expense of energy threshold and cut efficiency. The hard cuts are also useful as they reduce the systematic uncertainties in sky-map reconstruction by reducing the numbers of background events, relative to the signal. They also give a narrower PSF than the standard cuts, as the higher intensity cut selects better reconstructed events. A set of Loose cuts have been also optimised to give the maximum significance for a strong source, similar to the Crab, and a $\Gamma $of 3.0. The lower intensity cut here reduces the energy threshold of the analysis relative to the standard cut, and the fraction events passing the cuts is higher. When conducting source searches, the standard cuts are always used unless there is an a priori reason to expect a very hard or very soft spectrum from the source.

For analysis of large extended sources the cut on \ensuremath {\theta ^{2}} is usually set to be larger than the extension of the source, so that effectively all \ensuremath {\gamma }-rays from the source pass this cut. In order to demonstrate the effect of this strategy, a version of the standard cuts is described, with the \ensuremath {\theta ^{2}} cut set to a much larger value. These are referred to as extended cuts for the purposes of this paper. It should be emphasised that only the standard selection cuts are used in searches for point sources, extended source searches are carried out using an a priori \ensuremath {\theta ^{2}} cut suited to the source size, and trials are taken into account when testing multiple source extensions.

Figure 7 shows the distributions of MRSW, after standard selection cuts and background subtraction (see following section), for real data and simulations (with the same mean zenith angle). The cut on MRSW is not applied for this plot. The standard MRSW selection cuts are indicated, it can be seen that the cuts select \ensuremath {\gamma }-ray-like events. The small shift between the data and Monte Carlo simulations seen in this plot is due to differences in the optical efficiency; simulations with reduced efficiency (as in Fig. 4) agree well with the data. It can be seen that this shift has a negligible effect on the efficiency of the scaled parameter cuts.

5 Signal extraction and background estimation

When estimating the flux of \ensuremath {\gamma }-rays from a particular direction in the sky it is necessary to estimate the background level, due to non \ensuremath {\gamma }-ray events with directions reconstructed close to the source direction. The significance of the excess after background subtraction is then determined using the likelihood ratio method described by Li & Ma (1983). For the purpose of background estimation the distribution of background events is usually assumed to be azimuthally symmetric within the camera field of view. However, zenith angle dependent effects or variations in the NSB level across the field of view may introduce non-radial variations in the background level. A radial profile of the relative rate of background events passing shape cuts (the background acceptance) in the field of view is shown in Fig. 8 (dashed line). For comparison, a number of test observations (duration 30 min each) were made at a range of offsets from 0 \ensuremath {^{\circ }} to 2.5 \ensuremath {^{\circ }} on the Crab nebula. It can be seen that the relative rate of excess \ensuremath {\gamma }-ray events passing cuts (points) follow the background acceptance closely out to 1.5 \ensuremath {^{\circ }} offset in the camera. The Monte Carlo predicted \ensuremath {\gamma }-ray rate for this zenith angle is also shown (solid line), this is described in Sect. 6.3.


  \begin{figure} \par\includegraphics[width=7cm,clip]{5351f8.eps} \end{figure} Figure 8: A comparison of the predicted relative \ensuremath {\gamma }-ray rates (from MC simulations) and those measured from data (from 4-telescope test observations on the Crab nebula), as a function of off-axis angle in the field of view. Also plotted is the relative radial acceptance for background events passing selection cuts. This agrees well with the relative \ensuremath {\gamma }-ray acceptance out to 1.5 \ensuremath {^{\circ }}.
Open with DEXTER

When estimating the background, first the reconstructed shower direction for each $\gamma$-like event (i.e. an event that passed the shape cuts) is filled in a two dimensional histogram (so-called raw sky-map). The on signal for a given point in the sky is determined by selecting events within a circle around that point with radius  $\theta_{\rm cut}$. Two techniques are used to derive estimates of the background level within this region of the field of view, and are described below.

5.1 Reflected background model

The simplest background estimation technique uses the signal at a single position in the raw sky-map, offset in the opposite direction relative to the centre of the field of view, to estimate the background. This technique is used in the standard wobble observation mode (described in Sect. 3). However, it suffers from relatively low statistics in the measurement of the background level due to the choice of a single reflected background position, as well as possible systematic effects caused by local inhomogeneities at the background position.

The generalised reflected background technique, which is also suitable for wobble mode, uses a number of background regions equidistant from the observation position, as illustrated in Fig. 9. The combined events from these positions are used to estimate the background at the on position, scaled by the relative area of the on and off regions. In the case of a larger on integration region the number of background regions is reduced to eliminate overlapping. The normalisation, $\alpha$, is the ratio of the solid angles of the on and off regions. As the off positions are the same distance from the centre of the field of view as the on position, correction for the relative radial background acceptance of the detector is not required. However, this method cannot be used for positions closer to the centre of the field of view than the radius of the on region, as the background positions would overlap with the source position. As all of the data described here is taken in wobble mode, this method is used in this analysis for flux and spectral measurements.

The reflected-background method may also be susceptible to systematic effects caused by non-radial variations in the acceptance, especially for large offset positions in the field of view. Non-radial effects are strongest for data with only two telescopes, where the trigger efficiency can vary with the azimuth direction of the shower impact point on the ground relative to the telescope positions. This could be corrected by applying a small non-radial correction term across the field of view. This correction has not been applied in the analysis described here. However, in order to minimise systematic effects due to non-radial acceptance variations, the direction of the wobble offset is generally alternated run by run on either side of the target position.


  \begin{figure} \par\includegraphics[width=6.5cm,clip]{5351f9.eps} \end{figure} Figure 9: Schematic illustrating the background regions described in the text. The observation position of the telescopes is marked by a cross, while the target position is marked by an X. The on region surrounding the target position is marked by a cross hatched circle. The ring-background region is filled by horizontal lines, while the reflected-background regions, at constant offset from the observation position, are filled by diagonal lines. The two background regions have equal area in this case.
Open with DEXTER

5.2 Ring-background model


  \begin{figure} \par\subfigure[]{\includegraphics[width=6cm,clip]{5351f10a.eps} }... ...mm} \subfigure[]{\includegraphics[width=5.5cm,clip]{5351f10b.eps} } \end{figure} Figure 10: a) Uncorrelated 2-d plot of \ensuremath {\gamma }-ray excess from the Crab nebula, the reconstructed sky-map in RADec J2000 coordinates around the source position is shown. (b) The distribution of significance per bin in an uncorrelated significance sky map for the same data. The open circles denote the distribution for points in the map within 0.4 \ensuremath {^{\circ }} of the Crab position, while the filled circles represent points further from the source. The Gaussian function, fitted to the second distribution, has a mean of $0.04 \pm 0.006\ensuremath{_{{\rm \tiny{stat}}}} $ and a standard deviation of $0.98 \pm 0.004\ensuremath{_{{\rm \tiny{stat}}}} $. A mean of 0.0 and standard deviation of 1.0 is expected for an unbiased significance distribution of the background.
Open with DEXTER

The ring-background technique determines the background for each position in the field of view using the background rate contained in a ring around that position (Pühlhofer et al. 2003). The internal and external radii of the ring are here chosen such that the ratio of the areas of the off to on regions is close to 7, which makes for a convenient compromise between area within the ring and distance from the on position. The inner ring radius is chosen to be significantly larger than the on region, in order to avoid signal leakage into the off region. The normalisation ($\alpha$) is given by the area ratio modified by a weight factor to account for the radial background acceptance in the camera. The ring around the on position is illustrated in Fig. 9. When estimating the background for a test position close to a known source like the Crab nebula, the source position is cut out of the background ring in order to avoid signal pollution in the off region for the test position. This method has the advantage of allowing background estimation for all positions in the field of view. However, since the number of events at positions surrounding the source is used to estimate the background in the direction of the source, it is most suitable for sources with an small angular extent relative to the field of view of the detector. Figure 10a shows an excess map of the sky in the vicinity of the Crab nebula, after background subtraction. Figure 10b shows the distribution of significance of the excesses in each bin in the sky map. It can be seen that the significance is distributed normally in the off source regions of the map (filled circles), while the region close to the Crab nebula (open circles) shows a significant excess.

The ring-background method is less suitable for spectral analysis of sources than the reflected method as the background acceptance may not be constant as a function of energy, thus the background level may not be correctly estimated for the entire energy range of the spectral analysis. Thus this method is not used for the main spectral and flux analysis in this paper, however results are derived using this method for comparison purposes.

6 Energy reconstruction and effective areas

6.1 Energy reconstruction

The energy of the primary particle of a  \ensuremath {\gamma }-ray shower is estimated for each telescope as a function of the image amplitude and impact parameter using a lookup table. The lookup table contains the mean energy for Monte Carlo \ensuremath {\gamma }-ray simulations as a function of total image amplitude and the simulated true impact parameter. As for the scaled parameters, the lookup tables are created for a number of zenith angles and the resulting energy is estimated by linear interpolation in  $\cos~ ({Z})$, and averaged over the triggered telescopes for each event. Events with relative error in the reconstructed impact parameter greater than 25% are not used in the lookup table creation, in order to reject poorly reconstructed events, which may bias the lookup table. Events with a distance greater than 2 \ensuremath {^{\circ }} from the centre of the field of view are also rejected.

In the case where no estimate is available in the lookup table for a particular event, due to a lack of Monte Carlo statistics at that combination of impact distance and total amplitude, an alternative lookup is used with coarser impact distance binning. This occurs on average for 0.3% of events. In the case where the optical correction is applied, the energy is estimated using the corrected image intensity, as described in Sect. 3.4. The energy estimated for each telescope is averaged to give the mean energy for an event: ${E}_{\rm est} = \left(\sum_{{\rm tel}}{{E}_{{\rm tel}}}\right) / {N}_{{\rm tel}}$.

6.2 Energy resolution


  \begin{figure} \par\includegraphics[width=7cm,clip]{5351f11.eps} \end{figure} Figure 11: Relative bias ( $({E}_{\rm reco}-{E}_{\rm true})/{E}_{\rm true}$) in estimation of energy as a function of energy for three zenith angles. The vertical lines represent the safe energy thresholds for spectral analysis at each zenith angle.
Open with DEXTER

The error in the reconstructed energy ( $\Delta_{{E}}$) for a particular simulated \ensuremath {\gamma }-ray event with true energy  ${E}_{\rm true}$and reconstructed energy  ${E}_{\rm reco}$ is defined as $\Delta_{{E}} = ({E}_{\rm reco}-{E}_{\rm true})/{E}_{\rm true}$. The mean value of  $\Delta_{{E}}$ is shown as a function of  ${E}_{\rm true}$ in Fig. 11. For energies close to the threshold, there is a bias due to a selection effect, whereby events with energies reconstructed with too high a value are selected. In order to make an accurate energy spectrum it is necessary to define the useful energy range, so as to avoid the region of large energy bias. First the lowest energy bin in the bias histogram with a bias of less than 10% is found. The useful lower energy threshold is the maximum energy of this bin plus 10%. This safe energy threshold is also indicated in Fig. 11 for each zenith angle. This energy threshold for the analysis is increased to take account of the shift in the energy scale when the optical efficiency correction is applied. It can be seen that the energy bias above the safe threshold does not depend on the zenith angle, up to energies in excess of 60 TeV.

The distribution of  $\Delta_{{E}}$, for \ensuremath {\gamma }-rays simulated with a power law spectrum with $\Gamma = 2.6$, at 50 \ensuremath {^{\circ }} zenith angle, is shown in Fig. 12 for the standard analysis. The energy resolution for a particular energy range is defined as the width of this distribution. Events in this plot are selected above the safe threshold of 0.44 TeV in order to avoid the effect of the energy bias. The energy resolution defines the optimum binning for spectral reconstruction, as well as defining the minimum energy width of any resolvable spectral structure. It is possible to improve the energy resolution slightly by selecting only those events with higher telescope multiplicities and smaller impact parameters, at the expense of reduced event statistics.


  \begin{figure} \par\includegraphics[width=7cm,clip]{5351f12.eps} \end{figure} Figure 12: The distribution of the relative error in the reconstructed energy per event for Monte Carlo simulated \ensuremath {\gamma }-rays with a power law energy distribution (above 440 GeV) with $\Gamma = 2.6$ at 50 \ensuremath {^{\circ }} zenith angle. The root mean square (rms) width of this distribution is 16%. The width of the fitted Gaussian distribution is 14%.
Open with DEXTER

   
6.3 Effective areas

The \ensuremath {\gamma }-ray flux from a source is estimated from the number of excess events passing the selection cuts for a particular data set using the effective area of the instrument. The effective area is a function of the zenith angle and offset of the source from the pointing direction, the energy of the event and the particular selection cuts used. The effective area is modeled from Monte Carlo simulations by counting the fraction of simulated events which trigger the detector and pass the selection cuts. This effective area has been estimated in two ways: as a function of the Monte Carlo energy of the simulated events ( \ensuremath{{A}_{\rm true}}) and as a function of the reconstructed energy ( \ensuremath{{A}_{\rm reco}}). While  \ensuremath{{A}_{\rm true}} does not depend on the energy spectrum of the simulated \ensuremath {\gamma }-rays, the finite energy resolution makes \ensuremath{{A}_{\rm reco}} sensitive to this; the effective areas are usually estimated assuming a power law distribution of $\Gamma = 2.0$. However, when estimating a flux using events binned in reconstructed energy, it is correct to use \ensuremath{{A}_{\rm reco}} to estimate the effective area for each bin, when estimating the integrated effective area over the whole energy range one may use \ensuremath{{A}_{\rm true}}. In order to avoid a bias in the spectral reconstruction when using \ensuremath{{A}_{\rm true}}, it is necessary for the energy spectrum of the simulations to match that of the data, this is discussed further below. The effective areas as a function of true energy and reconstructed energy (for the standard selection cuts) are shown in Fig. 13 for three zenith angles. As the effective area of the telescope system depends strongly on the zenith angle of the observations, it is determined for a range of angles and the value for a particular energy and zenith angle is determined by linear interpolation in $\log~(E)$and  $\cos~ ({Z})$ (Aharonian et al. 1999a).


  \begin{figure} \par\mbox{\subfigure[]{\includegraphics[width=8cm,clip]{5351f13a.... ... \subfigure[]{\includegraphics[width=8.15cm,clip]{5351f13b.eps} }} \end{figure} Figure 13: The effective collecting area of the full HESS array versus energy a) as a function of true Monte Carlo energy and b) as a function of reconstructed energy for observations at zenith angles of  $20\ensuremath {^{\circ }} $, $45\ensuremath {^{\circ }} $, and  $65\ensuremath {^{\circ }} $. The vertical lines denote the safe energy thresholds for each zenith angle, increasing in zenith angle from left to right.
Open with DEXTER

In order to simplify the application of the optical correction discussed in Sect. 3.4, this correction is not applied in estimation of the effective area for each event. Since the distribution of Cherenkov light scales with the shower energy to a good approximation, the detection probability and hence the effective area depends primarily on the amount of light arriving at the camera, and not on the absolute energy of the  \ensuremath {\gamma }-ray. Thus it is not necessary to recalculate the effective area lookup tables when the actual optical efficiency changes by a small amount. For larger changes this correction breaks down due to effects of the system trigger on events near the energy threshold. The selection cuts are made using the image intensity without optical correction, and the corrected intensity is only applied in the energy estimation. This method has been tested on Monte Carlo simulated sources with reduced optical efficiency, and it was verified that the correct flux is reconstructed.

Instead of using the safe energy threshold as introduced above, the energy threshold for a set of observations has also been commonly defined as the peak in the differential rate vs. energy curve (Konopelko et al. 1999). This is formed by folding the effective area curve, as plotted in Fig. 13a with the expected \ensuremath {\gamma }-ray flux from the source. This energy threshold is generally slightly lower than the safe threshold defined above, which is designed to ensure an accurate spectral reconstruction. There may even be a significant \ensuremath {\gamma }-ray signal below the safe threshold. The vertical lines in Fig. 13 define the safe energy thresholds for each zenith angle. Figure 14 shows the predicted peak-rate energy threshold and \ensuremath {\gamma }-ray rate for a Crab-like source, based on simulations and projection from the Crab flux as measured by the HEGRA collaboration (Aharonian et al. 2000). Table 3 gives the pre-cut energy threshold as a function of zenith angle, along with the equivalent after the various selection cuts. It should be noted that the energy threshold, by this definition, depends on the spectrum of the source.

The effective area also varies with the position of the source in the field of view of the instrument. As larger energy showers are preferentially detected at higher impact distances, they appear closer to the edge of the field of view. Thus truncation of images for sources closer to the edge of the field of view tends to reject events at higher energies. Monte Carlo simulations are made at a range of source offset positions in order to make effective area curves and the effective area is interpolated for a particular observation. Table 3 also shows the event rates before selection cuts and after each of the selection cuts described above, as a function of zenith angle for a source similar in flux to the Crab nebula. These predictions are based on the effective areas estimated for a source offset by 0.5 \ensuremath {^{\circ }} from the observation position. The cut efficiencies in each case are also shown, with the (peak rate) energy thresholds.

7 Flux and spectral measurements

The data from observations of the Crab have been analysed using the technique described, individually by data set as outlined in Table 1, and combined. The analysis has been carried out for each of the sets of selection cuts as outlined in Table 2, using the reflected background method with 5 background regions. Table 4 outlines the numbers of events passing cuts in the on and off regions, as well as the background normalisation ($\alpha$), the number of excess events, the significance ($\sigma $) of the excess, the rate of \ensuremath {\gamma }-rays passing cuts and the \ensuremath{\sigma/\!\sqrt{t~({\rm h}^{-1})}}. Data Sets I-III are combined to give a total result. For example, the mean rate of \ensuremath {\gamma }-rays for the standard selection cuts is 6.0  \ensuremath{\gamma~{\rm min}^{-1}} with a significance of 27  \ensuremath{\sigma/\!\sqrt{t~({\rm h}^{-1})}}. It can be seen that the rate of events passing cuts is strongly dependent on the zenith angle of the observations, as well as on the selection cuts used. Comparisons between Tables 4 and 3 show that the effective area estimation correctly reproduces the \ensuremath {\gamma }-ray rate for the various data sets (within statistical errors). The mean rate for data set III, at a mean zenith angle of 54 \ensuremath {^{\circ }}, is predicted to be 4.8  \ensuremath{\gamma~{\rm min}^{-1}}, while the measured value is $(4.9 \pm 0.1\ensuremath{_{{\rm\tiny {stat}}}} ) \ensuremath{\gamma~{\rm min}^{-1}} $. The loose cuts at 55 \ensuremath {^{\circ }} zenith angle keep 53% of the possible \ensuremath {\gamma }-rays, making them more suitable for spectral studies of strong sources.


  \begin{figure} \par\includegraphics[width=7cm,clip]{5351f14.eps} \end{figure} Figure 14: The peak-rate energy threshold before and after selection cuts versus zenith angle, for three sets of selection cuts as described in Table 2. The safe threshold for spectral analysis in each case is slightly higher.
Open with DEXTER

   
7.1 Run by run flux measurements

For the purpose of producing a light curve of the \ensuremath {\gamma }-ray flux from a source, the integrated flux above the threshold energy is calculated for each time period ( $t_{\rm start}$ to  $t_{\rm stop}$), assuming a particular spectral form for the source, such as a power law with photon index $\Gamma $ and flux normalisation I0 in units of  \ensuremath{{\rm cm}^{-2}~{\rm s}^{-1}~{\rm TeV}^{-1}}. The excess number of events seen from a source ($\delta$) is given by the following:

 \begin{displaymath}% \delta = \int_{0}^{{E}_{\rm c}} \int_{t_{\rm start}}^{t_{\r... ...a} \ \ensuremath{{A}_{\rm true}} (E,Z(t))\ {\rm d}t\ {\rm d}E. \end{displaymath} (3)

Here \ensuremath{{A}_{\rm true}} is the effective area as a function of zenith angle Z and true energy E from Monte Carlo simulations. The flux normalisation can then be calculated by integrating the effective area up to some upper cutoff energy  ${E}_{\rm c}$ and over the integration time. The value of  ${E}_{\rm c}$ is imposed by the range of the Monte Carlo simulations and is normally above 100 TeV. The integral flux above threshold is usually quoted above the threshold energy for the observations, or alternatively above 1 TeV. Here the latter convention is used for simplicity.

The fluxes measured for each data set included in this analysis, along with the mean fluxes for data sets I-III are given in Table 4. The mean flux is also given for the various sets of selection cuts described in Table 2. The rms variation in the mean flux per data set is 8%, while the rms variation in the run-by-run fluxes is 14.9%, the typical statistical error on a single run is 5% for moderate zenith angles and offsets. Figure 15 shows light-curves for data sets I-III, while the distribution of run-by-run flux is shown in Fig. 16, for all of the data.

It can be seen that the long-term variations in the run-by-run flux after the correction for changes in the detector optical efficiency are small compared to the short-term variations, mainly due to atmospheric effects. The mean flux in data sets II and III, taken in January 2004, is 10% higher than that seen in data set I, which was taken in October of 2003. This difference is smaller than the rms spread of either data set, and can be explained by differences in the atmospheric conditions between the two periods. No correction is made to the run-by-run flux for short term variations in the atmospheric conditions. Such corrections are under study and will be the focus of a future paper.

Table 3: \ensuremath {\gamma }-ray rate predictions from simulations for the standard, hard and loose selection cuts. This table is valid for a source with an energy spectrum similar to the Crab, for observations at an offset of 0.5 \ensuremath {^{\circ }} (the usual observing mode). The cut selection efficiencies and (peak rate) energy thresholds in each case are also given.

Table 4: Events passing cuts in on and off regions for the Crab, listed by data set along with excesses, significance and \ensuremath {\gamma }-ray rates. Various selection cuts described in Table 2 are also compared for data sets I-III combined. The results using the ring-background model are given (denoted as Ring). The integrated flux from the Crab above 1 TeV is shown also, as described in Sect. 7.1. The \ensuremath {\chi ^{2}}/ d.o.f. for a fit to a constant flux for the data set is given, as is the percentage run-by-run rms deviation in the flux.


  \begin{figure} \par\subfigure[]{\includegraphics[height=7.9cm,width=4.3cm,clip]{... ...[]{\includegraphics[height=7.9cm,width=4.3cm,clip]{5351f15c.eps} } \end{figure} Figure 15: The run-by-run light curve of the integral flux above 1 TeV for a) data set I b) data set II c) data set III. All efficiency corrections as discussed in the text have been applied to these data.
Open with DEXTER

Systematic errors in the flux due to selection cuts and effective area estimates have been studied by applying the various selection cuts described to the data and measuring the integral flux. The results are outlined in Table 4. The rms in the measured fluxes is 15%. A possible systematic error in the integral flux measurement due to the background estimation method has been tested by applying the ring-background method, as described above to calculate a flux. It can be seen that the reconstructed flux differs only slightly with the two methods.

7.2 Energy spectrum

The energy spectrum of the Crab has been measured using the data described here. The method used for deriving the energy spectrum is similar to those described by Mohanty et al. (1998), Aharonian et al. (1999b) and Aharonian et al. (2004b). An energy spectrum is fit for each data set and for the combined data.

The bin size for the energy spectrum is set depending on the overall significance of the signal. The maximum possible energy bin is defined by the simulations, which extend to above 400 TeV at 50 \ensuremath {^{\circ }} zenith angle. Only energies above the safe threshold as defined above are used for spectral determination.

In each energy bin i above the minimum energy the differential flux is calculated by summing over the on source events  $N_{\rm on}$, weighted by the inverse of the effective area ( \ensuremath{{A}_{\rm reco}}) as a function of the reconstructed energy of each event. The normalised sum of the weighted off events  $N_{\rm off}$ is then subtracted. The difference is weighted by the live-time for that bin (T) and the bin width ( $\Delta {E}_{i}$):

 \begin{displaymath}% \frac{{\rm d}F_{i}}{{\rm d}E} = ({T}~\Delta {E}_{i})^{-1}.... ...\sum_{k=0}^{N_{\rm off}}{({A}_{\rm reco_{\it k}})^{-1}}\right) \end{displaymath} (4)

In the case where runs are combined with varying zenith angles, and thus varying useful energy thresholds, the live-time is calculated for each energy bin separately. The result is scaled by the appropriate live-time to give the \ensuremath {\gamma }-ray flux in each energy bin. The error on the flux in each energy bin is estimated using standard error propagation. A spectral energy function (for example a power law distribution) is fit to the flux points using the least-squares method. The maximum energy for spectral fitting is chosen so as to have a significance in that bin greater than 2$\sigma $.


  \begin{figure} \par\includegraphics[width=7cm,clip]{5351f16.eps} \end{figure} Figure 16: Distribution of run by run fluxes for data sets I-III. The fitted Gaussian distribution has a mean of $(2.21 \pm 0.06\ensuremath{_{{\rm \tiny{stat}}}} )$ $\times $ $10^{-11}~\ensuremath{{\rm cm}^{-2}~{\rm s}^{-1}} $ and a $\sigma $ of $(3.58 \pm 0.6\ensuremath{_{{\rm \tiny{stat}}}} )$ $\times $ $10^{-12}~\ensuremath{{\rm cm}^{-2}~{\rm s}^{-1}} $.
Open with DEXTER


  \begin{figure} \par\includegraphics[width=7cm,clip]{5351f17.eps} \end{figure} Figure 17: The ratio of the reconstructed effective area (estimated assuming an photon index of 2.0) to the true effective area per energy bin, for true photon indices from 1.1 to 3.2, based on Monte Carlo simulations at 45 \ensuremath {^{\circ }} zenith angle.
Open with DEXTER

Since the estimation of the effective area as a function of reconstructed energy, \ensuremath{{A}_{\rm reco}}, depends on an assumed spectral slope, it is strictly correct to adjust \ensuremath{{A}_{\rm reco}} for the fitted spectrum and then re-fit, repeating until the fit converges (Aharonian et al. 1999b). Figure 17 shows the bias as a function of energy introduced by using an assumed photon index of 2.0 for the effective area estimation, given true spectra with indices ranging from 1.1 to 3.2. It can be seen that at 440 GeV, the differential flux for a source with an intrinsic photon index of 2.6 is overestimated by 5%, while the differential flux given a true photon index of 1.5 is underestimated by 4%. For energies well above threshold the bias is less than 5% for a wide range of photon indices. Thus the effect of this correction on the Crab spectrum is small and was neglected for this analysis.

The number of excess events and significance is given for each energy bin along with the differential flux in Table 5. Only statistical errors are given here. A significant signal is seen in every energy bin from the threshold energy of 440 GeV up to 20 TeV, and a marginal signal is seen at the maximum bin of mean energy 30.5 TeV. The spectrum is shown in Fig. 18. The fit of a power law function to the combined data with the standard analysis cuts yields $\Gamma = 2.63$ $\pm$ $0.02\ensuremath{_{{\rm\tiny {stat}}}} $ and differential flux normalisation at 1 TeV $I_{0} = (3.45 \pm 0.05\ensuremath{_{{\rm\tiny {stat}}}} )$ $\times $ $10^{-11}~\ensuremath{{\rm cm}^{-2}~{\rm s}^{-1}~{\rm TeV}^{-1}} $ (PL in Table 6).

In the combined spectral fit there is evidence for a steepening of the energy spectrum; a fit of a power law with an exponential cutoff: $\frac{{\rm d}N}{{\rm d}E} = {I}_{0}\ ({E}/1~{\rm TeV})^{(-\Gamma)} \exp~(-{E}/{E}_{\rm c})$, gives a differential flux normalisation at 1 TeV ${I}_{0} = (3.76 \pm 0.07\ensuremath{_{{\rm\tiny {stat}}}} ) \times 10^{-11}~\ensuremath{{\rm cm}^{-2}~{\rm s}^{-1}~{\rm TeV}^{-1}} $, with $\Gamma = 2.39$ $\pm$ $0.03\ensuremath{_{{\rm\tiny {stat}}}} $ and a cutoff energy ${E}_{\rm c} = (14.3 \pm 2.1\ensuremath{_{{\rm\tiny {stat}}}} )$ TeV. The \ensuremath {\chi ^{2}} for this fit is 15.9 with nine degrees of freedom. This compares with a  \ensuremath {\chi ^{2}} of 104 with ten degrees of freedom for the straight power law fit, thus the fit including an exponential cutoff is clearly favoured. A broken power law fit (BPL in Table 6):

\begin{eqnarray*}\frac{{\rm d}N}{{\rm d}E} = {I}_{0}\ \left(\frac{{E}}{{E}_{\rm ... ...m c}}\right)^{1/{\rm S}}\right)^{{\rm S}(\Gamma_{1}-\Gamma_{2})} \end{eqnarray*}


gives a differential flux normalisation ${I}_{0} = (3.43 \pm 0.07\ensuremath{_{{\rm\tiny {stat}}}} )$ $\times $ $10^{-11}~\ensuremath{{\rm cm}^{-2}~{\rm s}^{-1}~{\rm TeV}^{-1}} $, with $\Gamma = 2.51$ $\pm$ $0.02\ensuremath{_{{\rm\tiny {stat}}}} $ below the break energy of ${E}_{\rm c} = (7.0 \pm 0.1\ensuremath{_{{\rm\tiny {stat}}}} )$ TeV. The photon index above the break is $\Gamma = 3.3$ $\pm$ $1.5\ensuremath{_{{\rm\tiny {stat}}}} $, and the \ensuremath {\chi ^{2}}/d.o.f. of the broken power law fit is 28.6/8. This fit includes a term (S) for the width of the transition region, which is fixed to 0.3.

Table 5: Flux measurements for each energy bin in the combined spectral fit on data sets I-III, as plotted in Fig. 18b. The flux errors are error-propagated 68% Feldman-Cousins confidence intervals (Feldman & Cousins 1998).


  \begin{figure} \par\subfigure[]{\includegraphics[width=7.25cm,clip]{5351f18a.eps... ... \subfigure[]{\includegraphics[width=7.45cm,clip]{5351f18b.eps} } \end{figure} Figure 18: a) Energy spectra for data sets I (filled circle), II (open circles) and III (filled triangle), fit residuals to the common power law fit are also shown. The dashed line indicates the best fit power law spectrum, while the solid line denotes the fit including an exponential cutoff. b) Combined average energy spectrum for data sets I-III, fit residuals to the combined power law fit are shown.
Open with DEXTER

Table 6: Flux and spectral measurements of the Crab, divided up by data set as outlined in Table 1, for a power law fit with an exponential cutoff. Results for the various selection cuts described in Table 2 are also compared. The results for a power-law fit (PL) and for a broken power-law fit (BPL) are also given. The spectral fit estimated using the ring background model is given (ring), as is that estimated without the optical efficiency correction (uncorr.). Only statistical results are shown in the table. Similar measurements from other experiments are given for comparison, the Whipple results is taken from Mohanty et al. (1998), the CAT results from Masterson et al. (1999) and the HEGRA results from Aharonian et al. (2004b).

The flux and spectral measurements for the separate data sets are summarised in Table 6. The rms spread of the photon index from data set to data set is 0.04. The rms spread in the integral flux, calculated from the fitted spectrum, is 15%, with a statistical error of 2% on the integrated flux for the combined data. Figure 18a shows the energy spectral points superimposed for the three data sets used in the spectral analysis, the residuals about the combined fitted spectrum are shown underneath. Figure 18b shows the energy spectral fit, along with residuals, for the combined data sets I-III. The spectral fits for the various selection cuts are also included in Table 6, as is the fit for the ring background model. For comparison the fit spectrum calculated without the optical efficiency correction is also included.


  \begin{figure} \par\includegraphics[width=7.05cm,clip]{5351f19.eps} \end{figure} Figure 19: Contour plot of the \ensuremath {\chi ^{2}} fit error as a function of the power law parameters F0 and $\Gamma $ for data sets I (filled circle), II (open circle) and III (filled triangle, dotted lines). The 68%, 95% 99.9% error contours are shown, and The best fit spectral parameters are marked in each case.
Open with DEXTER

Figure 19 shows a 2-d plot of the fitted photon index $\Gamma $ against the flux normalisation for each data set analysed, here a simple power law fit has been made for simplicity. The error contours, estimated using the least-squares method, are also shown. It can be seen that the three data sets are compatible at the $2\sigma$ level; data set II includes large zenith angle data, and is fitted with a slightly softer spectrum, caused by the higher energy threshold of these observations and the significant curvature seen in the spectral measurements.

The combined data sets I-III have also been analysed with the various selection cuts described in Table 2, and a spectrum fitted. The rms spread in the photon index between the various analyses is 0.08, however this includes the effect of the very different energy threshold for the hard selection cuts, which may give rise a softer photon index if the source spectrum is intrinsically curved. The rms spread of the reconstructed integral flux is 8%, which indicates that the reconstructed flux is not strongly dependent on the details of the analysis method. The use of the ring-background method results in a flux and spectral slope similar to that reconstructed using the standard reflection method.

7.3 Estimation of systematic errors

The systematic error on the absolute flux is estimated from the various independent contributing factors, as discussed in Sect. 2.2. These errors are summarised in Table 7. The total estimated systematic error, after correction for degradation in the optical efficiency, is 20%. The sources of systematic error include uncertainties due to the shower interaction model and the atmospheric model used in the Monte Carlo simulations. Also included is the estimated uncertainty in the flux due to the effect of missing pixels, which has been conservatively estimated at 5%, and the effect of uncertainty in the live-time measurement, which is less than 1%. The uncertainty in flux due to selection cuts is estimated from the rms of the flux and spectral slope measurements detailed in Table 6, as is the uncertainty due to the background model. The run-by-run rms over the entire set of data is 15%, this is thought to mainly be due to variations in the atmosphere, and thus is included in the systematic error as an independent factor. The rms of the spectral estimations for the various datasets in Table 6 is used to estimate the uncertainty in the spectral slope, which is 0.1.

8 Conclusions

A strong signal has been detected from the Crab nebula during the commissioning phase of the HESS instrument and with the complete instrument. An energy spectrum has been measured, with a differential spectrum described by a power law with slope $\Gamma = 2.39$ $\pm$ $0.03\ensuremath{_{{\rm\tiny {stat}}}}\pm 0.09\ensuremath{_{{\rm\tiny {sys}}}} $ and an exponential cutoff at $(14.3 \pm 2.1\ensuremath{_{{\rm\tiny {stat}}}}\pm 2.8\ensuremath{_{{\rm\tiny {sys}}}} )$ TeV. The integral flux above 1 TeV is $(2.26 \pm 0.08\ensuremath{_{{\rm\tiny {stat}}}}\pm 0.45\ensuremath{_{{\rm\tiny {sys}}}} ) \times 10^{-11}~\ensuremath{{\rm cm}^{-2}~{\rm s}^{-1}} $.

Table 7: Summary table showing the various estimated contributions to the systematic flux error.


  \begin{figure} \par\includegraphics[width=7cm,clip]{5351f20.eps} \end{figure} Figure 20: Comparison between the spectral measurements of CAT (open squares), Whipple (open circles) and HEGRA (open triangles) and the results if this study (filled circles). The residuals are shown for each spectrum relative to a power law fit to the HESS data (dotted line).
Open with DEXTER

Marginal steepening in the spectrum measured on the Crab nebula has been previously claimed by Aharonian et al. (2004b) in studies of the Crab with the HEGRA experiment. Figure 20 compares the Crab spectrum from this study with measurements by HEGRA, CAT and Whipple. Acceptable agreement is seen up to 10 TeV between the experiments, although the CAT result gives a somewhat steeper spectrum; the rms variation in integral flux between the four experiments is 22%. Above 10 TeV the energy spectrum as seen by HESS steepens significantly, in particular compared to HEGRA.


  \begin{figure} \par\includegraphics[width=7cm,clip]{5351f21.eps} \end{figure} Figure 21: Sensitivity of the HESS array, expressed as the amount of time required to detect a signal at the $5\sigma $ level, as a function of the flux of the source, for a source of similar spectral slope to the Crab nebula (for observations at 20 \ensuremath {^{\circ }} zenith angle, 0.5 \ensuremath {^{\circ }} offset from the source). A minimum of 10 excess events is also required for very short observations. Shown are curves for the standard, hard and loose selection cuts.
Open with DEXTER

The softening seen in the Crab spectrum at high energies is consistent with models of inverse Compton emission due to a population of electrons extending up to PeV energies. Due to the high magnetic field in the Crab nebula, the dominant target photon field for emission is probably created by synchrotron emission from the same electron population (Hillas et al. 1998). More detailed models of \ensuremath {\gamma }-ray emission in the Crab nebula are discussed by de Jager & Harding (1992); Atoyan & Aharonian (1996).

Given the agreement between the Monte Carlo simulations and the data, one can use the simulations to predict the time required to detect a point source of a certain strength as a function of zenith angle. In Fig. 21 we show the time as a function of signal strength required for a $5\sigma $ detection at 20 \ensuremath {^{\circ }} zenith angle, for the selection cuts described here. The HESS array is capable of detecting a point source with a flux of 1% of the Crab nebula in 25 h, or alternatively detecting a source of similar strength to the Crab in 30 s. The sensitivity for extended sources decreases approximately linearly with the source extension. This sensitivity is unprecedented in the field of VHE astrophysics and opens a new window for sensitive and precise measurements of VHE \ensuremath {\gamma }-ray sources.

Acknowledgements
The support of the Namibian authorities and of the University of Namibia in facilitating the construction and operation of HESS is gratefully acknowledged, as is the support by the German Ministry for Education and Research (BMBF), the Max Planck Society, the French Ministry for Research, the CNRS-IN2P3 and the Astroparticle Interdisciplinary Programme of the CNRS, the UK Particle Physics and Astronomy Research Council (PPARC), the IPNP of the Charles University, the South African Department of Science and Technology and National Research Foundation, and by the University of Namibia. We appreciate the excellent work of the technical support staff in Berlin, Durham, Hamburg, Heidelberg, Palaiseau, Paris, Saclay, and in Namibia in the construction and operation of the equipment.

References

 
Copyright ESO 2006