VHE Gamma-Ray Observation of the Crab Nebula and its Pulsar with the MAGIC telescope

We report about very high energy (VHE) gamma-ray observations of the Crab Nebula with the MAGIC telescope. The gamma-ray flux from the nebula was measured between 60 GeV and 9 TeV. The energy spectrum can be described with a curved power law dF/dE=f0 (E/300 GeV)^(a+b*log10(E/300 GeV)) with a flux normalization f0 of(6.0+-0.2)*10^-10 1/(cm^2 s TeV), a=-2.31+-0.06 and b=-0.26+-0.07. The position of the IC-peak is determined at (77+-47) GeV. Within the observation time and the experimental resolution of the telescope, the gamma-ray emission is steady and pointlike. The emission's center of gravity coincides with the position of the pulsar. Pulsed gamma-ray emission from the pulsar could not be detected. We constrain the cutoff energy of the spectrum to be less than 27 GeV, assuming that the differential energy spectrum has an exponential cutoff. For a super-exponential shape, the cutoff energy can be as high as 60 GeV.


INTRODUCTION
The Crab Nebula is the remnant of a supernova explosion that occurred in 1054 A.D. (e.g. Collins et al. 1999, and references therein) at a distance of ∼ 2 kpc. It is one of the best studied non-thermal celestial objects in almost all wavelength bands of the electromagnetic spectrum from 10 −5 eV (radio) to nearly 10 14 eV (γ-rays). There is little doubt that the engine of the nebula is the pulsar PSR B0531+21 (hereafter Crab pulsar) in its center.
In very high energy (VHE) γ-ray astronomy 23 the Crab Nebula was first detected with large significance at TeV energies by the pioneering Whipple telescope (Weekes et al. 1989). The nebula turned out to be the strongest source of steady VHE γ-ray emission in the Galaxy. It is therefore used as the standard "calibration candle" for ground-based γ-ray experiments. Apart from testing the performance of γ-ray instruments another aim is to increase the measurement precision of the Crab Nebula flux and the energy range covered. These continuous efforts provide insights necessary for the understanding of the very details of the emission mechanisms 23 We refer to VHE as γ-rays with energies between 10 GeV and 100 TeV of VHE γ-rays in the Crab Nebula and pulsar. Important questions remain to be answered concerning the emission mechanisms of the nebula and of the pulsar at GeV energies and of the nebula around PeV energies. Since its discovery, detailed studies of the VHE emission energy spectrum, ranging from several hundred GeV up to 80 TeV, have been carried out (e.g. by Akerlof et al. 1990;Vacanti et al. 1991;Baillon et al. 1991;Goret et al. 1993;Konopelko et al. 1996;Tanimori et al. 1998;Hillas et al. 1998;Amenomori et al. 1999;Majumdar et al. 2002;Aharonian et al. 2004Aharonian et al. , 2006. Between 10 GeV and ∼ 200 GeV, observations are sparse. A few results are provided by converted solar concentrator arrays that use the wavefront sampling technique (Oser et al. 2001;de Naurois et al. 2002;Arqueros et al. 2002). However, the wavefront sampling suffers from relatively large uncertainties in the calculation of the effective area and of the energy and suffers from poor γ/hadron discrimination, which make it difficult to perform differential flux measurements.
A good explanation of the nebular dynamics and the observed energy spectrum below GeV energies can be obtained with the magneto-hydrodynamic model suggested first by Rees & Gunn (1974) and developed further by Kennel & Coroniti (1984a,b). In this framework, the pulsar provides a continuous flow of charged particles (pulsar wind) with Lorentz factors of 10 6 − 10 7 . A standing reverse shock forms where the wind ram pressure balances the total pressure of the nebula. Wind particles accelerate in the shock to ultra-relativistic energies and subsequently lose their energy by synchrotron emission. The presence of synchrotron emission up to a few hundred of MeV in conjunction with the observed γ-ray spectrum at TeV energies shows that particle acceleration takes place up to energies of ∼ 10 15−16 eV. From the total luminosity of the synchrotron emission, it can be inferred that about 10% of the pulsar's energy loss rate is converted into the kinetic energy of particles.
Above 1 GeV the dominant source of γ-ray emission is most likely inverse Compton (IC) scattering of synchrotron photons by the synchrotron emitting electrons in the shocked wind region (Synchrotron Self Compton model (SSC)). (Gould 1965;de Jager & Harding 1992;Weekes et al. 1989). To explain the observed VHE γray spectrum, several other seed photon fields are also believed to contribute to the inverse Compton scattering, namely, far-infrared excess, cosmic microwave background and mm-photons (e.g. Aharonian et al. 2004).
Although the IC-mechanism gives a good description of the observed energy spectrum between 500 GeV and about ten TeV, other processes may contribute in part to the VHE γ-ray emission. It is likely that a significant fraction of the mechanical energy lost by the pulsar is taken away by a hadronic component in the wind. Following interactions of this component with the interstellar medium, VHE γ-rays are emitted by decaying π 0 s, which modify the energy spectrum at TeV energies and beyond (Atoyan & Aharonian 1996;Bednarek & Protheroe 1997;Bednarek & Bartosik 2003;Amato et al. 2003). Atoyan & Aharonian (1996) discuss the possibility of an "amplified" bremsstrahlungs flux at GeV-energies, which could account for the discrepancy between the measured GeV γ-ray flux and predictions within the SSC-framework (de Jager et al. 1996). If this is true, one should observe, in good approximation, a power-law spectrum between 100 GeV and 10 TeV with a spectral index 2.5-2.7. Another mechanism to be mentioned is IC-scattering of relativistic electrons in the unshocked pulsar wind. If the target photons are emitted by the pulsar, a pulsed component could extend to γray energies of several 100 GeV (Bogovalov & Aharonian 2000). An independent measurement in the intervening region between 60 GeV to 400 GeV would constrain further the parameters of various models. The MAGIC imaging air Cerenkov telescope (IACT) has a low energy trigger threshold (∼ 50 GeV) and is currently the only experiment capable of exploring this energy regime.
The spatially resolved morphology of the nebula is of a complex nature. Its size in optical bandwidths is about 4 ′ × 6 ′ . Due to synchrotron losses, the high energy electrons will have shorter cooling times and only the lower energy electrons will reach out further into the nebula. Thus the effective source size is expected to shrink with increasing energy of the radiation. The radio emission is expected to extend up to and beyond the filaments optically visible, whereas X-ray and multi-TeV gamma rays should be produced in the vicinity of the shock. On the other hand, the expected source size would increase if the presence of an ionic component is established. In a special study the HEGRA collaboration concluded that the RMS size of the VHE γ-ray emission region is < 1.5 ′ for energies > 1 TeV . A similar study at energies below 1 TeV has not been performed up to now.
The Crab pulsar is a source of pulsed radiation and has been detected up to GeV energies. The Crab pulsar has a period of 33 ms and was first discovered in radio (Staelin & Reifenstein 1968) and shortly afterwards as the first pulsar in the optical domain (Cocke et al. 1969). Since then, pulsed emission from the Crab pulsar has been detected at all accessible energies up to γ-rays (see Thompson et al. 1999, for a compilation of the broad band emission). EGRET detected the Crab pulsar in γ-rays up to energies of 10 GeV with a hint of a cutoff in the energy spectrum in the highest energy bin at 6 GeV (Nolan et al. 1993).
The high energy emission from the pulsar is assumed to be due to curvature and synchrotron radiation from relativistic charged particles which are forced to move along magnetic field lines inside the magnetosphere of the pulsar. The question of where the particles are being accelerated is the subject of ongoing theoretical activities. In the two most popular models, the production of electrons and positrons and their acceleration takes place either above the polar cap of the neutron star (e.g. Harding et al. 1978;Daugherty & Harding 1982) or in outer gaps in between the null surface and the light cylinder of the magnetosphere (e.g. Cheng et al. 1986b,a;Chiang & Romani 1992). We should not omit the slot-gap model, which places the acceleration zone at the outer rim of the polar cap (Arons 1983;Muslimov & Harding 2003). These models differ in the predicted shape and cutoff of the energy spectrum at the highest energies. A measurement of the turnover in the spectrum would shed light on the possible sites of the particle acceleration and can constrain models.
In this paper we report about the observation of the Crab Nebula and pulsar with the MAGIC telescope between October 2005 and December 2005. After describing the MAGIC telescope and the performed observations in Section 2 we present the analysis chain in Section 4. In Section 5 we present our results on the steady emission comprising the differential flux between 60 GeV and 9 TeV, a study of the morphology of the emission region and the measured integral flux > 200 GeV. In Section 5.4 we present the results of our search for pulsed emission and close with a discussion of our results in Section 6.

THE MAGIC TELESCOPE
The MAGIC (Major Atmospheric Gamma Imaging Cherenkov) telescope; see Lorenz (2004), is located on the Canary Island La Palma (2200 m asl, 28.45 • N, 17.54 • W). MAGIC is currently the largest single dish imaging atmospheric Cherenkov telescope (IACT). It has a 17 m diameter tessellated reflector consisting of 964 (0.5 × 0.5) m 2 diamond-milled aluminium mirrors, which are grouped onto support panels in units of four. Depending on the elevation of the telescope the position of every panel is adjusted by computer-controlled actuators of the so-called Active Mirror Control, thus providing optimal focusing. About 80% of the light from a point source is focussed within a radius of 0.05 • in the focal plane. The MAGIC telescope is focused to a distance of 10 km-the most likely location of the shower maximum for 50 GeV γ-ray air showers at small zenith angles.
The faint Cherenkov light flashes produced in air showers are recorded by a camera comprising 577 photomultiplier tubes (PMTs). The inner part of the camera (radius ∼ 1.1 • ) is equipped with 397 PMT (type ET 9116A; 1 ′′ ∅ from Electron Tubes Ltd (ET)) with a diameter of 0.1 • each . The outer part of the camera is equipped with 180 PMT (type ET 9117A; 1.5 ′′ ∅) of a diameter of 0.2 • . The central PMT is modified for optical pulsar studies (Lucarelli et al. 2005). Hollow hexagonal nonimaging light concentrators (often called light catchers) are placed in front of all photomultipliers to compensate for the dead space between them and in order to shrink the observation solid angle to the reflector. The entrance window of the PMTs is coated with a diffuse lacquer doped with a wavelength shifter (WLS) (Paneque et al. 2003). The combination of the hemispherically shaped PMT, the light catcher, the diffuse coating and the WLS results in a 15%-25% higher quantum efficiency compared to flat window PMTs. For protection purposes (humidity, dust) a thin entrance window made of plexiglas (type UG-218, with a UV cutoff around 290 nm) is placed in front of the camera.
The PMTs have 6 dynodes and operate at a gain of roughly 30,000 to slow down aging and damage from high currents by light during observations close to the galactic plane and during moonshine. After amplification, the fast analogue signals are converted to optical signals and transported by 160 m long optical fibers to the counting house. There the optical signal is converted back and split. Part of the signal is routed to the trigger. The current configuration of the MAGIC camera has a trigger region of 2.0 degrees in diameter . This provides a γ-ray trigger collection area of the order of 10 5 m 2 (at 200 GeV for a source close to zenith). Presently, the trigger energy range spans from 50-60 GeV (at small zenith angles) up to tens of TeV. An event is triggered if the signals in each of 4 neighbouring pixels exceed a threshold of ∼ 7 photoelectrons within a coincidence time window of 6 ns.
Before being digitized by an 8-bit, 300 MSamples/s FADC system, each signal is shaped with a time constant of about 6 ns. The FADC continuously writes the digitized amplitude information into a ring buffer. In case of a trigger the digitization stops and the corresponding part of the ring buffer is written onto a disk. The dead time introduced by the readout is 25 µs. In order to expand the dynamic range to ∼ 1000, the signal of every PMT is split into two branches, differing by a factor of ten in gain. The higher gain branch is read out for a 50 ns time interval. When the signal amplitude exceeds a preset threshold, the delayed lower gain is routed to the same FADC channel and recorded in the following 50 ns. Otherwise the signal from the high gain branch continues to be recorded and is used to determine the pedestal offset of each PMT channel.
The accuracy in reconstructing the direction of incoming γ-rays on an event by event basis, hereafter γ-ray point spread function or γ-PSF, is about 0.1 degrees, depending on the energy. With the information provided by a starguider camera, mispointing is corrected to an absolute precision of about 1 ′ . A γ-ray source with an absolute intensity of ∼ 2% of the Crab Nebula and similar spectrum can be detected with MAGIC within 50 hours at energies > 200 GeV on a significance level of 5 σ.

OBSERVATIONS AND DATA SELECTION
Observations of the Crab Nebula with MAGIC are conducted on a regular basis, as a means to monitor the performance of the telescope. In this report we restrict ourselves to the analysis of data obtained in the first observation cycle of the MAGIC telescope between October and December 2005. The observations were performed in the so-called ON/OFF mode. The telescope was pointed towards the Crab pulsar (ON) for about 16 hours. An OFF source position, a sky region where no γ-ray source is known, was observed in the same range of zenith angles as the ON source. For the background estimation we used OFF data collected for over 19 hours.
One of the main objectives of this analysis was to explore the lower energy range of accessible γ-ray energies. The energy threshold of IACTs depends strongly on the zenith angle of observation. Restriction to events with low zenith angles provides the lowest possible energy threshold. Therefore, we select events with zenith angles 20 • . For any given night the data affected by technical problems or fluctuations in the data rate in excess of 10 % were rejected. The atmospheric conditions were judged from the nightly averaged and publicly available atmo-spheric extinction coefficients from the nearby Carlsberg Meridian telescope (CMT 2007). Within the selected nights, the atmospheric light transmission changed by less than 5%. The nights with data that survived all the selection criteria, together with the corresponding observation times, trigger rates and zenith angle range are listed in Table 1. The selected data sample comprises 14 nights amounting to a total ON observation time of 955 minutes (∼ 16 hours).

DATA ANALYSIS
The data analysis was carried out using the standard MAGIC analysis and reconstruction software (Bretz & Wagner 2003). After removing faulty and unstable camera channels, which amount to 3-5% of the total number of PMTs, the signal amplitudes, extracted with the digital filtering method (Albert 2006), were converted to photoelectrons (phe) by using the F-factor method (Mirzoyan & Lorenz 1997). Using calibration events (Gaug et al. 2005) recorded interlaced to normal events, the conversion factors were updated every 10 seconds. There is a 10% uncertainty in the calibration that directly propagates to the uncertainty of the event energy scale. Time offsets between pixels are also corrected with a precision of better than 1 ns.
After calibrating the data, pixels with faulty reconstructed signal amplitudes and times were rejected and the corresponding amplitude and time information was interpolated from the signals and times of neighbouring pixels. Before image parametrization, a tail-cut cleaning of the image was performed, requiring signals higher than a pre-defined absolute amplitude level and time coincidences (3.3 ns) with neighbouring channels, thus effectively suppressing pixels containing only a signal from the night sky. For most of the analysis the minimum required pixel content is 6 phe for so-called core pixels and 4 phe for boundary pixels. For the morphology studies the minimum pixel contents were raised to 10 phe (core) and 8 phe (boundary), respectively, which improves the angular resolution, albeit at the expense of an increased analysis threshold.
Every cleaned event was parameterized by a principal component analysis, commonly referred to as Hillas parameterization (Hillas 1985). The parametrization was later used to separate between γ-ray event candidates and background event candidates. The Hillasparameters Dist, Length, Width and also Alpha are illustrated for a recorded shower image in Figure 1. Another useful parameter is the Size of a shower image, the intensity of the image after image cleaning in units of recorded photoelectrons. Note that Size depends on the applied tail cuts. Size is a good estimate of the primary particle energy, provided the shower impact distance to the telescope principal axis is below ∼ 120 m. An event pre-selection was performed by discarding event candidates affected by noise and pick-up (e.g. car flashes) and event candidates with a low number of pixels (after tail-cuts typically a minimum number of 5 core pixels was requested). In addition, an image-Size of at least ∼ 100 phe was requested. Figure 2 shows the energy distribution of simulated γ-ray events with Size > 100 phe. The distribution peaks at an energy of 75 GeV.  A Monte Carlo simulation properly describing data is a necessary requisite for a ground-based γ-ray experiment. In Figure 3 the image parameter Width is shown for simulated γ-rays and γ-rays extracted from data. The four figures are for consecutive bins in Size covering the entire range of analyzed γ-ray energies. For an unbiased comparison, loose cuts have been applied in the γ/hadronseparation (explained in the next section). The γ-ray excess was obtained by subtracting the scaled distribu- tion of the OFF-data sample from the distribution of the ON-data sample. The scaling factor was found by normalizing the |Alpha|-distributions of both samples between |Alpha|=30 • and 70 • degrees. The comparison was done by selecting events with small |Alpha| (typically less than 10 • ). A small |Alpha| is expected for γ-rays from the Crab Nebula, as will be explained later. The Monte Carlo simulated distributions are compatible with the distributions extracted from data in all four Size-bins.

Gamma Hadron Separation
Only a small fraction between 10 −3 and 10 −4 of the recorded data are γ-ray showers. The major fraction of the recorded events are due to cosmic-rays of hadronic origin. This unwanted background has to be suppressed offline. For the purpose of γ/hadron separation we applied a multivariate method, the Random Forest (Breiman 2001;Bock et al. 2004), which uses the Hillas-parameters to compute the Hadronness of an event. The Hadronness of an event quantifies the probability of an event to be "γ-like" or "hadron-like". The Random Forest is trained with Monte Carlo simulated γray events and either with simulated hadronic cosmic-ray events or, as in the present case, with background events recorded by MAGIC. In this study, we used the image pa-rameters Size, Dist, Width and Length, as well as the third moment along the major axis of a shower image and a parameter describing the Concentration of a shower image, to train the Random Forest. Figure 4 shows the Hadronness for Monte Carlo simulated γ-ray showers (red) and for recorded background (blue) as a function of Size. A clear separation between both populations is visible for Size 300 phe, corresponding to γ-ray energies 150 GeV. Below 300 phe both populations start to overlap until, at ∼ 200 phe (∼ 100 GeV), no more separation is possible. The lower panel in the figure shows the maximum quality factor Q (ǫ γ / √ ǫ B ) for each size interval obtained for an optimized Hadronness-cut. ǫ γ is the fraction of retained γ-ray events and ǫ B the fraction of retained background events. The corresponding Hadronness-cut is shown by the stars in the upper panel of the figure. In addition the Random Forest method is used to estimate the energy of each event (Albert et al. 2007). An energy resolution of ∼ 20% is achieved for events with energies > 200 GeV. The energy resolution reduces to ∼ 30% around 70 GeV.

Signal Extraction
After γ/hadron-separation and energy estimation, the γ-ray signal is extracted from an |Alpha|-distribution. Alpha is the angle between the major axis of the  recorded shower and the vector connecting its center of gravity (CoG) with the source position in the camera plane (cf. Figure 1). Shower images of γ-rays from the source point with their major axes toward the source position in the camera and appear as an excess at small values in the |Alpha|-distribution. Figure 5 shows two |Alpha|-distributions from the data (ON-data: black; OFF-data: red). The left panel in Figure 5 shows the |Alpha|-distribution of events with estimated energies between 80 GeV and 100 GeV; the right panel shows the distribution of events with reconstructed energies above 200 GeV. In both cases an excess of γ-ray events is clearly visible. However, the significance of the γ-ray signal at lower energies is considerably reduced compared to the higher energies (cf. Figure 5), due to the degradation of the background suppression towards lower energy.
At the lowest energies Alpha is currently the only mean by which it is possible to separate γ-rays and background events. This is illustrated in Figure 6, which shows the quality factor as a function of Size separate for an optimized Alpha-cut and optimized Hadronness-cut as well as the combination of both cuts. Below Size 250 phe a cut in Hadronness does not improve γ/hadron-separation and reduces only statistics. On the other hand, with Alpha a quality factor of two is still possible.

Sensitivity
The highest integral sensitivity for γ-ray emission from the Crab Nebula is obtained above ∼ 250 GeV. The integral sensitivities for energies above 250 GeV obtained on a daily basis are listed in Table 2. The sensitivities are calculated in units of σ LiMa / √ hours, the significance σ LiMa is calculated using equation 17 from Li & Ma (1983). In the last column of the table, the sensitivity is expressed as the minimum flux, normalized to the Crab flux that can be measured with 5 σ significance 24 in a 50-hour observation. The day-by-day sensitivities vary by about 10%, indicating a stable telescope performance throughout the observations. The energy dependence of the sensitivity was studied by calculating the integral sensitivity for several analysis thresholds. The Size-dependent Hadronness-cut (black solid line in Figure 4) was used for γ/hadron separation and the Alpha-cut was tightened with increasing energy based on studies from MC-simulations. The integral sensitivity for a 50-hours observation is 13% Crab at 75 GeV and improves continuously with increasing the analysis threshold to about 2.2% above ∼ 300 GeV. Note that above 1.6 TeV the background is estimated from only 2 events, which results in an uncertainty of more than 70% on the integral sensitivity.  the use of Monte Carlo simulation based on many linked processes using either physical models (example: the simulation of electromagnetic showers) or calibrating them in separate measurements (example: the spectral mirror reflectivity). In some cases one has to make reasonable guesses (example: photoelectron collection efficiency in the PMT front-end volume). The calibration of individual effects suffers partly from cross-correlations, which are not always well understood. Currently, the best approach is to estimate the systematic uncertainties (commonly called systematic errors) of the various parameters separately and to combine them to a global systematic error. Here, we follow the general practice of adding these individual errors in quadrature although this will result in a slight underestimate of the total systematic error. Table 4 lists the dominant systematic error contributions (≥ 2%) for the spectral parameters (flux, slopes, cut-offs). The systematic errors influence the spectrum in different ways. Some (Nr 1-7,10,11) result in an uncertainty of the energy scale and thus can enter as a large factor in the flux at a given energy in case of steep spectra; others (Nr 8,9,(12)(13)(14)(16)(17)(18) linearly influence the flux normalization for the spectrum.
The most critical contributions to the systematic error come from the uncertainties in the conversion of photons to measurable photoelectrons (combined under item 7), the so-called photon detection efficiency (PDE). The PDE is a combination of many small effects such as the reflectivity variation of the light catchers, tolerances in the light catcher geometry, angular effects on the PMT surface, non-uniformity of the diffuse lacquer coating, the QE-spread and cathode non-uniformity of the PMTs, the photoelectron collection efficiency in the PMT front-end volume and gain variations of the first dynode. Also some contribution of the signal transmission to the DAQ is included.
Fortunately, the PDE can normally be measured with a light source uniformly illuminating the camera with short blue or UV light pulses. Obviously, the light pulser itself introduces some systematic errors such as in the absolute light flux determination, small deviations from uniformity in the illumination, some (small) temperature drift and amplitude jitter. Also, the used method of determining the number of detected photoelectrons, the above mentioned so-called F-factor method introduces some uncertainty.  (Garczarczyk 2006). The measurements cover the Crab observation period from Oct. to Dec. 2005 Another rather big uncertainty is the effective reflectivity (defined as the light from a source at infinity being focussed onto the area of a pixel). Comparing the measured brightness of a star and its image back-reflected by a high quality diffuse reflector in the camera plane allows one to carry out a routine reflectivity measurement, see Figure 7, with an uncertainty of about 7%. A similarly large error contribution was estimated for the event reconstruction. Again, many small effects contribute to the reconstruction losses or to the wrong assignment of events. In contrast to the procedure to limit the uncertainty as in the example of the PDE, no simple method to cross-check the error range of the reconstruction efficiency is possible and a reasonable guess had to be made.

Saturation and nonlinearities of electronic and opto-electronic components
Effects which influence the slope of reconstructed γ-ray energy spectra (Class C effects in Table 4) are mostly dominating at the lowest and highest accessible energies. The estimate of the systematic slope error is rather difficult. In case of a power law or moderately curved power law we conservatively estimate an uncertainty on the slope of 0.2. We note that measurements by current second generation telescopes of the spectral slope of the Crab Nebula agree better than 0.1 in the overlapping energy range.
In summary we obtain a systematic energy scale error of 16%, a systematic error of 11% on the flux normalization (without the energy scale error) and a systematic slope error of ±0.2 (which is a combination of error 14 and the other relevant class A errors averaged over the energy).

Differential Energy Spectrum of the Crab Nebula
By extracting the γ-ray signal in each bin of the reconstructed energy E rec , a spectrum N i of γ-rays in each E rec -bin i can be constructed. The reconstructed energy is subject to a bias. Before determining a differential γray flux in true energy E true bins, the spectrum N i has, therefore, to be converted into a spectrum M j of γ-rays in bins of E true . This is done by applying an unfolding procedure with regularization (Anykeyev et al. 1991). An essential input for the unfolding procedure is the migration matrix, which describes the migration of events from bin i in E rec into a bin j of E true . The migration matrix is determined from MC simulated γ-ray showers. The unfolding is done independently for different regularization schemes (Tikhonov & Arsenin 1979;Bertero 1989;Schmelling 1994). Figure 8 shows one distribution of excess events from the Crab Nebula after unfolding by the method of (Bertero 1989). The integral rate of excess events is 0.4 Hz. The differences between the unfolded points M j obtained with the different regularization schemes are used to estimate a systematic error due to the unfolding. Figure 10 shows the differential γ-ray flux, which was obtained with the regularization scheme proposed by Bertero (1989) and by normalizing the unfolded spectrum M j , to the effective collection area (Figure 9), the effective observation time and the bin width of E true (given by the horizontal bars at each flux point in the figure). The average differential flux for each energy bin is presented in Table 5.
The influence of different choices for tail-cuts, Hadronness-cuts, Dist-cuts and core-pixel cuts on the measured flux is indicated in Figure 10 by the shaded region and quoted as systematic uncertainty in Table 5. Due to analysis uncertainties, the band broadens at low energies, mostly because of limited γ/hadron discrimination power. It broadens at the highest energies due to low event statistics.
The energy spectrum is parameterized with both a power-law and a curved power-law ansatz. The fit takes into account correlations between the spectral points that are introduced by the unfolding procedure. A correlated fit with a power-law dF dE provides a flux normalization f 0 of (5.7 ± 0.2 stat ) × 10 −10 cm −2 s −1 TeV −1 and a spectral index Γ of −2.48 ± 0.03 stat ± 0.2 syst . The χ 2 of the fit is 24 for 8 degrees of freedom, which disfavors a pure power-law description of the spectrum. The energy spectrum is better described by a curved power-law ansatz dF dE = f 0 (E/300 GeV) (a+b log 10 (E/300 GeV)) (2)  (Bertero 1989). The results of a fit of the spectrum with a power law and a broken power law are also shown. The lower panel shows the relative residuals between the fit and the data points. See text for further discussion.

TABLE 5
Mean energy and differential flux of the spectral points shown in Figure 10. The systematic errors are derived from different applied cuts and unfolding procedures.
where F i,j is the differential flux measured at E i,j . The four derived spectral indices at ∼ 150 GeV, ∼ 300 GeV, ∼ 1 TeV and ∼ 2.5 TeV shown in Figure 12 indicate a clear softening of the spectrum with increasing energy. The spectral index Γ ′ was also derived from the aforementioned results of the curved power-law fit, Γ ′ = a + 2 b log (E/300 GeV) (5) and is shown by the solid black curve and the ±1σconfidence band by the dashed black curves. A systematic uncertainty on the slope can cause an additional vertical shift of the measurement by ±0.2. Within uncertainties, the measured spectral index varies in good agreement with predictions by Aharonian et al. (2004) (blue curve), who, in addition to the IC-scattering on synchrotron photons, included several other soft photon fields such as mm-photons, CMB and far IR photons from dust and stars.
The predicted GeV γ-ray emission has a peak in the SED-representation (see Figure 11). If one assumes that the energy spectrum around the peak can be described with a curved power-law, the position of the peak can be determined from the measurement of the spectral index obtained from the result of the curved power-law fit. A necessary condition for the peak in the SED is that the spectral index Γ ′ is −2. With this condition the peak is determined at 77 ± 47 stat +107 −46 syst GeV (triangle in Figure  12).

VHE γ-Ray Light Curve of the Nebula Emission
In the VHE γ-ray astronomy community it is assumed that the Crab is a constant and stable γ-ray source and can therefore be used as a standard candle. However, with more sensitive measurements it is necessary to check the stability of the γ-ray source. Below we present a time-resolved measurement of the VHE γ-ray flux, i.e. the light curve for the Crab Nebula. Depending on the source strength and the analysis threshold, the time intervals can be as short as a few minutes.
We calculated light curves in bins of ten minutes from events with estimated energies above 200 GeV. The light curves of all 14 selected nights are shown in Figure 13. Note that the same loose cuts are used for the γ/hadronseparation as for the calculation of the energy spectrum, which reduces the sensitivity of the measurement. The probability that the light curve is described by a constant flux level is > 10% in all nights except the first night where the probability of the fit is 0.8%. The av- erage statistical uncertainty of each flux measurement is ∼ 20%. Figure 14 shows the average flux of each night. The dashed line in the figure denotes the average flux from all nights and the shaded region shows the statistical error in the flux. The average integral mean flux F >200 GeV is F >200 GeV = (1.96 ± 0.05 stat ) × 10 −10 cm −2 s −1 . (6) There is a probability of 67% that the measured daily flux values are compatible with a constant flux. We can, therefore, conclude that the reconstructed flux of the Crab Nebula, within statistical uncertainties, was constant over the entire observation period.

Morphology of the γ-Ray emitting Region
The morphology of the γ-ray emission was studied by generating sky-maps in three uncorrelated bins of Size. The reconstruction of the origin of a γ-ray event with a single telescope is possible with the so-called Dispmethod (Fomin et al. 1994;Lessard et al. 2001). For the studies presented here we used the following parameterization for Disp: where a and b are second order polynomials found by fitting Monte Carlo simulated γ-ray showers (Domingo-Santamaría et al. 2005). Strong tail cuts of 10 and 8 photoelectrons were used in the image cleaning for core-and boundary-pixels respectively and a tight Hadronness-cut < 0.1 was applied resulting in improved angular resolution. The reconstructed event origins were corrected for possible mispointing by using the information from the starguider camera. 2D-histograms with bin sizes of (0.057 • × 0.057 • ) were filled with the corrected event origins (events with energies < 500 GeV). A four times finer binning was chosen for the sky-map filled by events with energies above 500 GeV. Figure 15 shows the background-subtracted sky-maps of excess events from the Crab Nebula for γ-ray energies ∼ 160 GeV, ∼ 250 GeV and > 500 GeV.

Center of Gravity of the γ-Ray Emission
The CoG of the γ-ray emission was derived from the sky-maps of the excess events shown in Figure 15 by fitting them with a 2D-Gaussian of the form where F res is introduced to account for a possible constant offset of the background subtracted sky-map. In this representation σ defines the 39% containment radius of the observed γ-ray emission. Here we assume that the distribution of excess events is rotationally symmetric, i.e. σ x = σ y = σ. It is further assumed that σ is the convolution of the response of the detector σ psf (point spread function) and the apparent size of the γ-ray emission region σ src , i.e. σ 2 = σ 2 psf + σ 2 src . Note that here it is assumed that the γ-ray emission region follows the shape described by Equation 8, which in reality is not necessarily the case.
The CoGs obtained from the fittedx andȳ are listed in Table 6 and shown in Figure 16 superimposed on the composite image of optical, IR and X-ray observations of the Crab Nebula. The three measured CoGs are compatible among each other and coincide with the position of the pulsar. Note that the systematic uncertainty of the position is ∼ 1 ′ .

TABLE 6
Center of gravity of the γ-ray emission of the Crab nebula obtained for different energies. In the first column the applied Size-cut is stated. The second column shows the corresponding range of γ-ray energies covered (peak value and full width at half maximum of the distribution of Monte Carlo γ-ray events). The last two columns give the fitted position of the CoG and the statistical uncertainty.

Extension of the γ-Ray Emission Region
The extension of the γ-ray emission region was studied by comparing the width of the excess events distribution with that obtained for a simulated γ-ray pointsource. The simulated distributions were verified by comparing them with data from Mrk421, an extragalactic γray source that can be considered a point source for our The position of the Crab pulsar is marked with a black star. The Chandra X-ray image is shown in light blue, the Hubble Space Telescope optical images are in green and dark blue, and the Spitzer Space Telescope's infrared image is in red (picture from Chandra 2006). Overlaid are the center of gravity of the γ-ray emission at different energies (simple cross > 500 GeV; dot ∼ 250 GeV; square ∼ 160 GeV). The error bars indicate the statistical uncertainty in the position of the CoG. The upper limits (95% confidence level) on the 39% containment radius of the γ-ray emission region that were derived from the θ 2 -distributions are indicated with circles (dashed > 500 GeV; solid ∼ 250 GeV).
purpose. The width of the γ-ray excess extracted from Mrk421 and the simulated width for a point source agree within statistical uncertainties.
In the following, the average position of the CoGs obtained from the three sky-maps is assumed as the γ-ray source position. The angular distance squared (θ 2 ) between the reconstructed origin and the assumed source position is calculated for every event. The background subtracted θ 2 -distributions obtained for the three energy ranges are shown in Figure 17. Data (black) and MC (blue) are compatible within statistical uncertainties in all three θ 2 -distributions.
An exponential function of the form describes the expected θ 2 -distribution, where a is a normalization and σ is the same as in equation 8. Values for σ 2 and σ 2 psf found by fitting the corresponding θ 2distributions with equation 9 are shown in Table 7. Upper limits on σ src were calculated with a confidence level of 95% following the procedure outlined in Yao et al. (2006) for one-sided confidence intervals and Gaussian errors. The results are presented in Table 7, and, in addition, for energies ∼ 250 GeV and > 500 GeV they are shown in Figure 16. The limits obtained for γ-ray energies above 500 GeV and about 250 GeV constrain the γ-ray emission to a region within the optical synchrotron nebula. Search for pulsed γ-Ray Emission Among the most challenging tasks of ground-based γray experiments is the detection of a pulsar. Several detectors have tried but failed. Currently MAGIC is the only ground-based detector with threshold settings below 100 GeV that is appropriate for a search of pulsed γ-ray emission from the Crab pulsar. For the data a periodicity analysis was performed after γ/hadron-separation and selection of events with small |Alpha|-value. The cuts were chosen by Monte Carlo simulations to optimize the sensitivity of the analysis. After event selection, the event times 25 were transformed to the barycenter of the solar system with the TEMPO timing package (Taylor et al. 2000). Then, the corresponding phase φ j of the Crab pulsar was calculated for each transformed arrival time t j :   Table 8). We tested for periodicity with the H-test (de Jager et al. 1989), the Pearson's χ 2 -test and a test from Gregory & Loredo (1992) that is based on Bayesian statistics.  The analysis chain was tested with optical observations of the Crab pulsar by the MAGIC telescope. Within this 12.5 hours' observation, every time the readout of MAGIC was triggered by a cosmic ray shower, the signal of the pixel in the center of the MAGIC camera was recorded by the MAGIC DAQ for 100 ns. Along with an average trigger rate of 200 Hz, the effective observation time was only about one second. Figure 18 shows the reconstructed optical light curve of the Crab pulsar with the familiar main-and inter-pulse. For better readability the light curve is shown twice. The position of the main-pulse is shifted with respect to the position of the main-pulse in radio by (−252 ± 64µs), which is in agreement with the contemporary measurement of Oosterbroek et al. (2006).

Search for pulsed Emission in Differential Bins of Energy
We searched for pulsed γ-ray emission in five bins of reconstructed energy between 60 GeV and 9 TeV. This search was motivated by a possible pulsed γ-ray component at TeV energies (Hirotani 2001(Hirotani , 2007. However, no signature of periodicity was found in any of the tested energy intervals. For each energy bin an upper limit on the number of excess events was calculated with a confidence level of 95% in two different ways, first from the result of the H-test as described by de Jager (1994) and second from the pulse phase profile. In the calculation of the limit from the result of the H-test, it is assumed that the duty cycle (full width at half maximum) of the pulsed γ-emission is 20%, similar to the duty cycle of the light curve measured by the EGRET detector above 100 MeV (Fierro et al. 1998). No assumption about the position of the emission in the pulse phase profile enters into the calculation. This additional constraint is applied, however, when the upper limit is directly derived from the pulse phase profile. As signal regions we chose the phase intervals where EGRET had observed pulsed emission above 100 MeV, i.e. −0.06 . . . 0.04 and 0.32 . . . 0.43 (shaded region in Figure 20). The background was estimated from the remaining phase intervals. Having defined the signal and background regions in this way, the upper limit on the number of excess events was obtained by the method of Rolke et al. (2005), assuming a systematic uncertainty of 30% on the collection area. Because of the additional constraint made about the position of the expected pulsed emission, the limits obtained from the pulse phase profile are on average about a factor of two better than the limits obtained from the result of the H-test.
The upper limits derived from the pulse phase profiles were converted into flux limits. The collection area was calculated assuming a photon index of 2.6 for the γ-ray spectrum. The flux limits are shown in Figure 19 together with the upper limit on the cutoff energy, which is derived in the following section.

Upper Limit on the Cutoff Energy of the pulsed Emission
Apart from the search for pulsed emission in bins of reconstructed energy, we performed a periodicity analysis, Whipple (Lessard 2000) Outer Gap (Hirotani 2007) EGRET (Fierro 1998) and 27 GeV exponential cutoff EGRET (Fierro 1998) Fig. 19.-Upper limits on the pulsed gamma ray flux from the Crab pulsar; upper limits in differential bins of energy are given by the blue points. The upper limit on the cutoff energy of the pulsed emission is indicated by the dashed line. The analysis threshold to derive the upper limit on the cutoff energy is indicated by the red arrow.

MeV
this time selecting events with Size< 300 phe (γ-ray energies 180 GeV) and applying the same optimized Size dependent Hadronness-cuts and Alpha-cuts as above. Compared to the previously described analysis, this one is optimized for a search of pulsed emission close to the threshold of the experiment. The analysis threshold, defined as the peak of the energy-distribution of simulated γ-ray showers, is 60 GeV. Figure 20 shows the pulse phase profile obtained for the selected events. The result of a Pearson's χ 2 -test is 13.1 with 10 degrees of freedom, corresponding to a significance of 1.2 σ for periodic emission. The result of the H-test is 3.9, which is equivalent to a significance of 1.3 σ. The test by Gregory & Loredo (1992) results in a probability of 4.1 · 10 −4 that pulsed emission is present in the data. These tests do not make an assumption about the position of the pulsed emission in the pulse phase profile. However, an excess is evident in the pulse phase profile at the position of the inter-pulse in the same phase range where EGRET detected pulsed emission above 100 MeV. If the EGRET phase regions (shaded regions in the pulse phase profile) are defined as the signal region, and the remaining phase intervals as background region, the significance of the excess of 2.9 σ can be calculated. Note that the pulse phase profile is shown twice for better visibility. When increasing the upper Size-cut, the number of excess events in the phase interval between 0.32 and 0.43 monotonically increases and saturates for upper Sizecuts of 200 phe and above. This behavior is expected for a γ-ray source with a power-law energy spectrum that is attenuated by an exponential cutoff.
The observed excess is not sufficient to claim the detection of a pulsed signal, therefore, upper limits on the number of excess events were calculated with a confidence level of 95% (see Table 9). Note that because of the observed excess, the upper limit from the pulse phase profile is larger than the limit obtained from the H-test.
Using the different limits on the number of pulsed excess events we constrain, in the following, the cutoff energy of the pulsar spectrum under the assumption that the break in the energy spectrum can be described with an exponential cutoff. In the procedure we use the parametrization F (E) <10 GeV of the measured pulsar spectrum below 10 GeV (Fierro et al. 1998), extended with an exponential cutoff: The spectrum with a given cutoff energy E Cutoff is convoluted with the effective collection area after cuts. The collection area is derived from MC simulations, assuming the same γ-ray spectrum. The number of expected excess events for the assumed cutoff energy is obtained by multiplying the convoluted spectrum with the observation time. In an iterative algorithm E Cutoff is changed until the number of expected excess events matches the upper limit on the number of excess events. In this way we derive an upper limit on the cutoff energy of 27 GeV from the result of the H-test and 34 GeV from the limit obtained from the pulse phase profile respectively. Differential and integral upper limits on the flux were calculated and are shown in Table 9.

DISCUSSION
In this paper we report on the most detailed study to date of VHE γ-ray emission of the Crab Nebula below 500 GeV. This study includes • a measurement of the differential energy spectrum down to 60 GeV • a determination of the inverse Compton peak of the VHE γ-ray emission • the search for an extended source morphology • the calculation of a light curve of the VHE γ-ray emission from the nebula at 200 GeV • a search for pulsed emission from the Crab pulsar in differential bins of energy and in an optimized low energy analysis Most of the aforementioned studies were done in this energy region for the first time, they were possible only because the imaging air shower Cherenkov technique was used. The wavefront sampling technique, which, up to now, was the only experimental technique used in this energy domain, did, at best, allow one to arrange an integral flux measurement and to search for pulsed emission. The performance of MAGIC is superior even in the energy domain below 200 GeV, where a progressive degradation of the γ/hadron-separation power is observed.
In this context studies for improving the suppression of background events by exploiting the intrinsic time structure of the recorded PMT-signals are ongoing. Further improvement is expected with the second MAGIC telescope currently under construction. The measured energy spectrum of the Crab Nebula in Figure 10 extends over two decades in energy and five decades in flux. The spectral shape deviates from a pure power-law behavior and is, within experimental uncertainties, in agreement with a curved powerlaw. The spectrum clearly softens with increasing energy. The observation supports the generally accepted picture that the steady emission above a few tens of GeV and up to the highest measured energies can be described within the framework of the SSC-model (Gould 1965;de Jager & Harding 1992;Aharonian et al. 2004). Also the peak position of the inverse Compton emission of the tested predictions is in agreement with the determined one (77 ± 47 stat +107 −46 syst GeV). At GeV energies EGRET observed a γ-ray flux, which was a factor five above the flux predicted by the SSCmechanism (de Jager et al. 1996). Atoyan & Aharonian (1996) explain this γ-ray excess by an additional γ-ray component from bremsstrahlung of electrons that are partially captured in filaments of the nebula. Such an extra component can significantly change the spectral slope at several hundred GeV compared to a pure ICscenario (cf. blue and red curves in Figure 12), and results in an almost pure power-law behavior of the energy spectrum between ∼ 100 GeV and 10 TeV (constant Γ). At 300 GeV, where the measurement is most sensitive, the measured slope (black curve and data point) is considerably harder than predicted by Atoyan & Aharonian (1996). It is, therefore, unlikely that the γ-ray excess at GeV energies can be explained by bremsstrahlung as proposed. Later predictions by Aharonian & Atoyan (1998) that also include the mentioned bremsstrahlung mechanism, are in agreement with the presented measurement (red curve). However, all above mentioned predictions agree with the measurement if the measurement is shifted by 0.2 to negative values, which is within the range of the systematic uncertainty of the measurement. In the prediction by Aharonian & Atoyan (1998) also a γ-ray component from π 0 -decay is included, which results in a considerable harder spectrum above a few TeV compared to the pure IC-scenario (cf. Figure 12). However, given the limited statistics above 1 TeV one cannot exclude any such prediction from the measurement.
Studies about the morphology of the γ-ray emitting region of the Crab Nebula have been performed by Aharonian et al. (2000Aharonian et al. ( , 2004 for γ-ray energies above 1 TeV. In both cases it was found that within the resolution of the experiment the emission region is pointlike. They placed an upper limit on the source size of ∼ 2 ′ at energies between 1 and 10 TeV. In the VHE domain, the morphology of the emission region has not yet been studied at energies below 1 TeV. With the resolution of MAGIC it was possible to constrain the as the optical synchrotron nebula (see Figure 16). The upper limit on the size of the emission region is ∼ 2 ′ , which is about four times larger than the predicted size of the inverse Compton surface brightness for γ-ray energies below 500 GeV (de Jager & Harding 1992).
X-ray observations indicate variabilities in the acceleration and cooling times of electrons on timescale of months (e.g. Hester et al. 2002). However, variations in γ-rays could not be detected so far. The sensitivity of MAGIC allowed us to study the variability of the γ-ray emission above 200 GeV on timescales as short as a few minutes up to months. We measured a flux that is within statistics compatible with steady emission. During the 2.5 · 10 −11 7.9 · 10 −11 2σ Differential flux limit at 60 GeV [cm −2 s −1 GeV −1 ] 4.5 · 10 −12 8.9 · 10 −12 Peak Energy MC [GeV] 60 observation the stability of the integral flux was better than 10% on all tested time scales. In a search for pulsed VHEγ-ray emission with MC optimized cuts in Hadronness and Alpha an excess was found in the pulse phase profile at the same position were EGRET detected pulsed emission above 100 MeV. The significance of the excess is 2.9 σ if the phase regions where EGRET detected pulsed emission were chosen as signal regions and the remaining phase intervals are considered as background regions. The similarity of the distribution of excess events in the EGRET > 5 GeV and MAGIC data and the monotonic increase of the number of excess events with increasing upper Size-cut are strong indications that the observed excess is not a random fluctuation. However, a detection can not be claimed as the conventional 5 σ detection limit has not been reached.
With the result of the H-test an upper limit on the cutoff energy of 27 GeV was derived, assuming that the power-law spectrum of the pulsar at GeV energies is attenuated by an exponential cutoff. However, if the cutoff of the spectrum has a super-exponential shape a cutoff energy almost as high as the analysis threshold (∼ 60 GeV) cannot be excluded.
Also, no pulsed emission was detected for energies above 100 GeV, which could have its origin in ICupscattering of IR photons. Despite earlier claims of a strong component (Hirotani & Shibata 2001), latest models (Hirotani 2007) seem to disfavor a pulsed TeV component from the Crab pulsar due to dominant γγ absorption processes. In future, detailed spectroscopic studies of the pulsed emission by e.g. GLAST and ground-based experiments with lower thresholds and higher sensitivities like e.g. MAGIC II (under construction) and CTA (projected) will hopefully resolve the long standing question of the origin of the pulsed emission.

ACKNOWLEDGMENTS
We are grateful for discussions with Kouichi Hirotani. We also would like to thank the IAC for the excellent working conditions at the ORM in La Palma. The support of the German BMBF and MPG, the Italian INFN, the Spanish CICYT, ETH research grant TH 34/04 3, and the Polish MNiI grant 1P03D01028 are gratefully acknowledged.