CRAB NEBULA: FIVE-YEAR OBSERVATION WITH ARGO-YBJ

B. Bartoli; P. Bernardini; X. J. Bi; P. Branchini; A. Budano; P. Camarri; Z. Cao; R. Cardarelli; S. Catalanotti; S. Z. Chen; T. L. Chen; P. Creti; S. W. Cui; B. Z. Dai; A. D'Amone; Danzengluobu; I. De Mitri; B. D'Ettorre Piazzoli; T. Di Girolamo; G. Di Sciascio; C. F. Feng; Zhaoyang Feng; Zhenyong Feng; Q. B. Gou; Y. Q. Guo; H. H. He; Haibing Hu; Hongbo Hu; M. Iacovacci; R. Iuppa; H. Y. Jia; Labaciren; H. J. Li; G. Liguori; C. Liu; J. Liu; M. Y. Liu; H. Lu; L. L. Ma; X. H. Ma; G. Mancarella; S. M. Mari; G. Marsella; D. Martello; S. Mastroianni; P. Montini; C. C. Ning; M. Panareo; L. Perrone; P. Pistilli; F. Ruggieri; P. Salvini; R. Santonico; P. R. Shen; X. D. Sheng; F. Shi; A. Surdo; Y. H. Tan; P. Vallania; S. Vernetto; C. Vigorito; H. Wang; C. Y. Wu; H. R. Wu; L. Xue; Q. Y. Yang; X. C. Yang; Z. G. Yao; A. F. Yuan; M. Zha; H. M. Zhang; L. Zhang; X. Y. Zhang; Y. Zhang; J. Zhao; Zhaxiciren; Zhaxisangzhu; X. X. Zhou; F. R. Zhu; Q. Q. Zhu; G. Zizzi; (The ARGO-YBJ Collaboration); E. Striani

doi:10.1088/0004-637X/798/2/119

1. INTRODUCTION

The Crab Nebula is the remnant of a supernova exploded in 1054 A.D. at a distance of ∼2 kpc. It contains a 33 ms pulsar that powers a wind of relativistic particles. The interactions of these particles with the remnant gas, photons, and magnetic field produce a non-thermal radiation extending from radio waves to TeV gamma-rays (Bühler & Blanford 2014, and references therein). Most of the emission is generally attributed to synchrotron radiation of relativistic electrons and positrons. The spectral energy distribution (SED) peaks between optical and X-ray frequencies. A second component arises above ∼400 MeV, interpreted as Inverse Compton (IC) of the same electrons scattering off synchrotron photons and CMB photons.

The Crab Nebula is one of the most luminous sources of very high energy (VHE) gamma-rays in the sky and the first source to be detected at TeV energies (Weekes et al. 1989). Thanks to its high flux and apparent stability, it is considered a reference source in gamma-ray astronomy. Detected by many experiments, both Cherenkov telescopes (Aharonian et al. 2004, 2006; Albert et al. 2008) and air shower arrays (Amenomori et al. 2009; Abdo et al. 2012), the Crab Nebula was often used to check detector performance, including sensitivity, pointing accuracy, and angular resolution.

In 2010 September the AGILE satellite unexpectedly detected a strong flare from the direction of the Crab Nebula at energies above 100 MeV. It lasted two days, with a maximum flux three times higher than the average value (Tavani et al. 2011), later confirmed by Fermi (Abdo et al. 2011). From then on, Fermi and AGILE reported some more flares, characterized by a rapid increase and decay of the flux, typically lasting a few days. The most impressive occurred in 2011 April, when the observed flux was ∼10 times higher than usual (Buehler et al. 2012). The measured SED shows a new spectral component emerging during flares, peaking at high energies (up to hundreds of MeVs in the 2011 April flare), attributed to a synchrotron emission of a population of electrons accelerated up to energies of 10¹⁵ eV. The Fermi-LAT data also show that these sharp emission peaks are superimposed to long-lasting smoother modulations with timescales of weeks or months (Striani et al. 2013). The observed flux variations are attributed to the nebula, since the pulsar emission was found to be stable within 20% (Buehler et al. 2012). However the origin of this activity is still unclear. In this scenario, observations at higher energies could provide precious information to help understand the mechanisms responsible for this behavior.

A preliminary analysis of the data recorded by the air shower detector ARGO-YBJ during the flares, showed an increase of the Crab flux at TeV energies with a moderate statistical significance, in three out of four flares (Aielli et al. 2010b; Vernetto 2013). These results have not been confirmed by Imaging Atmospheric Cherenkov Telescopes (IACTs) because the Moon light hampered the observations. However, sporadic and short measurements carried out during the first part of the 2010 September flare by MAGIC and VERITAS show no evidence for a flux variability (Mariotti 2010; Ong 2010). The observations by VERITAS and HESS during a flare in 2013 March (when ARGO-YBJ was already switched off), report a counting rate consistent with the steady flux (Aliu et al. 2013; Abramowski et al. 2013).

In this paper we present a detailed analysis of the ARGO-YBJ data, carried out with a better reconstruction of the shower arrival direction, obtained by applying quality cuts on the events. The study concerns not only the flaring episodes, but the whole Crab Nebula data set, consisting of more than 5 yr of observation. The ARGO-YBJ layout and operation mode are presented in Section 2, with a particular attention to the performance in gamma-ray astronomy. In Section 3, the analysis technique to extract the gamma-ray signal is outlined, followed by the results obtained with Crab Nebula data. Section 4 reports the energy spectrum evaluation and discusses systematic errors. In Section 5 the analysis of the time behavior of the Crab Nebula signal is presented, with a search for possible flares and rate variations on different timescales. A time correlation analysis with the Fermi-LAT data at energy E > 100 MeV during ∼4.5 yr is also reported. Finally, Section 6 contains a summary of the results and concluding remarks.

2. THE ARGO-YBJ EXPERIMENT

The ARGO-YBJ is a "full coverage" air shower detector located at the Yangbajing Cosmic Ray Laboratory (Tibet, P.R. China, longitude 90 fdg 5 east, latitude 30 fdg 1 north) at an altitude of 4300 m above sea level, devoted to gamma-ray astronomy at energies above ∼300 GeV and cosmic-ray studies at energies above ∼1 TeV.

During its lifetime, from 2007 November to 2013 February, ARGO-YBJ monitored the gamma-ray sky with an integrated sensitivity ranging from 0.24 to ∼1 Crab Units (Bartoli et al. 2013a) and studied in detail the emission of the most luminous gamma-ray sources at energies above 300 GeV, namely the Crab Nebula, MGRO J1908+06 (Bartoli et al. 2012a), HESS J1841-055 (Bartoli et al. 2013b), the Cygnus Region (Bartoli et al. 2012b, 2014), and the blazars Mrk401 (Bartoli et al. 2011a; Aielli et al. 2010a) and Mrk501 (Bartoli et al. 2012c).

The detector consists of a ∼74 × 78 m² carpet made of a single layer of Resistive Plate Chambers (RPCs) with ∼92% of active area, surrounded by a partially instrumented (∼20%) area up to ∼100 × 110 m². The apparatus has a modular structure, with the basic data acquisition element being a cluster (5.7 × 7.6 m²) made of 12 RPCs (2.85 × 1.23 m²). Each RPC is read by 80 strips of 6.75 × 61.8 cm² (the spatial pixels), logically organized in 10 independent pads of 55.6 × 61.8 cm² which are individually acquired and represent the time pixels of the detector (Aielli et al. 2006). To extend the dynamical range up to PeV energies, each RPC is equipped with two large pads (139 × 123 cm²) to collect the total charge developed by the particles hitting the detector (Aielli et al. 2012). The full experiment is made of 153 clusters (18360 pads), for a total active surface of ∼6600 m².

ARGO-YBJ operated in two independent acquisition modes: the shower mode and the scaler mode (Aielli et al. 2008). In this analysis we refer to the data recorded from the digital read-out in shower mode. In this mode, an electronic logic was implemented to build an inclusive trigger, based on a time correlation between the pad signals, depending on their relative distance. In this way, all showers with a number of fired pads N_pad ⩾ N_trig in the central carpet in a time window of 420 ns generated the trigger. This trigger worked with high efficiency down to N_trig = 20, keeping the rate of random coincidences negligible (Aloisio et al. 2004).

The time of each fired pad in a window of 2 μs around the trigger time and its location were recorded. To calibrate the 18360 pads in time, a software procedure has been developed based on the Characteristic Plane method (He et al. 2007) that uses the secondary particles of large vertical showers as calibration beams. It iteratively reduces the differences between the time measurements and the time fit of the shower front (Aielli et al. 2009).

The full detector was in stable data taking with the trigger condition N_trig = 20 and an average duty cycle ∼86%. The trigger rate was ∼3.5 kHz with a dead time of 4%.

The detector performance and capabilities in gamma-ray astronomy have been studied and improved through Monte Carlo simulations describing the shower development in the atmosphere by using the CORSIKA code (Heck et al. 1998) and the detector response with a code based on the GEANT package (GEANT 1993).

2.1. Field of View

One of the distinctive features of air shower arrays is the large field of view (FOV), in principle including the entire overhead sky. Gamma-ray sources cross the FOV with different paths according to their declinations. The sensitivity is not uniform in the field of view, however. Given a photon flux, the atmospheric absorption reduces the rate of showers for increasing zenith angles. The cosmic-ray background also decreases, but more slowly, and the combination of the two rates determines the trend of the sensitivity as a function of the zenith angle. Figure 1 shows the event rate in ARGO-YBJ expected from a Crab-like source as a function of the zenith angle θ, normalized to the rate at θ = 0°, compared to the background rate. In the same figure the dependence of the detector sensitivity on θ is also reported. According to simulations, the sensitivity at θ = 30° (45°) is reduced by a factor ∼2 (∼10) with respect to the sensitivity at θ = 0°.

**Figure 1.** Sensitivity in the detector field of view. Solid line: zenith angle dependence of the sensitivity to a Crab-like source. The sensitivity is normalized to a zenith angle θ = 0°. Dotted line: zenith angle dependence of the gamma-ray event rate. Dashed line: zenith angle dependence of the background rate; both rates are normalized to θ = 0°.
Download figure:
Standard image High-resolution image

The capability to detect a given source depends on its path in the field of view (determined by the source declination), and in particular, on the amount of time that the source lies at different zenith angles. The maximum significance is for a declination δ_max = λ, where λ = 30 fdg 1 is the latitude of the detector. Given a Crab-like source, the sensitivity decreases by less than 10% for declinations |δ−δ_max| < 10°, while it is reduced by a factor ∼2 for declinations |δ−δ_max| ∼ 30°. The declination dependence is slightly stronger (weaker) for sources with softer (harder) spectra with respect to the Crab Nebula (Bartoli et al. 2013a).

At the ARGO-YBJ site, the Crab Nebula (declination δ = 22 fdg 01) culminates at a zenith angle θ = 8 fdg 1 and lies at zenith angles θ < 45° for 6.6 hr per sidereal day. In general, following a source for a longer time per day increases the signal significance, because of the increasing statistics. But, since the signal to background ratio decreases at large zenith angles, there is a maximum zenith angle beyond which the significance begins to reduce. According to simulations, the maximum zenith angle for the Crab Nebula is ∼45°.

2.2. Angular Resolution

The sensitivity needed to observe a gamma-ray source is related to the angular resolution, which determines the amount of cosmic-ray background. We evaluate the shower arrival direction by fitting the shower front with a conical shape centered on the shower core position, to take into account the time delay of secondary particles with respect to a flat front, a delay that increases with the distance from the core. We set this delay to 0.1 ns m⁻¹ (Aielli et al. 2009).

The high granularity of the detector allows the study of the shower profile in great detail, and the accurate determination of the core position by fitting the lateral density distribution with a Nishimura–Kamata–Greisen-like function. According to simulations, the core position error depends on the number of hit pads N_pad and the core distance from the detector center. For gamma-ray induced showers with a core distance less than 50 m, the average core position error is less than 8 (2) m for N_pad ⩾ 100 (1000).

The point-spread function (PSF) also depends on N_pad, and for a given N_pad value, it worsens as the shower core distance from the detector center increases. The angular resolution for showers induced by cosmic-rays has been checked by studying the Moon shadow, observed by ARGO-YBJ with a statistical significance of ∼9 standard deviations per month. The shape of the shadow cast by the Moon on the cosmic-ray flux provides a measurement of the detector PSF. This measurement has been found to be in excellent agreement with expectations, confirming the reliability of the simulation procedure (Bartoli et al. 2011b).

The PSF for gamma-ray showers is narrower than the cosmic-ray one by ∼30%–40%, due to the better defined time profile of the showers. To improve the angular resolution for gamma-ray astronomy studies, quality cuts have been implemented by rejecting the events with a core distance larger than a given value D_cut (depending on N_pad) and with an average time spread of the particles with respect to the fitted shower front exceeding 9 ns (Bartoli et al. 2013a). The values of D_cut are given in Table 1. The fraction of gamma-rays passing the selection cuts depends on N_pad averaging ∼80%, whereas the fraction of surviving background events is ∼76% for N_pad < 100 and ∼50% for N_pad ⩾ 100. The selection also acts as a mild gamma/hadron discrimination for events with N_pad ⩾ 100 (the sensitivity increases by a factor ∼1.1).

Table 1. Characteristics of Crab Nebula Simulated Events

N_pad	D_cut^a	Core Position	R₃₉^c	Median Energy
N_pad	(m)	Error^b (m)	(deg)	(TeV)
20–39	No limits	37	1.88	0.34
40–59	No limits	28	1.50	0.53
60–99	90	12	1.04	0.79
100–199	70	6.8	0.70	1.3
200–299	60	4.2	0.50	2.1
300–499	60	3.3	0.41	3.1
500–999	40	2.3	0.32	4.8
1000–1999	30	1.6	0.24	8.1
⩾2000	30	1.0	0.19	17.7

Notes. ^aMaximum distance of the shower core from the detector center, beyond which the events are rejected. ^bDistance between the true and reconstructed cores containing 68% of the events. ^cAngular resolution, defined as the 39% containment radius.

Download table as: ASCII Typeset image

The arrival directions of the selected showers are also corrected for the systematic error due to the partial sampling of the shower front when the core is close to the edge of the detector (Eckmann 1991). This systematic error is related to the angle between the vector "shower core-detector center" and the shower arrival direction. For events with N_pad ⩾ 100, for which the core position is determined with more accuracy, the error can be considerably reduced.

These selections and corrections shrink the PSF by a factor ranging from ∼1.1 for events with N_pad = 20–39, up to ∼2, for N_pad ⩾ 1000. The PSFs obtained by simulating the Crab Nebula along its daily path up to θ = 45° are shown in Figure 2 for different intervals of N_pad.

**Figure 2.** Angular resolution for different N_pad intervals, according to simulations. The curves represent the fraction of events beyond the angular distance d from the source, as a function of d.
Download figure:
Standard image High-resolution image

To describe the PSFs analytically, for small values of N_pad that cannot be simply fitted by a two-dimensional Gaussian function, the simulated distributions have been fitted with a linear combination of two Gaussians. In general, when the PSF is described by a single Gaussian (F(r) = 1/(2πσ²) exp (−r²/σ²), where r is the angular distance from the source position), the value of the root mean square σ is commonly defined as the "angular resolution." In this case, the fraction of events within 1σ is 39%. For our PSFs, the value of the 39% containment radius R₃₉ ranges from 0 fdg 19 for N_pad ⩾ 2000 to 1 fdg 9 for N_pad = 20–39. Table 1 reports the values of R₃₉ for different N_pad intervals, together with the core position error, after quality cuts, as obtained by simulating the source during the daily path in the ARGO-YBJ field of view.

2.3. Energy Measurement

The number of hit pads N_pad is the observable related to the primary energy that is used to infer the source spectrum. In general, the number of particles at ground level is not a very accurate estimator of the primary energy of the single event, due to the large fluctuations in the shower development in the atmosphere. Moreover, for a given shower, the number of particles detected in a finite area detector like ARGO-YBJ depends on the position of the shower core with respect to the detector center; for small showers this is especially poorly determined.

The relation between N_pad and the primary gamma-ray energy of showers surviving the selection cuts is illustrated in Figure 3, where the corresponding primary energy distributions for different N_pad intervals are reported, as obtained by simulating a Crab-like source with a power law spectrum with index −2.63. The distributions are broad, with extended overlapping regions, spanning over more than one order of magnitude for small N_pad values. The median energies for different N_pad intervals are given in Table 1. They range from 340 GeV for events with N_pad = 20–39, to ∼18 TeV for N_pad ⩾ 2000.

Since the variable N_pad does not allow the accurate measurement of the primary energy of a single event, the energy spectrum is evaluated by studying the global distribution of N_pad. The observed distribution is compared to a set of simulated ones obtained with different test spectra to determine the spectrum that better reproduces the data.

3. THE CRAB NEBULA SIGNAL

The data set used for this analysis contains all the events recorded from 2007 November to 2013 February, with N_pad ⩾ 20. The total on-source time is 1.12 × 10⁴ hr.

For each source transit, the events are used to fill a set of nine 12° × 12° sky maps centered on the Crab Nebula position, with a bin size of 0 fdg 1 × 0 fdg 1 in right ascension and declination ("event maps"). Each map corresponds to a defined N_pad interval: 20–39, 40–59, 60–99, 100–199, 200–299, 300–499, 500–999, 1000–1999 and N_pad ⩾ 2000.

To extract the excess of gamma-rays, the cosmic-ray background has to be estimated and subtracted. Using the time swapping method (Alexandreas et al. 1993), the shower data recorded in a time interval Δt = 2–3 hr are used to evaluate the "background maps," i.e., the expected number of cosmic-ray events in any location of the map for the given time interval. This method assumes that during the interval Δt the shape of the distribution of the arrival directions of cosmic-rays in local coordinates does not change, while the overall rate could change due to atmospheric and detector effects. The value of the time interval Δt is less than a few hours to minimize the systematic effects due to the environmental parameters variations that could change the distribution of the arrival directions.

The time swapping method is a sort of "simulation" based on real data: for each detected event, n_f "fake" events (with n_f = 10) are generated by replacing the original arrival time with new ones, randomly selected from an event buffer that spans the time Δt of data taking. By changing the time, the fake events maintain the same declination of the original event, but have a different right ascension. With these events, a new sky map (background map) is built, with statistics n_f times larger than the event map in order to reduce the fluctuations. To avoid the inclusion of the source events in the background evaluation, the showers inside a circular region around the source (with a radius related to the PSF and depending on N_pad) are excluded from the time swapping procedure. A correction on the number of swaps is applied to take into account the rejected events in the source region (Fleysher et al. 2004).

Event and background maps are then smoothed according to the PSF corresponding to each N_pad interval. Finally, the smoothed background maps are subtracted from the smoothed event maps, obtaining the "excess maps," where for every bin the statistical significance of the excess is calculated as:

$\begin{equation*} n_{\sigma } = \frac{N_E-N_B}{\sqrt{\delta N_E^2+ \delta N_B^2}} \end{equation*}$

with N_E = Σ_i n_i w_i and N_B = Σ_i b_i w_i/n_f. In these expressions n_i and b_i are the number of events of the ith bin of the event map and background map, respectively; w_i is a normalized weight, proportional to the value of the PSF at the angular distance of the ith bin; and n_f is the number of swappings. The sum is over all the bins inside a radius R_max, chosen to contain the signal events and depending on the PSF. Since the number of events per bin is large, the fluctuations follow the Gaussian statistics, hence the errors on N_E and N_B are: δN_E = $\scriptstyle\sqrt{\Sigma _i n_i w_i^2}$ and δN_B = $\scriptstyle \sqrt{\Sigma _i b_i w_i^2/n_f^2}$ .

The number of gamma-ray events from the source is:

$\begin{equation*} N_{\gamma } = \frac{N_E-N_B}{2 \pi \int _0^{R_{{\rm max}}} w(r)^2 r dr} \end{equation*}$

where w(r) is the weight used in the smoothing procedure, calculated at the angular distance r from the source position.

When adding all data, an excess consistent with the Crab Nebula position is observed in each of the nine maps, with a total statistical significance of 21.1 standard deviations. The number of excess events is ∼3.3 × 10⁵, corresponding to 189 ± 16 day⁻¹, where a "day" means a source transit. Table 2 gives the signal significance for each map and the corresponding event rates measured from the source. For comparison, the background rates measured inside an angular window of 1° radius around the source are given. Figure 4 shows the total significance map.

**Figure 4.** Significance map of the Crab Nebula region, after smoothing and background subtraction.
Download figure:
Standard image High-resolution image

Table 2. Summary of Crab Nebula Data

N_pad	Photon Rate	Background Rate^a	Significance	E_med	Differential Flux
N_pad	(events day⁻¹)	(events day⁻¹)	(s.d.)	(TeV)	(photons cm⁻² s⁻¹ TeV⁻¹)
20–39	56.7 ± 12.2	1.3 × 10⁴	4.6	0.34	(6.23 ± 1.34) × 10⁻¹⁰
40–59	75.4 ± 9.2	1.1 × 10⁴	8.2	0.53	(1.80 ± 0.21) × 10⁻¹⁰
60–99	34.7 ± 3.9	4.2 × 10³	9.0	0.79	(5.92 ± 0.66) × 10⁻¹¹
100–199	15.4 ± 1.7	1.9 × 10³	8.9	1.30	(1.37 ± 0.15) × 10⁻¹¹
200–299	5.23 ± 0.61	4.9 × 10²	8.5	2.1	(5.30 ± 0.63) × 10⁻¹²
300–499	3.51 ± 0.44	3.8 × 10²	8.0	3.1	(1.75 ± 0.22) × 10⁻¹²
500–999	2.07 ± 0.27	2.4 × 10²	7.6	4.8	(5.62 ± 0.74) × 10⁻¹³
1000–1999	0.50 ± 0.13	87.4	3.8	8.1	(1.00 ± 0.26) × 10⁻¹³
⩾2000	0.23 ± 0.07	34.2	3.5	17.7	(1.87 ± 0.54) × 10⁻¹⁴

Note. ^aAverage background rate within an angular distance of 1° from the source.

Download table as: ASCII Typeset image

Finally, the gamma-ray signal can be used to check the detector angular resolution since the Crab Nebula's angular size is small compared to the detector PSF. Figure 5 shows the distribution of the arrival directions of the excess showers with respect to the source position, for N_pad ⩾ 20, 100, and 500, compared to simulations. The agreement is excellent.

**Figure 5.** Comparison of experimental and simulated PSFs: event rate per solid angle as a function of the angular distance from the source position, for events with N_pad ⩾ 20, 100, and 500. The solid lines are the results of the simulation.
Download figure:
Standard image High-resolution image

4. ENERGY SPECTRUM

The energy spectrum is evaluated by comparing the number of events detected from the Crab Nebula in the previously defined N_pad intervals to the expected number given by a simulation assuming a set of test spectra. We consider the power law spectrum:

$\begin{eqnarray*} &&\frac{dN}{dE}(I_0,\alpha) = I_0 \left (\frac{E}{2 \; \mathrm{TeV}}\right)^{-\alpha } \end{eqnarray*}$

where the flux normalization I₀ and slope α are the parameters to be estimated with the fitting procedure.

The fit is made by minimizing the value of χ², evaluated for any couple of parameters as:

$\begin{eqnarray*} &&\chi ^2(I_0,\alpha) = \sum _{j=1,9} \frac{\left(N^j_{\gamma } - N^j_{{\rm MC}}(I_0,\alpha)\right)^2}{\left(\delta N^j_{\gamma }\right)^2 + \left(\delta N^j_{MC}\right)^2} \end{eqnarray*}$

where $N^j_{\gamma }$ and $N^j_{{\rm MC}}$ are the number of events detected and expected, respectively, in the jth N_pad interval.

The best-fit parameters obtained are I₀ = (5.2 ± 0.2) × 10⁻¹² photons cm⁻² s⁻¹ TeV⁻¹, and α = −2.63 ± 0.05, with χ² = 5.8 for 7 degrees of freedom, corresponding to a p-value p = 0.56. The integral flux above 1 TeV is 1.97 × 10⁻¹¹ photons cm⁻² s⁻¹. The flux at 1 TeV obtained in this work is 7% higher than that reported in a previous ARGO-YBJ paper (Bartoli et al. 2013a). The difference is due to the correction of the event rates applied in this work to reduce environmental and detector effects on the trigger rate, as described in Section 5.3.

Figure 6 shows the obtained spectrum compared with the results of other experiments. The energy of each point is the gamma-ray median energy for the corresponding N_pad interval. The values of energies and differential fluxes are given in Table 2. The spectrum is consistent with a constant slope from ∼300 GeV to ∼20 TeV and agrees rather well with the measurement by HEGRA and MAGIC, whereas the HESS and Milagro fluxes are about 20% higher in the ∼1–10 TeV energy range.

The data are less clear concerning a possible energy cutoff at higher energies. MAGIC (Albert et al. 2008) and HESS (Aharonian et al. 2006) show a steepening below 20 TeV, while the HEGRA spectrum is harder and continues with a slight softening up to ∼75 TeV (Aharonian et al. 2004). A possible cutoff is also observed by Milagro at ∼30 TeV (Abdo et al. 2012). The limited statistics of our data at high energy does not allow us to draw any conclusion about the spectral properties above 20 TeV. Selecting events with N_pad ⩾ 3000 (whose median energy is 26 TeV assuming a power law spectrum with index α = −2.63) the statistical significance of the signal is 0.75.

When fitting the data with a power law spectrum with an exponential cutoff:

$\begin{eqnarray*} &&\frac{dN}{dE}(I_0,\alpha,E_{{\rm cut}}) = I_0 \left (\frac{E}{2 \; \mathrm{TeV}}\right)^{-\alpha } \; \exp \; (-E/E_{{\rm cut}}) \end{eqnarray*}$

the obtained p-value is always smaller than without a cutoff, for any value of E_cut. For E_cut = 14.3 TeV (the best-fit value obtained by HESS) the p-value is 0.13. We found that the p-value is larger than 10% for any value of E_cut > 12 TeV, indicating that the presence of a cutoff above ∼10 TeV cannot be excluded, even if our data seems more consistent with a pure power law.

4.1. Estimation of Systematic Errors

The previous results can be affected by systematic errors of different origin. In the following, we discuss the possible sources of systematics, evaluating their effects both on the flux normalization and the spectral slope.

1.
Energy scale. In our measurement, the number of hit pads N_pad is used as an estimator of the primary energy. The relation between the primary energy and N_pad is given by Monte Carlo simulations. Possible uncertainties and simplifications in the simulation procedure (both in the shower development and the detector response) could produce an incorrect N_pad value and consequently an error in the energy scale.The energy scale reliability has been checked using the Moon shadow. Due to the geomagnetic field, cosmic-rays are deflected according to their energy, and the shadow that the Moon casts on the cosmic-ray flux is shifted with respect to the Moon's position by an amount depending on the energy. The westward shift of the shadow has been measured for different N_pad intervals and compared to simulations. From the analysis of the Moon data, we found that the total absolute energy scale error is less than 13% in the proton energy range ∼1–30 TeV (Bartoli et al. 2011b). This estimate includes the uncertainties of the cosmic-ray elemental composition and the hadronic interaction model.From this result, given a gamma-ray spectrum with index α = −2.63, the corresponding systematic error in the flux normalization would be less than 22%.
2.
Pointing error. Fitting the angular distribution of gamma-rays around the Crab Nebula position we found that the pointing error is less than 01. A pointing error affects the measured gamma-ray flux, since the number of photons is obtained by a smoothing procedure weighting the events with a PSF centered at the source nominal position. An incorrect position would produce a loss of signal. Since the PSF is narrower for events with large N_pad, the loss is larger at high multiplicities, and generates a steepening of the spectrum. According to our simulation, a pointing error of 01 would produce a loss of signal ranging from 0.1% for N_pad = 20–39 to 6.0% for N_pad > 2000. As a consequence, the spectral index would increase by 0.01 and the flux normalization would decrease by 2%.
3.
Background evaluation. Our measurement is based on a very precise evaluation of the background. As explained in Section 3, the number of gamma-rays is given by the difference between the number of events detected in the event map (that contains the source events plus the cosmic-ray background) and the number of background events estimated with the time swapping method. Since the ratio between the number of gamma-rays and the number of background events is very small (ranging from ∼3 × 10⁻⁴ for N_pad = 20–39 up to ∼4 × 10⁻² for N_pad > 300), even a small systematic error in the background evaluation could produce a big error in the source flux.Possible sources of systematics are: (a) the presence of cosmic-ray excess regions due to the medium scale anisotropy, as those reported in (Bartoli et al. 2013c), close to the source, (b) changes in atmospheric conditions able to modify the background distribution in local coordinates in less than 2–3 hr, and (c) the detector malfunctioning.Such effects, when present, could generate extended regions in the signal map with evident excesses or deficits, in some cases involving the whole map. Instead, an accurate evaluation of the background produces a map with all the bin contents consistent with zero, except at the position of real sources.Concerning the medium scale anisotropy, we adopted a particular procedure to correct the background systematics in the sky regions coincident or adjacent to cosmic-ray excesses (Bartoli et al. 2013a). This correction is not necessary in the Crab Nebula region.Concerning points (b) and (c), it has to be specified that the maps are built with data sets of 2–3 hr and individually checked. When a map shows significant anomalies, the corresponding data set is rejected, so only "good maps" are combined to build to the "total maps."To test the background reliability of the total maps, we use the regions that are not involved in the Crab Nebula emission, i.e., the bins with an angular distance from the source larger than a minimum value, depending on the PSF. From these "out-source" regions we expect no significant excess, since they do not include any other known gamma-ray source with a flux above the ARGO-YBJ sensitivity. For any of the nine maps, we have evaluated the distribution of the excesses in the out-source region bins (before smoothing). We found that all the distributions are well described by Gauss functions with mean values consistent with zero and rms consistent with unit.Adding all the nine maps together, the total number of events detected from the out-source regions is 1.18 × 10⁹. This value differs from the corresponding estimated background by −9.3 × 10³ events, corresponding to −0.3 standard deviations. Since there is no significant excess or deficit of events, we can calculate the upper limit of the systematic error in the out-source region. We found that the relative error in the background value is less than 3.7 × 10⁻⁵ at 90% confidence level.We can reasonably assume that a similar systematic error involves the region of the map containing the source signal, and that all nine maps have comparable systematic errors. Based on these assumptions, we can evaluate the effects of such an error on the signal, which are obviously more relevant for the maps in which the signal-to-background ratio is smaller. We found that the error in the photon number is <13% for N_pad = 20–39, <1% for N_pad = 100–199, and <0.01% for N_pad ⩾ 1000.According to these values, the corresponding systematic error in the spectrum flux normalization would be less than 2%, and the error in the spectral index would be less than 0.05.
4.
Event rate variations. Studying the rate of cosmic-ray showers over 5 yr, we observed variations on timescales from hours to months up to 10% with respect to the mean value. These variations are mostly due to: (a) variation of atmospheric pressure and temperature that modify the showers' propagation in the atmosphere, (b) variation of the detector efficiency due to changes of the local temperature and pressure, and (c) aging of the detector.

Gamma-rays are assumed to be subject to similar variations. To study the stability of the Crab Nebula flux, we corrected the rate of the events observed from the source using the cosmic-ray rate as a normalization factor (see Section 5.3). However, an absolute normalization cannot be performed. The Monte Carlo simulations refer to a fixed atmospheric condition and a given detection efficiency that cannot exactly reproduce the average effect over several years of different conditions. Considering the amount of the observed rate variation, a reasonable estimation indicates a possible systematic error in the flux smaller than 4%.

Total systematic error. Adding all these contributions linearly, we conservatively estimate the total systematic error to be less than 30% for the flux normalization and 0.06 for the spectral index.

5. CRAB NEBULA LIGHT CURVE

To study the stability of the Crab Nebula emission, we consider the events with N_pad ⩾ 40. Our total data set consists of 1851 days, with an average observation time of 6.0 hr per day. The average rate of events with N_pad ⩾ 40 is 137 ± 10 day⁻¹.

Figure 7 shows the observed rate for events with N_pad ⩾ 40, 100, and 500, as a function of the Julian date, in bins of 200 days. The median energies corresponding to these N_pad thresholds are 0.76, 1.8, and 5.1 TeV, respectively. The signal appears stable during 5 yr for any threshold within the statistical fluctuations. Assuming a constant rate, the obtained χ² are 15.7, 3.27, and 5.17 (with 9 d.o.f.) for N_pad ⩾ 40, 100, and 500. The corresponding p-values are 0.073, 0.95, and 0.82, respectively.

A 6 yr monitoring of the Crab Nebula was previously performed by the Tibet-III air shower array from 1999 to 2005, at energies ∼3 TeV, with a sensitivity 3–4 times lower than that of ARGO-YBJ, reporting a yearly flux consistent with a steady emission (Amenomori et al. 2009).

5.1. Search for Flares

To make a "blind" search for rate variations in short timescales, we consider all the time intervals of duration Δt ranging from 1 to 15 days, starting from every observation day. This time range has been chosen on the basis of the duration of the flares observed in the GeV energy region. For each interval, we compare the observed rate of Crab events with the average rate and evaluate the significance of the excess as: σ_i = (R_i − R_m)/δ(R_i − R_m), where R_i is the counting rate in the i-th interval, R_m is the average counting rate, and δ(R_i − R_m) is the statistical error of the difference R_i − R_m. Note that the values of σ_i are not independent, since the time intervals overlap.

Figure 8 shows the distributions of σ_i for Δt = 1 day and Δt = 2–15 days, for N_pad ⩾ 40. The total number of intervals is 1851 for Δt = 1 and 25911 for Δt = 2–15 days. The distributions can be fitted by a Gauss function with mean value m = −0.04 ± 0.03 and rms = 1.05 ± 0.02 for Δt = 1 day, and m = −0.06 ± 0.01 and rms = 1.061 ± 0.005 for Δt = 2–15 days. The root mean square values indicate rate variations slightly larger than what would be expected by statistical fluctuations. However no significant excess is observed for any of the considered time intervals.

Given the ARGO-YBJ sensitivity, a flare would produce a five standard-deviation signal (pre-trial) if the flux exceeds the average value by a factor f ∼ 10/ $\sqrt{\Delta t (\mathrm{days})}$ .

5.2. Correlation with Fermi-LAT Data

To reduce the number of trials in the search for possible flares, we can limit our analysis to the days in which a flare was observed by satellite instruments at lower energies. We consider the Fermi-LAT daily light curve at energy E > 100 MeV from 2008 August to 2013 February, obtained through the analysis of the scientific Fermi data publicly available at the Fermi Science Support Center.²²

The first panel of Figure 9 shows the daily light curve, representing the sum of the nebula and pulsar fluxes. The average flux is (2.66 ± 0.01) × 10⁻⁶ photons cm⁻² s⁻¹. Also, excluding the days with flares, the rate is variable, with modulations on timescales of weeks and months.

First, we consider the three largest Fermi flares, which occurred in 2009 February, 2010 September, and 2011 April (Abdo et al. 2011; Buehler et al. 2012). To define the time boundaries and duration of these flares, we select the days in which the Fermi flux is higher than 4 × 10⁻⁶ photons cm⁻² s⁻¹. The dates and the duration of the three flares are given in Table 3. The counting rates from the Crab Nebula measured by ARGO-YBJ with events with N_pad ⩾ 40 during the flares are compared with the average rate of 137 ± 10 day⁻¹. In all cases, the rates are slightly higher than the average value, but consistent with it within statistical errors (see Table 3). Summing the three flares the average rate is 205 ± 91 day⁻¹. Table 3 also shows the results concerning the events with N_pad ⩾ 100 and 500. No significant excess is present in this case, either.

Table 3. ARGO-YBJ Results during the Three Largest Crab Nebula Flares Observed by Fermi-LAT

	Δt	Max. Fermi Flux	ARGO-YBJ	ARGO Rate	ARGO Rate	ARGO Rate
	(days)	(photons cm⁻² s⁻¹)	Observation	N_pad ⩾ 40	N_pad ⩾ 100	N_pad ⩾ 500
	(days)	(photons cm⁻² s⁻¹)	Time (hr)	(photons day⁻¹)	(photons day⁻¹)	(photons day⁻¹)
Flare 1^a	8	6.3 ± 0.8 × 10⁻⁶	49.6	142 ± 151	21 ± 28	2.5 ± 4.6
Flare 2^b	5	6.4 ± 0.8 × 10⁻⁶	31.5	265 ± 190	58 ± 36	−3.8 ± 5.7
Flare 3^c	9	19.8 ± 0.8 × 10⁻⁶	58.0	228 ± 144	51 ± 27	2.9 ± 4.4
Sum of three flares	22		139	205 ± 91	41 ± 17	1.3 ± 2.8

All ARGO data				137 ± 10	27 ± 2	2.8 ± 0.3

Notes. ^aStart time MJD 54864 (2009 February 2). ^bStart time MJD 55457 (2010 Septemebr 18). ^cStart time MJD 55662 (2011 April 11).

Download table as: ASCII Typeset image

Our preliminary analysis reported in (Aielli et al. 2010b) showed a four standard-deviation excess observed in the time interval from 2010 September 17 to 22, from a direction consistent with the Crab Nebula. However, removing the contribution of the steady flux and taking into account the number of trials, the post-trial significance was about two standard deviations. A further excess with a similar post-trial statistical significance was observed during the 2011 April flare (Vernetto 2013). In the present work, based on a better shower reconstruction and the event selection described in Section 2.2, the significance of the Crab Nebula signal integrated over 5 yr increases by about 15% with respect to the old analysis, but the signal observed during the Fermi flares decreases. The flux measured during the flares appears slightly higher than what would be expected from the steady emission, but consistent with it within one standard deviation. Both our previous analysis and the current one hint at a possible flux enhancement during the flares, but the reduced significance prevents us from drawing a definitive conclusion.

To extend the search for flares to the whole observation time, and not limit the analysis to the largest flares, we selected the Fermi data according to the measured daily flux and checked the corresponding ARGO-YBJ event rate. Table 4 reports the ARGO-YBJ rates for different levels of the Fermi flux and different N_pad thresholds. The rates are consistent with the average rate for any Fermi flux level. In particular, the ARGO-YBJ rate (for N_pad ⩾ 40) measured in the 62 days in which the Fermi flux exceeds 4 × 10⁻⁶ photons cm⁻² s⁻¹, is 190 ± 55 day⁻¹, i.e., 1.4 ± 0.4 times higher than the average rate.

Table 4. ARGO-YBJ Photon Rate for Different Flux Levels Measured by Fermi

Fermi Flux	Number of Days	ARGO Rate	ARGO Rate	ARGO Rate
(10⁻⁶ photons cm⁻² s⁻¹)		N_pad ⩾ 40	N_pad ⩾ 100	N_pad ⩾ 500
(10⁻⁶ photons cm⁻² s⁻¹)		(photons day⁻¹)	(photons day⁻¹)	(photons day⁻¹)
<2.0	175	197 ± 34	26 ± 6	3.1 ± 1.0
2.0–3.0	915	118 ± 15	27 ± 3	2.8 ± 0.4
3.0–4.0	435	148 ± 21	28 ± 4	3.0 ± 4.3
4.0–5.0	46	188 ± 64	33 ± 12	1.9 ± 2.0
>5.0	16	198 ± 107	54 ± 20	0.2 ± 3.2

All ARGO data		137 ± 10	27 ± 2	2.8 ± 0.3

Download table as: ASCII Typeset image

Finally, to study a possible correlation on timescales of months or years, we compare the light curves of the two detectors over the common observing time (∼4.5 yr), dividing the data into bins of 200 days. The bin width is chosen in order to have a significant signal in the ARGO-YBJ data (about seven s.d.).

Since the flux measured by Fermi is the sum of the nebula and pulsar contributions, and since the pulsar flux F_P averaged over the pulsation period is also stable during flares (F_P = (2.04 ± 0.01) × 10⁻⁶ photons cm⁻² s⁻¹ for E >100 MeV (Buehler et al. 2012)), the flux of the pulsar has been subtracted. The obtained average nebula flux is (6.2 ± 0.1) × 10⁻⁷ photons cm⁻² s⁻¹. The nebula flux shows variations up to ∼30% of the average flux, with χ² = 126, for 8 d.o.f. (see the second panel of Figure 9). The large variations are not only due to flares. In the same figure, the dashed curve shows the flux obtained excluding the 62 "flaring days." In this case, the average value is (5.6 ± 0.1) × 10⁻⁷ photons cm⁻² s⁻¹, with χ² = 80.

The lower panel of Figure 9 shows the corresponding ARGO-YBJ data for N_pad ⩾ 40. The average rate is 139.3 ± 10.6 events day⁻¹ (χ² = 14.2 for 8 d.o.f., p-value p = 0.077). Even if the ARGO-YBJ rate variations are consistent with statistical fluctuations, the Fermi and ARGO-YBJ data seems to follow a similar trend. The ARGO-YBJ rate appears higher in the "hot" Fermi periods. The dashed curve is obtained after the exclusion of the flaring days.

Figure 10 shows the ARGO-YBJ percentage rate variation with respect to the mean value (ΔF_ARGO) as a function of the corresponding variation of the Fermi rate (ΔF_Fermi), for the nine bins of the light curve. The Pearson correlation coefficient between the two data sets is r = 0.56 ± 0.22. The quoted error for r is the root mean square of the distribution of the correlation coefficients obtained by simulating the fluctuations of the counting rates of each bin, according to their statistical errors. Fitting the nine points with the function ΔF_ARGO = a ΔF_Fermi + b, the values of the best-fit parameters are a = 0.88 ± 0.37 and b = 0.018 ± 0.079, with χ² = 8.3 for 7 d.o.f. Discarding the 62 "flaring" days, the correlation coefficient becomes r = 0.45 ± 0.23 and the parameters of the linear fit are a = 0.96 ± 0.45 and b = 0.018 ± 0.082, with χ² = 10.4.

The same analysis has been performed using a different bin width, ranging from 10 to 450 days. The corresponding correlation coefficient steadily increases from r = 0.10 ± 0.06 (10 days) to r = 0.59 ± 0.23 (450 days). It has to be noted, however, that when using a small bin width, the ARGO-YBJ signal is not significant enough to search for a correlation unless the flux variations are very large. In 10 days, for example, the average ARGO-YBJ signal is 137 ± 135 events day⁻¹. The statistical fluctuations would hide a possible flux variation, unless the flux becomes more than a factor of 4–5 higher than the average.

The above results refer to events with N_pad ⩾ 40. The correlation coefficient is lower when selecting more energetic events: using a bin width of 200 days, for N_pad ⩾ 100, r = 0.19 ± 0.31; and for N_pad ⩾ 500, r = 0.46 ± 0.28.

5.3. Stability of the ARGO-YBJ Data

When studying the time evolution of a signal over several years, a discussion on the possible causes of detector instabilities is mandatory, to exclude systematic effects that could produce artificial rate variations. Since the measured number of events from the source N_S = N_E − N_B is the difference between the number of events N_E detected in the source map and the number of background events N_B estimated with the time swapping method, one must separately analyze the stability of the different contributions.

1.
A loss of signal events N_S could be produced by variations of the pointing accuracy. Studying the Moon shadow month by month, we have verified that the pointing is stable within 0.1 deg (Bartoli et al. 2011b). Given the moderate angular resolution for events with N_pad ⩾ 40, such a value could produce signal fluctuations of less than 2%.
2.
A worsening of the detector angular resolution (due to an increase of the time resolution of RPCs occurring at particularly low temperatures) could produce a loss of signal events N_S. A broadening of the PSF would also cause a decrease of the Moon shadow signal. That, however, is found to be stable within statistical fluctuations.
3.
Atmospheric pressure and temperature variations can affect the RPC detection efficiency, which can also be altered by some RPC not working properly or by aging effects.
4.
Pressure and temperature produce changes in the shower rate on the order of a few percent due to the different conditions in which the showers propagate in the atmosphere.The two latter effects modify N_S, N_E, and N_B by about the same factor (neglecting the different behavior of cosmic-ray and gamma-ray showers, which in this contest can be considered a second-order effect). This allows the use of N_B to correct the Crab rate, by multiplying the Crab rate observed in a given time interval by the correction factor f_c = B_m/B, where B_m is the average background rate and B is the background rate in that interval. The light curve in Figure 7 has been corrected according to this method, with f_c ranging from 0.91 to 1.07.
5.
Further possible systematics could be an incorrect evaluation of the background N_B. In Section 4 we evaluated the accuracy of the background for the total source signal. To check the accuracy of the background along the years, we can use the same out-source regions previously defined. For events with N_pad ⩾ 40, the out-source light curve in 200-day bins has a mean value of −7.9 ± 19.0 events day⁻¹ and a χ² = 10.2 for 9 d.o.f., corresponding to a p-value p = 0.67. According to these results, the background is stable and should not introduce any systematic effect on the rate of the Crab signal.

6. SUMMARY AND CONCLUSIONS

The ARGO-YBJ events recorded over 5 yr have been analyzed to evaluate the Crab Nebula spectrum and study the temporal behavior of the gamma-ray emission. Using the events with N_pad ⩾ 20, the statistical significance of the gamma-ray signal reaches more than 21 standard deviations, and the observed photon rate is 189 ± 16 day⁻¹. The event angular distributions around the source are well described by the PSFs obtained by simulations.

The source spectrum extends over nearly two decades in energy and five decades in flux. The spectral shape is consistent with a power law behavior in the range of 0.3–20 TeV with a spectral index −2.63 ± 0.05. An exponential cutoff would be consistent with our data in case of a cutoff energy higher than 12 TeV at 90% confidence level.

The study of the Crab Nebula light curve has been carried out to check the stability of the flux over years and to search for possible flares on the timescale of days. All the known sources of rate instabilities have been examined and the effects corrected.

Concerning flares, a blind search for flux increases of a duration between 1 and 15 days shows no significant excess. The average rate of events with N_pad ⩾ 40 measured by ARGO-YBJ during the three most powerful flares detected by Fermi-LAT (in 2009 February, 2010 September, and 2011 April) is 205 ± 91 day⁻¹, which is consistent with the average value of 137 ± 10 day⁻¹.

The 5 yr ARGO-YBJ light curve with a binning of 200 days is consistent with a constant flux with a probability of 0.07. A correlation analysis with the corresponding Fermi-LAT data gives a Pearson correlation coefficient r = 0.56 ± 0.22. The small statistical significance of these results does not allow the claim for a flux variability correlated with the observations at lower energies. If such a correlation was due to a real astrophysical phenomenon, the found regression coefficient a = 0.88 ± 0.37 would imply a similar percentage variation in Fermi and ARGO-YBJ rates, suggesting a similar behavior of the gamma-ray emission at energies ∼100 MeV and ∼1 TeV.

So far, no variation of the Crab Nebula flux at TeV energies has been reported by any detector. Assuming the flares observed by AGILE and Fermi due to synchrotron radiation from a population of electrons accelerated up to 10¹⁵ eV, the inverse Compton emission associated with this population would occur in the Klein-Nishina regime and would produce gamma-rays of energy approximately equal to that of the electrons. Such a flux would not be detectable by any of the existing gamma-ray experiments. With these assumptions, a TeV excess could hardly be interpreted as IC emission associated with the synchrotron radiation observed at lower energies, and would require a completely new interpretation.

This work is supported in China by NSFC (No. 10120130794), the Chinese Ministry of Science and Technology, the Chinese Academy of Science, the Key Laboratory of Particle Astrophysics, CAS, and in Italy by the Istituto Nazionale di Fisica Nucleare (INFN).

We also acknowledge the essential support of W. Y. Chen, G. Yang, X. F. Yuan, C. Y. Zhao, R. Assiro, B. Biondo, S. Bricola, F. Budano, A. Corvaglia, B. D'Aquino, R. Esposito, A. Innocente, A. Mangano, E. Pastori, C. Pinto, E. Reali, F. Taurino, and A. Zerbini in the installation, debugging, and maintenance of the detector.

CRAB NEBULA: FIVE-YEAR OBSERVATION WITH ARGO-YBJ

Article metrics

Permissions

Author e-mails

Author affiliations

Dates

ABSTRACT

1. INTRODUCTION