EDP Sciences
Free Access
Issue
A&A
Volume 541, May 2012
Article Number A13
Number of page(s) 11
Section Astrophysical processes
DOI https://doi.org/10.1051/0004-6361/201218846
Published online 19 April 2012

© ESO, 2012

      

1. Introduction

W51 is a massive molecular complex located at the tangential point (l = 49°) of the Sagitarius arm of the Galaxy, at a distance of  ~5.5 kpc (Sato et al. 2010). As seen in radio continuum images, three main components are identified: the star-forming regions W51A and W51B and, attached to the south-eastern boundary of W51B, the supernova remnant (SNR) W51C. The estimated age of this SNR is 30 kyrs (Koo et al. 1995a). Evidence of interaction between W51C and W51B is provided by several observations. Most crucial of them are the existence of two 1720 MHz OH masers (Green et al. 1997) and the detection of about 103 solar masses of atomic gas at a velocity shifted between 20 and 120 km s-1 with respect to its ambient medium (Koo & Moon 1997a). The high-velocity atomic gas exhibits a counterpart in high density molecular gas clumps (Koo & Moon 1997b) sharing the same location and velocity shift. Koo & Moon showed that the shocked gas is displayed in a thin layer in the interface between the SNR shell, as delimited by the X-ray image from ROSAT and the unshocked molecular gas. This can be taken as the existence of a J-type shock penetrating the dense gas in a particular region of W51B (Koo & Moon 1997b), whereas in the location of the 1720 MHz OH masers the shock should be continuous (C-type). Moreover, recent measurements (Ceccarelli et al. 2011) showed over-ionization of the gas in W51B in certain locations close to W51C coinciding with the shocked gas. They conclude this excess in ionization implies the existence of an intense flow of freshly accelerated cosmic rays (CRs) that, through proton-proton collisions, ionize the hydrogen in the adjacent cloud. However,  ~0.2 degrees South-East to the shocked gas region, a hard X-ray source CXO J192318.5+140305 is detected. This object was first resolved by ASCA (Koo et al. 2002) and later confirmed by Chandra (Koo et al. 2005). Its X-ray spectrum, together with its morphology, suggests that it is a possible pulsar wind nebula (PWN) associated with the SNR. Therefore, the presence of CXO J192318.5+140305 plays a role in the interpretation of the gamma-ray emission from the W51 region. For these reasons, W51C represents an interesting case for the study of the acceleration of particles to very high energies (VHE) and their interaction with the interstellar medium.

An extended source of gamma rays was first detected by the H.E.S.S. telescopes with an integral flux above 1 TeV of about 3% that of the Crab Nebula (Fiasson et al. 2009). However, the presented morphological and spectral information was not enough to attribute the origin of the emission to any particular object in the field of view. Also, the Large Area Telescope (LAT) on board the Fermi satellite detected an extended source between 200 MeV and 50 GeV coincident with the H.E.S.S. source (Abdo et al. 2009a). Moreover, the reanalysis of the archival MILAGRO data after the release of the first Fermi catalog revealed a 3.4σ excess with median energy of 10 TeV coincident with the Fermi/LAT source (Abdo et al. 2009b). At radio wavelengths, synchrotron radiation on ambient magnetic field explains the emission detected from W51C. At higher energies, there are several processes that yield emission of gamma rays: inverse Compton scattering of electrons on seed photons (cosmic microwave background, starlight), non-thermal bremsstrahlung of electrons on charged target, and decay of neutral pions created in flight from a proton-nucleon collision. The modeling done by (Abdo et al. 2009a) of the spectral energy distribution (SED) of W51C disfavors leptonic models and suggests a hadronic origin for the emission. For the hadronic channel, two main (non-exclusive) mechanisms are to be considered: molecular cloud illumination by cosmic rays that escaped the accelerating shock (Gabici et al. 2009; Ohira et al. 2011) or emission from clouds that are being overtaken by the SNR blast wave (Uchiyama et al. 2010; Fang & Zhang 2010). It is well known that a 10% of the energy released by the supernova explosions in the Galaxy can account for the energy budget of the CR spectrum up to energies close to the knee (~1015 eV). Nevertheless, the evidence that SNRs can accelerate particles up to such high energies is still missing. Since W51C is one of the most luminous Galactic sources at Fermi/LAT energies, observation of gamma rays up to several TeV would have serious implications regarding the SNR contribution to the Galactic CRs: such an observation would show that SNRs are not only capable to provide a sufficient flux, but could also shed light on the question of the maximum energy of CR’s still achievable in such a medium age SNR.

However, the object from which the gamma rays originate has not yet been identified within the W51 field, and the gamma-ray spectrum has so far been precisely measured only up to some tens of GeV. In what follows, we report observations with the MAGIC telescopes, which will help to address some of the remaining questions on the gamma-ray source in the W51 region, both regarding its precise location and the physical processes needed to explain the observations. In Sect. 2 we describe the observations that we performed; in Sect. 3 we show the observed morphology and spectral properties; and, finally, in Sect. 4 we apply a theoretical framework that can explain the detected gamma-ray emission.

2. Observations

MAGIC consists of two 17 m diameter Imaging Atmospheric Cherenkov Telescopes (IACT) located at the Roque de los Muchachos observatory, on La Palma island, Spain (28°46′ N, 17°53′ W), at the height of 2200 m a.s.l. The stereo observations provide a sensitivity1 of 0.8% of the Crab Nebula flux at energies  >300 GeV, see Aleksić et al. (2012). MAGIC has the lowest trigger threshold of all operating IACTs, enabling it to observe gamma rays between 50 GeV and several tens of TeV. MAGIC observed W51 in 2010 and 2011. In the first period of observations between May 17 and August 19, 2010 about 31 h effective time remained after quality cuts. Between May 3 and June 13, 2011 additional 22 h effective time of good quality data were taken, resulting in a total amount of 53 h effective dark time and covering a zenith angle range from 14 to 35 degrees. The observations were carried out in the so-called wobble mode around the center of the Fermi/LAT source W51C (RA = 19.385 h,Dec = 14.19°). All data were taken in stereoscopic mode, recording only events which triggered both telescopes. To minimize systematic effects in the exposure and to optimize the coverage for an unknown extension of the emission a total of six pointing positions (npoint = 6), were used. In all pointing positions the wobble distance (offset from the central position) was 0.4°, as it is regularly done in MAGIC observations.

The analysis of the data was performed using the MARS analysis framework (Moralejo et al. 2009) including the latest standard routines for stereoscopic analysis (Lombardi et al. 2011). After calibrating the signal and cleaning the images of the two telescopes individually, the two images of each stereo event are combined. The arrival direction is determined from the combination of the individual telescope information. To suppress the background, a global variable dubbed hadronness is determined by using the so-called random forest method (Albert et al. 2008). The energy of individual events is estimated using look-up tables generated from gamma-ray Monte-Carlo events. For a detailed description of the complete analysis chain described above see Aleksić et al. (2012). The gamma-ray signal is estimated by comparing the spatial distribution of gamma-like events around the assumed source position (ON region) with respect to those recorded in signal-free (OFF) regions. The total signal of the source is evaluated using a cut on the squared angular distance between reconstructed gamma-ray direction and source position of θ2 < 0.07. For each pointing position, the ON sample is compared to an OFF sample obtained from the combination of the npoint − 1 OFF regions observed at the same focal plane coordinates but from the complementary pointing positions. Four of the pointing positions have an observation time of the order of  ~12 h each. Therefore, three background samples per pointing can be averaged. The remaining two positions have an observation time around 2 h each and, in this case, the background was estimated from one sample only. This method ensures a maximum usage of symmetrical OFF positions without introducing big scaling factors due to differences in the observation time. The significance of the excess is determined from the combined θ2 distribution of all individual pointing positions using Eq. (17) in Li & Ma (1983, hereafter Li&Ma).

3. Results

3.1. Detection

Figure 1 shows the relative flux map2 above an energy threshold of 150 GeV around the center of the observations. The angular resolution of MAGIC for this analysis is 0.085° defined as one sigma of a Gaussian distribution, see Aleksić et al. (2012) for details. The map was smeared with a two-dimensional Gaussian kernel with a sigma equivalent to that of angular resolution3. Contours represent isocurves of test statistics (TS) evaluated from the excess of gamma-like events over a background model. This test statistic is Li&Ma significance, applied on a smoothed and modeled background estimation. Its null hypothesis distribution mostly resembles a Gaussian function, but in general can have a somewhat different shape or width. The signal region is defined within 0.265° radius around the Fermi/LAT position. This radius is selected in order to include the emission observed in the relative flux map. We compute an excess of 1371.7 ± 122.5 events inside the signal region, yielding a statistical significance of 11.4 standard deviations. The centroid of the emission (black dot in Fig. 1, statistical errors are represented by the ellipse) has been derived by fitting a 2 dimensional Gaussian function to the map, prior to the smearing. As the centroid we find:

This deviates by 0.04° from the position reported by Fermi/LAT, marked as the center of the sky map (green cross) (see Fig. 1).

thumbnail Fig. 1

Relative flux (excess/background) map above 150 GeV around W51. Overlaid are contour levels from test statistics starting at 3 and increasing by one per contour. The map was smoothed with a Gaussian kernel of 0.085°. The green cross represents the center of the observations, while the green dashed circle represents the integration area. The black dot is the determined position of the centroid with the statistical uncertainties shown by the surrounding black ellipse. The region of shocked atomic and molecular gas (Koo & Moon 1997b,a) is represented by the red dashed ellipse. The blue diamond shows the position of the possible PWN CXO J192318.5+140305. In the left lower corner the Gaussian sigma of a point-like source (PSF) after the applied smearing is shown.

Open with DEXTER

To determine the extension of the source we computed the distribution of the squared angular distance θ2 between the arrival direction of the gamma-like events and the centroid of the MAGIC source (see Fig. 2), both for the integration area represented in Fig. 1 and for a combination of signal-free regions from where we estimate the background.

thumbnail Fig. 2

θ2 distribution of the excess events towards the centroid of the emission determined from Fig. 1, showing a clear and extended signal. The excess has been fitted by an exponential (blue curve) to determine the extension. For comparison the shape of a point-like source with the same excess determined from Monte-Carlo simulations is shown (red curve). The energy threshold of this analysis is 150 GeV.

Open with DEXTER

We then fit the difference between ON and OFF θ2 distributions using an exponential function (corresponding to a Gaussian-shaped source). For illustration, the shape of a point source with the same excess was calculated from Monte-Carlo simulations and is shown as comparison (red curve) to the fit to the data (blue curve). After correcting for the angular resolution (0.085 degrees  >150 GeV) of the instrument the intrinsic extension of the source is determined to be: 0.12 ± 0.02stat ± 0.02syst degrees.

thumbnail Fig. 3

Differential energy spectrum of W51 obtained by MAGIC. The red points represent the differential flux points after unfolding. The red line represents a power law fit to the data. The error bars represent the statistical errors. For comparison, the dotted line represents the spectrum of the Crab Nebula as shown in Aleksić et al. (2012).

Open with DEXTER

3.2. Spectrum

We extracted the energy spectrum of the gamma-ray emission. The effective area was estimated using a Monte-Carlo data set with photons simulated uniformly on a ring of 0.15 to 0.55° distance to the camera center. This accounts for variations of the acceptance across the area of the source. The effect of using this ring Monte-Carlo compared to standard point-like ones turns out to lie well within the statistical uncertainties. The spectrum needs to be unfolded in order to take into account the finite energy resolution and the energy bias of the instrument (Albert et al. 2007). The spectrum shown in Fig. 3 starts at 75 GeV and is well described (χ2/NDF = 5.26/6) by a simple power law of the form: (1)with a photon index of Γ = 2.58 ± 0.07stat ± 0.22syst, and a normalization constant at 1 TeV of N0 = (9.7 ± 1.0stat) × 10-13   cm-2   s-1   TeV-1. This is the first time that the differential energy spectrum at VHE is published. The energy threshold of MAGIC allows us to almost connect the spectrum to the Fermi/LAT points (Abdo et al. 2009a). The systematic error on the flux normalization is 15%, which includes the systematic uncertainties of the effective area (11%) and the background calculation. In addition, the systematic uncertainty in the energy scale is estimated to be 17% at low (~100 GeV) and 15% at medium (~250 GeV) energies. The integrated flux above 1 TeV is equivalent to  ~3% of the flux of the Crab Nebula above the same energy, and therefore agrees with the previous flux measurement by the H.E.S.S collaboration (Fiasson et al. 2009). The spectral index measured by MAGIC agrees well with the one measured by Fermi/LAT above 10 GeV (Paneque et al. 2011) of Γ = 2.50 ± 0.18stat. The emission from W51 can be described by a single power law between 10 GeV and 5.5 TeV.

3.3. Detailed morphology

MAGIC reaches its best sensitivity in the energy range from  ~300 to  ~1000 GeV. At energies of 300 GeV the angular resolution of MAGIC is 0.075° and it improves until reaching the saturation value of 0.054° at energies above 1 TeV. We investigate sky maps in two energy ranges. The first map covers the estimated energy range from 300 to 1000 GeV, and the second the energies above 1000 GeV. Both maps were smeared with a Gaussian kernel of a width equal to the angular resolution of the instrument in each energy range.

thumbnail Fig. 4

Relative flux maps: from 300 GeV to 1000 GeV (top) and  >1000 GeV (bottom). On the left hand side the MAGIC data are combined with the 13CO (J = 1–0) intensity maps from the Galactic Ring Survey (see http://www.bu.edu/galacticring/new_index.html) integrated between 63 and 72 km s-1 shown as green countours. On the right hand side the green contours represent the 21 cm radio continuum emission is shown from (Koo & Moon 1997a). In all maps the blue diamond represents the position of CXO J192318.5+140305 and the black cross the position of the OH maser emission (Koo et al. 2005; Green et al. 1997). The red dashed ellipse represents the region of shocked atomic and molecular gas (Koo & Moon 1997b,a). The 3 counts contour above 1 GeV determined by Fermi/LAT is displayed by the pink contour. In each picture the Gaussian sigma of a point-like source (PSF) after the applied smearing is shown. The color scale (blue to red) represents the relative flux as measured with MAGIC. In addition the TS contours (cyan) are shown starting at 3 and increasing by one per contour.

Open with DEXTER

In Fig. 4 (top panels) the relative flux map between 300 and 1000 GeV is shown. The overall shape of the emission appears to be elongated showing a tail towards the lower left. The maximum of the emission coincides with the shocked-gas region, represented by the red dashed circle, where the lack of molecular material at the systemic velocity is clear (top left panel). The determined centroid and extension agree within statistical errors with those found above 150 GeV.

Above 1000 GeV (Fig. 4, bottom panels) the centroid and the extension of the emission are in agreement with those obtained at lower energies. The South-Eastern tail of the source, evident in the 300 to 1000 GeV map, becomes a prominent feature coincident with the possible PWN CXO J192318.5+140305 at energies above 1 TeV. However, the main part of the emission is still coincident with the shocked gas region.

While the centroid of the emission is consistent with the position of the shocked gas, we see a tail towards the PWN candidate. We note that, in any case, the VHE emission does not strictly follow the SNR shell (as seen from the 21 cm continuum emission represented by green contours in the right panels), nor does it follow the molecular gas with the velocity expected due to Galactic rotation, as traced by the 13CO (green contours, left panels). The tail seen towards the PWN rises the question of a possible substructure in the emission.

3.3.1. Projections

In order to investigate the source for underlying structures, we project the unsmeared excess distribution of the source along a line. The line is 2° long divided in 40 bins with 0.05° width. The orientation of the line is defined by the position of the PWN candidate and the centroid of the shocked clouds identified by Koo & Moon (RA = 19.380 h,Dec = 14.19°). Events within a distance of 2 Gaussian sigma of the instrumental PSF to the line were projected. Since the angular resolution is energy dependent, the width of the projected rectangle is 0.3° and 0.216° for the energy ranges from 300 to 1000 GeV and above 1000 GeV, respectively. OFF events were estimated from the background model. The number of projected excess events is not the same as in the spectral calculation, where we used a circular region of 0.265° radius around the center of the observations. Therefore, the projected excess does not allow for direct determination of the fluxes from specific regions of the map. The projection has been carried out in both energy ranges independently on the unsmeared excess distribution and is shown in Fig. 5.

thumbnail Fig. 5

Projection of the excess inside the marked box in both differential sky maps: 300 GeV to 1000 GeV (top) and above 1000 GeV (bottom) along the line connecting the PWN and the shocked-gas region described in Koo & Moon (1997a,b). The projection is done with the unsmeared distribution. The excess is fitted with one (black) and two (red) Gaussian curves. The positions of the shocked gas and the PWN are marked with red arrows. On the right-hand side a sketch of the skymaps in both energy ranges is shown to illustrate the projected areas, as well as the position of the cloud and the PWN, respectively. The box has a length of 1° and a width of 4 Gaussian sigma of the instrumental PSF. The sky maps show the smeared excess (for comparison with Fig. 4) with the black contour representing the 3 TS contour.

Open with DEXTER

We fit the projection alternatively using one and two Gaussian functions. χ2/d.o.f. values are 28/17 (one Gaussian) and 18/14 (two Gaussians) for the medium-energy range and 16/17 (one Gaussian) versus 12/14 (two Gaussians) for the high-energy events. The data are very well described with the two Gaussian functions, where the centroid of the individual functions coincides within statistical errors with the position of the shocked gas and the PWN. The tail-like feature towards the possible PWN is more peaked in the energy range above 1000 GeV.

The statistics are not sufficient to clearly discriminate between an extended source of Gaussian excess, an extended source of a more complicated shape, or two individual sources. However, the fact that there is no region of dense gas close to the PWN makes it difficult to explain the enhancement of TeV emission in this area under the assumption of uniform CR density. A possible scenario of two emission regions could manifest in different spectral behaviours.

3.3.2. Energy spectra of individual regions

To quantify the results obtained from the projections we investigated in more detail the spectral properties of the detected signal, we concentrated on two individual regions within the source and analyzed them separately. One was defined to cover the shocked cloud region with centroid at RA = 19.380 h,Dec = 14.19°; this will be called the cloud region. The second one was defined by the position of CXO J192318.5+140305 and will be called the PWN region. To avoid contamination from the surrounding emission, and their possible spread due to the worse angular resolution at lower energies, we use an integration radius of 0.1°. We compared the same analysis on data of the Crab Nebula and find that such a region contains at least 70% of the excess from a point-like source above 300 GeV. For an easier comparison, the integration radii were chosen to be the same. The distance between the chosen positions is 0.19°. There is an area of overlap of 1.7% compared with the integration area of each region, therefore they can be treated as independent. The combined areas of both regions represent about 29% of the area used to determine the overall spectrum.

The small distance between the regions and a very similar average distance to the camera center allow us to assume the same acceptance of gamma-like events for both regions, at least within 5%.

For each individual region we determined the amount of excess events above three different energies, and calculated the contribution to the overall emission. The resulting values are shown in Table 1. Excesses used to calculate these ratios show a significance of at least 2.9σ.

Table 1

Number of excess events determined for the PWN-region and the cloud-region and their contribution to the overall emission.

The excess contribution arising from the cloud region is about 30% and shows no dependence on energy. We performed a spectral analysis of a point source for the cloud region above 350 GeV. The emission can be well described by a pure power law with a flux normalization constant at 1 TeV of Ncloud = (4.3 ± 0.9stat) × 10-13   cm-2 s-1 TeV-1. The integrated flux above 350 GeV is equivalent to 1.2% of the flux of the Crab Nebula. The spectral index of the cloud emission is −2.4 ± 0.5stat and agrees within statistical uncertainties with the spectral index of the overall emission.

Assuming a point-like emission, the flux from the PWN region above 350 GeV is equivalent to 0.7% of the flux of the Crab Nebula representing about 20% of the overall observed emission. The emission between 350 GeV and 2 TeV can be well described by one single power law with a spectral index of −2.5 ± 0.6stat and a flux normalization at 1 TeV of NPWN = (2.3 ± 0.8stat) × 10-13   cm-2 s-1 TeV-1.

The excess contribution of each of the regions shows no dependence on energy, suggesting no intrinsic morphological changes in the energy ranges investigated here. This is in agreement with the spectra, with the differential maps, and with the projections of the excess distribution. We want note that the number of excess events within the PWN region and the cloud region (Table 1) agrees within statistical errors with the projected excess (Fig. 5) found within  ± 0.1 degree from the PWN and the cloud positions, respectively. By looking at the skymaps (Fig. 5) only, the emission around the PWN seems to be more intense above 1 TeV. However this can be explained by the worse angular resolution at lower energies and by a much higher signal-to-noise ratio at the higher energies.

4. Discussion

Before modeling the multi-wavelength emission we address shortly the possibility that the PWN alone is the source of all the emission. Then we qualitatively discuss the possibility of a contribution of the PWN to the overall emission and justify the approach using a one-zone model (i.e. one homogeneous emission region filled with one particle distribution per species) to investigate the processes underlying the emission.

First, we want to recall that we have found no spectral or morphological energy dependence. In order to assess whether the VHE emission can be originated only by the PWN candidate, we consider the estimated rate of rotational energy loss [Ėentity!#xa0!] = 1.5 × 1036 erg s-1 estimated by Koo et al. (2005) with the empirical relation from Seward & Wang (1988). Given the observed luminosity of the order of  ~1036 erg s-1 reported in Abdo et al. (2009a), it seems unlikely that the PWN alone is the source of all gamma-ray emission since, it would require an extremely high efficiency in the conversion from rotational energy into gamma rays.

Second, we consider the possibility of a two-zone model. The PWN region can account for the 20% of the gamma-ray emission confined in a point-like source, however the brightest part is the cloud region. This scenario would require an cĖonversion into gamma rays of the order of 10%, in agreement with the generally accepted value.

With the current statistics and resolution, it cannot be established if there is a spectral difference between the cloud and the PWN region, but the contribution of the PWN is in any case small. For the reasons above, the simplest approach is to assume one overall particle distribution underlying the emission we observe. This assumption introduces an error in the flux normalization of about 20% in case part of the emission originates from the PWN candidate; this uncertainty lies within the statistical and systematic errors of the MAGIC measurement.

4.1. Model description

We model the SNR as a sphere homogeneously filled with hydrogen, helium and electrons, with respective average number densities nH, nHe and ne. For the relative abundances of helium we assume the cosmic abundance ratio nHe = 0.1nH. For the electron ratio we assume full ionization of the medium, such that ne = 1.2nH. The magnetic field B is assumed to be homogeneous inside the sphere; Koo et al. (2010) derived an upper limit for it of B < 150 μG, but Brogan et al. (2000) measured a local magnetic field as high as 1.5–1.9 mG towards the maser sites.

The geometric model of the SNR and molecular cloud interaction region as proposed by Koo & Moon (1997a), describes a scenario in which the spherical blast wave of the supernova explosion interacts with part of the cylindrical molecular cloud contained inside the SNR volume. Carpenter & Sanders (1998) estimated the total mass of the molecular cloud to be mcloud = 1.9 × 105M. From the radio measurements in Moon & Koo (1994) the angular extent of the partial radio shell of the SNR is known to be θ ≈ 30′. We see a clear displacement between the morphology presented here and the center of the spherical extended SNR as seen in thermal X-ray emission (Koo et al. 1995b). The maximum of the emission is located at the interaction region of remnant and the molecular cloud. Therefore we conclude that the size of the remnant is not physically related to the size of the VHE emission region. We adopt the intrisnic extension determined in this work to determine the radius of a spherical emission zone. Assuming a distance to W51C of 5.5 kpc, as measured by Sato et al. (2010) and Moisés et al. (2011), the radius of the sphere is estimated to be 24pc.

The explosion energy of the SNR has been estimated in Koo et al. (1995b) as ESN ≈ 3.6 × 1051erg, using both a Sedov and an evaporative model to derive the parameters of the SNR. We will compare this value with the one obtained from the integral of our initial spectra, after we fix the normalisation constants; we will determine how much of the initial explosion energy of the supernova has been converted into particles (We, Wp). The different parameters of the supernova, of the SNR, and of the molecular cloud are summarized in Table 2.

Table 2

Parameters of the W51C supernova, supernova remnant and molecular cloud.

We model the spectral energy distribution folding input spectra of accelerated particles with cross sections of processes yielding photons; this includes synchrotron radiation, inverse Compton scattering (IC), non-thermal bremsstrahlung and π0 decay (Blumenthal & Gould 1970; Baring et al. 1999; Kelner et al. 2006).

For IC, we consider three seed photon fields: the cosmic microwave background (kTCMB = 2.3 × 10-4 eV, uCMB = 0.26eVcm-3), infrared (kTIR = 3 × 10-3eV,uIR = 0.90eVcm-3) and optical (kTOPT = 0.25eV,uOPT = 0.84eVcm-3), with temperatures and energy densities for the infrared and optical components adopted from Abdo et al. (2009a). Bremsstrahlung is computed on a target of electrons and ions. For the π0 production cross section, we use the parametrization of Kelner & Aharonian (2006) with a constant nuclear enhancement factor of 1.85 (Mori 2009).

The multi-wavelength data considered here include radio continuum measurements (Moon & Koo 1994), high-energy observations by the Fermi/LAT (Abdo et al. 2009a) and the new VHE data taken with MAGIC, presented in this paper. Included is also one data point by MILAGRO (Abdo et al. 2009b). Note that the lowest energy radio data point may be affected by free-free absorption, see Moon & Koo (1994) or Copetti & Schmidt (1991), which we do not consider here. However, this single point does not affect the fitting of the radio data. The radio measurements in Moon & Koo (1994) indicate a spectral index of αr ≈ −0.26 (as defined by Sν ∝ ναr). This can be attributed to electrons emitting synchrotron radiation and fixes the initial power-law index of the electron spectrum to s ≈ 1.5. We adopt this value both for electrons and protons.

As an upper limit for the non-thermal X-ray emission we consider the integrated thermal X-ray flux of the whole remnant as measured by ROSAT (Koo et al. 1995b) converted into a differential flux in the sub-keV range. We use the thermal emission observed by Chandra from CXOJ192318.5+140305 as an upper limit to the non-thermal emission of the possible PWN. The MILAGRO measurement has a significance of 3.4σ, was derived assuming a gamma-spectrum  ∝ E-2.6 without a cut-off and is given at an energy of 35 TeV. For details see Abdo et al. (2009b).

We consider seperate scenarios in which one of the following emission processe dominates over the others, pion decay, inverse Compton, or Bremsstrahlung. The models discussed here are obtained using as equilibrium particle spectra a broken power law with an exponential cut-off, both for electrons and protons, of the form: (2)The spectral index changes here from s to s + Δs at an energy Ebr with a smooth transition. The exponential cut-off at Ecut,e,p reflects the roll-off of the particle spectrum near the maximum energy, arising from the acceleration and confinement mechanism, as well as energy losses.

The break energy Ebr is fixed from the Fermi/LAT data, while the new MAGIC data allow us to fix the spectral break Δs. A spectral break in the particle spectrum at these energies is traditionally thought to be inconsistent with both standard or non-linear diffusive shock acceleration theory, see Malkov & O’C Drury (2001) and references therein. However, Malkov et al. (2011) have recently proposed a mechanism which can also explain a spectral break in the cosmic ray spectrum of Δs = 1 by strong ion-neutral collisions in the surroundings of a SNR, leading to a weakening in the confinement of the accelerated particles. The spectral break that we have derived here is Δs = 1.2, not far off this prediction, giving a hint that this mechanism might be responsible for the observed break. Note also that other authors have proposed scenarios in which the CR spectrum, and consequently the gamma-ray spectrum, can show one or more spectral breaks, for example due to finite-size acceleration or emission region (Ohira et al. 2011) or energy dependent diffusion of run-away CRs from the remnant (Gabici et al. 2009; Aharonian & Atoyan 1996).

The luminosity of W51C in the energy range 0.25 GeV–5.0 TeV, which is roughly the energy range of the Fermi and MAGIC data, is Lγ ≈ 1 × 1036 erg s-1, assuming a distance of 5.5 kpc, which is one of the highest compared with other SNRs.

4.2. Adjustment of model parameters

First we consider the case where the emission is dominated by leptonic emission mechanisms. We find the same problems already reported by Abdo et al. (2009a), namely that it cannot reproduce the radio and gamma-ray data simultaneously. Furthermore these models need an unusually high electron to proton ratio of the order of one.

When we model the emission with pion decay as the dominant process, both radio and gamma-ray emission can be reasonably reproduced, as shown in Fig. 6. A hadronic scenario is particularly interesting, as the shock-cloud interaction naturally favors a CR-matter interaction mechanism. Moreover, the parameters used in this model, see Table 3 are a reasonable description of the interstellar medium around W51.

thumbnail Fig. 6

Model of the multi-wavelength SED in the hadronic-dominated scenario. The dashes with error bars are 21 cm radio continuum, circles represent Fermi/LAT data, squares are the data obtained in this work and the star represent the MILAGRO data point. The upper limit in the X-ray regime is obtained from ROSAT data as discussed in the text. The details of the scenario are discussed in the text.

Open with DEXTER

Table 3

Parameters used in the modeling of the multi-wavelength spectral energy distribution for the hadronic scenario.

Compared to the hadronic model suggested in the work of Abdo et al. (2009a), the main difference is the index of the partcile distribution after the break. The spectrum after the break is more precisely determined by the data presented here. The index we obtain is harder, allowing for the explanation of all the gamma-ray data up to the end of the MAGIC spectrum.

A detailed view of the high energy and VHE region is shown in Fig. 7. It shows that the index above the break is clearly determined by the data presented here. In addition, the hadronic model by Abdo et al. (2009a) is displayed. In addition to the good aggrement between the model and the data, the plot shows that the results presented here clearly improve the determination of the underlying particle distribution. In this scenario a cut-off energy of Ecut,p ≥ 100 TeV is needed to fit the MAGIC data, indicating the existence of protons at least to this energy.

thumbnail Fig. 7

Detailed view of the hadronic model in the high energy and VHE region. For comparison to the hadronic model by Abdo et al. (2009a), shown as double dotted line. This models ends at the highest energy shown in that publication. The major difference between our model and that of Abdo et al. (2009a) is the harder particle spectrum above  ~100 GeV, which is now precisely constrained by the measurements presented here.

Open with DEXTER

The precise cut-off energy of the electron spectrum Ecut, e is not well constrained, since the synchrotron peak is not resolved. Therefore, the energy Ecut, e used here is only a lower limit, as enforced with the radio data. However, a 1 TeV electron in a magnetic field of 50 μG has a lifetime of about 4700 years, determined by synchrotron losses. This value is much lower than the age of the remnant, suggesting that for such high energies the electron spectrum should develop a break, with the consequent spectral steepening at higher energies. Assuming constant electron injection over time, the electron spectrum steepens at 100 GeV by a factor 1/E. This yields a very similar gamma-ray emission as in the hadronic model presented here, even for a higher value of Ecut, e.

4.3. Physical outcome of the models

We discuss what general conclusions can be drawn from the model which fit the data: the hadronic scenario.

The volume-averaged hydrogen density is obtained as a parameter of the fit. From that, we compute the volume filling factor f, which is the fraction of the mass of the clumpy molecular cloud that is contained inside the SNR interaction volume (defined as volume of the emission zone) as . Here mcloud is the total mass of the molecular cloud, V is the volume of the radiation sphere and , are the masses of a hydrogen and helium atom, respectively. This would imply that around 11% of the mass of the molecular cloud is contained in the emission volume and is interacting with the SNR. This value is consistent with the filling factors of around 8–20% for other SNRs interacting with molecular clouds, obtained by other authors (Uchiyama et al. 2010).

The total amount of kinetic energy in electrons and protons is about 16% of the explosion energy of the supernova. This fraction is just slightly higher than the value normally assumed, of around 10%, of the explosion energy converted into CRs to maintain the observed flux of Galactic CRs (Hillas 2005). The proton to electron ratio is not far from value observed at earth of Kp/Ke ≈ 50, see for example Simpson (1983).

Since the hadronic gamma-ray emission is proportional to the product of the kinetic energy in protons and the density of the medium, this parameters are striclty correlated. Assuming that the complete mass of the molecular cloud acts as target material (f = 1), this would imply a density of n = 100 cm-3. Therefore the lower limit of the energy in relativistic protons is about 1.6% of the explosion energy of the supernova. We note that such a scenario would need either a higher magnetic field (B ~150   μG) or a much lower electron to proton ratio (Ke/Kp ~ 1/800) to still reproduce the broadband emission. In addition, the morphology presented in this work shows that only a fraction of the molecular cloud is emitting VHE gamma emission (see Fig. 4). Therefore we conclude that the amount of kinetic energy in protons is clearly above this lower limit and in the order of 10–20%.

In the scenario investigated here all of the gamma-ray emission was attributed to π0 decay. It was not possible to model the broad-band emission with a purely leptonic scenario. The radio data could not be fitted and the model parameters were not physically reasonable (too low density nH, too high energy content We in electrons, too low magnetic field B). However, that could also point to problems in the modeling, especially to oversimplifications concerning the homogeneity of the medium and of the magnetic field.

We conclude that the Fermi/LAT data and the MAGIC data can be explained in terms of hadronic interactions of high-energy protons with the molecular cloud and subsequent decay of neutral pions. With the current data it is not possible to decide what process causes the hint of emission observed by MILAGRO which, if confirmed at this flux level, would require the introduction of an additional component at the highest energies.

4.4. Discussion on the acceleration process

Following the result of the modeling we assume the observed gamma-ray emission to be of hadronic origin. As mentioned in Sect. 1, there are two main possible scenarios: a cloud illuminated by runaway CRs or acceleration of CRs in the shock wave propagating through the cloud.

In the first case, CRs escaping the SNR will homogeneously fill a sphere with a radius where D is the diffusion coefficient and t is the time since particles are diffusing (Gabici et al. 2010). For a distance of 5.5 kpc and 10 TeV protons, responsible for gamma-ray emission of 1 TeV, the radius of this sphere would be about 350 pc, assuming the average Galactic CR diffusion coefficient at 10 TeV to be  ~3 × 1029   cm2 s. Here we assumed that the high-energy particles escape the SNR early enough such that the diffusion time can be approximated to be the age of the SNR. The distance between the maximum of the emission measured by MAGIC above 1 TeV and the assumed center of the SNR (RA = 19.384 h, Dec = 14.11°) is about 8 pc. The distance to other parts of the SNR/cloud complex W51C/B is of similar order. This implies that the complete cloud should be uniformly illuminated by CRs. As can be seen in Fig. 4 we do not detect the complete W51B/C complex at energies above 1 TeV (lower right skymap): parts of the outer regions, both on the side towards the SNR and on the opposite side, do not emit gamma radiation. In the scenario of runaway CRs, we would expect diffusion from the SNR to W51A (northern region in the 21 cm emission); no significant emission from W51A is detected. However, the distance of the regions A and B is measured with an error of the order of hundreds of parsecs, which means that the relative distance between the two could be high enough to explain the lack of diffusion from one to the other. The scenario of runaway CRs also can not explain the incomplete illumination of W51B/C, especially towards the outer regions.

Concerning the acceleration of CRs in the shocked cloud scenario, the gamma radiation should be originated very close to the acceleration site of the radiating particles due to the high density of the surrounding medium. This is in agreement with the morphology described in this work. The unusually high ionization reported by Ceccarelli et al. (2011) close to the maximum VHE emission region indicates the presence of freshly accelerated low-energy protons. The missing emission towards the edges of the cloud could be explained with a lower diffusion coefficient in the shocked cloud region, or with a shielding effect, either of which is possible in a surrounding medium of high density.

Both the morphology at TeV energies and the measured high ionization are hints for an ongoing acceleration. This suggests that the particle distribution, whose gamma emission we observe, may represent the source spectrum of cosmic rays currently being produced in W51. However, the differentiation between ongoing acceleration of particles in the shocked region or reacceleration of already existing CRs, like in the crushed cloud scenario (Uchiyama et al. 2010), in the same region is not obvious and is not addressed in this work.

5. Conclusions

MAGIC has performed a deep observation of a complex Galactic field containing the star-forming regions W51A and W51B, the SNR W51C and the possible PWN CXO J192318.5+140305. As a result of this observation, emission of gamma rays above 150 GeV has been detected with 11σ statistical significance. The spectrum of this emission has been measured between 75 GeV and 4 TeV. Spectral points are well fitted with a power law with a photon index of 2.6, compatible with the Fermi/LAT measurement between 2 and 40 GeV. The spectrum measured by MAGIC allows for the first time a precise determination of the spectral slope of the underlying particle distribution above the spectral break measured at around a few GeV by Fermi/LAT.

The MAGIC source spatially coincides with those previously reported by H.E.S.S. and Fermi/LAT. We are able to restrict the emission region to the zone where W51C interacts with W51B and, in particular, to the region where shocked gas is observed. This clearly pinpoints the origin of the emission to the interaction between the remnant and the molecular cloud.

Non-thermal X-ray emission which could help to trace the relativistic electron distribution was found only from a compact region around the position of the possible PWN CXO J192318.5+140305 (Koo et al. 2005). The MAGIC source exhibits a morphological feature extending towards CXO J192318.5+140305, more prominent in the image at higher energies.

The projection of the gamma-like events on the line connecting the putative PWN and the centroid of the shocked clouds shows a hint of an underlying distribution that may be described as the sum of two Gaussian functions. However, the existence of two independent, resolved sources cannot be statistically established. We thus investigate the contribution to the total excess of two regions of 0.1° radius centered on the cloud region and the PWN region. We find that they contribute about 30% and 20% of the total emission, respectively, and the contribution is not energy dependent within the uncertainties. Spectra of the individual regions above 350 GeV could be obtained, but do not allow for detailed conclusions due to the weak individual fluxes. Given the small possible contribution of the PWN candidate in the energies investigated in this work, it is very unlikely that the main conclusion drawn here will be significantly affected even if the PWN contribution can be established.

MAGIC observations determine the VHE spectral energy distribution of W51 over more than one order of magnitude in energy. We have produced a physically plausible model of the emission of the SNR by considering a spherical geometry and uniform distribution of the ambient material. We note that this system is clearly anisotropic (as seen in the multi-wavelength data), and more detailed modeling may achieve a better description of the source. We find that the VHE emission from W51C cannot be explained by any of the considered leptonic models. The emission is best described when neutral pion decay is the dominant gamma-ray production mechanism. In the proposed model, the SNR has converted about 16% of the explosion energy into kinetic energy for proton acceleration and the emission zone engulfs a 10% of a molecular cloud of 105 solar masses, which provides the target material. In this scenario, protons are required to reach at least an energy of the order of 100 TeV to produce the observed emission.

The morphology of the source cannot be explained by CRs diffusing from the SNR to the cloud. It can instead be qualitatively explained with VHE gamma-ray emission being produced at the acceleration site of CRs. This involves ongoing acceleration of CRs or re-acceleration of already existing CRs at the shocked cloud region. Given the high luminosity of this source and its plausible hadronic origin, we conclude that W51C is a prime candidate cosmic ray source in the Galaxy.

Finally, we want to give a short outlook and address a few issues connected with W51C. The detection of neutrinos from this source would be the final proof about the hadronic nature of the emission. But, according to the calculations by Yuan et al. (2011), the chances for detection are low. However, also an extension of the high-energy gamma emission towards lower energies, as performed for example in Giuliani et al. (2011), may also provide more clues to the nature of particle acceleration in this region. To reveal the morphology and the possible emission of the PWN, more data at energies above 1 TeV are necessary. Extension of the spectrum towards higher energies would constrain the maximum achievable energy in the system and might shed light on the meaning of the MILAGRO measurement, which cannot be accommodated in the theoretical framework proposed here.


1

Sensitivity is defined here as the minimal integral flux to reach 5σ excess in 50 h of observations, assuming a spectral index like that of the Crab Nebula.

2

Relative flux means excess events over background events. This quantity accounts for acceptance differences between different parts of the camera.

3

The PSF shown in all skymaps is the sum in quadrature of the instrumental angular resolution and the applied smearing.

Acknowledgments

We would like to thank the anonymous referee as well as the Associate Editor M. Walmsley for fruitful comments and suggestions. We would like to thank the Instituto de Astrofísica de Canarias for the excellent working conditions at the Observatorio del Roque de los Muchachos in La Palma. The support of the German BMBF and MPG, the Italian INFN, the Swiss National Fund SNF, and the Spanish MICINN is gratefully acknowledged. This work was also supported by the Marie Curie program, by the CPAN CSD2007-00042 and MultiDark CSD2009-00064 projects of the Spanish Consolider-Ingenio 2010 programme, by grant DO02-353 of the Bulgarian NSF, by grant 127740 of the Academy of Finland, by the YIP of the Helmholtz Gemeinschaft, by the DFG Cluster of Excellence “Origin and Structure of the Universe”, by the DFG Collaborative Research Centers SFB823/C4 and SFB876/C3, and by the Polish MNiSzW grant 745/N-HESS-MAGIC/2010/0. This publication makes use of molecular line data from the Boston University-FCRAO Galactic Ring Survey (GRS). The GRS is a joint project of Boston University and Five College Radio Astronomy Observatory, funded by the National Science Foundation under grants AST-9800334, AST-0098562, AST-0100793, AST-0228993, & AST-0507657.

References

All Tables

Table 1

Number of excess events determined for the PWN-region and the cloud-region and their contribution to the overall emission.

Table 2

Parameters of the W51C supernova, supernova remnant and molecular cloud.

Table 3

Parameters used in the modeling of the multi-wavelength spectral energy distribution for the hadronic scenario.

All Figures

thumbnail Fig. 1

Relative flux (excess/background) map above 150 GeV around W51. Overlaid are contour levels from test statistics starting at 3 and increasing by one per contour. The map was smoothed with a Gaussian kernel of 0.085°. The green cross represents the center of the observations, while the green dashed circle represents the integration area. The black dot is the determined position of the centroid with the statistical uncertainties shown by the surrounding black ellipse. The region of shocked atomic and molecular gas (Koo & Moon 1997b,a) is represented by the red dashed ellipse. The blue diamond shows the position of the possible PWN CXO J192318.5+140305. In the left lower corner the Gaussian sigma of a point-like source (PSF) after the applied smearing is shown.

Open with DEXTER
In the text
thumbnail Fig. 2

θ2 distribution of the excess events towards the centroid of the emission determined from Fig. 1, showing a clear and extended signal. The excess has been fitted by an exponential (blue curve) to determine the extension. For comparison the shape of a point-like source with the same excess determined from Monte-Carlo simulations is shown (red curve). The energy threshold of this analysis is 150 GeV.

Open with DEXTER
In the text
thumbnail Fig. 3

Differential energy spectrum of W51 obtained by MAGIC. The red points represent the differential flux points after unfolding. The red line represents a power law fit to the data. The error bars represent the statistical errors. For comparison, the dotted line represents the spectrum of the Crab Nebula as shown in Aleksić et al. (2012).

Open with DEXTER
In the text
thumbnail Fig. 4

Relative flux maps: from 300 GeV to 1000 GeV (top) and  >1000 GeV (bottom). On the left hand side the MAGIC data are combined with the 13CO (J = 1–0) intensity maps from the Galactic Ring Survey (see http://www.bu.edu/galacticring/new_index.html) integrated between 63 and 72 km s-1 shown as green countours. On the right hand side the green contours represent the 21 cm radio continuum emission is shown from (Koo & Moon 1997a). In all maps the blue diamond represents the position of CXO J192318.5+140305 and the black cross the position of the OH maser emission (Koo et al. 2005; Green et al. 1997). The red dashed ellipse represents the region of shocked atomic and molecular gas (Koo & Moon 1997b,a). The 3 counts contour above 1 GeV determined by Fermi/LAT is displayed by the pink contour. In each picture the Gaussian sigma of a point-like source (PSF) after the applied smearing is shown. The color scale (blue to red) represents the relative flux as measured with MAGIC. In addition the TS contours (cyan) are shown starting at 3 and increasing by one per contour.

Open with DEXTER
In the text
thumbnail Fig. 5

Projection of the excess inside the marked box in both differential sky maps: 300 GeV to 1000 GeV (top) and above 1000 GeV (bottom) along the line connecting the PWN and the shocked-gas region described in Koo & Moon (1997a,b). The projection is done with the unsmeared distribution. The excess is fitted with one (black) and two (red) Gaussian curves. The positions of the shocked gas and the PWN are marked with red arrows. On the right-hand side a sketch of the skymaps in both energy ranges is shown to illustrate the projected areas, as well as the position of the cloud and the PWN, respectively. The box has a length of 1° and a width of 4 Gaussian sigma of the instrumental PSF. The sky maps show the smeared excess (for comparison with Fig. 4) with the black contour representing the 3 TS contour.

Open with DEXTER
In the text
thumbnail Fig. 6

Model of the multi-wavelength SED in the hadronic-dominated scenario. The dashes with error bars are 21 cm radio continuum, circles represent Fermi/LAT data, squares are the data obtained in this work and the star represent the MILAGRO data point. The upper limit in the X-ray regime is obtained from ROSAT data as discussed in the text. The details of the scenario are discussed in the text.

Open with DEXTER
In the text
thumbnail Fig. 7

Detailed view of the hadronic model in the high energy and VHE region. For comparison to the hadronic model by Abdo et al. (2009a), shown as double dotted line. This models ends at the highest energy shown in that publication. The major difference between our model and that of Abdo et al. (2009a) is the harder particle spectrum above  ~100 GeV, which is now precisely constrained by the measurements presented here.

Open with DEXTER
In the text

Current usage metrics show cumulative count of Article Views (full-text article views including HTML views, PDF and ePub downloads, according to the available data) and Abstracts Views on Vision4Press platform.

Data correspond to usage on the plateform after 2015. The current usage metrics is available 48-96 hours after online publication and is updated daily on week days.

Initial download of the metrics may take a while.