GeV Gamma-Ray Emission and Molecular Clouds toward Supernova Remnant G35.6–0.4 and the TeV Source HESS J1858+020

We note that there are two 4FGL-DR2 catalog sources 4FGL J1857.6+0212 and 4FGL J1858.3+0209 toward the SNR-HIIR complex G35.6−0.4 region. They are first treated as backgrounds to check the residual emission around our target. We use the fit method to refit the spectral parameters of the sources within 3° from the ROI center with a significance above 4σ and the normalization parameters of the two diffuse background components. Then, we fix all parameters except for the normalization of the Galactic diffuse emission to their best-fit values and generate the residual test-statistic (TS) map which is shown in the left panel of Figure 2. Here, the TS value is defined as TS = 2log(L₁/L₀), in which L₀ is the maximum likelihood of the null hypothesis and L₁ is the maximum likelihood with a putative source located in this pixel. As can be seen in Figure 2(a), there is still strong residual emission to the west of the SNR-HIIR complex. We thus add three point sources with simple power-law spectra located at the peak pixels with TS > 20, the positions of which are listed in Table 1, to model the residuals and repeat the above procedures. After this, we find that there is still some excess to the south of the complex. We then try adding this access as the fourth point source with the power-law spectrum. Compared to the three-point-source model, the likelihood of the four-point-source model can be increased by about 15, resulting in TS_model = 2log (L_4ps/L_3ps) ≈ 30. So we use the four point sources to model the residual emission. Finally, the background-subtracted TS map is displayed in Figure 2(b).

Zoom In Zoom Out Reset image size

Next, we study the γ-ray morphology toward SNR-HIIR complex G35.6−0.4 in detail. After excluding the 4FGL-DR2 catalog sources 4FGL J1857.6+0212 and 4FGL J1858.3+0209 from the source model, the TS maps in 0.2–500 GeV and 8–500 GeV are shown in Figures 3(a) and (b), respectively. As can be seen, there is little emission above 8 GeV in the SNR region (represented by the bigger radio contours), while the emission at such energies concentrates on the H ii region (represented by the smaller radio contours), suggesting that there are likely two sources. Meanwhile, we use the likelihood ratio to test three spatial models: one point source, one uniform disk, and two point sources. The model with the largest likelihood value will be preferred. To do so, we first use one point source with a LogParabola (LogP) spectrum, and then explore its best-fit position and extension to obtain the likelihood values of L_1ps for the one-point-source model and L_disk for the uniform disk model. For the two-point-source model, we use the templates of 4FGL J1857.6+0212 and 4FGL J1858.3+0209 and refit their best-fit positions, resulting in SrcA (R.A. = 284418, decl. = 2192, 1σ error = 0016) and SrcB (R.A. = 284582, decl. = 2142, 1σ error = 0042) for 4FGL J1857.6+0212 and 4FGL J1858.3+0209, respectively. The extension of the two sources is also explored, giving TS_ext = 2.4 and 4.0 for SrcA and SrcB, respectively. Finally, according to the likelihood values for the three spatial models, we obtain ${\mathrm{TS}}_{\mathrm{disk}}=2\mathrm{log}({L}_{\mathrm{disk}}/{L}_{1\mathrm{ps}})=13$ and TS_2ps = 30 for the uniform disk and two-point-source models, respectively. Considering the TS distribution in the different energy ranges and the likelihood ratio test for the three spatial models, the two-point-source model for the SNR-HIIR complex G35.6−0.4 is preferred and used in the following analysis. The TS values of SrcA and SrcB in the two-point-source model are fitted to be 761.2 and 55.3, corresponding to the significance of 27.1σ and 6.7σ, respectively.

Zoom In Zoom Out Reset image size

In Figure 3, the TS maps for SrcA and SrcB are also presented in the second-row panels. Our best-fit positions with 1σ uncertainty are displayed with the red pluses and circles, respectively. For comparison, the 4FGL-DR2 positions and the best-fit ones in Cui et al. (2021) are shown with green pluses and green diamonds, respectively. As can be seen, the new best-fit positions for both sources are very close to their original ones listed in the catalog, but are different from those in Cui et al. (2021), in particular for SrcA. To check this difference, we fit the positions by just using 5–500 GeV data and obtain similar results to those of our above treatment in 0.2–500 GeV, ruling out the effect of different energy ranges used in the two studies. Further, we repeat the analysis in Cui et al. (2021) by replacing the IRFs "P8R3_SOURCE_V3_v1" with "P8R3_SOURCE_V2_v1" and using 8 yr catalog for the background, and reproduce the results in Cui et al. (2021). Thus, the discrepancy in the best-fit position for SrcA is verified to be caused by using different IRFs and catalogs.

In the 4FGL-DR2 catalog, the spectral types of SrcA and SrcB are LogP and PowerLaw (PL), respectively. To study the spectral properties of both sources in the whole energy range of 0.2–500 GeV, other spectral types including ExpCutoffPowerLaw (ECPL) and BrokenPowerLaw (BPL) will also be explored. We first change the spectral type of SrcA to PL, ECPL, and BPL, respectively, to find the best spectral type. At the same time, the spectral type of SrcB remains to be PL. After this, we keep the best choice for SrcA and change the spectral type of SrcB to LogP, ECPL, and BPL, respectively, to find the best spectral formula for SrcB. The formulae of these spectra are listed in Table 2. The spectral type is favored if it has the largest TS value defined as ${\mathrm{TS}}_{\mathrm{model}}=-2\mathrm{log}({L}_{\mathrm{PL}}/{L}_{\mathrm{model}})$. As shown in Table 3, a LogP spectrum is preferred for SrcA. But for SrcB, there is no obvious difference between the four spectral types and the simplest PL form is adopted. Finally, we obtain the fluxes in 0.2–500 GeV as 4.2 × 10⁻¹¹ erg cm⁻² s⁻¹ (with Γ = 2.66 ± 0.07 and β = 0.44 ± 0.06 for E₀ = 1.57 GeV) and 1.1 × 10⁻¹¹ erg cm⁻² s⁻¹ (with Γ = 2.19 ± 0.09) for SrcA and SrcB, respectively.

Table 2. Formulae for γ-Ray Spectra

Name	Formula	Free Parameters
PL	${dN}/{dE}={N}_{0}{\left(E/{E}_{0}\right)}^{-{\rm{\Gamma }}}$	N0, Γ
ECPL	${dN}/{dE}={N}_{0}{\left(E/{E}_{0}\right)}^{-{\rm{\Gamma }}}\exp (-E/{E}_{\mathrm{cut}})$	N0, Γ, Ecut
LogP	${dN}/{dE}={N}_{0}{\left(E/{E}_{0}\right)}^{-{\rm{\Gamma }}-\beta \mathrm{log}(E/{E}_{0})}$	N0, Γ, β
BPL	${dN}/{dE}={N}_{0}\left\{\begin{array}{l}{\left(E/{E}_{{\rm{b}}}\right)}^{-{{\rm{\Gamma }}}_{1}},E\leqslant {E}_{{\rm{b}}}\\ {\left(E/{E}_{{\rm{b}}}\right)}^{-{{\rm{\Gamma }}}_{2}},E\geqslant {E}_{{\rm{b}}}\end{array}\right.$	N0, Eb, Γ1, Γ2

Download table as:  ASCIITypeset image

The spectral energy distributions (SEDs) within 0.2–500 GeV of SrcA and SrcB are generated by using the maximum likelihood analysis in seven logarithmically spaced energy bins. During the fitting process, the free parameters only include the normalization parameters of the sources with the significance 4σ within 5° from the ROI center as well as the Galactic and isotropic diffuse background components, while all the other parameters are fixed to their best-fit values from the above analysis in the whole energy (0.2–500 GeV) range. In the energy bins where the TS value of SrcA or SrcB is smaller than 4, we calculate the 95%-confidence level upper limit of its flux. The results are displayed in Figure 9.

As shown in Figure 3(c), SrcA is located within the shell of SNR G35.6−0.4 and is projectively close to two pulsars, PSR J1857+0212 and PSR J1857+0210. We first discuss the possibility of the pulsar's origin. According to the dispersion measures, PSR J1857+0212 and PSR J1857+0210are at distances 7.98 kpc (Han et al. 2006) and 15.4 kpc (Morris et al. 2002), respectively. But the distances are revised to 6.0 and 7.3 kpc in the ATNF pulsar catalog (Manchester et al. 2005), respectively, based on the new electron-density model developed by Yao et al. (2017). Taking 6.0 kpc (the smaller one) for example, SrcA would have a luminosity of $\sim 2\times {10}^{35}{d}_{6\,\mathrm{kpc}}^{-2}\ \mathrm{erg}\ {{\rm{s}}}^{-1}$ in the energy range 0.2–500 GeV, where d_6kpc is the distance in the unit of 6 kpc, which is significantly higher than the spin-down luminosity of 2.2 × 10³⁴ erg s⁻¹ for PSR J1857+0212 (Manchester et al. 2005; Han et al. 2006) and 2.2 × 10³³ erg s⁻¹ for PSR J1857+0212 (Morris et al. 2002; Manchester et al. 2005). Thus the pulsar's origin for SrcA can be ruled out unless the true distance is smaller than 2 kpc for PSR J1857+0212 and 0.7 kpc for PSR J1857+0210.

SNR G35.6−0.4 may be responsible for the γ-ray emission of SrcA, and we provide estimates for the emission in cases of hadronic and leptonic interaction separately. We simply assume that the particles accelerated in the SNR have a power-law form with a high-energy cutoff:

$$\begin{eqnarray}&&{{dN}}_{i}/{dE}={A}_{i}{\left(E/{E}_{0}\right)}^{-{\alpha }_{i}}\exp (-E/{E}_{i,c})\end{eqnarray} \tag{ 1 }$$

where, i = e, p, α_i , and E_i,c are the power-law index and high-energy cutoff, respectively. The normalization A_i is determined by the total energy (W_i ) in particles with energy above 1 GeV. Since the SNR is very likely to be associated with the ∼+50 km s⁻¹ MC at a distance of ∼10.5 kpc as revealed above, the γ-ray emission arising from the SNR-accelerated particles' hadronic interaction is first considered. In such a hadronic scenario, SrcA's spectrum can be well fit (as shown in the left panel of Figure 7) and we obtain α_p ≈ 3.0 and n₀ W_p ≈ 8 × 10⁵¹ erg cm⁻³, where n₀ is the number density of the target gas for proton–proton hadronic interaction. The cutoff energy cannot be constrained by the current data and is fixed as 3 PeV in our calculation. With a mass of 9 × 10³ M_⊙ and an average atomic hydrogen density n₀ ∼ 140 cm⁻³ for the ∼+50 km s⁻¹ filamentary molecular gas, the energy budget in protons W_p is about 6 × 10⁴⁹ erg, which is acceptable and reasonable for the SNR scenario. The somewhat large index may be attributed to the strong ion-neutral collisions (Malkov et al. 2011) or the escaped process (e.g., Aharonian & Atoyan 1996; Li & Chen 2012). In addition, the luminosity of SrcA in 1–100 GeV is 2 × 10³⁶ erg s⁻¹, which is consistent with the known γ-ray-bright interacting SNRs (Liu et al. 2015; Acero et al. 2016). Thus the hadronic scenario in which the energetic protons are from SNR G35.6−0.4 is a plausible explanation for SrcA.

For the leptonic case in which γ-rays are produced via the inverse Compton (IC) process, due to the lack of constraint on the index by the GeV data, we fix α_e = 2.0 based on the radio index α_r = 0.47 (Green 2009). To explain the data (see the model curve in orange in the left panel of Figure 9), E_e,c ≈ 80 GeV and W_e ≈ 7 × 10⁵⁰ erg are obtained when the seed photons include the IR emission with a temperature of 35 K and an energy density of 0.6 eV cm⁻³ estimated from the interstellar radiation field (ISRF) model (Porter et al. 2006; Shibata et al. 2011) and the cosmic microwave background. Considering the large amount of energy in electrons, which is comparable to the canonical supernova explosion energy, the leptonic process seems not proper to explain the γ-ray emission of SrcA unless there is an unusually strong IR emission.

Zoom In Zoom Out Reset image size

As shown in Figure 3(d), the GeV source SrcB is spatially coincident with the unidentified TeV source HESS J1858+020. Moreover, the GeV spectrum of SrcB can be well connected with the TeV spectrum of HESS J1858+020 by a simple power-law function with an index of ∼2.2 (see the right panel in Figure 9). Thus, SrcB is very likely the GeV counterpart of TeV source HESS J1858+020. Meanwhile, both sources are spatially coincident with the MCs at a velocity around the ∼+53 to +57 km s⁻¹ which have a mass of 1.3 × 10³ M_⊙ and an average atomic hydrogen number density of 600 cm⁻³. As suggested in Paron & Giacani (2010) and this work (Section 3.2), this molecular arc is likely associated with HIIR G35.6−0.5. Both the association with MCs and the soft GeV–TeV spectra (index >2) suggest the GeV–TeV γ-ray emissions likely have a hadronic origin. One possible scenario for SrcB and HESS J1858+020 is that the energetic protons accelerated in HIIR G35.6−0.5 bombard the MC to produce the hadronic γ-ray emission. By adopting a distance of 3.4 kpc, α_p ≈ 2.15 and W_p ≈ 1.7 × 10⁴⁷ erg are obtained. To explain the TeV data, which seem not to have an obvious cutoff, the cutoff energy of protons E_p,c should reach the PeV energy (see the right panel of Figure 9). This implies that there may be a potential PeV accelerator in HIIR G35.6−0.5. This can be tested by the LHAASO experiment in the future. As one of the PeVatron candidates, the lower limit of the proton cutoff energy ∼30 TeV for HESS J1858+020 is obtained by only using the TeV data (Spengler 2020). In fact, the index and the cutoff energy are strongly degenerated. We note that the index constrained in Spengler (2020) for this source is obviously smaller than 2.0. With the help of the GeV data, the index is fitted as ∼2.2 in this study, resulting in a larger cutoff energy.

HIIRs vary from ∼0.1 pc to a hundred pc in size, are generally related to the active star formations, and are considered to be particle accelerators. It has been proposed that particles can be accelerated by the colliding winds of binaries (e.g., Eichler & Usov 1993) and/or the young massive clusters (Aharonian et al. 2019), which are the ionizing sources of HIIRs. This is supported by the TeV γ-ray observations toward η Car (H.E.S.S. Collaboration et al. 2020) for the colliding winds and the Cygnus region (Aliu et al. 2014; Bartoli et al. 2014; Abeysekara et al. 2021; Amenomori et al. 2021) for the massive star clusters. In HIIR G35.6−0.5, although there are no massive stars detected as yet, Paron & Giacani (2010) found six young stellar object (YSO) candidates toward the molecular clump at ∼+53 km s⁻¹ (namely the southern clump in their paper) and concluded that there is an active star formation region. Further observations found at least one evolved YSO (i.e., IRS1) is embedded in the clump, although no signature of outflows was confirmed (Paron et al. 2011). By analyzing the Chandra data, Paredes et al. (2014) found seven X-ray sources, and suggested that sources X1–4, almost projectively distributed on the shell of the HIIR (see the right panel of Figure 1), might be embedded protostars and that source X5, close to the "southern" (in the Galactic coordinate system) clump, might be coincident with the star formation region. Thus, it could not be excluded that there may be some nondetected massive stars embedded in the MC.

Alternatively, we also consider the lepto-hadronic hybrid case in which the GeV emissions are from the hadronic process and the TeV γ-rays are generated by the IC process. We assume the electrons and protons have the same power-law index and high-energy cutoff. For the seed photons in the IC process, the IR component of ISRF is also included and is estimated as 35 K and 0.6 eV cm⁻³ by using the similar method for SrcA above. To explain the data (see the purple curve in the right panel of Figure 9), α_e = α_p ≈ 2.3 and E_c,e = E_c,p ≈ 200 TeV are obtained. Electrons will suffer the synchrotron radiation loss during the acceleration, giving a cooling-limited maximum energy ${E}_{\max ,\mathrm{cool}}=35{u}_{3}/\sqrt{{\eta }_{g}{B}_{1}}$ TeV (e.g., Zirakashvili & Aharonian 2007; Ohira et al. 2012), where u₃ is the shock velocity in unit of 1000 km s⁻¹, η_g the gyrofactor, and B₁ the magnetic field strength in unit of 10 μG. For the wind velocity u₃ = 1 and the Bohm limit η_g = 1, it requires B < 0.3 μG to boost the energy of electrons up to 200 TeV. This magnetic strength is as weak as an order of magnitude lower than the mean value for Galactic ISM, which means that it is hard to accelerate electrons to the 100 TeV band for the standard interstellar magnetic field. Thus, the hybrid model seems rather unlikely according to the fitted results.