Observations of the Crab pulsar between 25 GeV and 100 GeV with the MAGIC I telescope

We report on the observation of $\gamma$-rays above 25\,GeV from the Crab pulsar (PSR B0532+21) using the MAGIC I telescope. Two data sets from observations during the winter period 2007/2008 and 2008/2009 are used. In order to discuss the spectral shape from 100\,MeV to 100\,GeV, one year of public {\it Fermi} Large Area Telescope ({\it Fermi}-LAT) data are also analyzed to complement the MAGIC data. The extrapolation of the exponential cutoff spectrum determined with the Fermi-LAT data is inconsistent with MAGIC measurements, which requires a modification of the standard pulsar emission models. In the energy region between 25 and 100\,GeV, the emission in the P1 phase (from -0.06 to 0.04, location of the main pulse) and the P2 phase (from 0.32 to 0.43, location of the interpulse) can be described by power laws with spectral indices of $-3.1 \pm 1.0_{stat} \pm 0.3_{syst}$ and $-3.5 \pm 0.5_{stat} \pm 0.3_{syst}$, respectively. Assuming an asymmetric Lorentzian for the pulse shape, the peak positions of the main pulse and the interpulse are estimated to be at phases $-0.009 \pm 0.007$ and $0.393 \pm 0.009$, while the full widths at half maximum are $0.025 \pm 0.008$ and $0.053 \pm 0.015$, respectively.


INTRODUCTION
The pulsar B0531+21, also commonly known as the Crab pulsar, is the compact object left over after a historic supernova explosion that occurred in the year 1054 AD. Its energetic pulsar wind creates a pulsar wind nebula, the Crab Nebula, an emitter of strong and steady radiation. The pulsar and pulsar wind nebula have been observed and studied in almost the entire accessible electromagnetic spectrum from about 10 −5 eV (radio emission) to nearly 100 TeV (very high energy γ -rays; e.g., Aharonian et al. 2004). The nebular emission is commonly used as a standard candle for astronomy in various energy ranges. Recently, γ -ray flares from the Crab Nebula were discovered in the GeV range (Tavani et al. 2011;Abdo et al. 2011) and a hint of increased flux in the TeV range during a GeV flare was also reported (ATel 2921).
The Crab pulsar and several other pulsars are amongst the brightest known sources at 1 GeV. However, a spectral steepening made their detection above 10 GeV elusive despite numerous efforts (e.g., Aharonian et al. 2007;de Naurois et al. 2002;Lessard et al. 2000). The energy thresholds of imaging atmospheric Cherenkov telescopes (IACTs) were, in general, too high, while the γ -ray collection area of satellite-borne detectors was too small to detect pulsars above 10 GeV.
On the other hand, a precise measurement of the energy spectrum at and above the steepening is an important test for the standard pulsar models, such as the polar cap (PC), outer gap (OG), and slot gap (SG) models.
In the PC model, emission takes place within a few neutron star (NS) radii above a PC surface (Arons & Scharlemann 1979;Daugherty & Harding 1982, 1996. There, high-energy γ -rays above ∼1 GeV should be absorbed by a strong magnetic field (Baring 2004), which results in a very sharp cutoff (so-called super-exponential cutoff) in the energy spectrum.
Extending the original idea by Arons (1983), Muslimov & Harding (2003, 2004a, 2004b and Dyks et al. (2004) investigated the possibility of high-energy emission along the flaring field lines at high altitudes. This type of emission, an SG emission, can be observed at all viewing angle and for most cases emission from the two poles can be observed. The SG model predicts an exponential cutoff above 1 GeV (Harding et al. 2008). However, such geometrically thin emission models reproduce less than 20% of the observed γ -ray fluxes of the Crab and Vela pulsars (Hirotani 2008).
Seeking a different possibility for high-altitude emissions, Cheng et al. (1986aCheng et al. ( , 1986b proposed the OG model, hypothesizing that the emission zone is located in higher altitudes, beyond the so-called null-charge surface. Subsequently, Romani 29 Now at Ecole polytechnique fédérale de Lausanne (EPFL), Lausanne, Switzerland. 30  & Yadigaroglu (1995) and Romani (1996) developed the caustic model of OG emissions. A three-dimensional version of such a geometrical model of OG emissions was investigated by Cheng et al. (2000). Recently, Romani & Watters (2010) presented an atlas of pulse properties and proposed a method to discriminate different emission models from the geometrical point of view.
It is noteworthy that all the existing OG or SG models predict that the highest-energy photons are emitted via curvature radiation and that an exponential cutoff appears around 10 GeV in the spectrum of the Crab pulsar (e.g., Tang et al. 2008).
The MAGIC I telescope in its standard trigger mode has the worldwide lowest threshold of all currently operating IACTs, around 60 GeV. A previous study of MAGIC data above 60 GeV revealed a 2.9σ excess from the Crab pulsar (Albert et al. 2008c). Following up on this hint, we investigated an alternative trigger concept, the sum trigger (see Section 2.2), which lowered the energy threshold of MAGIC to about 25 GeV. Using this new trigger, we observed the Crab pulsar between 2007 October and 2008 February and detected high-energy γ -ray emission from the Crab pulsar with a significance of 6.4σ (Aliu et al. 2008). This detection suggests the distance of the emission region from the stellar surface to be larger than 6.2 ± 0.2 stat ± 0.3 syst times the stellar radius, which ruled out the PC model as a viable explanation of the observed emission. This initial detection was briefly reported in Aliu et al. (2008). In winter 2008/2009, the Crab pulsar was observed again with MAGIC using the sum trigger.
In 2008 August, the new satellite-borne γ -ray detector with 1 m 2 collection area, the Fermi Large Area Telescope (Fermi-LAT), became operational and measured the spectra of γ -ray pulsars up to a few tens of GeV (Abdo et al. 2010c). All the energy spectra are consistent with a power law with an exponential cutoff, though statistical uncertainties above 10 GeV are rather large. 34 The cutoff energies are typically between 1 GeV and 4 GeV. These Fermi-LAT measurements also disfavor the PC model and support the OG and the SG model (Abdo et al. 2010c).
However, the cutoff energy of the Crab pulsar determined with the Fermi-LAT under the exponential cutoff assumption is ∼6 GeV, an unlikely value for the signal above 25 GeV detected by MAGIC. In order to verify the exponential cutoff spectrum, a precise comparison of the energy spectra measured by the two instruments is needed. The recent detection of the Crab pulsar above 100 GeV by the VERITAS Collaboration has shown that indeed the energy spectrum above the break is not consistent with an exponential cutoff but that it is better described by a broken power law (Aliu et al. 2011). It is, however, not clear whether the spectrum continues as a power law after the break or there is another component above 100 GeV in addition to the exponential cutoff spectrum because of missing flux measurements in the intermediate energy range from 25 GeV to 100 GeV.
The main objectives of this paper are the evaluation of the exponential cutoff spectrum of the Crab pulsar with the MAGIC data and the presentation of its energy spectrum between 25 GeV and 100 GeV. We also give details of the MAGIC observations, the data selection, the analysis, and physics results. We report, for the first time, separate energy spectra and a pulse profile analysis above 25 GeV for both the main pulse and the interpulse. The large majority of the results presented in this paper are extracted from the PhD thesis of Takayuki Saito (Saito 2010). The paper has the following structure. After describing the MAGIC telescope and the sum trigger in Section 2, we present the observation and details of the data processing in Section 3. The detection of pulsed emission is described in Section 4. Based on the MAGIC detection and the Fermi-LAT measurements, the evaluation of the exponential cutoff assumption is reported in Section 5. Energy spectra for the main pulse, the interpulse, and the summed pulsed emission are presented in Section 6, followed by a discussion of the pulse profile in Section 7. We conclude the paper in Section 8 together with a theoretical interpretation of the spectrum and an outlook on what can be expected in the near future.

System Overview
The MAGIC I telescope is a new generation IACT located on the Canary island of La Palma (27. • 8 N, 17. • 8 W, 2225 m asl). Some of the Cherenkov photons emitted by charged particles in an air shower are collected by a parabolic reflector with a diameter of 17 m and focused onto a fine pixelized camera, providing an image of the air shower. The camera comprises 577 photomultipliers (PMTs) and has a field of view of ∼3. • 6 diameter.
The fast analog PMT signals are transported via optical fibers to the counting house where the signals are processed by the trigger system and recorded by the data acquisition (DAQ) system. The trigger is normally derived from the pixels in the innermost camera area (the trigger area) of around 1 • radius (325 pixels). Each signal from the pixels in the trigger area is amplified and split into two signals equal in amplitude. The signals are routed to the trigger system and to the DAQ system, respectively. Signals from the non-trigger pixels enter the DAQ system directly. The trigger criteria of the standard digital trigger (hereafter called standard trigger) system are applied in two steps: each optical signal from the trigger area is converted to an electric one and examined by a discriminator with a computercontrolled threshold level; the threshold level is typically 6-7 photoelectrons (PhEs). The digital signals are then processed by a topological pattern logic, which searches for a close-packed group of four compact-next-neighbor pixels firing within a time window of ∼ 6 ns.
In the standard trigger concept, only signals above the preset threshold in four compact-next-neighbor pixels can generate a trigger, while signals below the threshold or pixels not situated closely together cannot contribute to the decision. This deficiency is particularly pronounced in the case of shower images of a light content close to the threshold, i.e., in the interesting energy region below 60 GeV. For this reason, a new trigger system called sum trigger has been developed to explore the energy region down to ∼25 GeV. Details of the new sum trigger will be presented in the Section 2.2 below.
The signals entering the DAQ system are recorded by flash analog-to-digital converters (FADCs) with a sampling rate of 2 GHz. For each event, 50 FADC slices are recorded for all the pixels. The details of the DAQ system are described in Goebel et al. (2007). For the pulsar study, the central pixel of the camera was modified to record the optical flux from the object under study, i.e., to measure the optical pulsations of the Crab pulsar. The details of the central pixel system can be found in Lucarelli et al. (2008). The telescope tracked the Crab position with a typical precision of 0. • 02. In addition, we regularly recorded calibration and pedestal events with a frequency of 25 Hz each. Further details on the telescope can be found in Baixeras et al. (2004).

The Sum Trigger
As mentioned above, the standard trigger scheme is not very efficient below 60 GeV because even at 25 GeV the image covers well over four pixels and the signals show a wide spread in amplitude. In order to improve the trigger efficiency just above threshold, a new technique was developed, the so-called sum-trigger method. The main feature of the sum trigger is the summation of the analog pixel signals from a wider camera area, so-called patches, followed by a discrimination of this summed-up signal. There are 24 partially overlapping patches in an annulus with inner and outer camera radii of ∼0. • 25 and ∼0. • 8, respectively. Each patch comprises 18 pixels. The threshold level for the summed signal from a patch is an amplitude of 27 PhEs.
This sum-trigger concept has some clear advantages compared to the standard trigger. The summation of the analog signals allows any pixel signal in the patch to contribute to the trigger, even if its amplitude is below the pixel threshold of the standard trigger. The concept, however, has the disadvantage of being quite sensitive to accidental triggers from afterpulses. Afterpulses are caused by PhEs hitting the dynodes and sometimes releasing ions from adsorbed water or adsorbed gases. These ions are back-accelerated by the relatively high voltage between the photocathode and the first dynode, hit the photocathode, and liberate many secondary electrons. The afterpulse amplitude spectrum follows basically an exponential distribution, which drops significantly slower than the Poissonian night sky light distribution and completely dominates the rate of signals above 5-6 PhEs. In a single patch, i.e., from the sum of 18 pixels, the rate of afterpulses above 27 PhEs was found to be around 20-30 kHz, which is far beyond the MAGIC DAQ rate limit of 1 kHz. To suppress this undesirable background, individual pixel signals were, before summing, clipped at an amplitude of ∼6 PhEs, thus making the trigger insensitive to large afterpulses of individual PMTs. The sum-trigger area (the annulus with inner and outer radii of ∼0. • 25 and ∼0. • 8), the patch size (18 pixels), the threshold level (27 PhEs in amplitude after sum), and the clipping level (∼6 PhEs) were optimized by measurements and detailed Monte Carlo (MC) studies (Rissi 2009). The chosen settings for the sum trigger result in a trigger threshold of 25 GeV for a γ -ray source with an index of −2.6 as shown in the top panel of Figure 1 (according to the convention in ground-based γ -ray astronomy the threshold is defined as the peak of the reconstructed differential energy spectrum). The sum trigger also improved significantly the collection area for low-energy showers when compared to the area of the standard trigger. At 20 GeV the collection area of the sum trigger is 10 times larger and at 60 GeV is still twice as large when compared to that of the standard trigger, as shown in the bottom panel of Figure 1. A more detailed description of the sum trigger is presented in Rissi (2009) and Rissi et al. (2009).

Observation
The first observation of the Crab pulsar with the sum trigger started on 2007 October 21 and extended up to 2009 February 3. In total the Crab pulsar was observed for 48 hr in winter 2007/2008 (Aliu et al. 2008) and for 78 hr in winter 2008/2009. In the first year, all observations were restricted to zenith angles below 20 • where the air mass between the showers and the telescope is the lowest possible, i.e., the atmospheric transmission for Cherenkov light is highest. In this zenith angle range the correlation between energy and observed number of PhEs is almost independent of the zenith angle and the trigger threshold is nearly constant as a function of the zenith angle. In the second year, some of the observations were done at zenith angles above 20 • . These data are not used in the following analysis. In the second campaign, five sub-patches malfunctioned. The losses are estimated and corrected using MC simulations.

Data Processing
In the calibration process, the conversion factors from the FADC counts to the number of PhEs and the relative timing offsets of all pixels are computed using the calibration and pedestal events. The details of the procedure can be found in Albert et al. (2008a). After the calibration, an image cleaning is performed in order to remove pixels which do not contain a useful Cherenkov photon signal (e.g., pixels only containing FADC pedestal, night sky background (NSB) photons, and afterpulses). The standard procedure for the image cleaning can be found in Aliu et al. (2009). Since this study aims for the lowest possible threshold, a more sophisticated method of image cleaning was used in this analysis. At first the algorithm searches for the core pixels of the shower image. The definitions of the core pixels are as follows.
1. If two neighboring pixels have more than 4.7 PhEs each and the time difference is less than 0.8 ns, these two pixels are core pixels. 2. If three neighboring pixels have more than 2.7 PhEs each and the arrival times of all three are within 0.8 ns, these three pixels are core pixels. 3. If four neighboring pixels have more than 2.0 PhEs and the arrival times of all four are within 1.5 ns, these four pixels are core pixels.
After the core search, boundary pixels of the image are selected. For a pixel to be defined as a boundary pixel, the following three conditions must be fulfilled.
1. The pixel must be a neighbor to at least one of the core pixels. 2. The pixel must have more than 1.4 PhEs. 3. The time difference to at least one of the neighboring core pixels must be less than 1.0 ns.
The charge and timing information of pixels which are neither "core" nor "boundary" is discarded. In this way, accidental trigger events (half of the recorded data) are efficiently rejected and the contamination of, e.g., NSB photons to the shower image can be mostly suppressed. Further details of this method can be found in Shayduk et al. (2005).
After the image cleaning, the conventional image parameters are calculated in the standard way (Hillas 1985;Aliu et al. 2009).

Data Pre-selection
Only data taken under stable atmospheric and hardware conditions were used in the analysis. Selection criteria include, for instance, the performance of the mirror focusing and reflection, the cloud coverage, and the stability of the event rate after image cleaning. A detailed description of the pre-selection criteria can be found in Saito (2010). As mentioned before, only data with zenith angles below 20 • were used in the analysis. After applying all criteria 25 hr remained from the winter 2007/2008 period and 34 hr from the winter 2008/2009 period.

Event Selection
A cleaned image of a 25 GeV γ -ray on average has ∼8 pixels, which is barely sufficient to perform a moment analysis to obtain the conventional Hillas image parameterization of the shower. However, the only effective way to separate γ -rays from the background below 100 GeV exploits the orientation of the images in the camera plane. γ -rays from the source have images that point with their major axis toward the source location. Background events (mainly hadron showers and muon arcs/rings), on the other hand, have images that are randomly oriented. The parameter that describes the orientation of an image in the camera is called ALPHA, which is the angle between the major axis of an image and the line connecting the center of gravity (COG) of that image with the location of the source in the camera. The ALPHA distribution of the γ -ray MC events as a function of SIZE is illustrated in Figure 2, where SIZE is the total number of PhEs in the image. Red stars indicate the cut values which maximize the so-called quality factor Q defined as where γ and BG are the survival efficiencies of the γ -ray events and hadron background events, respectively. In the analysis we used an ALPHA cut depending on SIZE, which was derived by fitting a numerical function a(log 10 (SIZE) + b) c to the best cut values found in the individual SIZE bins (stars in Figure 2). γ , BG , and Q for the SIZEdependent ALPHA cut used are shown in Figure 3.
In addition, the so-called spark-like events, which are created by the discharge of charge accumulated at the glass envelope of some PMTs, are rejected. These events are efficiently identified by applying the condition 4.0 × log 10 (CONC) > 1.5 − log 10 (SIZE), where CONC is defined as the sum of the charges in the two pixels with the highest content divided by SIZE, which indicates how strongly the charge is concentrated in a small region.

Energy Estimates
The energy of each event is estimated by means of the Random Forest method. After the training with MC γ -ray events, the  Random Forest assigns the most probable energy to each event by using several image parameters in a comprehensive manner. The details of the method can be found in Albert et al. (2008b). In this analysis SIZE, LENGTH, DIST, and the zenith angle were used for the Random Forest energy estimation. LENGTH is the second moment of the charge distribution along the major axis of the shower image, while DIST is the distance between the source position in the camera and the COG of the image. WIDTH (the second moment along the minor axis of the image) and CONC, which are normally included, were not used because it was found that they do not contribute to the energy estimate in the very low energy regime.

Pulse-phase Calculation
Each event is marked with a time stamp, which gives the time when the event was triggered. The time stamps are derived from a GPS-controlled rubidium clock and have an absolute accuracy of less than 1 μs. In order to compensate the varying propagation times of the γ -rays within the solar system, which are mainly due to Earth's movement around the Sun, the recorded times are transformed to the barycenter of the solar system. The barycentric correction was done with the software package TEMPO [24]. The rotation frequency ν 0 , its time derivativeν 0 , and the barycentric times of the main pulse peak of the Crab pulsar are monitored in radio at 610 MHz by the Jodrell Bank radio telescope (Lyne et al. 1993) and the values of the parameters are published every month. 35 Based on these, the rotational phase of each event is computed from the barycentric times with the following formula: The second and higher derivative terms of this Taylor series are negligible on the scale of one month.

Optical Pulsation
In order to assure that the timing system of MAGIC and the pulse-phase calculation worked properly, the optical pulsation of the Crab pulsar was checked first. The optical pulsation was measured with the central pixel, which was modified for this purpose to be sensitive and to digitize the light flux variation from the source. Every time a shower event was triggered, the signal of the central pixel was recorded by the DAQ for 25 ns.
The phase distribution (hereafter pulse profile) of the central pixel data is shown in Figure 4. Two peaks are clearly visible at the expected phases. Phase 0 corresponds to the main peak position in radio at 610 MHz. A delay of ∼0.01 in phase can be seen with respect to the radio main peak position, which is known and consistent with other measurements (see, e.g., Oosterbroek et al. 2008).

Pulsation above 25 GeV
The pulse profile of the γ -ray events detected with MAGIC is shown in Figure 5. Events with SIZE below 25 PhEs and above 500 PhEs were discarded. Note that every event is shown twice (three times for the first bin) in order to generate a pulse profile that spans the phase region from −45/44 (−1.0227) to 45/44 (1.0227). The bin width is 1/22 in phase, which corresponds to about 1.5 ms. An excess is evident in the profile at the position of the main pulse and interpulse of the pulsar. Following the often-used convention (Fierro et al. 1998) of P1 (main pulse phase from −0.06 to 0.04) and P2 (interpulse phases from 0.32 to 0.43), the numbers of excess events in P1 and P2 are 6200 ± 1400 (4.3σ ) and 11300 ± 1500 (7.4σ ), respectively. By summing up P1 and P2, the excess corresponds to 7.5σ . The background level was estimated using the so-called off pulse (OP) phases (0.52-0.88; Fierro et al. 1998). Above 25 GeV, the flux of P2 is nearly twice that of P1. The width of the main pulse is significantly smaller than the conventional P1 phase interval. By defining the signal phase of the main pulse to be −1/44 to 1/44 in phase (the bin with the largest number of events), the excess is 6400 ± 970 corresponding to 7.0σ . No significant emission between the main pulse and the interpulse was detected. A detailed study of the pulse profile is given in Section 7.
In order to verify the soundness of the signal, further tests were made. First, the phase distribution of the events which are rejected by the ALPHA cut is examined. The results are shown in Figure 6. As expected, the distribution is consistent with statistical fluctuations of the background without any signal. Also, the growth of the number of excess events as a function of the number of background events is checked. As one can see in Figure 7, the excess grows linearly, assuring that the signal is constantly detected.

Variability Study
The linear growth of the excess shown in Figure 7 implies a constant flux of the pulsed signal. Nevertheless, we also applied the χ 2 method to test for a possible yearly variability. The number of excess events as a function of SIZE is compared between the two years in Figure 8. The difference of observation time and the effect of the malfunctioning of the sum-trigger sub-patches are corrected for the second year data. Using MC simulations, the sub-patch malfunction effects on the acceptance were estimated to be 21%, 17%, 11%, and 7% for the SIZE ranges of 25-50, 50-100, 100-200, and 200-400 PhEs, respectively. The χ 2 values for the comparison of the two years are 1.0 and 3.1 for P1 and P2, respectively, for 4 degrees of freedom. No significant yearly variability was found in the flux. We also studied a possible variability of the pulse profile. Figure 9 shows the pulse profiles for the first (winter 2007/2008) and the second (winter 2008/2009) year. The two profiles are compared with each other with a χ 2 test from phase −0.0682 to 0.432 (−3/44 to 19/44), which is roughly from the beginning of P1 to the end of P2. The χ 2 obtained is 5.0 for 10 degrees of freedom. Therefore, no significant yearly variability of the pulse profile was found.

EVALUATION OF THE EXPONENTIAL CUTOFF SPECTRUM SUGGESTED FROM THE FERMI-LAT DATA
It is very important for the verification of the standard OG and SG models to check if the energy spectrum follows an exponential cutoff. All energy spectra measured by Fermi-LAT up to a few tens of GeV are indeed consistent with the OG/SG model. In this section, we evaluate the exponential cutoff hypothesis based on the measurements performed by Fermi-LAT and MAGIC. performed a customized analysis of the specific phase intervals in order to properly compare the Fermi-LAT and MAGIC data. One year of Fermi-LAT data taken from 2008 August 4 to 2009 August 3 are analyzed. Events with an energy between 100 MeV and 300 GeV and with an arrival direction within a radius of 20 degrees around the Crab pulsar were downloaded from the Fermi Science Support Center. 36 Only events with the highest probability of being photons, those in the diffuse class, were used in this analysis. Events with imperfect spacecraft information and events taken when the satellite was in the South Atlantic Anomaly were rejected. In addition, a cut on the maximum zenith angle (<105 • ) was applied to reduce the contamination from the Earth-albedo γ -rays, which are produced by cosmic rays interacting with the upper atmosphere. The pulse-phase assignment to each event was carried out by the Fermi-LAT analysis tool with the monthly ephemeris information from Jodrell Bank. 37 For the detector response function, "P6_V3_Diffuse" is used, while "isotropic_iem_v02.txt" and "gll_iem_v02.fit" are used for the extragalactic and galactic diffuse emission models. In order to estimate the contribution of the Crab Nebula, the pulse-phase interval between 0.52 and 0.87 is used. 38 The contamination from nearby bright sources such as Geminga and IC 443 is also taken into account in the calculation of the spectrum.   Table 1 and Figure 12). Dotted lines are the expectation in the case where the energy scale of MAGIC is 23% higher than that of Fermi-LAT due to systematic uncertainties of both instruments.

Analysis of Public Fermi-LAT Data
(A color version of this figure is available in the online journal.) unbinned likelihood spectral analysis assuming a power-law spectrum with an exponential cutoff gives F 1 = (2.32 ± 0.05 stat ) × 10 −10 cm −2 s −1 MeV −1 , Γ 1 = 1.99 ± 0.02 stat , and E c = 6.1 ± 0.5 stat GeV as best-fit parameters, consistent with the values in Abdo et al. (2010a). Results for the total pulse are shown by the solid black line in Figure 10, together with the pulse profile above 100 MeV. The green curve in the figure is the spectrum given in Abdo et al. (2010a). The points are obtained by applying the same likelihood method in the limited energy intervals, assuming a power-law spectrum within each interval. The same analysis was applied to P1, P2, and the sum of the two (P1 + P2). The results are shown in Figure 12 and the best-fit parameters are summarized in Table 1.

Statistical Evaluation of the Difference between the Extrapolated Fermi-LAT Spectrum and the MAGIC Data
For MAGIC, the energy resolution below 50 GeV is about 40%. In addition, near the trigger threshold, the energy is overestimated because only events which deposit more Cherenkov photons onto the MAGIC mirrors are selectively triggered. These effects can be corrected with MC simulations, but the correction introduces additional systematic uncertainties. In order to minimize these uncertainties, we adopt the following method. Assuming that the exponential cutoff spectrum determined by the Fermi-LAT below a few tens of GeV is valid up to 2 TeV, we calculate the expected SIZE distribution in the MAGIC data. Taking into account the trigger threshold of 27 PhEs, the events with SIZE below 30 PhEs are not used to avoid a possible mismatch between MC and real data near the threshold. Then, we compute the χ 2 value between the expected and measured distributions. Statistical errors of E c measured by Fermi-LAT are taken into account as an error of the expected distribution.
The results are shown in Figure 11. The χ 2 values are 54.2, 15.8, and 42.3 for P1 + P2, P1, and P2, respectively. The number of degrees of freedom is 3. The exclusion probabilities correspond to 6.7σ , 3.0σ , and 5.8σ . It should be noted that the possible energy scale shift between the two instruments is not taken into account here.
The systematic uncertainty in the energy scale of the Fermi-LAT is estimated to be less than 7% above 1 GeV (Abdo et al. 2009) while that of MAGIC is estimated to be 16% (Albert et al. 2008c). We performed the same statistical test with an increased Fermi-LAT energy scale of 23%, corresponding to the linear sum of the systematic errors in both instruments. The results are shown as dotted lines in Figure 11. Though the discrepancies become smaller, the χ 2 values are 42.3, 12.6, and 30.0 with 3 degrees of freedom for P1 + P2, P1, and P2, respectively. The exclusion probabilities correspond to 5.8σ , 2.5σ , and 4.7σ . Even with the systematic uncertainties taken into account, the inconsistency between the extrapolated Fermi-LAT spectrum and the observations by MAGIC is significant.

Notes.
a Obtained by fitting Equation (4) to Fermi-LAT data. b Obtained by fitting Equation (5) to MAGIC data.

Spectra of P1, P2, and P1 + P2
The energy spectra of the Crab pulsar were computed based on the detected excess events found in P1 and P2, using the standard MAGIC software. The energy resolution and the trigger bias effect were corrected by an unfolding procedure which includes the Tikhonov regularization method (Tikhonov & Goncharsky 1987). The results are shown in Figure 12. The combined spectrum of P1 + P2 is consistent with a power law, which can be described using the following formula: F 30 = (14.9 ± 2.9 stat ± 9.6 syst ) × 10 −9 cm −2 s −1 TeV −1 and Γ 2 = 3.4 ± 0.5 stat ± 0.3 syst were obtained as best-fit parameters. Also, the spectrum seems to connect smoothly to the VERITAS measurements above 100 GeV (Aliu et al. 2011). The individual spectra of P1 and P2 were also calculated using the same data set. Results are also shown in Figure 12. The best-fit parameters are summarized in Table 1. Figure 12 clearly shows the deviation of the MAGIC spectra with respect to the extrapolation of the exponential cutoff spectra determined by Fermi-LAT, which is consistent with our statistical analysis in the previous section.

Combined Fit above 5 GeV
To get a better estimate of the power-law index for the higher energies, Fermi-LAT data points above 5 GeV and MAGIC data points are combined and fitted by a power law: It should be mentioned that Fermi-LAT points are obtained using the likelihood analysis for each energy interval assuming a power law, while for the fit each point was assumed to follow Gaussian statistics with the standard deviation being the error obtained by the likelihood results. Though this does not statistically correspond to the exact likelihood function of the problem, it is a good and appropriate approximation. The results are shown in Table 2.

Peak Phase and Pulse Width
We examine the peak phase and the pulse width in the MAGIC energy range assuming a pulse shape a priori and fitting it to the measured profile. The functions used are a Gaussian and a Lorentzian where μ corresponds to the peak phase, while σ can be translated into a pulse width. The full width at half-maximum (FWHM) is equal to 2.35×σ and 2×σ for the Gaussian and the Lorentzian, respectively. In order to study the asymmetry of the pulses, asymmetric Gaussian and Lorentzian functions are also assumed with σ being different below and above x = μ.
The results are shown in Figure 13 and Table 3. The peak phase of the main pulse is compatible with 0.0 (defined by the radio peak) for all fitted parameterizations. The Fermi-LAT Collaboration reported that the pulse shape is well modeled by an asymmetric Lorentzian and the peak phase above 100 MeV is −0.0085±0.0005 (Abdo et al. 2010a). The MAGIC result under the asymmetric Lorentzian assumption is also consistent with this. The peak phase of the interpulse depends on the assumption of the shape, while it is approximately 0.39 ± 0.01, which is also consistent with the value above 100 MeV (0.398 ± 0.003; Abdo et al. 2010a). The FWHM of the main pulse is approximately 0.03±0.01 independently of the assumed shape, while that of the interpulse is 0.07 ± 0.01 and 0.05 ± 0.01 for the assumptions of a Gaussian shaped and a Lorentzian shaped pulse, respectively. The main peak is narrower than the interpulse. The asymmetric assumptions imply that, for the main pulse, the rising edge is steeper than the falling edge, while the opposite is true for the interpulse. In Figure 14, the half-widths σ for the rising and falling edges of both the main pulse and the interpulse are compared with the values reported in Abdo et al. (2010a). The rising half of the main pulse become narrower as the energy increases. Though the uncertainty is larger, a similar tendency is also visible in the rising and falling halves of the interpulse, while no such energy dependence is visible in the falling half of the main pulse.

Other Emission Components
The AGILE Collaboration reported a possible third peak at phase between 0.65 and 0.8 above 100 MeV with a significance  of 3.7σ (Pellizzoni et al. 2009). A hint of a third peak is also seen in the Fermi-LAT data at phase ∼0.74 only above 10 GeV with significance of 2.3σ (Abdo et al. 2010a). They coincide with the radio peak observed between 4.7 and 8.4 GHz. In the MAGIC data above 25 GeV, a similar peak is seen at phase ∼0.73. Defining the signal phases as 0.72-0.75 and the background control phases as the OP phases (0.52-0.88) excluding the signal phases, 1600 ± 700 excess events were found, corresponding to 2.2σ . This pre-trial significance is too low to claim a detection and it is within the range of expected fluctuation of the background. An emission between the main pulse and the interpulse, i.e., a so-called bridge emission, is seen in some energy bands. With the MAGIC data above 25 GeV, it is not visible though the statistical uncertainty is large. Defining the signal phases as 0.04-0.32, i.e., from the end of P1 to the beginning of P2, and using the OP phases as the background estimate, 3200 ± 2800 excess events were found, corresponding to 1.1σ . The upper limit on the number of excess events with the 95% confidence level is 8800, which corresponds to half of the flux of P1 + P2.

Summary of Findings
The findings of this study can be summarized as follows.
1. 59 hr of MAGIC observations of the Crab pulsar during the winters 2007/2008 and 2008/2009 resulted in the detection of 6200 ± 1400 and 11300 ± 1500 excess events from P1 and P2, respectively. The flux of P2 is a factor ∼2 larger than that of P1 above 25 GeV. 2. No yearly variability in the pulse profile or in the flux was found. 3. The flux measured with MAGIC is significantly higher than the extrapolation of the exponential cutoff spectrum determined by Fermi-LAT. 4. The energy spectra extend up to at least 100 GeV and can be described by a power law between 25 GeV and 100 GeV. The power-law indices of P1, P2, and P1 + P2 are −3.4 ± 0.5 stat ± 0.3 syst , −3.1 ± 1.0 stat ± 0.3 syst , and −3.4 ± 0.5 stat ± 0.3 syst , respectively. The sensitivity of MAGIC I above 100 GeV is not sufficient to clarify if the spectrum continues with a power law or drops more rapidly. However, our spectrum of P1 + P2 and the VERITAS measured spectrum above 100 GeV seem to be a good extrapolation of each other. 5. Assuming an asymmetric Lorentzian for the pulse shape, the peak positions of the main pulse and the interpulse are estimated to be −0.0089 ± 0.0065 and 0.3929 ± 0.0087 in phase, while the FWHMs are 0.031 ± 0.010 and 0.053 ± 0.015. Compared with the Fermi-LAT measurements, the pulse widths are narrower in the MAGIC energy regime. 6. The bridge emission between P1 and P2 is weak. With the current sensitivity no signal was found. A potential third peak with a pre-trial significance of 2.2σ is seen at a similar position as in the AGILE data above 100 MeV and the Fermi-LAT data above 10 GeV, but it is consistent with the background fluctuation.
The spectrum of the Crab pulsar does not follow an exponential cutoff but, after the break, it continues as a power law. This is inconsistent with the OG and SG models in their simplest versions, where it is assumed that the emission above ∼1 GeV comes only from curvature radiation, leading to an exponential cutoff in the spectrum.
A theoretical interpretation of this deviation from the exponential cutoff is discussed in the next section.

Theoretical Interpretation of the Spectrum
In a pulsar magnetosphere, high-energy photons are emitted by the electrons and positrons that are accelerated by the magnetic-field-aligned electric field, E . To derive E , we must solve the inhomogeneous part of the Maxwell equations (Fawley et al. 1977;Scharlemann et al. 1978;Arons & Scharlemann 1979) ∇ where ρ denotes the real charge density, and ρ GJ ≡ −Ω · B/(2πc) the Goldreich-Julian charge density (Goldreich & Julian 1969;Mestel 1971); the angular-velocity vector Ω points in the direction of the NS spin axis with magnitude |Ω| = 2π/P , B refers to the magnetic field, and c is the speed of light. If ρ coincides with ρ GJ in the entire magnetosphere, E vanishes everywhere. However, if ρ deviates from ρ GJ in some region, it is inevitable for a non-vanishing E to arise around the region and the particle accelerator, or the so-called "gap," appears.
To predict the absolute luminosity of the gap, as well as any phase-averaged and phase-resolved emission properties, we must constrain E ≡ |E |, ρ, and the gap geometry in the three-dimensional magnetosphere. We solve the Poisson equation (9) together with the Boltzmann equation for electrons and positrons (e ± 's) and with the radiative transfer equation between 0.005 eV and 10 TeV. The position of the gap is solved within the free-boundary framework so that E may vanish on the boundaries, and turned out to distribute in the higher altitudes as a quantitative extension of previous OG models, which assume a vacuum (i.e., ρ = 0) in the gap. In the present paper, we propose a new, non-vacuum (i.e., ρ = 0) OG model, solving the distribution of ρ self-consistently from the e − -e + pair creation at each point (Beskin et al. 1992;Hirotani & Okamoto 1998;Hirotani & Shibata 1999a, 1999bHirotani et al. 1999;Takata et al. 2006).
The e ± 's created are polarized and accelerated by E in the gap to attain high Lorentz factors, ∼10 7.5 . Such ultra-relativistic e + 's and e − 's emit primary γ -rays via synchro-curvature and inverse Compton (IC) processes. For the IC process, the target photons are emitted from the cooling NS surface, from the heated PC surface, and from the magnetosphere in which pairs are created. The primary γ -rays that are emitted by the (inwardly accelerated) e − 's efficiently (nearly head-on) collide with the surface X-rays to materialize as the primary e ± 's within the gap and as the secondary e ± 's outside the gap. The secondary e ± 's efficiently lose their energy via the synchrotron process in the inner magnetosphere and cascade into tertiary and highergeneration e ± 's via two-photon and one-photon (i.e., magnetic) pair-creation processes. The primary γ -rays that are emitted by the (outwardly accelerated) e + 's via the curvature process collide with the magnetospheric X-rays to materialize as secondary e ± 's with initial Lorentz factors γ ∼ 10 3.5 outside the gap, while those emitted via the IC process collide with the magnetospheric IR-UV photons to materialize with γ ∼ 10 7 . The secondary pairs with γ ∼ 10 3.5 emit synchrotron emission below 10 MeV and synchrotron self-Compton (SSC) emission between 10 MeV and a few GeV (by upscattering the magnetospheric X-ray photons). The secondary pairs with γ ∼ 10 7 , on the other hand, emit synchrotron emission below 10 GeV and SSC emission between 10 GeV and 1 TeV (by upscattering the magnetospheric IR-UV photons). Such secondary SSC photons (between 10 GeV and 1 TeV) are efficiently absorbed colliding with the magnetospheric IR-UV photons to materialize as tertiary pairs with 10 4 < γ < 10 6 . In the high-energy end (i.e., 10 5 < γ < 10 6 ), the tertiary pairs upscatter magnetospheric IR-UV photons into the energy range between 1 GeV and 300 GeV. Therefore, if we focus on the emission component above 10 GeV, the photons are emitted by the secondary and tertiary pairs when they upscatter the magnetospheric IR-UV photons (see also Lyutikov et al. 2011 for an analytical discussion of this process using the multiplicity factor of highergeneration pairs). Note that the primary IC photons are emitted from a relatively inner part of the outer magnetosphere and hence are totally absorbed by the magnetospheric IR-UV photon field, the specific intensity of which is self-consistently solved together with E and the particle distribution functions at each point. Note also that the surface X-ray field little affects the pair creation or the IC process in the outer magnetosphere of the Crab pulsar. For the details of this self-consistent approach, see Hirotani (2006). In Figure 15, we present the spectral energy distribution of the total pulse component when the magnetic axis is inclined 60 • with respect to the rotational axis and when the observer's viewing angle is 75 • . The temperature of the cooling NS emission is assumed to be 70 eV. The distance to the pulsar is assumed to be 2 kpc. The thick, green solid line represents the spectrum to be observed (i.e., with magnetospheric absorption), while the thin black solid line indicates the photons emitted by the primary positrons (with Lorentz factor ∼10 7.5 inside the gap), and the thin dashed line indicates those emitted by the secondary pairs (with Lorentz factor ∼10 3.5 or 10 7 outside the gap). Interstellar absorption is not taken into account. The primary, unabsorbed IC component becomes prominent above 40 GeV; however, most of these photons are absorbed by two-photon pair production and reprocessed in lower energies as the secondary SSC component. In the secondary emission component, which is depicted by the dashed curve, there is a transition of the dominant component: the synchrotron component dominates the IC component (due to the SSC process) below a few GeV, whereas the latter dominates the former above this energy.
It should be noted that the pulsed emission between 25 GeV and 180 GeV is dominated by the SSC component emitted by the secondary and tertiary e ± 's. Since the higher-generation components, which are denoted by the dotted curves, are emitted from the higher altitudes (near the light cylinder), they are less efficiently absorbed, thereby appearing as pulsed flux above 20 GeV. The resultant spectrum (green, thick solid line) exhibits a power-law-like shape above 20 GeV, rather than an exponential cutoff. We therefore interpret that the detected γrays above 25 GeV are mainly emitted via the SSC process when the secondary and tertiary pairs upscatter the magnetospheric synchrotron IR-UV photons.
Although the present theoretical result is obtained by simultaneously solving the set of Maxwell and Boltzmann equations under appropriate boundary conditions, it does not rule out other possibilities such as the synchrotron emission (e.g., Chkheidze et al. 2011) or the IC emission (e.g., Bogovalov & Aharonian 2000) from the wind zone by ultra-relativistic particles, which may be accelerated by MHD interactions or by magnetic reconnection, for instance. However, the fuller study of other theoretical models lies outside the scope of the present paper.

Outlook
For further studies of the pulsar emission mechanisms, observations with a higher sensitivity are essential. The stereoscopic system comprising the two MAGIC telescopes has a 50 GeV threshold and its sensitivity above 100 GeV is nearly three times higher than that of MAGIC I (Aleksić et al. 2011a). The results of stereoscopic observations are given elsewhere (Aleksić et al. 2011b).