REVEALING W51C AS A COSMIC RAY SOURCE USING FERMI-LAT DATA

Supernova remnants (SNRs) are widely believed to be the sources of Galactic cosmic rays ($E\lt {10}^{15}{\rm{}}\;{\rm{eV}}$). Verification of this hypothesis is possible by studying the γ-ray emission from SNRs. Accelerated cosmic rays interact with surrounding matter and produce ${\pi }^{0}$mesons that subsequently decay into γ-rays, leading to a characteristic break in the γ-ray spectrum at about $E\;\sim$ 200–300 MeV due to the finite rest mass of the ${\pi }^{0}$. Competing leptonic γ-ray production mechanisms such as bremsstrahlung and inverse Compton (IC) emission require strong fine tuning to produce a similar feature.³ The pion bump is thus an unmistakable feature of the hadronic origin of γ-ray emission that is observable by the Fermi-Large Area Telescope (LAT). This feature was recently found in the two brightest SNRs detected by Fermi-LAT, IC 443 and W44, as reported in Ackermann et al. (2013). The third brightest SNR detected by Fermi-LAT, W51C, is another prime candidate to search for the pion-decay signature in the γ-ray emission because it is interacting with a molecular cloud (MC) that poses an excellent target for cosmic-ray interactions and subsequent ${\pi }^{0}$ decay.

W51C is an SNR identified by its radio shell and assumed to be 30 kyr old and at a distance of 5.5 kpc (Sato et al. 2010). The SNR belongs to the larger W51 complex composed of one of the largest star-forming regions in our Galaxy and divided into two parts denoted W51A and W51B. The SNR is located toward the south-eastern end of W51B and evidence for interaction between the SNR shell and the MC W51B is provided by the discovery of two OH (1720 MHz) masers by Green et al. (1997). Further evidence is provided by high-velocity atomic gas that shows velocity shifts between 20 and $120{{\rm{}}{\rm{km}}{\rm{s}}}^{-1}$ with respect to the ambient medium (Koo & Moon 1997). The OH masers and the high-velocity clouds are commonly interpreted as the result of shock waves penetrating into MCs. More recent measurements by Ceccarelli et al. (2011) find an overabundance of overionized gas in certain W51B locations close to W51C that they interpret as the result of ionization through by-products of freshly accelerated cosmic rays interacting with nucleons. All these measurements outline a compelling scenario in which the SNR shell of W51C is interacting with the MC W51B.

X-ray imaging of the region revealed a hard source denoted CXO J192318.5+140305, which was identified as a possible pulsar wind nebula (PWN) related to the SNR (Koo et al. 2002, 2005 ). Such a PWN could also power relativistic particle acceleration resulting in γ-ray emission.

W51C was detected in 11 months of data by the Fermi-LAT and showed spatially extended γ-ray emission at $E\gt 2{\rm{}}\;{\rm{GeV}}.$ The γ-ray spectrum showed curvature that required a break at a few GeV, and subsequent modeling by a hadronic spectrum required a broken power law for the proton spectrum with a break momentum of 10–15 GeV/c (Abdo et al. 2009). At very high energies ($E\gt 100{\rm{}}\;{\rm{GeV}},$ VHE) an energy-dependent morphology was revealed by the MAGIC Collaboration (Aleksić et al. 2012). No hint of spectral variation could be identified and the emission is attributed to the SNR W51C due to its morphology and multi-wavelength modeling of the spectral energy distribution (SED). The authors mention that a contribution of up to 20% of the W51C flux could come from the PWN candidate CXO J192318.5+140305. A final determination of the parent particle population that creates the γ-rays was not possible, though multi-wavelength interpretation of the SED of W51C favors a hadronic origin of the γ-ray emission.

Here we analyze five years of Fermi-LAT data and search for features that can identify the parent particle type in the most detailed γ-ray spectrum of W51C derived to date.

The LAT, the primary instrument of the Fermi Gamma-ray Space Telescope mission, is a pair-conversion telescope sensitive to γ-rays in an energy range from $20\;{\rm{MeV}}\;$ to $E\gt 300\;{\rm{GeV}}$. The reader is referred to Atwood et al. (2009) for a more detailed description and to Ackermann et al. (2012) for the on-orbit performance. We analyzed public Fermi-LAT data between MJD 54682.7 and MJD 56516.5 corresponding to about five years of Pass 7 Reprocessed data. The data are analyzed using the Fermi ScienceTools version v9r32p04.⁴ We select events with high probability of being γ-rays by choosing the SOURCE event class. In order to evade γ-ray contamination generated by cosmic rays hitting the Earth's atmosphere, we remove time intervals when the field of view of the LAT came too close to the Earth's limb (zenith angle $\lt 100^\circ$). We analyze data between 60 MeV and 300 GeV and use the P7REP_SOURCE_V15 instrument response functions (IRFs). To account for the interstellar γ-ray emission caused by cosmic rays interacting with gas or interstellar radiation in our Galaxy, the model gll_iem_v05_rev1.fit is used. This interstellar emission model (IEM) is our standard IEM provided by the Fermi-LAT Collaboration for point-source analysis. The isotropic γ-ray emission is accounted for by the model iso_source_v05.txt, which also includes any residual charged-particle background present in the Fermi-LAT data. We chose a $20^\circ \times 20^\circ$ region of interest (ROI) centered on W51C and performed a binned likelihood analysis with $0\buildrel{\circ}\over{.} 1$ bins and 30 bins logarithmically spaced in energy.

We construct a test statistic (TS) following Mattox et al. (1996) to evaluate the improvement of the likelihood fit to the ROI when adding a new source. In the case of one additional source with one additional free parameter we can define $\mathrm{TS}=2(L-{L}_{0})$ and the significance as $\sigma =\sqrt{\mathrm{TS}}$ where L₀ and L are the log-likelihood values without the source and with it respectively.

In order to study the effect of nearby sources on the spectral fit of W51C we use gtobssim³ to obtain Monte Carlo simulations of the W51C region. In our simulations we assume spectral values for all sources as given in the 3FGL catalog (Acero et al. 2015a). We vary all independent spectral parameters randomly by either $+2\sigma$ or $-2\sigma$ of the associated error. The $\pm 2\sigma$ is chosen to provide an estimation in the case of a strong deviation of the real spectral value from that in the 3FGL catalog and at the same time to require running only a few simulations. We run ten variations of the simulation so that different combinations of sources have positive and negative fluctuations. Afterwards we analyze these ten Monte Carlo simulation samples with varying free parameters for the likelihood fit of the background sources. Finally, we compare the fit results of W51C with the simulated values for each analysis. We find that we have only to leave free parameters of the five sources marked with cyan crosses in Figure 1 in the likelihood fit to obtain within statistical uncertainties the simulated parameter values of W51C. For sources with $\mathrm{TS}\lt 50$ in the 3FGL catalog we leave only the normalization free, and for more significant sources all other independent spectral parameters too. In addition to the sources in Figure 1 we also leave free the normalizations of the IEM and the isotropic emission template in our analysis.

Zoom In Zoom Out Reset image size

To estimate the systematic uncertainties on flux and spectral properties of W51C caused by our limited knowledge of the IEM we compare the standard IEM results with the results obtained with eight alternative models generated with GALPROP (Vladimirov et al. 2011) and afterwards refined by fitting to the Fermi-LAT data. Each of these models consists of eight map cubes (i.e., three-dimensional models of gamma-ray intensity as a function of position and energy) inferring the interstellar γ-ray emission from gas by tracing it with H i and CO maps. The H i and CO maps are divided into the contributions from four galactocentric rings. In addition to these contributions there are map cubes added for large residual structures in the Fermi-LAT γ-ray sky, namely Loop I (Casandjian et al. 2009), the Fermi bubbles (Su et al. 2010; Ackermann et al. 2014), and an IC emission map. For a more detailed explanation of the alternative models see Acero et al. (2015b). Inside the W51C ROI only the H i rings 2–4, the CO rings 2–4, the IC map, and the Loop I template give any background contribution and all other components are removed from our fits. In addition to these components an IEM-specific isotropic component is added.

In our IEM systematic uncertainty study the normalizations of the individual components of each IEM are free in the likelihood fit. The alternative IEMs provide more freedom to the fit by leaving the components individually adjustable instead of only one global normalization, as in the case of the standard IEM .

W51C was already reported to be an extended source for the Fermi-LAT in Abdo et al. (2009). We use the pointlike code (Kerr et al. 2011) to verify whether our data set, which is more than five times larger than that used in Abdo et al. (2009), would allow us to further constrain the morphology of W51C. Therefore we divide the γ-ray emission into energy bands from 60 MeV to 1 GeV, 1–5 GeV, and 5–300 GeV. The energy ranges are chosen to separate the individual spectral regimes present in the W51C γ-ray spectrum (see Section 4) and simultaneously include high photon statistics in each bin. No solid evidence for changing morphology between the tested energy bands is found. Slight differences in the position, geometry, and shape of the γ-ray excess that are seen between the individual energy bands have no impact on the obtained photon spectrum and thus are negligible for the spectral analysis of W51C. In our study W51C is modeled by a flat elliptical disk with five free parameters (R.A., decl., minor axis, major axis, rotation angle). In the best-position fit we use the best-fit spectral values obtained with the previously published template⁶ and afterwards verify that leaving the spectral parameters free in the position fit does not alter the spatial shape or position. A fit using a two-dimensional Gaussian instead of the elliptical disk does not result in an improvement of the overall description of the region and has only insignificant effects on the W51C fit results. Our best-fit morphology for the overall energy range is in good agreement with the previously published template of W51C and we use this publicly available template for all our studies.

Comparing the morphology of W51C seen by the Fermi-LAT with that obtained by MAGIC shows a spatial coincidence as seen in Figure 2. The MAGIC contours are not an exact match to the Fermi-LAT template but they completely overlap. The difference between the two shapes might be explained by the lower angular resolution of the Fermi-LAT, making it impossible as shown in Lande et al. (2012) to distinguish between very similar source shapes. We conclude from these comparisons that the changes in energy-dependent morphology seen by MAGIC cannot be detected at $E\lt 75{\rm{}}\;{\rm{GeV}}$ by Fermi-LAT.

Zoom In Zoom Out Reset image size

We fit the W51C flux in 20 logarithmic energy bins independently between 60 MeV and 300 GeV. The resulting SED is shown in Figure 3 and a clear low-energy break around 300 MeV is visible. The SEDs obtained using the alternative IEMs differ by more than the statistical errors only at energies below a few GeV. All of them indicate a low-energy break at a very similar energy to the standard IEM. Above about 1 GeV there are only very small differences between the SEDs derived with the alternative IEMs and the standard one, as expected since the diffuse emission from cosmic-ray interactions decreases more rapidly than the source spectrum.

Zoom In Zoom Out Reset image size

The significance of the break is obtained by fitting W51C using the likelihood approach between 60 MeV and 3 GeV once with a smoothly broken power law of the form $F{(E)={N}_{0}(E/{E}_{0})}^{{{\rm{\Gamma }}}_{1}}{[1+{(E/{E}_{{\rm{break}}})}^{({{\rm{\Gamma }}}_{1}-{{\rm{\Gamma }}}_{2})/\alpha }]}^{-\alpha },$ with E₀ = 200 MeV and $\alpha =0.1$, and once with a power law of the form $F{(E)={N}_{0}(E/{E}_{0})}^{-{\rm{\Gamma }}}.$ The improvement in likelihood when fitting a smoothly broken power law and assuming a nested model with two additional degrees of freedom yields a $9\sigma$ significance for a spectral break at ${E}_{{\rm{break}}}=290\pm 20\;{\rm{}}\;{\rm{MeV}}$ and change of spectral index from ${{\rm{\Gamma }}}_{1}={0.70}_{-0.28}^{+0.23}$ to ${{\rm{\Gamma }}}_{2}=2.15\pm 0.03.$ To estimate the systematic uncertainty of the break energy due to our limited knowledge of the Fermi-LAT IRFs the bracketing IRFs method is used (Ackermann et al. 2012). Here the worst-case scenario, in which the IRFs change maximally at the break energy, is assumed to provide conservative systematic uncertainties. Additionally we apply the same method using the eight alternative IEMs and obtain very similar results for ${E}_{{\rm{break}}}.$ The alternative IEM fits allow us to estimate the systemaitc uncertainty of the break energy by calculating the variance with respect to the standard IEM (Acero et al. 2015b). Furthermore, the smoothly broken power law is favored above a power law in all alternative models and yields at least a significance of $8.7\sigma$, rendering the break significant in all IEMs tested. The variation of the IEMs is currently our best method to estimate uncertainties associated with our limited knowledge of the IEM but cannot provide a complete coverage of all possible deviations. Taking into account the aforementioned uncertainties, we obtain a clear low-energy break in the W51C spectrum at ${E}_{{\rm{break}}}=290\pm 20({\rm{stat}})\pm 40({\rm{syst,}}\;{\rm{IEM}})\pm 30({\rm{syst,}}\;{\rm{IRF}})$ MeV. In addition we also test for the effect of energy dispersion that results in general in an overprediction of the lowest flux point in the SED by about $30\%.$ Similar to the study in Ackermann et al. (2013) we find that neglecting energy dispersion in our analysis leads to a lower significance of the detected low-energy break. The break in W51C shows very similar properties (${E}_{{\rm{break}}},$ ${{\rm{\Gamma }}}_{1},$ ${{\rm{\Gamma }}}_{2}$) to those of the SNRs W44 and IC 443 that were already identified convincingly as cosmic-ray accelerators (Ackermann et al. 2013).

Closer inspection of the W51C SED suggests that there is also a high-energy break in the spectrum. Indeed such a break is expected since the spectral index above the low-energy break, ${{\rm{\Gamma }}}_{2}=2.15\pm 0.03$, is not compatible with that obtained by the MAGIC Collaboration at energies beyond 75 GeV, ${\rm{\Gamma }}=2.58\pm 0.07$ (Aleksić et al. 2012). Applying the same likelihood approach as for the low-energy break but using the energy range 400 MeV to 300 GeV to compare between smoothly broken power law and simple power law, we find a break at ${E}_{{\rm{break,HE}}}={2.7}_{-0.8}^{+1.0}({\rm{stat}})$ GeV with $7.5\;\sigma$ significance. Systematic uncertainties due to the IEM are unimportant at these energies and, given the large statistical uncertainty, all systematic errors are negligible in comparison. The power-law spectral index changes from ${{\rm{\Gamma }}}_{1}={2.11}_{-0.05}^{+0.06}$ in the low-energy regime to ${{\rm{\Gamma }}}_{2}={2.52}_{-0.07}^{+0.06}$ above the break. The high-energy spectral index is in excellent agreement with the spectral index reported by the MAGIC Collaboration at energies above 75 GeV, and thus the MAGIC spectrum can be well explained by the same particle population as the Fermi emission above ${E}_{{\rm{break,HE}}}.$ In the next section the origin of the γ-ray emission is tested with hadronic and leptonic scenarios.

The obtained W51C SED is fitted by a proton-induced and alternatively by an electron-induced γ-ray spectrum. For the case of the electron-induced spectrum we consider IC- and bremsstrahlung-dominated γ-ray production. The γ-ray spectra for the leptonic cases are calculated according to the cross sections given in Blumenthal & Gould (1970) and Ellison et al. (1999) for electron–electron bremsstrahlung. The IC seed photons are provided by the cosmic microwave background (CMB), the emission of infrared light, mainly from dust, and starlight. Each of these seed photon fields is modeled by a blackbody spectrum with ${{kT}}_{{\rm{CMB}}}=2.3\times {10}^{-4}{\rm{}}\;{\rm{eV}},$ ${u}_{{\rm{CMB}}}=0.26\;{{\rm{}}{\rm{eV}}{\rm{cm}}}^{-3}$ for the CMB, ${{kT}}_{{\rm{IR}}}=3\times {10}^{-3}{\rm{}}\;{\rm{eV}},$ ${u}_{{\rm{IR}}}=0.90\;{{\rm{}}{\rm{eV}}{\rm{cm}}}^{-3}$ for infrared light, and ${{kT}}_{{\rm{Star}}}=0.25\;{\rm{}}\;{\rm{eV}},$ ${u}_{{\rm{Star}}}=0.84{{\rm{}}{\rm{eV}}{\rm{cm}}}^{-3}$ for starlight. The values are taken from Abdo et al. (2009). We also consider cooling effects due to synchrotron emission for electrons above 1 TeV, for which we assume a constant injection. Other cooling channels are unimportant for the assumed acceleration lifetime of the SNR (30 kyr) and are thus neglected. The proton–proton cross sections are taken from Kamae et al. (2006) and we multiply them by a factor of 1.85 to account for heavier element abundance as suggested by Mori et al. (2009). Secondary spectra from the decay of charged pions are computed as well but found to contribute insignificantly to the total γ-ray emission and hence are neglected.

For all tested processes the parent particle spectrum is modeled by a broken power law with an exponential cut-off in momentum space. A comparison with a simple power law with exponential cut-off yields much worse results for bremsstrahlung and IC processes and a slightly worse description in the case of the hadronic spectrum. We note that we use only statistical errors when fitting the SED.

The best-fit spectra are shown in Figure 4 and the only adequate description of the Fermi-LAT data is provided by the proton-induced emission spectrum. Including the MAGIC data in the fit yields again a valid description only with the proton spectrum with ${({\chi }^{2}/\mathrm{dof})}_{{\rm{pp}}}=26.3/22,$ whereas the ${\chi }^{2}$ values of the IC ${({\chi }^{2}/\mathrm{dof})}_{{\rm{IC}}}=484/23$ and bremsstrahlung spectra ${({\chi }^{2}/\mathrm{dof})}_{{\rm{Brems}}}=98.5/23$ indicate an insufficient description of the data. More complicated spectral shapes that allow for another low-energy break in the lepton spectrum might be able to explain the W51C spectrum as due to bremsstrahlung processes but again that would require strong fine tuning. Therefore we find the hadronic explanation of the γ-ray spectrum to be the most convincing model.

Zoom In Zoom Out Reset image size

We also fitted the SEDs obtained using the alternative IEMs to assess the modeling dependence of the dominant systematic uncertainties with the aforementioned processes. The proton–proton mechanism is always significantly favored above the bremsstrahlung description.

The proton spectrum has a rather soft spectral index below the break of ${{\rm{\Gamma }}}_{1}=2.48\pm 0.02$ compared to the ${\rm{\Gamma }}=2.0$ or harder predicted by diffusive shock acceleration (DSA). However, we note that fixing the spectral index to ${\rm{\Gamma }}=2.0$ results in a worse but perhaps still acceptable value ${\chi }^{2}/\mathrm{dof}=47.8/23$. Also the two other SNRs, IC 443 and W44, that are identified as proton accelerators show rather soft spectral indices (${\rm{\Gamma }}=2.36$) very similar to W51C.

We determine the break momentum to be ${p}_{{\rm{break}}}={81.3}_{-5.4}^{+5.8}\;{\rm{GeV}}/c$ and the high-energy spectral index to be ${{\rm{\Gamma }}}_{2}=2.73\pm 0.03.$ Our fit does not allow for a determination of the exponential cut-off; we find by profiling the cut-off momentum that any value between ${p}_{{\rm{cut}}}=30\;{\rm{TeV}}$ and ${p}_{{\rm{cut}}}=300\;{\rm{TeV}}$ yields ${\rm{\Delta }}{\chi }^{2}\lt 3.7.$

We establish a low-energy break at ${E}_{{\rm{break}}}=290\pm -20({\rm{stat}})\pm {40}_{{\rm{syst,IEM}}}\pm {30}_{{\rm{syst,IRF}}}$ MeV in the γ-ray spectrum of W51C. This break is associated with the energy threshold of ${\pi }^{0}$production and hence is strong evidence of the hadronic origin of the γ-ray emission. Furthermore, we establish a second break in the photon spectrum at ${E}_{{\rm{break,HE}}}={2.7}_{-0.8}^{1.0}({\rm{stat}})$ GeV that smoothly connects the Fermi-LAT spectrum and the spectrum obtained by the MAGIC Collaboration above 75 GeV. The best description of the combined Fermi-LAT and MAGIC SED is given by a proton–proton-induced γ-ray spectrum with a broken-power-law spectrum in momentum space of the parent particles.

These results make W51C the third definitely identified SNR accelerating cosmic rays. As with the other two SNRs, IC 443 and W44, W51C is a middle-aged SNR that is not expected to contain the highest-energy cosmic rays it accelerated in its youth. All three of these objects require a broken power-law momentum spectrum of the protons. The break in the momentum spectrum is not connected to the ${\pi }^{0}$break in the γ-ray spectrum but is required to describe the high-energy part of W51C's γ-ray spectrum. The nature of the break in the parent particle population is unknown, but since there is no energy-dependent morphology found in the Fermi-LAT energy range⁷ two independent proton particle populations seem unlikely to explain the break. According to models of cosmic-ray diffusion (Aharonian & Atoyan 1996; Gabici et al. 2009), the highest-energy cosmic rays start to escape from middle-aged SNRs, but this scenario results in an exponential cut-off in the parent particle spectrum and not a break. Such an exponential cut-off is not sufficient to model the W51C data or those of the other two SNRs. Another possible explanation for a break in the proton spectrum is a modification of the shock acceleration due to the interaction with the surrounding MC, maybe by neutral ion damping as suggested in Malkov et al. (2011) for the case of W44. Another possibility might be that emission is generated through reaccelerating existing cosmic rays in shocks in MCs triggered by SNR blast waves (Uchiyama et al. 2010). In this scenario a break could be observed as an overlay of different acceleration zones in the MC that have varying densities (Uchiyama et al. 2010). In the case of W51C there is evidence for interaction between the SNR shell and adjacent MCs from multi-wavelength data. Also in W44 there is evidence for interaction with surrounding MCs (Uchiyama et al. 2012). On the other hand IC 443 is most likely not directly interacting with an MC and there is some diffusion required before the cosmic rays can interact with the MCs in the vicinity of the SNR shock; hence models for escaping cosmic rays have been developed (Torres et al. 2010). The cause of the rather soft proton spectrum and the momentum break at a few to hundreds of ${\rm{GeV}}/c$ in these three SNRs cannot be conclusively explained at present. The differences in the environments of these three SNRs might explain the differences in the break momentum and the change in the spectral index, but our knowledge of the exact environmental parameters is limited. Interestingly, in all three cases an additional exponential cut-off in the γ-ray emission and in the parent particles is not detected. Therefore the highest-energy protons still accelerated by the SNR (or reaccelerated in the MCs) cannot be determined but must exceed about 30 TeV. Finding such high-energy protons in a relatively old SNR such as W51C is surprising, and future measurements with imaging Cherenkov telescopes, like the Cherenkov Telescope Array, will reveal the acceleration capabilities of W51C. Also more detailed morphology studies of W51C will help to provide further insight into what part of the MC–SNR interaction leads to modification of the cosmic-ray spectrum around the high-energy break and provide vital input to the development of theoretical models that will be able to explain the current data.

The Fermi-LAT Collaboration acknowledges generous ongoing support from a number of agencies and institutes that have supported both the development and the operation of the LAT as well as scientific data analysis. These include the National Aeronautics and Space Administration and the Department of Energy in the United States, the Commissariat à l'Energie Atomique and the Centre National de la Recherche Scientifique/Institut National de Physique Nucléaire et de Physique des Particules in France, the Agenzia Spaziale Italiana and the Istituto Nazionale di Fisica Nucleare in Italy, the Ministry of Education, Culture, Sports, Science and Technology (MEXT), High Energy Accelerator Research Organization (KEK), and Japan Aerospace Exploration Agency (JAXA) in Japan, and the K. A. Wallenberg Foundation, the Swedish Research Council, and the Swedish National Space Board in Sweden.

Additional support for science analysis during the operations phase is gratefully acknowledged from the Istituto Nazionale di Astrofisica in Italy and the Centre National d'Études Spatiales in France.

Facility: Fermi-(LAT).