© ESO 2021

1 Introduction

HESS J1857+026 was discovered by the H.E.S.S. observatory during a survey of the Galactic plane (Aharonian et al. 2008). The source was clearly identified as extended. These latter authors obtained a centroid for the TeV emission at (l, b) = (35.96, −0.06), semi-major and -minor axes of about 0.11° and 0.08°, respectively, and some emission projecting towards the north. At present, HESS J1857+026 has no clear counterpart in other spectral bands. It harbours the Vela-like radio pulsar PSR J1856+0245, located ~8′ offset from the centroid of the γ-ray emission. This latter is a 81 ms pulsar discovered by the Arecibo PALFA survey of the Galactic plane. The timing analysis in the radio domain shows that the pulsar has a spin-down age of ~ 21 kyr and a spin-down energy Ė ~ 4.6 × 10³⁶ erg s⁻¹ (Hessels et al. 2008). A distance of ~ 6.3 kpc was obtained using the dispersion measure (DM) of 624 cm⁻³ pc and the electron density model YMW16 (Yao et al. 2017) (we refer to the ATNF Pulsar Catalog; Manchester et al. 2005). Hessels et al. (2008) suggest that HESS J1857+026 could be a pulsar wind nebula (PWN) powered by PSR J1856+0245, based on the position and energetics of the pulsar, which translates into a γ-ray efficiency ${ϵ}_{\gamma }=L_{\gamma }/\stackrel{.}{E,}$ similar to other TeV PWN/pulsar associations. This scenario is supported by the presence of the ASCA X-ray source coincident with the radio pulsar. Regardless of the low number of counts, the authors derived a X-ray efficiency (ϵ_X = L_X∕Ė) similar to other Vela-like pulsars powering synchrotron nebulae. Nice et al. (2013) performed a detailed study of the X-ray source using data from the XMM-Newton and Chandra telescopes. These authors found a point source coincident with PSR J1856+0245 and no evidence of extended emission related to a PWN. The spectrum is likely non-thermal and highly absorbed (indicating a distant source) and would be produced by particles accelerated at the pulsar’s magnetosphere.

In the GeV band, HESS J1857+026 was first detected by Fermi-LAT as a point source with a centroid coincident with the H.E.S.S. position (Rousseau et al. 2012). No pulsed emission from PSR J1856+0245 was reported in the Fermi data. The broad band spectral analysis (from the X-ray to the TeV bands) reveals that the leptonic scenario is compatible with the pulsar’s age and produces an acceptable fit assuming a low magnetic field (~ 4 μG). The hadronic mechanism for γ-ray production also gives an adequate fit, but requires a very high injection energy. More recently, Ackermann et al. (2017) analysed 6 yr of observations with Fermi-LAT above 10 GeV and found evidence of extended emission towards HESS J1857+026. Even if the extended nature of the source is undoubted, the large extent derived from the spatial analysis (~ 0.6°), which is indeed larger than the TeV extension, could be affected by two nearby hot spots likely unrelated with HESS J1857+026.

The most detailed study of HESS J1857+026 in the TeV band was presented by MAGIC Collaboration (2014), using data obtained with the MAGIC telescope. These authors find that in the 0.3−1 TeV energy band the emission is well fitted with a 2D Gaussian function with a centroid coincident with that of H.E.S.S., while above 1 TeV the source is resolved into two distinct sources, labelled MAGIC J1857.2+0263 (hereafter MAG1) and MAGIC J1857.6+0297 (hereafter MAG2). The authors propose a two-source model for HESS J1857+026, in which the γ-ray appearance is produced by the superposition in the line of sight of two unrelated sources located at two different distances.

On the one hand, the pulsar PSR J1856+0245 appears projected over MAG1, which presents an elongated shape. The authors conclude that the pulsar is powering a PWN, which has the expected energy-dependent morphology in the very high-energy domain: the lower-energy electrons diffuse further from the pulsar before cooling through inverse Compton scattering, causing a smaller extension of the nebula for energies > 1 TeV. Therefore, according to MAGIC Collaboration (2014), MAG1 would be a TeV PWN located at ~ 9 kpc, as they used the out-of-date NE2001 electron density model of Cordes & Lazio (2002) to derive the DM distance of the pulsar. The elongated morphology of the PWN could be a consequence of the interaction with the reverse shock of the (so far undetected) supernova remnant (SNR) that was created together with the pulsar.

Regarding MAG2, which is resolved only above 1 TeV as a point source, its nature remains unclear. In the direction of this source, MAGIC Collaboration (2014) found molecular clouds, HII regions, high-density gas, and H₂O maser emission; they also report on the detection of an incomplete shell of CO in the 54− 56 km s⁻¹ velocity range that partially overlaps MAG2. The ultracompact HII region U36.40+0.02 (located at ~3.7 kpc) is probably related to this molecular structure. Based on the presence of HII regions, together with the detection of H₂O masers, massive stars and star forming activity are expected to be present in the field. These latter authors propose different scenarios to explain the nature of MAG2. TeV emission could arise from particles accelerated in the colliding winds of massive stars or in the shocks produced by the outflows of protostars. Alternatively, the TeV radiation could be produced by electrons accelerated by an unidentified PWN evolving inside the cavity blown by the wind of its progenitor star and/or by the HII region U36.40+0.02. As the peak of the TeV emission does not match the CO peaks, a leptonic mechanism for γ-ray production is favoured with respect to a hadronic scenario.

It is worth noting that the pulsar PSR J1857+0300 lies projected onto the TeV emission from HESS J1857+026, in the direction of MAG2. Discovered by the Arecibo telescope (Lyne et al. 2017), it is an old (τ_c ~ 4.6 × 10⁶ yr) and low-energy (Ė ~ 2.5 × 10³² erg s⁻¹) neutron star. From the DM of ~691 pc cm⁻³, a distance of ~6.7 kpc is obtained using the YMW16 electron density model.

So far, the nature of HESS J1857+026 remains unclear and the two-source scenario of MAGIC Collaboration (2014) needs to be further investigated. For this purpose, we obtained high-quality radio observations using the Karl G. Jansky Very Large Array¹ (VLA) covering thewhole extension of the source. We plan to detect the radio counterpart of the putative TeV PWN powered by PSR J1856+0245 and the shell of the associated SNR, allowing us to characterise the sources that may be accelerating particles and producing γ-rays. We also investigate the distribution of the neutral and molecular gas using data taken from public surveys to look for dense material that might correlate with the γ-ray emission and that could be acting as an efficient proton target. Additionally, we plan to find evidence in the interstellar medium (ISM) of explosive and/or expanding events in the region. This will allow us to test the most recent models of very-high-energy emission, which suggest that TeV-emitting particles are efficiently accelerated by the combined action of stellar winds and SNRs in the interior of superbubbles (SBs).

2 Data

2.1 New radio observations

HESS J1857+026 was observed with the VLA in the C configuration (project ID 18B-001) using a mosaicing technique with 14 different pointings following a hexagonal pattern. The data were taken at L-band using the wide-band 1 GHz receiver system centred at 1.5 GHz which consists of 16 spectral windows, each with a bandwidth of 64 MHz, spread into 64 channels. The source 3C 48 was used as a primary flux density and bandpass calibrator, while the phase was calibrated with J1822+0938. We also performed observations with the VLA in the C configuration towards a field of 0.4° × 0.4° in the direction of PSR J1856+0245 using the mosaicing technique with 14 different pointings. These observations were planned to detect the putative PWN associated with the pulsar. In this case, we used the wide-band 4− 8 GHz receiver system. We requested the 3-bit sampler with two sub-bands of 2 GHz, each comprising 16 spectral windows with a bandwidth of 128 MHz each, spread into 64 channels. The source 3C 48 was used as flux density and bandpass calibrator, and J1922+1530 as a phase calibrator. Table 1 summarises the observations.

All data were processed with the VLA CASA Calibration Pipeline. We applied extra flagging before imaging to improve the quality of the calibrated data. The final images were obtained with the TCLEAN task of the CASA software using the multi-frequency deconvolution algorithm for both spectral bands. At L band we found structures with different spatial scales, so we used a multiple scale algorithm in the deconvolution procedure. The final image at L band has a resolution of 17′′.38 × 15′′.51, an angle of 6°.87, and an effective noise of 0.7 rmmJy beam⁻¹. The final image at C band has a resolution of 3′′.55 × 2′′.62, an angle of − 50°.8, and an effective noise of 12 μJy beam⁻¹.

Fig. 1 Radio continuum emission at 1.5 GHz. The white crosses indicate the positions of the pulsars PSR J1856+0245 and PSR J1857+0300. The white square shows the region mapped in the C band with the VLA and the contours represent the TeV emission below 1 TeV extracted from H. E. S. S. Collaboration (2018) with values of 1.5, 2.0, 2.7, 3.5, and 4.5 in units of 10−13 ph cm−2 s−1. The angular resolution is 17.′′38 × 15.′′ 51, PA = 6°.87 and the noise is 0.7 mJy beam−1. The color scale is expressed in Jy beam−1.

2.2 The surrounding medium

The neutral hydrogen (HI) data were extracted from the Very Large Array Galactic Plane Survey (VGPS, Stil et al. 2006), which maps the 21 cm line emission with angular and spectral resolutions of 1′ and 1.3 km s⁻¹, respectively.We also studied the distribution of the molecular gas using the emission of the ¹³CO extracted from the Galactic Ring Survey (GRS, Jackson et al. 2006). The survey maps the Galactic ring in the ¹³CO J = 1−0 line with angular and spectral resolutions of 46′′ and 0.2 km s⁻¹, respectively. Besides, we used the ¹²CO J = 1−0 emission from the FOREST Unbiased Galactic Plane Imaging (FUGIN) survey² (Umemoto et al. 2017). The angular and velocity resolutions are 20′′ and 0.65 km s⁻¹, respectively.

3 Results

3.1 The new VLA images

The new radio image at 1.5 GHz is shown in false colour in Fig. 1 with contours of the γ-ray emission for energies <1 TeV. The positions of the pulsars PSR J1856+0245 and PSR J1857+0300 are indicated with white crosses. There is no evidence, down to the noise level of the image, of enhanced emission in the direction of PSR J1856+0245 that could represent a probable PWN, or of the surrounding radio shell associated with the host SNR.

As can be seen from Fig. 1, several extended sources are detected in the field. To the north, we detect the HII region G036.459−00.183, which has a systemic velocity of ~73 km s⁻¹ and a kinematic distance of ~8.9 kpc (refer to the WISE Catalog of Galactic HII Regions³ and references therein). To the south, we identify the SNR G35.6−0.4 and a complex of HII regions. This complex has a systemic velocity of ~ 52 km s⁻¹ and a kinematic distance of ~10.4 kpc (see the aforementioned catalogue of HII regions). In Fig. 2 we display a zoomed view of the area.

This is another case in which the new dedicated radio observations carried out with higher sensitivity than current surveys do not show any indication of diffuse emission around the pulsar. In Sect. 4, we discuss the possible reasons for this. In principle, we rule out that it is a consequence of the limitations of the interferometer when observing extended sources. On the one hand, the largest angular scale of the VLA in the C configuration at L band is ~ 8′, similar to the extent of the source MAG1 (10′ × 4′), and so we would expect no significant losses of flux from a putative PWN of a similar size. On the other hand, considering the large extent of HESS J1857+026 with respect to the primary beam (θ_PB ~ 17′ FWHM at L band), the mosaicing technique guarantees a constant sensitivity over the field, avoiding the attenuation of the primary beam.

The radio continuum emission at 6.0 GHz is shown in Fig. 3. The new high-sensitivity radio map does not reveal any extended emission around PSR J1856+0245. We detect a point source at the position of the pulsar, with a flux density S_{6.0 GHz} ~ 0.23 mJy. A zoomed-in view of this region is shown in the upper left corner of Fig. 3. The flux of the pulsar at 1.4 GHz was measured by Nice et al. (2013), who obtained S_{1.4 GHz} = 0.58 ± 0.02 mJy. We note that the flux is belowthe sensitivity of our 1.5 GHz observation (Table 1), explaining the non-detection of the pulsar in the new VLA image of Fig. 1. Nice et al. (2013) also derived a spectral index α = 0.4 ± 0.2 (S_ν ∝ ν^−α) by fitting the emission of their four 100 MHz sub-bands. We estimate the spectral index by combining the flux at 1.4 GHz from Nice et al. (2013) and the flux at 6.0 GHz from our observations. We obtain α = 0.64 ± 0.06 for PSR J1856+0245.

3.2 Neutral gas distribution

We investigated the neutral gas distribution towards HESS J1857+026, looking for: (i) high-density gas that correlates with the TeV emission and could be acting as an efficient proton ‘target’ for hadronic production of γ-rays, and (ii) the vestiges in the ISM of the action of expansive or explosive events, such as the explosion of the SN that could have created the pulsars PSR J1856+024 and PSR J1857+0300.

After inspecting the complete HI cube, we found a cavity-like structure surrounded by an incomplete shell of neutral gas in the velocity range between 81 and 100 km s⁻¹. Figure 4 shows the integration of the HI emission in this interval, together with contours of the TeV emission from HESS J1857+026, and the positions of the pulsars PSR J1856+0245 and PSR J1857+0300. Interestingly, the size and morphology of the cavity approximately match those of the TeV emission. Both of them display an elliptical shape elongated in the direction parallel to the Galactic plane. Adopting a systemic velocity of 91 km s⁻¹ for the HI structure, and using the Galactic rotation model of Persic et al. (1996) with the parameters of Reid et al. (2014), we obtained near and far distances of 5.5 and ~8.3 kpc, respectively. Taking into consideration the distance uncertainties of the kinematical method, the near distance is compatible with the DM distances of the pulsars PSR J1856+0245 (~ 6.3 kpc) and PSR J1857+0300 (~ 6.7 kpc).

Without any possibility of discriminating between the two kinematic distances, we assume as a working hypothesis that the HI cavity and shell were swept by the expansion of the SNR created together with the pulsar and/or by the action of the winds of its progenitor star. Then, following the association between the pulsars and the HI structure, we adopt 5.5 kpc as the most probable distance to the cavity of neutral gas.

To characterise this HI structure, we estimate the neutral gas mass integrating the HI column density between 81 and 100 km s⁻¹. Assuming that the emission is optically thin, the HI column density is given by N_HI = 1.8 × 10¹⁸ ∫ (ΔTdv) cm⁻², where ΔT is the mean brightness temperature of the integration region and dv the velocity interval of integration. The mass is therefore $$M=1.3\times {10}^{-3}{D}_{\text{kpc}}^2\Omega \Delta T\text{d}v M_{\odot },$$(1)

where D_kpc is the distance in kiloparsecs and Ω is the solid angle of the integration region.

We approximate the shape of the shell with an elliptical ring with inner and outer boundaries of 35.5′ × 18′ and 55.8′ × 35′, respectively.We note that the inner boundary defines the outer limit of the HI cavity. Then, the mean HI brightness temperature of the swept shell integrated between 81 and 100 km s⁻¹ is T_s = 1265 K. The emission from the shell must be corrected with the background emission, that is, the emission of the ISM prior to the formation of the shell. It follows that ΔT_s = T_s − T_bg, where T_bg is the background temperature. A reasonable assumption to obtain T_bg is to consider it equal to the temperature of the unperturbed ISM just outside the shell. From Fig. 4we note that the HI distribution around the shell is very inhomogeneous and it is difficult to choose an appropriate background region. Therefore, we followed the method of Suad et al. (2019) to derive T_bg, which assumes that all the mass evacuated from the cavity (M_cav) is accumulated in the shell (M_s) and that ΔT_cav = T_bg − T_cav. For the HI cavity, we derived a mean HI brightness temperature T_cav = 970 K. Therefore, using Eq. (1), the condition M_cav = M_s gives (T_s − T_bg)Ω_shell = (T_bg − T_cav)Ω_cav. The latter can be solved to obtain T_bg. Taking T_s = 1265 K, T_cav = 970 K, Ω_s = 3370 arcmin², and Ω_cav = 2010 arcmin², we obtain T_bg = 1155 K, and therefore ΔT_s = 110 K.

Using the mean HI brightness temperature of the shell, the column density of the neutral hydrogen is N_HI ~ 1.6 × 10²⁰ cm⁻² and the mass of the swept shell is M_s ~ 1.6 × 10⁴ M_⊙. To estimate the original ambient density n₀, we performed the integration of Eq. (1) over the cavity area using the background temperature, which is the temperature of the unperturbed ISM prior to the expansion of the HI shell. The obtained mass is distributed in the volume of the cavity, which is approximated by a sphere of radius $R_{\text{eff}}=\sqrt{ab}$, where a and b are the mean major and minor semi-axis of the shell. Then, for R_eff = 56 pc, we obtain an original ambient density n₀ ~ 7.5 cm⁻³. The kinetic energy is derived from $E_{\text{k}}=1/2M_{\text{s}}{v}_{\text{exp}}^2$ and the dynamical age from t_dyn = 0.6R_eff∕v_exp, where v_exp is the expansion velocity of the shell. As we show later, the cavity was probably created by stellar winds, so a factor 0.6 is included according to evolutionary models of interstellar bubbles (Weaver et al. 1977). We take v_exp = 10 km s⁻¹, which is half the velocity range where the HI shell is detected. This expansion velocity is a lower limit, because confusion from unrelated foreground or background emission hampers the detection of the caps of the expanding shell. We obtain E_k ~ 1.6 × 10⁴⁹ erg and t_dyn ~ 3.3 × 10⁶ yr. Generally, this method produces large uncertainties of the parameters (up to 40%) due to the error in the distance, the choice of the background level, and the integration limits. Despite the errors involved, the age of the HI shell is two orders of magnitude greater than the age of PSR J1856+0245. The parameters that characterise the shell are summarised in Table 2.

The derived morphological and physical characteristics of the HI shell are similar to those found in Galactic SBs (Suad et al. 2014, 2019). The origin of these structures is generally attributed to the action of the stellar winds of massive stars and their subsequent explosion as SNe. As neither associations of OB stars nor other types of early stars were found in the literature within the limits of the cavity, we analyse the types of stars that could have formed the SB. We compare the kinetic energy E_k of the SB with the mechanical energy E_w injected by the stellar wind into the ISM during the lifetime of the bubble. The latter is estimated from the relation E_w = ϵ_wL_wt_dyn, where ϵ_w is the fraction of the stellar wind that is converted into mechanical energy and L_w is the mechanical luminosity of the wind, defined as $L_{\text{w}}=1/2\dot{M}{v}_{\text{w}}^2$, where Ṁ is the mass-loss rate of the star and v_w is the terminal velocity of the wind. Taking typical values of the stellar wind of an O-type star (Ṁ = 2 × 10⁻⁶ M_⊙ yr⁻¹, v_w = 2000 km s⁻¹; Prinja et al. 1990; Mokiem et al. 2007), we obtain L_O-star ~ 2.5 × 10³⁶ erg s⁻¹. Considering a canonical ϵ_w = 0.2 (Weaver et al. 1977), the total energy released during the lifetime of the SB (t_dyn = 3.3 Myr) that is converted into mechanical energy is E_w ~ 5 × 10⁴⁹ erg, which exceeds the kinetic energy of the shell. This result supports a stellar wind origin for the SB, showing that a single massive star is energetic enough to evacuate the HI cavity and accumulate the material in the shell.

Fig. 2 Magnification of the southern region of the image obtained at 1.5 GHz. The circles represent the positions of the HII regions taken from the WISE Catalog of Galactic HII Region. The synthesized beam is represented in the upper left corner. The color scale is expressed in Jy beam−1.

Fig. 3 Radio continuum emission at 6.0 GHz. The resolution is 3.′′55 × 2.′′ 62 while the effective noise level is 12 μJy beam−1. The color scale is expressed in Jy beam−1.

Fig. 4 HI distribution towards HESS J1856+026, integrated inthe velocity range between 81 and 102 km s−1. The contours correspond to the TeV emission and the crosses mark the position of the pulsars PSR J1856+0245 and PSR J1857+0300. The color scale is expressed in arbitrary units.

3.3 Molecular gas distribution

We inspected the distribution of the ¹³CO J = 1–0 in a wide area surrounding the very-high-energy source and in the velocity interval where the HI shell spans. We found several molecular condensations between 78 and 90 km s⁻¹, which seem to form an incomplete shell around the brightest γ-ray emission (see Fig. 5). We have labelled five molecular components: cloud A is located to the east of the TeV peak; clouds B and C, to the north; cloud D is a compact feature to the southwest; and cloud E is a large feature partially overlapping the TeV emission.

To derive the physical properties of the molecular clouds, we consider that the gas is at local thermodynamic equilibrium (LTE). The ¹³CO J = 1–0 line is assumed to be optically thin, and we use it to derive the H₂ column density N(H₂). The excitation temperature T_ex was obtained from the brightness temperature $T_{\text{B}}^{\text{12CO}}$ of the ¹²CO J = 1–0 (taken from the FUGIN survey), which is assumed to be optically thick. It follows that $T_{\text{ex}}=5.5{[ln(1+5.5(T_{\text{B}}^{\text{12CO}}+0.82))]}^{-1}$, where all the temperatures are expressed in K. We obtained a mean $T_{\text{B}}^{\text{12CO}}~4\text{K}$ over clouds A, B, C, and D, which yields T_ex ~ 7.7 K. For cloud E, we obtain $T_{\text{B}}^{\text{12CO}}~6.5\text{K}$ and T_ex ~ 9.8 K. The ¹³CO column densities and H₂ masses were derived using Eqs. (1)–(3) of Petriella et al. (2012). We took a background temperature T_b = 2.7 K, a H₂ to ¹³CO relative abundance of 5 × 10⁵ (Ortega et al. 2017), and a helium abundance of 25%. The projected shapes of the clouds were approximated by ellipses. To estimate the number densities, we approximated the cloud volume by ellipsoids with a third axis equal to the mean of the ellipse semi-axis. We use a distance of 5.5 kpc, coincident with the distance of the HI cavity. The results are listed in Table 3.

4 Discussion

The source first identified as a possible origin of γ-ray emission from HESS J1857+026 is the pulsar PSR J1856+0245. However, the complex morphology in the TeV band suggests that it could be formed by two different astrophysical objects (Aharonian et al. 2008). On the one hand, MAGIC Collaboration (2014) suggest that MAG1 could be a PWN powered by PSR J1856+0245 and located at ~ 9 kpc. On the other hand, MAG2, which is only resolved above 1 TeV, appears as an extension of the most intense TeV emission from MAG1 and has no clear counterpart in other wavelengths, but it may be associated with a cavity of CO and a UC HII region located at ~ 3.7 kpc. This two-source scenario resembles that of VER J1907+062, whose complex morphology in the TeV band was explained by the superposition of two unrelated γ-ray sources located at two different distances (Crestan et al. 2021; Duvidovich et al. 2020).

Our analysis of the distribution and kinematics of the atomic gas revealed the existence of a SB towards HESS J1857+026 in the velocity range between 81 and 100 km s⁻¹. Superbubbles are formed by the injection of a large sum of energy from stellar winds and SN explosions and were proposed as alternative sources of Galactic cosmic rays, in contrast to isolated SNRs (see e.g. Parizot et al. 2004; Binns et al. 2007). Indeed, the collective action of massive stars and SNRs within SBs could energise protons up to PeV energies, making them the most promising pevatron candidates (Aharonian et al. 2019). In the context of the discovery of a SB coincident with HESS J1857+026, we analysed the possible contributions of all astrophysical objects to the γ-ray emission.

We obtained a kinematic distance of ~5.5 kpc for the SB, which is compatible, within the errors, with the DM distance of PSR J1856+0245 (~ 6.3 kpc). The TeV emission appears completely immersed within the HI cavity and both of them present elongated shapes in the direction parallel to the Galactic plane. The spatial and morphological match between the SB and the entire TeV emission from HESS J1857+026 suggests an association between them, arguing in favour of a single γ-ray source confined inside the HI cavity and located at a distance of ~5.5 kpc. Even if there is no unequivocal evidence of the physical association between HESS J1857+026 and the neutral gas structure, we hypothesise that the SB may have been created by the stellar winds of massive stars located near its centre, with the possible contribution of the SNRs formed after the explosion of at least one of these stars. Such an SN event is suggested by the presence of the pulsar PSR J1856+0245.

In our scenario, PSR J1856+0245 is located inside the SB and part of the γ-ray emission from HESS J1857+026 may correspond to a TeV PWN powered by the pulsar. The non-detection of a radio PWN resembles other energetic pulsars powering PWNe detected in the high- and very high-energy domains but not observed in the radio band. Examples of such pulsars are PSR J1907+0602 (Duvidovich et al. 2020), PSR J1826−1256 (Duvidovich et al. 2019), and PSR J1809−1917 (Castelletti et al. 2016); all of them, like PSR J1856+0245, have magnetic fields of a few 10¹² G, as commonly found in young pulsars.

Despite the non-detection of a radio PWN around PSR J1856+0245, we use the radio observations to put some constraints on its characteristics. A rough estimate of the flux density of the PWN at 1.4 GHz is calculated from $F_{1.4}[\text{mJy}]=2.15\times {10}^8{ϵ}_R{\dot{E}}_{37}/D_{\text{kpc}}^2$ (Gaensler et al. 2000), where ϵ_R is the efficiency of conversion of the pulsar power into radio emission, Ė₃₇ is the spin-down power of the pulsar inunits of 10³⁷ erg s⁻¹, and D_kpc the distance in kiloparsec (taken equal to the SB distance of 5.5 kpc). We obtain F_1.4 ~ 330 mJy. We assume that the radius of the radio PWN is ~0.1° (= 10 pc for a distance of 5.5 kpc), similar to the spatial extension of the TeV source MAG1. Taking a typical ϵ_R = 10⁻⁴ (Gaensler & Slane 2006), the expected surface brightness at 1.5 GHz is S_1.5 = F_1.5∕A ~ 0.6 mJy beam⁻¹. The latter was obtained after expressing the area A of the PWN in units of the beam at 1.5 GHz (Table 1) and taking F_1.5 ~ F_1.4, where F_1.5 is the flux density at 1.5 GHz. Under all the assumptions made, the radio PWN remains undetectable as its surface brightness falls below the sensitivity of the radio image (rms ~ 0.7 mJy beam⁻¹). For an evenlarger PWN, the radio emission is expected to be dimmer. On contrary, if the radio PWN is more compact, we expect a brighter radio nebula. For example, taking a radius of ~ 0.05° (= 5 pc), we get S_1.5 ~ 2.5 mJy beam⁻¹, which is three times the sensitivity of our data. In this case, the non-detection could be explained by the shift of the particle injection of the pulsar to higher energies, which translates into ϵ_R ~ 10⁻⁵ (Gaensler et al. 2000). The non-detection of the radio PWN around PSR J1856+0245 could therefore be a consequence of its large extent and/or an inefficient conversion of the pulsar power into radio emission.

The pulsar PSR J1857+0300, with a DM distanceof ~6.7 kpc, also appears projected over the SB. Its age (~4.6 Myr) is similar to the dynamical age of the SB (~3.3 Myr). We cannot therefore discard that the SN explosion that gave birth to PSR J1857+0300 has occurred within the SB at its early formation. Due to the kick velocity acquired after the SN event (which may be some hundreds of km s⁻¹), the pulsar might have overcome the SB and now is observed distant from it. The contribution of PSR J1857+0300 to the γ-ray emission from HESS J1857+026 is certainly negligible as its spin-down power (~2.5 × 10³² erg s⁻¹) is three orders of magnitude less than the TeV luminosity.

We have shown that there is molecular material projected into the HI cavity. This material seems to form an incomplete shell around the peak of the TeV emission. We can speculate two origins for this material. Either it could be preexisting molecular material that was reached by the expanding SB and was not swept away or evaporated by it, or it could be material from the parent molecular cloud that was swept away by the stellar wind of the central stars and now is detected as an expanding shell of molecular material. As Parizot et al. (2004) point out, the interaction of the stellar winds inside SBs with high-density clumps can amplify turbulence and MHD waves. In addition, dense gas could be acting as a target for hadronic interactions and subsequent production of γ-rays. The explosion that gave birth to PSR J1856+0245 has left a SNR which remains undetectable in the radio band, despite the improved quality of the new VLA observations. This remnant is expected to be evolving inside the SB and could be an efficient accelerator of cosmic rays. In a hadronic scenario, the integral flux of γ-ray photons with energies greater than E is given by $F_{\gamma }(>E)={10}^{-10}{f}_{\Gamma }{E}_{\text{TeV}}^{-\Gamma +1}A$ (Torres et al. 2003), where the factor f_Γ depends on the spectral index Γ, E_TeV = E∕TeV, and $A=\theta E_{51}{D}_{\text{kpc}}^2{n}_0\text{ph}{\text{cm}}^{-2}{\text{s}}^{-1}$, where θ is the efficiency of conversion of the SN energy into cosmic ray energy, E₅₁ is the mechanical energy of the SN explosion in units of 10⁵¹ erg, D_kpc is the distance in kiloparsecs, and n₀ the ambient density of the ISM in cm⁻³. Spectral data from HESS H1857+026 were taken from H. E. S. S. Collaboration (2018): Γ ~ 2.57 (for which, we took f_Γ = 0.19) and F_γ(>1 TeV) ~ 4 × 10⁻¹² ph cm⁻² s⁻¹. For a canonical SN explosion (E₅₁ = 1) at 5.5 kpc and θ = 0.3, an ambient density ~ 22 cm⁻³ is necessary to produce the γ-ray flux from HESS J1857+026. Therefore, based on the density derived from the ¹³CO emission (Table 3), all the clouds are dense enough to act as a target for hadronic emission associated with the SNR. If this mechanism contributes significantly to the TeV emission from HESS J1857+026, we would expect a coincidence between the distribution of the γ-ray emission and the molecular gas. Figure 5 shows that the ¹³CO overlaps with the γ-ray source, but we do not detect an enhancement of the TeV emission at the position of the molecular clouds. The hadronic mechanism associated with the putative SNR could therefore be contributing to the TeV emission from HESS J1857+026 but further evidence is necessary to prove this scenario.

Finally, the SB appears surrounded by the HII region G036.459−00.183 to the northand a complex of HII regions to the south (see Figs. 1 and 2), suggesting that the expanding bubble may have triggered the formation of a second generation of stars. There is evidence of star forming activity in the shells of SBs (Suad et al. 2016, 2012). This is not straightforward for the SB associated with HESS J1857+026. On the one hand, G036.459−00.183 is located at a distance of ~8.9 kpc, which is not compatible with the distance of the SB. On the other hand, the HII complex is probably located at ~ 10.4 kpc based on the distance estimate of G035.649−00.053 obtained from the association with molecular gas at ~52 km s⁻¹ (Anderson & Bania 2009). However, we note that the HII region G035.467+00.139 is not related with the complex. There is associated molecular gas at a velocity of ~ 77.2 km s⁻¹, which places this HII region at ~4.6 kpc, compatible with the distance of the SB. A detailed study of the star forming activity is required to determine whether or not the expanding HI shell is triggering the formation of new stars.

Fig. 5 13CO emission integrated between 78 and 90 km s−1 with TeV emission from HESS J1857+026 (red contours). The color scale is expressed in arbitrary units. The cyan ellipse indicates the inner boundary of the HI cavity, and the yellow cross the positions of the pulsars PSR J1856+0245 and PSR J1857+0300. The color scale is expressed in arbitrary units.

5 Conclusions

Following the discovery of a SB in the direction of HESS J1857+026, we analysed an alternative scenario to the two-source model proposed by MAGIC Collaboration (2014) to explain the very-high-energy emission.

The positional and morphological coincidence between the SB and the TeV emission from HESS J1857+026 suggests an association between them and argues in favour of a single γ-ray source rather than the separation into two different sources at different distances. SBs are created by the collective effects of stellar winds and SNRs and are considered cosmic-ray factories and sources of γ-ray emission in the GeV and TeV bands. We argue that the pulsar PSR J1857+0245 is associated with the SB and could also be contributing to the γ-ray emission. Indeed, this pulsar is energetic enough to power the entire TeV emission from HESS J1857+026. The SNR associated with this pulsar remains undetected, but due to its young age it is probably still expanding and accelerating particles inside the bubble. Part of the TeV emission from HESS J1857+026 could be of hadronic origin. Indeed, we detect dense molecular material probably associated with the SB that could be acting as a target for hadronic interactions of protons accelerated inside the bubble.

On the contrary, the old pulsar PSR J1857+0300 could also be associated with the SB, although its contribution to the γ-ray emission is negligible due to its low energetics. In addition, the associated SNR has already dissipated into the ambient gas, though it might have contributed in the past to the injection of energy within the SB.

We conclude that the γ-ray emission from HESS J1857+026 is produced in the interior of a SB. Very-high-energy photons are probably radiated by a population of relativistic particles accelerated in the stellar winds and/or SNRs within the expanding bubble. Furthermore, particles could also be accelerated in a PWN powered by PSR J1857+0245, which is located inside the SB. This study provides evidence that SBs are promising candidates for relevant sources of Galactic cosmic rays.

Acknowledgments

A.P. and E.G. are members of the Carrera del Investigador Científico of CONICET, Argentina. L.D. is postdoctoral fellow of CONICET, Argentina. This work was partially supported by Argentina grants awarded by UBA (UBACyT) and ANPCYT. The National Radio Astronomy Observatory is a facility of the National Science Foundation operated under cooperative agreement by Associated Universities, Inc. 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, and AST-0507657. This publication makes use of data from FUGIN, FOREST Unbiased Galactic plane Imaging survey with the Nobeyama 45-m telescope, a legacy project in the Nobeyama 45-m radio telescope.

References

All Tables

All Figures

Fig. 1 Radio continuum emission at 1.5 GHz. The white crosses indicate the positions of the pulsars PSR J1856+0245 and PSR J1857+0300. The white square shows the region mapped in the C band with the VLA and the contours represent the TeV emission below 1 TeV extracted from H. E. S. S. Collaboration (2018) with values of 1.5, 2.0, 2.7, 3.5, and 4.5 in units of 10−13 ph cm−2 s−1. The angular resolution is 17.′′38 × 15.′′ 51, PA = 6°.87 and the noise is 0.7 mJy beam−1. The color scale is expressed in Jy beam−1.

In the text

	Fig. 2 Magnification of the southern region of the image obtained at 1.5 GHz. The circles represent the positions of the HII regions taken from the WISE Catalog of Galactic HII Region. The synthesized beam is represented in the upper left corner. The color scale is expressed in Jy beam−1.
In the text

	Fig. 3 Radio continuum emission at 6.0 GHz. The resolution is 3.′′55 × 2.′′ 62 while the effective noise level is 12 μJy beam−1. The color scale is expressed in Jy beam−1.
In the text

	Fig. 4 HI distribution towards HESS J1856+026, integrated inthe velocity range between 81 and 102 km s−1. The contours correspond to the TeV emission and the crosses mark the position of the pulsars PSR J1856+0245 and PSR J1857+0300. The color scale is expressed in arbitrary units.
In the text

	Fig. 5 13CO emission integrated between 78 and 90 km s−1 with TeV emission from HESS J1857+026 (red contours). The color scale is expressed in arbitrary units. The cyan ellipse indicates the inner boundary of the HI cavity, and the yellow cross the positions of the pulsars PSR J1856+0245 and PSR J1857+0300. The color scale is expressed in arbitrary units.
In the text