© ESO, 2014

1. Introduction

One of the longstanding goals of very-high-energy (VHE) gamma-ray astronomy is to trace the particle populations responsible for producing TeV photons and, in so doing, to search for the sources of cosmic rays in our Galaxy. VHE gamma rays may be produced through hadronic interactions via the π⁰ decay channel and by energetic electrons through inverse compton (IC) scattering on soft photon fields or via nonthermal Bremsstralung. These parent particles are thought to be energized by shock acceleration, which very likely happens in astrophysical objects such as supernova remnants (SNRs), pulsar wind nebulae (PWNe) and compact binaries. More recently, some Galactic TeV gamma-ray sources have been linked to regions of massive star formation, although these associations are complicated by the expected presence of additional objects of the aforementioned types. Finally, nearly one third of all Galactic sources, many of which were discovered through the H.E.S.S. galactic plane survey (Carrigan et al. 2013) remain unidentified and could provide new insight into these extreme astrophysical environments.

HESS J1857+026 was discovered as a source without clear associations at other wavelengths during the H.E.S.S. survey of the inner Galaxy (Aharonian et al. 2008). The source had an intrinsic extension of (0.11 ± 0.08_stat)° × (0.08 ± 0.03_stat)°, with an inclination (measured counter-clockwise with respect to the RA-axis) of (−3 ± 49_stat)°. A significant tail-like structure toward the north was also seen, which hinted at either a larger overall extension or a weaker but distinct northern source. The spectrum of HESS J1857+026, which covered the range of 0.8−45 TeV, was well described by a single power law with a spectral index of 2.39 ± 0.08_stat.

Subsequent searches at radio wavelengths have led to the discovery of a nearby energetic pulsar, PSR J1856+0245, with a spin period of 81 ms and characteristic spin-down age of 21 kyr (Hessels et al. 2009). Given a pulsar dispersion measure (DM) of 622 ± 2 cm^-3 pc, the distance was estimated to be ~9 kpc based on the NE2001 model for the Galactic distribution of free electrons (Cordes & Lazio 2002). The uncertainty on this distance is not well determined and the value may vary by factors of 2−3 (Hessels et al. 2009), which in turn leads to large uncertainties on the estimated size and energetics of the candidate PWN.

HESS J1857+026 was also detected at GeV energies (Neronov & Semikoz 2010; Paneque et al. 2011; Rousseau et al. 2012), although no significant pulsation from the pulsar was found in 3 years of Fermi-LAT data. Rousseau et al. (2012) extracted a spectrum assuming a point source, and fit both leptonic and hadronic models to a broadband spectral energy distribution (SED) that used the original H.E.S.S. results and preliminary results from MAGIC (Klepser et al. 2011). Their multiwavelength SED also included an upper limit on the extended X-ray synchrotron flux based on archival Chandra X-ray data, the analysis of which was presented by Nice et al. (2013). The authors found no evidence of an X-ray PWN in the region around PSR J1856+0245 and derived a 3σ upper limit on the unabsorbed 1−10 keV flux that translated into a limit on the luminosity of 5 × 10³² erg s^-1 for an assumed distance of 9 kpc. The extraction region used to derive this limit was an annulus centered on the pulsar with inner and outer radii of 2″ and 15″, respectively. This extended X-ray flux upper limit implies a weak magnetic field in the case of leptonic scenarios, with best-fit values in the models of Rousseau et al. (2012) being in the order of a few μG.

We present our follow-up study of the complex morphology of HESS J1857+026. Section 2 describes the observations and data analysis performed using MAGIC. We summarize the MAGIC results in Sect. 3. In Sect. 4 we describe our multiwavelength study of this region, which makes use of archival data in the radio and infrared bands. In Sect. 5 we present our physical interpretation of the VHE gamma-ray and multiwavelength results. We summarize our conclusions in Sect. 6.

2. Analysis with MAGIC

MAGIC is a system of two Imaging Atmospheric Cherenkov Telescopes (IACTs) located at the Observatorio Roque de los Muchachos on the Canary Island of La Palma, Spain, which has been operating in stereoscopic mode since 2009. Prior to this, MAGIC comprised a single IACT that began operating in 2004. MAGIC achieves a low trigger threshold of 50−60 GeV, while its sensitivity to point-like sources is 0.8 Crab units (C.U.) above 290 GeV in 50 h (Aleksic et al. 2012).

We observed HESS J1857+026 during July-October 2010, taking 29 hours of data, for zenith angles between 25°−36°. The data were taken in wobble mode by pointing the telescopes at four different pairs of pointing directions symmetric with respect to the source position, in order to obtain a more uniform and flat exposure. Two of these pairs were chosen to be at 0.4° distance from the direction RA: 18^h57^m27^s and Dec: 02°42′60′′, while the other two were 0.5° away.

The results presented here were obtained using the MARS analysis framework (Moralejo et al. 2009) and using the image “sum-cleaning” algorithm (Lombardi et al. 2011), which was found to improve the performance at lower energies, close to 100 GeV. The calibrated photomultiplier signals were cleaned to reduce the effect of the night sky background light, and the resulting image shape and timing information were combined between the two telescopes for each stereoscopic event. The gamma/hadron separation and the event direction reconstruction both made use of the random forest method (Albert et al. 2008). The estimated energy of a given event was determined from the brightness of the shower images, its reconstructed impact parameter and its incidence direction, using Monte Carlo (MC) filled look-up tables. Given that HESS J1857+026 is an extended source our analysis used random forests and look-up tables that were trained using diffuse MC gamma-ray events.

The determination of excess events in the signal region accounts for geometric exposure inhomogeneities by computing the background from the corresponding wobble partner data set and extracting it at the same relative camera plane coordinates as the signal region (further details in Krause 2013). Our skymapping procedure models the background directly from the data and in bins of azimuth, thus taking azimuthal dependencies of the off-axis exposure into account (Lombardi et al. 2011). To compute the flux, we estimated the average effective area using a diffuse MC gamma-ray data sample selected as an annulus with radii 0.25° and 0.65°, in order to account for variations in acceptance across the extent of the source. The resulting spectra were cross-checked in terms of the normalization of the estimated background inside the signal region by event numbers or by effective on-time with respect to the off-source regions. Furthermore, the unfolded spectrum was checked against four different spectral unfolding algorithms (Albert et al. 2007), including the so-called forward folding method, which fits the assumed spectrum by folding it with the energy response matrix and comparing the resulting distribution with the measured distribution of excess events.

Fig. 1 Spectral energy distributions of HESS J1857+026 measured by MAGIC (this work), H.E.S.S (Aharonian et al. 2008) and Fermi-LAT (Neronov & Semikoz 2010; Paneque et al. 2011; Rousseau et al. 2012). The MAGIC data are unfolded to correct for migration and energy biasing effects (Albert et al. 2007), therefore the errors are not independent.

Fig. 2 MAGIC gamma-ray flux map in arbitrary units (a.u.) for events with estimated energy 300 GeV <Eest< 1 TeV (left) and Eest> 1 TeV (right). The median energies of these maps are estimated to be around 480 GeV and 1.7 TeV, respectively. The gamma-ray flux in arbitrary units is calculated as the number of smeared excess events in units of the background flux within 0.1°. Overlaid are test statistic (TS) value contours in steps of 1, starting at 3. They roughly correspond to Gaussian significances. Also shown in the bottom-left corner is the instrumental point spread function (PSF) after the applied smearing. The Fermi-LAT source position is that determined by Rousseau et al. (2012). The markers labeled HESS J1857+026 and HESS J1858+020 are centroid positions determined by Aharonian et al. (2008).

Fig. 3 13CO(J = 1 → 0) line emission integrated intensity (in units of K km s-1) in the vicinity of Hii region U36.40+0.02 (VLSR = 53.3 km s-1) and of the cataloged clouds G036.59-00.06 (VLSR = 53.49 km s-1) and G036.74-00.16 (VLSR = 55.19 km s-1), spanning three ranges of VLSR: 51.21−53.56 km s-1 (left), 53.76−56.11 km s-1 (middle), 57.38−59.72 km s-1 (right). The cyan contours indicate the MAGIC > 1 TeV TS levels in steps of 1, starting at 3. The circle-square green markers indicate Hii regions in the catalog of Anderson et al. (2009). The dashed green circle indicates the position of a possible cavity in the molecular gas.

3. Spectral and morphology results

The MAGIC detection and spectrum of HESS J1857+026 were obtained using a large circular signal extraction region with a radius of 0.4°, chosen to cover the entire extent of the source. Coupling this with a set of loose background rejection cuts with an MC gamma-ray efficiency of 95%, we detected the source above s150 GeV with a significance of 12σ, using Eq. (17) of Li & Ma (1983). The spectrum measured by MAGIC was fitted through forward folding and is well described (χ²/ d.o.f. = 9.0/11) by a single power law of the form dN/ dE = N₀(E/1TeV)^−Γ with index Γ = 2.16 ± 0.07_stat ± 0.15_sys and differential flux at 1 TeV of N₀ = (5.37 ± 0.44_stat ± 1.5_sys) × 10^-12 (TeV^-1,cm^-2 s^-1). The MAGIC spectrum is shown in Fig. 1 together with Fermi-LAT (Neronov & Semikoz 2010; Paneque et al. 2011; Rousseau et al. 2012) and H.E.S.S. (Aharonian et al. 2008) results. Our spectral points have been unfolded to correct for migration and energy biasing effects (Albert et al. 2007), therefore the errors are not independent. The measurement made by MAGIC connects those by H.E.S.S. and Fermi-LAT, providing a continuous coverage of the spectral turnover close to 100 GeV.

The gamma-ray flux sky maps in Fig. 2 show the energy-dependent morphology of HESS J1857+026. Based on MC gamma rays weighted using a spectral index of 2.3, the median energies of the low and high energy maps are around 480 GeV and 1.7 TeV, respectively. The source HESS J1858+020 is not investigated here, given our relatively low exposure at its angular distance from the MAGIC pointing positions. In the estimated energy range of 0.3−1 TeV, we have an instrument point spread function (PSF) of 0.079° and applied a smearing kernel of 0.077° such that the total PSF is 0.11°, while above 1 TeV, our instrument PSF of 0.062° and smearing kernel of 0.05° resulted in a total PSF of 0.08°. Each of the aforementioned numbers corresponds to the 1σ value of a symmetric 2D Gaussian function. We fitted a symmetric 2D Gaussian function to the 0.3−1 TeV gamma-ray flux map and obtained a centroid position of RA: and Dec: which is compatible with the H.E.S.S. centroid position. We also measured an intrinsic source extension of (0.20 ± 0.03_stat ± 0.02_sys)° from the fit, after removing the combined effect of the instrument PSF and the applied smearing.

Above 1 TeV we found that the VHE emission is due to two spatially distinct statistically significant components, which we denote MAGIC J1857.2+0263 and MAGIC J1857.6+0297, as indicated in Fig. 2. Using integration radii of 0.21° for MAGIC J1857.2+0263 and 0.14° for MAGIC J1857.6+0297, we computed significances of 6.7σ and 6.0σ respectively, by applying Eq. (17) of Li & Ma (1983). These radii were chosen to be as large as possible without having any overlap. In order to measure the positions and intrinsic extensions of the two components, we applied an iterative procedure in which we fitted one peak at a time while including the adjacent previously fitted 2D Gaussian function into the background. The elongated source MAGIC J1857.2+0263 was fitted by a noncircular 2D Gaussian which yielded a centroid position RA: and Dec: and intrinsic extensions of (0.17 ± 0.03_stat ± 0.02_sys)° and (0.06 ± 0.03_stat ± 0.02_sys)° along its major and minor axes, respectively, with the major axis having an inclination of (37 ± 6_stat)° counter-clockwise with respect to the RA-axis. MAGIC J1857.6+0297 was fitted by a circular 2D Gaussian which gave a mean position: RA: and Dec: and was compatible with a point source.

We also tried to extract spectra around MAGIC J1857.2+0263 and MAGIC J1857.6+0297 using integration radii of 0.21° and 0.14°, respectively. The limited event statistics above 1 TeV and the inability to resolve the two peaks at lower energies made it difficult to extract spatially-resolved spectra. We estimate that the respective differential fluxes of MAGIC J1857.2+0263 and MAGIC J1857.6+0297 are s45% and s20% of the total differential flux at 1 TeV, which was measured from a circular region with of radius of 0.4° and centered on RA: 18^h57^m27^s and Dec: 02°42′60′′. We also determined a spectral index of Γ = 2.2 ± 0.1_stat for MAGIC J1857.2+0263, which is compatible with that obtained for the entire region.

4. Multiwavelength view

While we expect MAGIC J1857.2+0263 to be associated with PSR J1856+0245, the origin of the VHE emission from MAGIC J1857.6+0297 is unidentified and we analyzed archival multiwavelength data to look for possible counterparts. At radio wavelengths we used the 21 cm continuum and atomic hydrogen (Hi) line emission data from the Very Large Array (VLA) Galactic Plane Survey (VGPS; Stil et al. 2006). At mm wavelengths we used the ¹³CO(J = 1 → 0) molecular line emission data from the Galactic ring survey (GRS; Jackson et al. 2006) to trace the molecular hydrogen (H₂) in this region. Finally, we examined 8 μm images from the Spitzer Galactic Legacy Infrared Midplane Survey Extraordinaire (GLIMPSE) Legacy Project (Benjamin et al. 2003). Further details regarding these data sets and our analysis are given in Appendix A.

From the GRS data cube we computed the average ¹³CO(J = 1 → 0) emission spectrum over Galactic longitudes of 36.0° to 37.0° and Galactic latitudes of −0.5° to 0.5°, finding two prominent peaks at kinematic local standard of rest (LSR) velocity (V_LSR)s55 km s^-1 and s80 km s^-1. We associate these peaks with three molecular gas clouds in the catalog of Roman-Duval et al. (2009): G036.59-00.06, G036.74-00.16 and G036.49-00.16, with kinematic distances of 3.50 kpc, 3.62 kpc and 8.38 kpc, respectively. Most of the emission from G036.49-00.16 is outside the TS = 3 contour of MAGIC J1857.6+0297, which disfavors it as a counterpart. Figure 3 shows the molecular line emission toward MAGIC J1857.6+0297 for a range of velocities that span those of G036.59-00.06 and G036.74-00.16. The range 51.21−53.56 km s^-1 shows mainly the emission from G036.59-00.06, while at larger kinematic distances the range of 53.76−56.11 km s^-1 shows emission from both G036.74-00.16 and G036.59-00.06.

Fig. 4 Left: 13CO(J = 1 → 0) line emission (in units of K) of a single velocity channel at 55.36 km s-1. Middle: GLIMPSE 8μm emission (in units of MJysr-1). Right: VGPS 21 cm continuum emission (in units of K). The cyan contours indicate the MAGIC >1 TeV TS levels in steps of 1, starting at 3. The green circle-square markers indicate Hii regions in the catalog of Anderson et al. (2009), while the white circle marker indicates the Hii region HRDS G036.347+0.048 (Bania et al. 2012). The larger dashed circle marks the position of a possible cavity in the molecular gas.

Star formation takes place in the dense cores of massive molecular clouds, and this leads to a natural connection between molecular clouds and Hii regions. The molecular gas properties of Hii regions were explored by Anderson et al. (2009), and their kinematic distances were determined by Anderson & Bania (2009) using the Galactic rotation curve of McClure-Griffiths & Dickey (2007). We searched these catalogs for Hii regions close to the direction of MAGIC J1857.6+0297, finding two candidate counterparts. The first is the compact Hii region C36.46-0.18, has a radio recombination line velocity of V_LSR = 72.7 ± 2.7 km s^-1 (Lockman 1989) and a fitted molecular gas velocity of V_LSR = 75.99 km s^-1, which corresponds to a distance of 9.0 kpc (Anderson & Bania 2009). Hence it is likely associated with the molecular cloud G036.49-00.16 (V_LSR = 76.87 km s^-1) and we consider it unrelated to the VHE emission. The second candidate is the ultra-compact Hii region U36.40+0.02. Observations of H110α emission from this region determined a velocity of V_LSR = 53.3 ± 1.4 km s^-1 and a corresponding near kinematic distance of kpc (Watson et al. 2003), assuming the Galactic rotation curve of Brand & Blitz (1993). H₂CO absorption line measurements by Watson et al. (2003) indicated the presence of molecular clouds at V_LSR = 52.5 ± 0.2 km s^-1 and at V_LSR = 57.9 ± 0.2 km s^-1. U36.40+0.02 has a fitted molecular gas velocity of V_LSR = 52.42 km s^-1 with corresponding near kinematic distance of 3.3 kpc (Anderson & Bania 2009), and hence is likely to be associated with the cloud G036.59-00.06 (V_LSR = 53.49 km s^-1).

The velocity ranges of 51.21−53.56 km s^-1 and 57.38−59.72 km s^-1 in Fig. 3 were selected to illustrate the two prominent ¹³CO emission peaks in the molecular gas complex surrounding U36.40+0.02. Between them, in the range of 53.76−56.11 km s^-1, we discovered an incomplete shell-like structure coincident with the direction of MAGIC J1857.6+0297. Furthermore, the Galactic longitude-velocity and latitude-velocity plots in Fig. B.1 show a perturbed spatial and velocity structure in the side of the cloud G036.59-00.06, in a region that roughly corresponds to the directions and velocity channels over which this feature is visible. In particular, the latitude-velocity signature indicates that the partial shell spans s3 km s^-1 in velocity space. Figure 4 shows a single GRS velocity channel at 55.36 km s^-1 that provides a close-up view of what could be a gas cavity or wind-blown bubble.

In Fig. 4 we also show the 8 μm emission, which primarily traces emission from polycyclic aromatic hydrocarbons (PAHs) heated by UV starlight. These carbon molecules are found in photodissociation regions (PDRs), which interface the ionized matter in Hii regions and the molecular gas of their parent clouds. The 8 μm emission from U36.40+0.02, which was originally called IRAS 18551+0302 (Beichman et al. 1988), shows a possible diffuse component that extends toward the center of the proposed gas cavity. A second, more compact, component corresponds to the newly discovered Hii region HRDS G036.347+0.048 (Bania et al. 2012). While a kinematic distance has not been determined for this source, its hydrogen radio recombination line (Hnα) spectrum yielded a V_LSR of 76.9 km s^-1, so it is not physically related to U36.40+0.02. A survey of IRAS sources using the high-density gas tracer CS(2 → 1) provided a velocity of V_LSR = 58.0 km s^-1 for IRAS 18551+0302 (Bronfman et al. 1996). This measurement likely traces the gas clump seen in the V_LSR range of 57.38−59.72 km s^-1 of Fig. 3. Ammonia emission, which traces high density (s10⁴ cm^-3) gas, was also detected by the Red MSX Source (RMS) survey (Urquhart et al. 2011), in which the source, MSX6C G036.4057+00.0230, had a NH₃ (1,1) velocity of V_LSR = 57.77 km s^-1. The authors also detected H₂O maser emission, which indicates on-going star formation. This was found over the velocity range 47.2−53.1 km s^-1 with the peak emission at V_LSR = 48.4 km s^-1, which is just below the V_LSR range of 51.21−53.56 km s^-1 shown in Fig. 3.

Fig. 5 Hi line emission integrated intensity (in units of K km s-1) in the vicinity of Hii region U36.40+0.02 (VLSR = 53.3 km s-1) for three ranges of VLSR: 51.72−53.36 km s-1 (left), 54.19−55.83 km s-1 (middle), 57.49−59.13 km s-1 (right). The red contours indicate 13CO(J = 1 → 0) integrated intensity levels of 1 and 2 K km s-1 from the corresponding velocity intervals in Fig. 3. The cyan contours indicate the MAGIC > 1 TeV TS levels in steps of 1, starting at 3. The circle-square green markers indicate Hii regions in the catalog of Anderson et al. (2009). The dashed yellow circle indicates the position of a possible cavity in the molecular gas.

In Fig. 5 we display the Hi 21 cm line emission from the VGPS over three velocity ranges that approximately correspond to those of Fig. 3. The most compact Hi absorption features, produced against strong background continuum emission (Kuchar & Bania 1990), can be either galactic (e.g. Hii regions) or extragalactic in origin. The velocity ranges 51.72−53.36 km s^-1 and 54.19−55.83 km s^-1 exhibit absorption from cold neutral atomic hydrogen that is associated with a molecular cloud absorbing against more distant warm 21 cm emitting gas (e.g. see Roman-Duval et al. 2009). Cold, optically thick Hi can form in the dense parts of molecular clouds through the dissociation of H₂ by penetrating cosmic rays. Goldsmith & Li (2005) found that its resulting self-absorption features tend to correlate with peaks in the ¹³CO(J = 1 → 0) emission. We note, however, that the dark region located within the yellow dashed circle shown in Fig. 5 does not spatially coincide with any bright ¹³CO(J = 1 → 0) emission in the same velocity range, nor with any strong background continuum source (see Fig. 4). This feature may be explained as self-absorption by cold atomic gas components surrounding the molecular cloud. These are expected to form via the photodissociation of H₂ by interstellar UV photons (Jackson et al. 2002). Such optically thick Hi may also trace H₂ that is not detected in ¹³CO, or alternatively, it may point to a strong shock induced depletion of the dust grains, which provide photoelectric heating and facilitate H₂ formation (see Gibson 2010, and references therein). Finally, the feature could also be a real void in the atomic gas. We note that the possible presence of OB stars would lead to a natural connection between such a cavity and the surrounding ¹³CO shell-like feature.

5. Discussion

The broadband SED of HESS J1857+026 in Fig. 1 exhibits a strong turnover at energies close to 100 GeV. This may be attributed to the IC peak expected from a leptonic scenario, although Rousseau et al. (2012) showed that the SED could also be modeled in terms of gamma-ray emission with a hadronic origin. Here we consider whether the proposed PWN scenario can adequately explain the morphology of HESS J1857+026 observed with HESS, MAGIC and Fermi-LAT.

While we find that the source is significantly extended with respect to our PSF in the energy range of 0.3−1 TeV, no significant extension was found by Rousseau et al. (2012) using Fermi-LAT data for energies above 10 GeV. The MAGIC 0.3−1 TeV intrinsic extension of 0.20 ± 0.03_stat however, is comparable to the size of the instrument PSF (68% containment radius) of Fermi-LAT above 10 GeV, and this may explain the nondetection of an intrinsic extension for this source using Fermi-LAT. Indeed Acero et al. (2013) recently specified a 99% confidence level upper limit of 0.28° on the extension of this source. Comparing with H.E.S.S., the intrinsic extension measure of (0.11 ± 0.08_stat)° × (0.08 ± 0.03_stat)° in their energy range is consistent with our 0.3−1 TeV extension measurement. Furthermore, if we only attribute the PWN emission above 1 TeV to MAGIC J1857.2+0263, its size of (0.17 ± 0.03_stat)° × (0.06 ± 0.03_stat)° is in good agreement with the H.E.S.S. result, although it is not clear to what extent the north tail-like emission, which is seen in the H.E.S.S. skymap and may correspond to MAGIC J1857.6+0297, influenced their 2D Gaussian fit.

By itself, the MAGIC energy-dependent morphology would favor a leptonic PWN scenario in which lower energy electrons can diffuse out to larger distances from the pulsar, but only if we consider MAGIC J1857.2+0263 in the >1 TeV map of Fig. 2 as an independent source that becomes less compact at 0.3−1 TeV. This implies that the fainter MAGIC J1857.6+0297 emission, which we only resolve above 1 TeV, is unrelated to PSR J1856+0245. This is because if MAGIC J1857.6+0297 was due to the PWN electrons, its distance from the pulsar would be s43 × (d/ 9 kpc) pc, where d is the distance to PSR J1856+0245, and this would make the entire emission region more extended above 1 TeV than below. Moreover, if MAGIC J1857.6+0297 did indeed trace the same PWN as MAGIC J1857.2+0263, the apparent hole between these emission regions makes it less likely that both originate from the IC scattering of ambient photons on the same electron population.

While we do not consider any substructure within MAGIC J1857.2+0263 to be significant, we note that PSR J1856+0245 is located close to one end of its major axis. Using the observed angular elongation of MAGIC J1857.2+0263 we can estimate the physical length of the PWN emission region to be s40 × (d/ 9 kpc) pc. The VHE gamma-ray emission of a PWN traces a relic electron population accumulated over its lifetime. If we consider the characteristic spin-down age of PSR J1856+0245 (21 kyr) and assume a typical velocity of 400−500 km s^-1, the offset resulting from its proper motion would be at most s10 × (d/ 9 kpc) pc. This maximum value corresponds to the extreme case in which the velocity direction is perpendicular to the line of sight. However, as outlined by Gaensler & Slane (2006), the morphology of middle-aged PWN can be strongly influenced by the complex interaction between the PWN and the reverse shock of its SNR. This can distort the PWN morphology and cause an offset in the position of the pulsar with respect to the center of VHE emission. Therefore the observed morphology of MAGIC J1857.2+0263 may indicate that such evolutionary effects have played some role in its elongated appearance.

The nature of the VHE emission from MAGIC J1857.6+ 0297 remains a mystery, although our targeted multiwavelength study provides some clues to its origin. The shell-like morphology of the ¹³CO channel maps seen in Fig. 3, coupled with the corresponding signatures of perturbed gas in Fig. B.1, points to the presence of a cavity or wind-blown bubble of radius s8 × (d/ 3.7 kpc) pc that coincides with the direction of the VHE gamma-ray emission peak. Here we consider possible scenarios in which the proposed cavity and the nearby ultra-compact Hii region U36.40+0.02, which is embedded in this molecular gas complex, are linked to MAGIC J1857.6+0297. That the peak of the TeV emission does not spatially correspond to the ¹³CO emission peaks in any of the velocity ranges shown in Fig. 3 would argue against a hadronic origin, where molecular hydrogen acts as target matter for the accelerated hadrons. The emission could therefore have a leptonic origin, where accelerated electrons may IC scatter on the IR to UV photon fields of U36.40+0.02. We do not, however, exclude a hadronic origin due to the possible presence of cold, dense atomic gas within the proposed ¹³CO cavity, which could also act as target material for accelerated hadrons (e.g. see Fukui et al. 2012).

A leptonic scenario that may explain the radio data has the VHE emission originating from a PWN whose progenitor star was in an OB association formed within U36.40+0.02, the stellar winds of which created the cavity. As outlined in de Jager & Djannati-Ataï (2008), this could then facilitate an initial unimpeded expansion of the PWN, allowing for weaker magnetic fields that would result in lower levels of synchrotron cooling for the VHE gamma-ray emitting electrons. A similar scenario was also proposed by Bock & Gvaramadze (2002) to explain the possible association of the pulsar PSR B1706-44 with SNR G343.1-2.3, which are candidate counterparts of the extended source HESS J1708-443 (Abramowski et al. 2011). Bock & Gvaramadze (2002) suggested that in this case, the massive progenitor star escaped its molecular cloud and produced a wind-blown bubble in the less dense ISM together with a wind-driven cavity in the parent cloud. They also speculated that the subsequent SN blast wave then interacted on one side with the shell of the bubble and with the parent cloud on the other side. In the case of MAGIC J1857.6+0297, additional high-resolution multiwavelength data at radio wavelengths and in X-rays will be necessary to fully test such a scenario.

An alternative source of particle acceleration to explain the observed VHE emission could be the combined effect of the strong stellar winds themselves. For instance, VHE gamma rays may be produced by particles accelerated in the colliding winds of early-type (O, early B, Wolf-Rayet) star binary systems, although Reimer et al. (2006) showed that the spectra of such systems are expected to suffer from strong absorption above s50 GeV due to pair production in the surrounding photon fields. Given that U36.40+0.02 and the surrounding dense molecular gas host regions of on-going star formation, as evidenced by the H₂O maser emission, the VHE emission could instead stem from particles accelerated to multi-TeV energies at the shocks produced by outflows from massive protostars (Araudo et al. 2008; Bosch-Ramon et al. 2010). Since the kinetic luminosity of one such jet is expected to be around 10³⁶ erg s^-1, a collection of these objects would also be required to explain the VHE emission.

6. Conclusions

MAGIC has performed observations of the TeV PWN candidate source HESS J1857+026 in order to better understand its complex morphology. We have also extended the VHE spectrum of the source down to 100 GeV, thus bridging the gap between the earlier spectral measurements of H.E.S.S. and those obtained more recently with Fermi-LAT. Assuming that the observed VHE gamma-ray emission is produced by energetic electrons via the IC channel, their energy-dependent diffusion would dictate a more extended PWN about PSR J1856+0245 when going down in energy from the multi-TeV regime to the range covered by Fermi-LAT. Our study of the energy-dependent morphology has revealed for the first time that this is is not the case, with the source becoming more extended above 1 TeV. Based on this finding we have proposed a two-source scenario for the separated and statistically significant peaks in emission that we detected above 1 TeV, which we have called MAGIC J1857.2+0263 and MAGIC J1857.6+0297. We then interpret MAGIC J1857.2+0263 as the emission from the relic PWN of PSR J1856+0245. We expect that the elongated shape of the emission above 1 TeV as well as its offset position from the pulsar are signatures of complex PWN evolutionary processes, which involve an interaction with the SNR reverse shock, while the initial pulsar velocity could also play a role.

We have used archival multiwavelength data at radio and infra-red wavelengths to search for possible counterparts of MAGIC J1857.6+0297. We found a possible association with an ultra-compact Hii region located at 3.7 kpc and its neighboring molecular clouds. The morphology of molecular gas in this region suggests that a cavity or wind-blown bubble may be present just behind the Hii region and in the same direction as the peak of the VHE emission. Assuming that MAGIC J1857.6+0297 is indeed associated with this cloud complex, the overall lack of overlap with the most dense parts of the molecular clouds hints at a leptonic origin for the VHE gamma-ray emission. A hadronic origin is not excluded, however, due to the possible presence of an atomic gas cloud within the proposed cavity. We have outlined alternative scenarios that could explain the proposed link between the gamma-ray emission and this molecular cloud complex. They include a PWN whose progenitor star created a cavity in the cloud prior to its SN explosion, particle acceleration based on the combined effects of strong stellar winds from massive stars, or the acceleration of particles in the outflows of massive protostars. More targeted multiwavelength data will be needed, especially in X-rays, in order to characterize this interesting region and test whether it really is the source of the TeV gamma rays.

Acknowledgments

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 CPAN CSD2007-00042 and MultiDark CSD2009-00064 projects of the Spanish Consolider-Ingenio 2010 programme, by grant 127740 of the Academy of Finland, by the DFG Cluster of Excellence “Origin and Structure of the Universe”, by the Croatian Science Foundation Project 09/176, by the DFG Collaborative Research Centers SFB823/C4 and SFB876/C3, and by the Polish MNiSzW grant 745/N-HESS-MAGIC/2010/0.

References

Appendix A: Multiwavelength data

Here we provide a technical description of the multiwavelength data used in Sect. 4. The VGPS Hi data (Stil et al. 2006) consists of continuum images with angular resolution of 1′ and an rms noise of 0.3 K, while the Hi line images also have an angular resolution of 1′, a spectral resolution of 1.56 km s^-1 and an rms noise of 0.3 K per channel (0.824 km s^-1). For the ¹³CO(J = 1 → 0) data, the GRS survey Jackson et al. (2006) employed a grid spacing of 22′′, and had an angular resolution of 46′′ and a spectral resolution of 0.212 km s^-1, with rms noise of s0.13 K. Finally, the 8μm GLIMPSE data (Benjamin et al. 2003) has arcsecond-scale resolution.

Our analysis of the radio data cubes used the Miriad (Sault et al. 1995) package, while the images in this section were generated using the Karma (Gooch 1996) and DS9 (Joye & Mandel 2003) software packages.

Appendix B: Avergae ¹³CO(J = 1 → 0) emission

Figure B.1 shows Galactic longitude-velocity and latitude-velocity plots of ¹³CO emission in the direction of a possible gas cavity that is discussed in Sect. 4.

Fig. B.1 Galactic latitude-velocity (left) and longitude-velocity (right) plots showing 13CO(J = 1 → 0) line emission averaged over the range of longitudes: 36.21°−36.45° (left) and latitudes: −0.06°−+0.21° (right). The green rectangles indicate the approximate location of a possible gas cavity. Their angular positions correspond to the extent of the dashed circle in Fig. 3 while in velocity they span the VLSR range of 53.23−56.21 km s-1. The white contours indicate levels of 0.3 K.

All Figures

	Fig. 1 Spectral energy distributions of HESS J1857+026 measured by MAGIC (this work), H.E.S.S (Aharonian et al. 2008) and Fermi-LAT (Neronov & Semikoz 2010; Paneque et al. 2011; Rousseau et al. 2012). The MAGIC data are unfolded to correct for migration and energy biasing effects (Albert et al. 2007), therefore the errors are not independent.
In the text

Fig. 2 MAGIC gamma-ray flux map in arbitrary units (a.u.) for events with estimated energy 300 GeV <Eest< 1 TeV (left) and Eest> 1 TeV (right). The median energies of these maps are estimated to be around 480 GeV and 1.7 TeV, respectively. The gamma-ray flux in arbitrary units is calculated as the number of smeared excess events in units of the background flux within 0.1°. Overlaid are test statistic (TS) value contours in steps of 1, starting at 3. They roughly correspond to Gaussian significances. Also shown in the bottom-left corner is the instrumental point spread function (PSF) after the applied smearing. The Fermi-LAT source position is that determined by Rousseau et al. (2012). The markers labeled HESS J1857+026 and HESS J1858+020 are centroid positions determined by Aharonian et al. (2008).

In the text

Fig. 3 13CO(J = 1 → 0) line emission integrated intensity (in units of K km s-1) in the vicinity of Hii region U36.40+0.02 (VLSR = 53.3 km s-1) and of the cataloged clouds G036.59-00.06 (VLSR = 53.49 km s-1) and G036.74-00.16 (VLSR = 55.19 km s-1), spanning three ranges of VLSR: 51.21−53.56 km s-1 (left), 53.76−56.11 km s-1 (middle), 57.38−59.72 km s-1 (right). The cyan contours indicate the MAGIC > 1 TeV TS levels in steps of 1, starting at 3. The circle-square green markers indicate Hii regions in the catalog of Anderson et al. (2009). The dashed green circle indicates the position of a possible cavity in the molecular gas.

In the text

Fig. 4 Left: 13CO(J = 1 → 0) line emission (in units of K) of a single velocity channel at 55.36 km s-1. Middle: GLIMPSE 8μm emission (in units of MJysr-1). Right: VGPS 21 cm continuum emission (in units of K). The cyan contours indicate the MAGIC >1 TeV TS levels in steps of 1, starting at 3. The green circle-square markers indicate Hii regions in the catalog of Anderson et al. (2009), while the white circle marker indicates the Hii region HRDS G036.347+0.048 (Bania et al. 2012). The larger dashed circle marks the position of a possible cavity in the molecular gas.

In the text

Fig. 5 Hi line emission integrated intensity (in units of K km s-1) in the vicinity of Hii region U36.40+0.02 (VLSR = 53.3 km s-1) for three ranges of VLSR: 51.72−53.36 km s-1 (left), 54.19−55.83 km s-1 (middle), 57.49−59.13 km s-1 (right). The red contours indicate 13CO(J = 1 → 0) integrated intensity levels of 1 and 2 K km s-1 from the corresponding velocity intervals in Fig. 3. The cyan contours indicate the MAGIC > 1 TeV TS levels in steps of 1, starting at 3. The circle-square green markers indicate Hii regions in the catalog of Anderson et al. (2009). The dashed yellow circle indicates the position of a possible cavity in the molecular gas.

In the text

Fig. B.1 Galactic latitude-velocity (left) and longitude-velocity (right) plots showing 13CO(J = 1 → 0) line emission averaged over the range of longitudes: 36.21°−36.45° (left) and latitudes: −0.06°−+0.21° (right). The green rectangles indicate the approximate location of a possible gas cavity. Their angular positions correspond to the extent of the dashed circle in Fig. 3 while in velocity they span the VLSR range of 53.23−56.21 km s-1. The white contours indicate levels of 0.3 K.

In the text