FERMI LAT DISCOVERY OF EXTENDED GAMMA-RAY EMISSION IN THE DIRECTION OF SUPERNOVA REMNANT W51C

The origin of cosmic rays remains unsolved. The idea that supernovae produce cosmic rays (e.g., Hayakawa 1956) has been developed both in theoretical and observational aspects, so that galactic cosmic rays are widely considered to be accelerated in the expanding shocks of supernova remnants (SNRs) through a diffusive shock acceleration process (see, e.g., Blandford & Eichler 1987). This conjecture has been strengthened by recent observations of young SNRs in synchrotron X-rays and very high energy (VHE) gamma-rays (see, e.g., Reynolds 2008; Aharonian et al. 2008). High efficiency for converting the kinetic energy released by supernova explosions into the energy of relativistic protons and nuclei is a key requirement for explaining the galactic cosmic rays; it has yet to be confirmed observationally.

Enhanced π⁰-decay emission expected from SNRs that are interacting with molecular clouds could provide direct evidence of the nuclear component of cosmic rays being accelerated at supernova shocks (Aharonian et al. 1994). To identify the π⁰-decay gamma-rays as evidence of the nuclear cosmic rays, observations in the GeV domain are crucial because the π⁰-decay spectrum has a characteristic bump around 70 MeV. Previous measurements of GeV gamma-rays with EGRET onboard the Compton Gamma-Ray Observatory found some gamma-ray sources near radio-bright SNRs (Esposito et al. 1996). However, the possible origins of the EGRET sources, namely SNR shell contributions or pulsar associations, remained unclear, mainly due to poor localization.

The advent of the Large Area Telescope (LAT) onboard the Fermi Gamma-ray Space Telescope provides a new opportunity to study the gamma-ray emission from SNRs at GeV energies. An initial source list based on the first three months of Fermi observations (Abdo et al. 2009c) includes 0FGL J1923.0+1411, which is spatially coincident with SNR W51C. Even with the three months of data, the detection was at ∼23σ level. There are no EGRET counterpart(s) to the LAT source (Hartman et al. 1999); this is the first SNR discovered by Fermi as confirmed in this Letter.

W51C (G49.2 − 0.7) is a radio-bright SNR at a distance of D ≃ 6 kpc with an estimated age of ∼3 × 10⁴ yr (Koo et al. 1995). It has an elliptical shape with an extent of 50' × 38'. The radio continuum map exhibits thick shell-like structures (Subrahmanyan & Goss 1995). The bulk of the X-ray emission comes from thermal plasma with a temperature of kT ∼ 0.3 keV (Koo et al. 2005). A molecular cloud–shock interaction is known, as evidenced by observations of shocked atomic and molecular gases (Koo & Moon 1997a, 1997b). Most recently, the H.E.S.S. Collaboration has announced the detection of extended VHE gamma-ray emission coincident with W51C (Fiasson et al. 2009). Also, the Milagro Collaboration has reported a possible excess of multi-TeV gamma-rays in this direction (Abdo et al. 2009a).

In this Letter, we report the analysis results for the LAT source coincident with SNR W51C, using data accumulated over the first year of Fermi's operation. The Fermi observations of SNR W51C permit a refined study of cosmic-ray acceleration. Specifically, the LAT data suggest that π⁰-decay emission is the dominant contribution to the gamma-ray signal.

The Fermi Gamma-ray Space Telescope was launched on 2008 June 11 by a Delta II Heavy launch vehicle. The LAT onboard Fermi is a pair-conversion gamma-ray detector capable of measuring gamma-rays in a very wide range of energy from 20 MeV up to 300 GeV. The LAT tracks the electron and positron resulting from pair conversion of an incident gamma-ray in thin high-Z foils, and measures the energy deposition due to the subsequent electromagnetic shower that develops in the calorimeter. The effective area is ∼8000 cm² above 1 GeV (on-axis) and the per-photon 68% containment radius is ∼08 at 1 GeV. The point-spread function (PSF) depends largely on photon energy and improves at higher energies. The tracker of the LAT is divided into two regions, front and back. The front region (first 12 planes) has thin converters to optimize the PSF while the back region (four planes after the front section) has thicker converters to enlarge the effective area. The angular resolution for the back events is approximately twice as broad than that for the front events. The details of the LAT and data processing are given in Atwood et al. (2009), and the on-orbit calibration is described in Abdo et al. (2009b).

The gamma-ray data acquired from 2008 August 5 to 2009 July 14 are analyzed. The diffuse class events as defined in Atwood et al. (2009) are chosen for gamma-ray analysis. A cut on earth zenith angles greater than 105° is applied to reduce the residual signal from earth albedo gamma-rays. The instrument response functions (IRFs) called "Pass 6 V3" are used (Rando et al. 2009).

The maximum likelihood technique is employed for spectral and spatial parameter estimation using gtlike, which is publicly available as part of Fermi Science Tools.⁵⁸ The likelihood is the product of the probability of observing the gamma-ray counts of each spatial and energy bin given the emission model, and parameter values are estimated by maximizing the likelihood of the data given the model (Mattox et al. 1996). The gamma-ray emission model includes individual sources at fixed coordinates, galactic diffuse emission (resulting from cosmic-ray interactions with interstellar gas and photons), and an isotropic component (extragalactic and residual background). The so-called mapcube file of gll_iem_v02.fit is used for modeling the galactic diffuse emission, together with the corresponding tabulated model for the isotropic diffuse emission. Other versions of the galactic diffuse models, generated by the GALPROP code (Strong et al. 2004), are also utilized to assess systematic error. The maximum likelihood analysis is performed inside a square region of 12° × 12° centered on W51C with a pixel size of 01, unless otherwise mentioned. Background point sources detected in six months data are included in the likelihood analysis with free normalization and power-law index, though none of them affect the results in this Letter.

3.1. Spatial Distribution

In Figure 1, the maps of photon counts in the 2–10 GeV band in the vicinity of SNR W51C are shown; the right panel is a close-up view of the left panel. Gamma-ray events that converted in the front section of the tracker are selected. A bright gamma-ray source is enclosed by the outer boundary of W51C. The average surface brightness inside the SNR boundary is about 2 and 5 times larger than neighboring regions on the galactic plane in the 0.5–2 GeV and 2–10 GeV bands, respectively. The gamma-ray distribution appears to be somewhat clumpy, suggesting the presence of sub-structures. Due to the limited statistics, we defer the investigation of possible sub-structures to a future publication.

Zoom In Zoom Out Reset image size

The spatial distribution of the gamma-ray source is significantly extended compared to a simulated point source that has the same spectrum. A one-dimensional profile in the right ascension direction of the counts map is compared with that expected for a point source in Figure 2. Though the width of the PSF above 5 GeV is known to deviate from the Monte Carlo simulations, it has negligible effects on the simulated point source and other results in this Letter. Assuming a two-dimensional gaussian distribution fixed at the W51C centroid (α, δ) = (290818, 14145), we calculate the extension parameter of the source as σ_ext = 022 ± 002 by varying σ_ext to find the best match with the data. Note however that a different assumption on the spatial distribution results in different quantification of the spatial extension.

Zoom In Zoom Out Reset image size

The extended nature of this LAT source cannot be ascribed to a superposition of two point-like sources. To quantify this, we define a grid of 60 positions inside the remnant as trial point-source positions, and calculate maximum likelihood values (L_2ps) for each possible combination by running gtlike for a box of 8° × 8°. On the other hand, assuming that the gamma-ray source has uniform surface brightness inside the SNR boundary, we obtain another maximum likelihood value (L_uni). The resulting values of the likelihood test statistics, −2ln(L_2ps/L_uni)>32, demonstrate that the uniform emission gives a much better result.

The extension of the source precludes magnetospheric radiation from a pulsar as a dominant gamma-ray source. However, a small fraction (∼10%) of the gamma-ray flux could be contributed by a pulsar. Our pulse search results in non-detection, which places a 5σ upper limit on the pulsed gamma-rays as ≃6 × 10⁻⁸ photons cm^-2 s⁻¹ above 100 MeV (see Abdo et al. 2009d).

3.2. Spectrum and Its Uncertainties

The gamma-ray spectrum of the W51C source is shown in Figure 3. It is obtained by performing likelihood analysis for each energy bin with gtlike. The lower energy bound is set at 0.2 GeV, below which the systematic uncertainties become too large due to the background diffuse emission (see below). The surface brightness of the W51C source is assumed to be uniform inside the SNR boundary. The normalization of the galactic diffuse emission model is left free at each energy bin. As evident from Figure 3, the gamma-ray spectrum cannot be described by a simple power law and steepens above a few GeV. A likelihood-ratio test indicates that a power-law hypothesis has a chance probability of 5 × 10⁻⁵ for obtaining the spectral data. The physical interpretation of the spectral energy distribution will be discussed in Section 4.

Zoom In Zoom Out Reset image size

The gamma-ray luminosity can be estimated as ≃1 × 10³⁶(D/6 kpc)² erg s^-1 in the LAT bandpass, making this object one of the most luminous gamma-ray sources in our Galaxy. Note that the distance to the remnant is well constrained. X-ray absorption (Koo et al. 1995) indicates that W51C is situated behind the W51 molecular cloud complex (a so-called 68 km s^-1 cloud in particular) whose distance is determined as 7.0 ± 1.5 kpc by the proper motion of W51 Main H₂O maser (Genzel et al. 1981). On the other hand, the angular extent and the X-ray measurements (Koo et al. 2002) require D ≲ 10 kpc. We adopt D = 6 kpc following previous publications (e.g., Koo et al. 2005).

It should be noted that in addition to the systematic errors commonly assigned to LAT spectral data,⁵⁹ the accuracy of the flux estimated for each energy bin of the W51C source is limited by the following errors. First, a possible effect of the uncertainty of underlying galactic diffuse emission is considered. The observed gamma-ray intensity of nearby source-free regions on the galactic plane is compared with the intensity expected from the galactic diffuse models. The difference, namely the local departure from the best-fit diffuse model, is found to be ≲6%. By changing the normalization of the galactic diffuse model artificially by ±6%, we estimate the possible systematic error to be 40% (0.2–0.4 GeV), 22% (0.4–0.8 GeV), and <10% (>0.8 GeV).

The fact that we do not know the true gamma-ray morphology of the W51C source introduces another error in our flux estimation. Since the gamma-ray source is point like below 1 GeV given the PSF, this uncertainty matters only above 1 GeV. We adopt a uniform surface brightness inside the SNR boundary (a flat elliptical template) as the spatial distribution of the source gamma-rays. Different spatial distributions such as a flat elliptical template reduced in size (scaled by 0.5) are tested to estimate the systematic error. Our conservative estimate is ≲20% in 1–6 GeV and ∼30% above 6 GeV as the systematic uncertainty attributable to the unknown shape of the source.

The extended gamma-ray emission positionally coincident with SNR W51C has been studied using the Fermi LAT. The gamma-ray spectrum presented in Figure 3 is not fitted by a simple power law, exhibiting a remarkable steepening. Here, we discuss the origin of the extended emission and the underlying particle spectra that give rise to the observed spectrum of photons.

The expanding shock waves driven by the supernova explosion are expected to be the sites of the acceleration of multi-GeV particles. To phenomenologically interpret the spectral curvature in the LAT+TeV bands, a broken power law is adopted for the momentum distribution of the radiating electrons/protons

$$\begin{equation} N_{{\rm e,p}}(p) = a_{{\rm e,p}} \left(\frac{p}{p_0}\right)^{-s} \left(1+ \left(\frac{p}{p_{\rm br}}\right)^{2} \right)^{-\Delta s/2}, \end{equation} \tag{ 1 }$$

where p₀ = 1 GeV c⁻¹. For simplicity, the indices and the break momentum are assumed to be identical for both accelerated protons and electrons. As we argue below, the break may reflect the character of magnetohydrodynamic turbulence. To account for the radio synchrotron index α ≃ 0.26 (Moon & Koo 1994), we adopt s = 1.5, though a steeper index (say, s = 1.7) could be reconciled with the radio observations within the uncertainty. The energetic particles are assumed to be uniformly distributed in the volume of V = (4π/3)R³_eff with an effective radius R_eff = 30 pc. The age and radius imply a shock velocity of v_sh ∼ 400 km s^-1 and E_SN/n₀ ∼ 1.6 × 10⁵² erg cm³, where E_SN and n₀ represent the explosion kinetic energy and the interstellar density into which the main blast wave is propagating, respectively. The X-ray data suggest n₀ ∼ 0.3 cm^-3. The radio images indicate that radio-emitting electrons are smoothly distributed in a thick shell. The model of Koo & Moon (1997a) suggests the presence of a molecular cloud of ∼1 × 10⁴ M_☉ engulfed by the blast wave. The engulfed cloud can act as target material for relativistic particles. The total (atomic and molecular) hydrogen mass contained in the volume is denoted by $M_{\rm H}= \bar{n}_{\rm H} m_{\rm p} V$. Note that $\bar{n}_{\rm H} \gg n_0$ can be expected.

Given the interaction with a molecular cloud, we first attribute the observed gamma-rays to the decay of π⁰ mesons produced in inelastic collisions between accelerated protons (and nuclei) and target gas (Figure 3). The gamma-ray spectrum of π⁰ decay is calculated based on Kamae et al. (2006) using a scaling factor of 1.85 for helium and heavy nuclei (Mori 2009). Note that the scaling factor assumes the local interstellar abundance for target material and the observed cosmic-ray composition. Contributions from bremsstrahlung and inverse Compton (IC) scattering by accelerated electrons are also shown in Figure 3. Electron–ion and electron–electron bremsstrahlung spectra are computed as in Baring et al. (1999). The interstellar radiation field for IC (see Table 1) is comprised of two diluted blackbody components (infrared and optical) and the cosmic microwave background (CMB). The infrared and optical components are adjusted to reproduce the interstellar radiation field in the GALPROP code (Porter et al. 2008). Cooling effects due to ionization and synchrotron (or IC) losses, which introduce cooling breaks in particle spectra in addition to p_br, are taken into account assuming constant particle injection over a period T₀ ∼ 3 × 10⁴ yr. The synchrotron cooling becomes important in the TeV band for leptonic models.

Table 1. Parameters of Multiwavelength Models

	Parameters					Energetics
Model	ae/ap	Δs	pbr (GeV c−1)	B (μG)	$\bar{n}_{\rm H}$ (cm−3)	Wp (1050 erg)	We (1050 erg)
(a) π0 decay	0.02	1.4	15	40	10	5.2	0.13
(b) Bremsstrahlung	1.0	1.4	5	15	10	0.54	0.87
(c) Inverse Compton	1.0	2.3	20	2	0.1	8.4	11

Notes. Seed photons for IC include the CMB (kT_CMB = 2.3 × 10⁻⁴ eV, U_CMB = 0.26 eV cm⁻³), infrared (kT_IR = 3 × 10⁻³ eV, U_IR = 0.90 eV cm⁻³), and optical (kT_opt = 0.25 eV, U_opt = 0.84 eV cm⁻³). The total energy content of radiating particles, W_e,p, is calculated for p > 10 MeV c⁻¹.

Download table as:  ASCIITypeset image

Figure 4(a) shows the radio+gamma-ray spectrum together with the radiation model that uses the parameters in Table 1. We adopt here M_H = 2.8 × 10⁴ M_☉ $(\bar{n}_{\rm H}=10\ {\rm cm^{-3}})$, which is somewhat larger than the value quoted above. The total energy of the high-energy protons amounts to W_p = 5.2 × 10⁵⁰ erg, which is inversely proportional to M_H, but insensitive to other parameters. The large luminosity of the LAT source can be explained naturally by a large number of accelerated protons and dense environments as expected for the case of W51C.

Zoom In Zoom Out Reset image size

The break feature in the proton spectrum, which is introduced on phenomenological grounds, is not expected in a shock acceleration zone if acceleration proceeds close to the Bohm limit (e.g., Baring et al. 1999). Protons can be accelerated up to ∼200 TeV via diffusive shock acceleration with a power-law (or even slightly concave) momentum distribution. A possible explanation for the falling proton spectrum above p_br = 10–30 GeV c⁻¹ (which depends on the choice of other parameters) could be the effects of damping of magnetohydrodynamic turbulence due to ion–neutral collisions. According to Ptuskin & Zirakashvili (2003), the maximum attainable momentum in the presence of wave dissipation by ion–neutral collisions is estimated as p_max(T₀) ∼ 30(n/1 cm^-3)^−1/2 GeV c⁻¹, using the SNR parameters described above. Here, n is the ambient neutral gas density. The shock precursor cannot confine accelerated protons with p > p_max(T₀) even if they were accelerated earlier. Therefore, energy-dependent escape of accelerated particles from the remnant (e.g., Aharonian & Atoyan 1996) may account for the steepening of the proton distribution. In any case, a break is needed to fit the observed spectral energy distribution (SED).

In contrast, leptonic models for the GeV gamma-ray production face difficulties (see Figures 4(b) and (c), and Table 1). While the chosen model parameters are not unique in their ability to fit the broadband spectrum, they are representative; no reasonable choices for them can easily account for the radio+gamma-ray data. Leptonic scenarios require a_e/a_p far in excess of cosmic-ray abundance ratios. It is difficult to reconcile the bremsstrahlung-dominated model with the observed radio synchrotron spectrum (Figure 4(b)). Also, the IC-dominated model (Figure 4(c)) requires an unrealistically large energy content in radiating electrons, W_e ≃ 1 × 10⁵¹ erg, a low magnetic field of B ≃ 2 μG (to suppress the radio synchrotron emission), and a low density of $\bar{n}_{\rm H} < 0.1\ {\rm cm^{-3}}$ (to reduce the bremsstrahlung component). Such a low density is problematic because even X-ray-emitting gas has a density of ∼1 cm^-3 (Koo et al. 2002). It is also difficult to account for the LAT signal by an IC process in a possible pulsar wind nebula seen inside W51C, CXO J192318.5+143035 (Koo et al. 2005), which has only a few parsec extent and much weaker radio emission compared with the SNR shell.

The most plausible explanation for the LAT source is therefore offered by π⁰-decay gamma-rays in a dense environment, spawned by efficient acceleration of protons and nuclei taking place or having taken place at the shocked shell of SNR W51C. The discovery of the GeV-scale gamma-rays presented in this Letter raises the intriguing possibility that Fermi has discerned accelerated ions in a middle-aged remnant that is interacting with dense clouds.

The Fermi LAT Collaboration acknowledges support from a number of agencies and institutes for both development and the operation of the LAT as well as scientific data analysis. These include NASA and DOE in the United States, CEA/Irfu and IN2P3/CNRS in France, ASI and INFN in Italy, MEXT, KEK, and JAXA in Japan, and the K. A. Wallenberg Foundation, the Swedish Research Council, and the National Space Board in Sweden. Additional support from INAF in Italy and CNES in France for science analysis during the operations phase is also gratefully acknowledged.