A SEARCH FOR ENHANCED VERY HIGH ENERGY GAMMA-RAY EMISSION FROM THE 2013 MARCH CRAB NEBULA FLARE

The Crab Nebula is one of the best-studied cosmic particle accelerators. Its distance of ∼2 kpc and absolute luminosity of 5 × 10³⁸ erg s⁻¹ allow the study of the nebula in great detail across the entire electromagnetic spectrum. From radio to GeV energies, the emission is consistent with synchrotron emission of relativistic electrons (Hester 2008). However, at higher energies, the dominant emission mechanism is thought to be inverse-Compton upscattering of low-energy photons by the same population of electrons (Gould 1965; Weekes et al. 1989; de Jager & Harding 1992; Aharonian et al. 2004).

The energy source powering the nebula is believed to be the Crab pulsar located at its center (Staelin & Reifenstein 1968). With the pulsar as the central engine, a self-consistent magnetohydrodynamic model can be developed that explains the main features of the nebula (Rees & Gunn 1974; Kennel & Coroniti 1984). The discovery of flaring episodes by the AGILE (Tavani et al. 2011) and Large Area Telescope on board the Fermi Gamma-ray Space Telescope (Fermi-LAT; Abdo et al. 2011) teams was unexpected in this framework. The Crab Nebula flux was seen to increase by more than a factor of ten in less than a day between 100 MeV and 1 GeV in the most extreme of these flares.

Determining the cause of these flares is a major experimental and theoretical challenge. The observed flaring timescales of 12 hr (Balbo et al. 2011) and 8 hr (Buehler et al. 2012) imply that the emission region is less than 3 × 10⁻⁴ pc in diameter. This size constraint coupled with the observation that the emitted isotropic power peaks at about 1% of the pulsar spin-down power argues in favor of an emission region that moves mildly relativistically (Buehler et al. 2012; Clausen-Brown & Lyutikov 2012; Lyutikov et al. 2012; Bednarek & Idec 2011). As no enhancement of the pulsed emission has been observed during flares, it has been concluded that the emission region likely resides outside the corotating magnetosphere (Buehler et al. 2012; Balbo et al. 2011).

The investigation of the origin of the flares is complicated because no correlated enhancements have been observed at other wavelengths to date (Balbo et al. 2011; Striani et al. 2011, 2013; Buehler et al. 2012). Multiwavelength campaigns have been executed every time a flare has been observed since the detection of the 2010 September flare (Tavani et al. 2011). Extensive simultaneous coverage over the entire synchrotron emission from radio to X-rays did not reveal correlated activity (Horns et al. 2010; Weisskopf 2012) that could have shed light on the location of the flares due to better angular resolution at these energies.

The non-detection of correlated activity favors a monoenergetic population of relativistic electrons as the origin of the observed flares. While multiwavelength coverage has been excellent in radio, optical, and X-rays, it has been rather sparse at energies above 100 GeV, i.e., in the inverse-Compton component. No enhancement of the TeV emission was reported by MAGIC or the Very Energetic Radiation Imaging Telescope Array System (VERITAS) during the 2010 September flare (Mariotti 2010; Ong 2010). The ARGO-YBJ Collaboration have reported enhanced signals with a median energy of 1 TeV from the direction of the Crab Nebula contemporaneous to GeV-band flares, although these enhancements did not reach the 5σ level (Aielli et al. 2010; Bartoli et al. 2012; Vernetto 2013).

The electrons responsible for the flares should also upscatter soft photons in the nebula to produce TeV photons, which enables constraining the dynamics of the electrons. In this paper we present the most sensitive observations at TeV energies performed during a flare of the Crab Nebula to date. These observations with VERITAS are discussed in the context of observations with Fermi-LAT.

2.1. VERITAS

VERITAS is an array of four 12 m diameter imaging atmospheric Cherenkov telescopes located at the base of Mt. Hopkins in southern Arizona, USA that observes very high energy (VHE; E > 100 GeV) gamma rays. Each telescope in the array has a reflector that is composed of 345 hexagonal mirror facets that focus light onto a 499-pixel photomultiplier tube camera at the focal plane with a field of view (FoV) of ∼35. The array operates in the energy range ∼0.1–30 TeV with an energy resolution of ∼15% at energies above 1 TeV and an angular resolution of 015 (Holder et al. 2008).

VERITAS observations of the Crab Nebula in its flaring state were triggered by an automated Fermi-LAT analysis pipeline at Barnard College-Columbia University (Errando et al. 2011) on 2013 March 2, two days prior to the ATel from the Fermi-LAT collaboration announcing the gamma-ray flare (Ojha et al. 2013). The VERITAS data during the flare are composed of 10 nights of observations in the period MJD 56353–56366 (2013 March 2–15, henceforth referred to as the flare time window, FTW). Observations of the Crab Nebula as part of the standard observing schedule from 2012 October 13 to 2013 April 2 excluding the FTW comprise a data set on the source in its non-flaring state, which is used as a baseline with which to compare the flare data.

All VERITAS Crab observations were taken in wobble mode with an offset of 05 from the source position alternately in each of the four cardinal directions, so that the background can be estimated from simultaneously gathered data, and systematic effects in the background estimation cancel out (Aharonian et al. 2001; Berge et al. 2007). Observations were conducted using the full four-telescope array in a range of zenith angles 12°–55°, giving a total of 10.3 hr of live time on the source during the FTW and 17.4 hr during the rest of the season. Two nights of flare observations (MJD 56353 and 56358) were conducted at large zenith angles, which has the effect of increasing the effective energy threshold of the array. Due to this dependence of the energy threshold, the low-energy threshold for the spectral analysis is set to a common value of 1 TeV.

The recorded images are first flat-fielded using information from nightly calibration runs taken with a pulsed UV LED (Hanna et al. 2010). The images are cleaned using a form of the picture/boundary method (Daniel 2008) and parameterized (Hillas 1985) to suppress the cosmic ray background. The shower direction is reconstructed from the data in each telescope, and a set of selection criteria is applied to reject background events (Konopelko 2001; Daniel 2008).

Energy spectra are calculated >1 TeV both for the FTW and the baseline observations and are shown in Figure 1. The spectra are parameterized as power laws of the form

$$\begin{equation} {dN \over dE} = N_0 \left({E \over {1\,{\rm TeV}}} \right)^{\gamma }. \end{equation} \tag{ 1 }$$

The baseline spectral fit gives a normalization of $N_0^{{\rm baseline}} = (3.48 \pm 0.14_{{\rm stat.}} \pm 1.08_{{\rm sys.}}) \times 10^{-7} \ {{\rm TeV}}^{-1} \ {\rm m}^{-2} \ {\rm s}^{-1}$ and γ^baseline = −2.65 ± 0.04_stat. ± 0.3_sys., with a χ² value of 16.6 with 12 degrees of freedom (dof). The FTW spectral fit gives a normalization of $N_0^{{\rm flare}} = (3.53 \pm 0.15_{{\rm stat.}} \pm 1.12_{{\rm sys.}}) \times 10^{-7} \ {\rm TeV}^{-1} \ {\rm m}^{-2} \ {\rm s}^{-1}$ and a spectral index γ^flare = −2.72 ± 0.05_sys. ± 0.3_sys., with a χ² value of 10.1 with 12 dof. The fit probabilities are 16% and 61%, respectively. These spectral parameters are mutually consistent, implying no change of the TeV flux during the FTW. The systematic uncertainties on the flux normalization and spectral index are expected to vary slowly with time, and a paper containing a proper treatment of these uncertainties is currently in preparation.

Zoom In Zoom Out Reset image size

2.2. Fermi-LAT

Fermi-LAT is a pair-conversion telescope sensitive to gamma-ray photons with energies between 20 MeV and 300 GeV. It has a wide FoV of ∼2.5 sr and surveys the entire sky every three hours. For a complete description of the instrument, see Atwood et al. (2009); Ackermann et al. (2012).

In order to extract spectral parameters of the Crab, the Fermi-LAT Science Tools v9r27p1 with P7V6 instrument response functions (IRFs) and the standard quality cuts described in Nolan et al. (2012) are used. Two years of "source"-class events with energies between 100 MeV and 300 GeV collected between MJD 54832 and 55562 within 20° of the Crab are processed with the maximum likelihood fitting routine. A model of the background is obtained in a binned likelihood analysis by fitting spectral models for all sources in the 2FGL catalog within 20° of the Crab in addition to the galactic and isotropic diffuse backgrounds (gal_2yearp7v6_v0.fits, iso_p7v6source.txt). Photon arrival times are barycentered with Tempo2 (Hobbs et al. 2006) using a publicly available Jodrell Bank radio ephemeris for the Crab pulsar (Lyne et al. 1993) to allow a selection of the off-pulse phase region 0.48–0.88. Under the assumption that emission from the pulsar is negligible in the off-pulse region, spectral parameters for the synchrotron and inverse-Compton components of the Crab Nebula are calculated. These parameters are fixed in the model to allow fitting of the pulsar spectral component after undoing the selection on pulsar phase.

The Crab Nebula synchrotron differential spectrum is parameterized as a power law (PowerLaw2 in the Fermi-LAT Science Tools) of the form

$$\begin{equation} {dN \over dE} = {F_0(\gamma + 1)E^{\gamma } \over ({300\,{\rm GeV}})^{\gamma + 1} - ({0.1\,{\rm GeV}})^{\gamma + 1}}. \end{equation} \tag{ 2 }$$

The fit of the quiescent state yields a synchrotron integral flux above 100 MeV of F₀ = (6.40 ± 0.11) × 10⁻⁷ cm⁻² s⁻¹ and photon index of γ = −3.69 ± 0.11, which are consistent with previously published results (e.g., Buehler et al. 2012).

A similar analysis is done for the FTW. Since Fermi-LAT carried out a targeted observation of the Crab during the flare, the recommended P7V6MC IRFs and pointed mode data selection criteria are used in this analysis.³⁰ The synchrotron integral flux above 100 MeV for the FTW is found to be (5.30 ± 0.13) × 10⁻⁶ cm⁻² s⁻¹ with a harder photon index of −3.10 ± 0.05. A combined spectral energy distribution (SED) showing the Fermi-LAT and VERITAS spectra is given in Figure 2.

Zoom In Zoom Out Reset image size

A test for variability in the VERITAS FTW light curve (shown in Figure 3) is performed by fitting the light curve with a constant flux. This fit gives a flux >1 TeV of (2.05 ± 0.07) × 10⁻⁷ m⁻² s⁻¹ with a χ² value of 19.1 with 9 dof (probability ∼2.4%). By fitting a light curve of data taken outside of the FTW, the Crab Nebula is detected with a baseline VHE flux >1 TeV of (2.10 ± 0.06) × 10⁻⁷ m⁻² s⁻¹ with a χ² value of 21.7 with 22 dof (probability ∼47.8%). The FTW flux is thus consistent with the baseline flux and with no statistically significant variability during the flare. An analysis of a subset of the data with energies extending down to ∼150 GeV was also conducted (shown in Figure 3), however no variability is revealed at these energies.

Zoom In Zoom Out Reset image size

In order to test for correlated Fermi-LAT and VERITAS (>1 TeV) flux variability in the light curves shown in Figure 3, a publicly available implementation of the z-transformed discrete correlation function (ZDCF) is employed (Alexander 1997, 2013). The ZDCF method requires a minimum of 12 observations in each light curve for a statistically valid analysis, so two nights of pre-flare VERITAS Crab Nebula observations taken on MJD 56331 and 56339 (February 8 and 16, respectively) are added before the cross-correlation is performed. The zero time-lag bin reported a ZDCF correlation coefficient of

$$\begin{equation} {\rm DCF} = -0.07 \pm 0.31 \end{equation} \tag{ 3 }$$

which is consistent with no correlation at zero lag. Results for all other time-lag bins are also consistent with no statistically significant correlation.

Relative flux changes during the FTW are calculated for VERITAS and Fermi-LAT. The i^th relative flux change $\Delta F_{{\rm rel.}}^{i}$ for both VERITAS and Fermi-LAT observations on the i^th night is computed as

$$\begin{equation} \Delta F_{{\rm rel.}}^{i} = \frac{ {F^i - \overline{F}}}{\overline{F}}. \end{equation} \tag{ 4 }$$

For VERITAS, Fⁱ is the average flux for one night. For Fermi-LAT, Fⁱ is the average flux in one 12 hr time bin centered on midnight Arizona time (MST, 0700 UTC). $\overline{F}$ is the average non-flare flux from the nebula. The VERITAS and Fermi-LAT relative flux changes for simultaneous observations are shown in Figure 4. Averaged over the simultaneous observations in the FTW, the relative flux changes are

$$\begin{eqnarray} &\overline{\Delta F^{{\rm VTS}}_{{\rm rel.}}} = -0.026 \pm 0.035 \quad {({\rm VERITAS}} > {1\,{\rm TeV})} \end{eqnarray} \tag{ 5 }$$

$$\begin{eqnarray} &\overline{\Delta F^{{\rm Fermi}}_{{\rm rel.}}} = 6.14 \pm 0.38 \quad {({\it Fermi}-{\rm LAT}} > {100\,{\rm MeV})}. \end{eqnarray} \tag{ 6 }$$

Zoom In Zoom Out Reset image size

From $\overline{\Delta F^{{\rm VTS}}_{{\rm rel.}}}$, a 95% confidence level upper limit (UL) is computed for an elevated VHE flux. Given the assumption of a positive and non-zero correlation of Fermi-LAT and VERITAS flux changes, a Bayesian prior is introduced in the limit calculation that is zero for negative relative flux changes and one elsewhere. This prior is equivalent to invoking the physical constraint that all of the VHE flux changes are at least zero. The UL is calculated over the Bayesian interval [0, x^up] such that

$$\begin{equation} {\frac{ \int _0^{x^{{\rm up}}}\ {\rm exp}\left({-\frac{(\overline{\Delta F^{{\rm VTS}}_{{\rm rel}}}-x)^2}{2\sigma ^2}}\right)dx}{ \int _0^{\infty }\ {\rm exp}\left({-\frac{(\overline{\Delta F^{{\rm VTS}}_{{\rm rel}}}-x^{\prime })^2}{2\sigma ^2}}\right)dx^{\prime }}} = 0.95 \ \ (x^{{\rm up}} > 0), \end{equation} \tag{ 7 }$$

where σ is the error on $\overline{\Delta F^{{\rm VTS}}_{{\rm rel.}}}$, and the 95% CL UL is given by x^up, which is obtained by solving the equation numerically. Limits are calculated for three different energy thresholds shown in Table 1.

Table 1. 95% CL Bayesian Upper Limits on the VHE Relative Flux Increase During the Flare Period for Three Energy Thresholds

Energy Band	$\overline{\Delta F^{{\rm VTS}}_{{\rm rel.}}} 95\%$ CL UL	95% CL Integral UL at Threshold
(TeV)		(TeV m−2 s−1)
>1	5.3%	8.7 × 10−9
>4	6.8%	5.9 × 10−9
>6	37.4%	2.7 × 10−8

Download table as:  ASCIITypeset image

By adopting the assumption that the relative flux change seen by VERITAS is linearly related to that seen by Fermi-LAT during the flare:

$$\begin{equation} \overline{\Delta F^{{\rm VTS}}_{{\rm rel.}}} = \alpha \overline{\Delta F^{{\rm {\it Fermi}}}_{{\rm rel.}}}, \end{equation} \tag{ 8 }$$

a constraint on the linear correlation factor α can be calculated, which can be used to test model predictions. Taking the ratio of the >1 TeV UL and the average Fermi-LAT relative flux change, we find that α < 8.6 × 10⁻³ (95% CL) for the average of the ten nights of simultaneous observations. The constraint on α is also computed night-by-night, though only MJD 56358 gives the slightly better constraint of α < 8.1 × 10⁻³(95% CL).

In this paper we present observations of the Crab Nebula with VERITAS and Fermi-LAT during the 2013 March flare. The light curve and reconstructed energy spectrum between 1 TeV and 10 TeV do not indicate any flux enhancement at TeV energies, while the flux above 100 MeV was six times elevated during our observations.

Earlier flares had very hard spectra with peak energy reaching up to _flare ≈ 500 MeV (Buehler et al. 2012), though in the present flare, a peak could not be resolved in the MeV–GeV spectrum leaving the electron spectrum unconstrained at lower energies. The synchrotron spectrum above 100 MeV is slightly harder than for the quiescent Crab, which may reflect a separate electron population and/or an increase in the magnetic-field strength in the emission zone that shifts a harder section of a curved synchrotron spectrum into the frequency band accessible with Fermi-LAT. Neglecting the weak modifications arising from the possibility of mildly relativistic bulk motion, we suggest that some excess electron acceleration took place.

From classical electrodynamics, the Lorentz factor of electrons that would emit 200 MeV synchrotron radiation is

$$\begin{equation} \gamma _\mathrm{sy}=3\times 10^9\, \left(\frac{B}{\mathrm{mG}}\right)^{-0.5} \ , \end{equation} \tag{ 9 }$$

and their energy-loss rate and life time are

$$\begin{eqnarray} \dot{E}_\mathrm{sy} &=& (8\times 10^{-3}\ \mathrm{}{\rm erg \,s}^{-1})\,\left(\frac{B}{\mathrm{mG}}\right),\nonumber\\ \tau _\mathrm{sy} &=& (3\times 10^5\ \mathrm{s})\,\left(\frac{B}{\mathrm{mG}}\right)^{-1.5}\ . \end{eqnarray} \tag{ 10 }$$

Assuming a magnetic field of 1 mG in the emission region, similar to that deduced in Bednarek & Idec (2011), the flare duration τ_sy is on the order of a few days, which is consistent with observed flares at a few hundred MeV. If the magnetic field were significantly stronger than 1 mG, the synchrotron lifetime would become very short compared to the flare duration, and so the electron population would need to be continuously replenished to sustain the flare. Thus, the main cause of the synchrotron flare was likely the injection of a large number of excess electrons at PeV energies.

Bednarek & Idec (2011) consider a model in which electrons are injected into the magnetic field of the pulsar wind zone and produce synchrotron gamma rays through acceleration in reconnection regions of the magnetic field. Assuming the electrons reach an equilibrium spectrum described by a differential power law with index between 3.0–3.6 and with a characteristic cut-off at γ = 3 × 10⁹ for flares, they suggest variability above ∼1 TeV of roughly 10% with more substantial changes above ∼10 TeV as a result of inverse-Compton scattering. However, inverse-Compton scattering of soft photons by electrons with Lorentz factors ∼10⁹ is heavily Klein-Nishina suppressed and would provide gamma rays in the PeV band, beyond the reach of VERITAS. Excess electrons with Lorentz factors of γ_IC ≃ 10⁷ may produce a flux enhancement at TeV energies, but the non-detection with VERITAS poses challenges for this model and thus constrains the number of electrons with Lorentz factors of γ_IC.

The number of electrons with Lorentz factors of ∼3 × 10⁹ can be estimated as

$$\begin{equation} N_\mathrm{e,sy} = \frac{L_\mathrm{sy}}{\dot{E}_\mathrm{sy}}\simeq 6\times 10^{37} \,\left(\frac{B}{\mathrm{mG}}\right)^{-1}\ , \end{equation} \tag{ 11 }$$

where L_sy is the synchrotron luminosity at 200 MeV. To calculate the number of electrons that may inverse-Compton scatter soft (infrared, IR) photons into the TeV band, we need to know the density of low-frequency radiation in the nebula. To this end we use L_soft ∼ 10³⁷ erg s⁻¹ (Marsden et al. 1984) as the pulsar wind nebula (PWN) luminosity in IR photons, _soft ∼ 0.1 eV as the photon energy, and d_PWN ≃ 1 pc as the characteristic size of the Crab Nebula. The density of IR photons is then

$$\begin{equation} n_\mathrm{soft} \simeq \frac{L_\mathrm{soft}}{4 \pi \, d_\mathrm{PWN}^2\, c \epsilon _{{\rm soft}}}\simeq 20\ \mathrm{cm^{-3}}. \end{equation} \tag{ 12 }$$

Using the UL on an extra flux component >1 TeV given in Table 1, we find that the inverse-Compton luminosity L_IC ≲ 4 × 10³² erg s⁻¹. The number of electrons that upscatter photons to TeV energies is given by

$$\begin{equation} N_\mathrm{e,IC} = \frac{L_\mathrm{IC}}{\sigma _\mathrm{T} n_\mathrm{soft} c \times \gamma m_e c^2}, \end{equation} \tag{ 13 }$$

where σ_T is the Thomson cross-section. Ignoring the moderate Klein–Nishina suppression (the kinematic parameter $4\,\epsilon _\mathrm{soft}\epsilon _\gamma / (m_e^2\,c^4)\simeq 10$), the UL derived on excess TeV gamma rays corresponds to at most

$$\begin{equation} N_\mathrm{e,IC} \left(\gamma \approx 10^7\right)\lesssim 10^{44}\ . \end{equation} \tag{ 14 }$$

Assuming for ease of exposition that the spectrum of excess electrons follows a power law, N_e(γ)∝γ^−s, the corresponding constraint on the spectral index is

$$\begin{equation} s\lesssim \frac{6.2+\log \left(\frac{B}{\mathrm{mG}}\right)}{2.5-\frac{1}{2}\log \left(\frac{B}{\mathrm{mG}}\right)}\ , \end{equation} \tag{ 15 }$$

which permits s ≃ 2.5 for the fiducial magnetic-field strength of 1 mG. This index is harder than that assumed by Bednarek & Idec (2011) and constrains the number of electrons that may be responsible for the Crab flare. Future observations with VERITAS or next-generation telescope arrays will likely provide more stringent constraints.

This research is supported by grants from the U.S. Department of Energy Office of Science, the U.S. National Science Foundation and the Smithsonian Institution, by NSERC in Canada, by Science Foundation Ireland (SFI 10/RFP/AST2748) and by STFC in the U.K. We acknowledge the excellent work of the technical support staff at the Fred Lawrence Whipple Observatory and at the collaborating institutions in the construction and operation of the instrument.

The authors are grateful to M. Mayer, E. Hays, and R. Buehler for a useful discussion of the Fermi-LAT data.

A. O'Faoláin de Bhróithe acknowledges the support of the Irish Research Council "Embark Initiative."

N.O. and G.R. gratefully acknowledge the support of a Fermi Guest Investigator grant.