A Study of GeV Gamma-ray Emission toward Supernova Remnant G51.26+0.11 and Its Molecular Environment

TL;DR

This study reanalyzes GeV gamma-ray emission from SNR G51.26+0.11, revealing a likely SNR-MC interaction and explaining the gamma-ray emission with a hadronic model based on 13.2 years of Fermi-LAT data.

Contribution

It provides the first detailed analysis of gamma-ray emission and molecular environment interaction for SNR G51.26+0.11, supporting a hadronic origin of the gamma rays.

Findings

Extended gamma-ray emission overlaps with SNR G51.26+0.11

Evidence of SNR-MC interaction through molecular line profiles

Gamma-ray emission explained by a hadronic model with proton spectral index ~2.25

Abstract

We reanalyze the Fermi-LAT GeV $\gamma$-ray emission in the region of supernova remnant (SNR) G51.26+0.11 and investigate its interstellar molecular environment with the CO-line data. At GeV energies, based on 13.2 years of Fermi-LAT data, the extended $\gamma$-ray emission observed in this region is resolved into a uniform-disk source ('Src A') with a significance of 19.5$\sigma$ and a point source (4FGL J1924.3+1628) with a significance of 4.2$\sigma$ in 0.2$-$500 GeV. With an angular radius of $\sim$ 0.17$^{\deg}$, 'Src A' overlaps with SNR G51.26+0.11 significantly in the line of sight. On the other hand, the morphological coincidence between the SNR and the $\sim$ +54 km s$^{-1}$ molecular clouds (MCs) together with the asymmetric or broad $^{12}$CO line profiles near the SNR boundary provides evidence for the very likely SNR-MC interaction. The SNR-MC interaction and the HI absorption features indicate that SNR G51.26+0.11 is located at a kinematic distance of 6.2 $\pm$ 0.5 kpc. Combined with the results from the multi-wavelength analysis, the $\gamma$-ray emission of the SNR ('Src A') can be naturally explained by a hadronic model with a soft power-law proton spectrum of index $\sim$ 2.25.

Tables

Table 1. Table 1: Best-fitting parameters in different spatial models ( > > 5 GeV).

Morphological model		R.A.(J2000) a	Dec(J2000) a	Extension ($\sigma$) b	$\Delta k$ c	$\Delta$AIC c
		(∘)	(∘)	(∘)
Single-source	1Gauss d	291.31$\pm 0.02$	16.30$\pm 0.02$	0.22${}_{-0.03}{}^{}{}_{}{}^{+0.04}$	0	0
Two-source	1Gauss	291.35$\pm 0.02$	16.28$\pm 0.02$	0.18${}_{-0.03}{}^{}{}_{}{}^{+0.04}$	2	-7.87
	1PS e	291.10	16.48	$-$
	1Disk	291.38$\pm 0.02$	16.27$\pm 0.02$	0.17${}_{-0.02}{}^{}{}_{}{}^{+0.02}$	2	-12.49
	1PS e	291.10	16.48	$-$

Table 2. Table 2: Formulae and free parameters for different γ 𝛾 \gamma -ray spectra models.

Name	Formula	Free parameters
PL	$dN/dE=N_0{\left(E/E_0\right)}^{-\Gamma }$	$N_0$, $\Gamma$
PLEC	$dN/dE=N_0{\left(E/E_0\right)}^{-\Gamma }\exp \left(-E/E_{\mathrm{cut}}\right)$	$N_0$, $\Gamma$, $E_{\mathrm{cut}}$
LogP	$dN/dE=N_0{\left(E/E_0\right)}^{-\Gamma -\beta \log \left(E/E_0\right)}$	$N_0$, $\Gamma$, $\beta$
BPL	$dN/dE=\{\begin{array}{cc}{N}_0{\left(E/E_0\right)}^{-{\Gamma }_1},E\le E_b& \\ N_0{\left(E/E_0\right)}^{-{\Gamma }_2},E>E_b& \end{array}$	$N_0$, ${\Gamma }_1$, ${\Gamma }_2$, $E_b$

Table 3. Table 3: Results from spectral analysis of ‘Src A’.

Spectral model	$\Delta k$ a	$\Delta$AIC a
PL	0	0
PLEC	1	20.54
LogP	1	2.00
BPL	2	2.51

Table 4. Table 4: Fitted parameters for region ‘E1’ and ‘E2’.

Region	Line	Centre	FWHM	$T_{\mathrm{peak}}$	$W$
		(km s-1)	(km s-1)	(K)	(K km s-1)
E1	12CO ($J=1-0$) a	52.8	3.0	6.5	66.4
		55.4	10.0	4.2
	13CO ($J=1-0$)	54.1	4.0	3.8	17.2
E2	12CO ($J=1-0$) a	53.9	3.3	2.8	39.5
		53.9	12.9	2.2
	13CO ($J=1-0$)	54.1	2.8	1.6	6.0

Table 5. Table 5: Properties of the molecular clouds ‘Cloud A’ and ‘Cloud B’.

MC	(R.A., Dec) a	$R$	$d$	$N$(H2)	$M$	$n_H$	$T_{\mathrm{ex}}^{}{}_{}{}^b$	$\tau$(13CO) c
	(∘)	(∘)	(kpc)	(1021 cm-2)	(10${}_{}{}^{4}d_{}^2{M}_{\odot }$)	(cm-3)	(K)
Cloud A	(291.38, 16.27)	0.17	6.2	2.8	6.4	70	13.4	0.31
Cloud B	(291.59, 16.38)	0.08	6.2	8.0	4.6	430	24.5	0.86

Equations

$$\mathrm{d}{N}_i/\mathrm{d}{E}_i\propto E^{-{\alpha }_i}\exp (-E_i/E_{i,\mathrm{cut}}),$$

Full text

††pubyear: 2023††pagerange: A Study of GeV Gamma-ray Emission toward Supernova Remnant G51.26+0.11 and Its Molecular Environment–A Study of GeV Gamma-ray Emission toward Supernova Remnant G51.26+0.11 and Its Molecular Environment

1 Introduction

The acceleration of cosmic rays (CRs) can be explored through observations of the $\gamma$-ray emission resulting from the decay of $\pi^{0}$ mesons produced via proton-proton collisions (Ackermann et al., 2013). Supernova remnants (SNRs) are a sort of the popular candidates for Galactic CR accelerators by virtue of their strong shocks, while molecular clouds (MCs) can provide target gas for proton-proton hadronic interaction. Therefore, studies of the $\gamma$-ray emission from SNR-MC associations are crucial to search for the signatures of proton acceleration, and multi-wavelength observations are necessary to reveal their interstellar environments. Dozens of GeV $\gamma$-ray sources associated with SNRs have been discovered with the Large Area Telescope (LAT) on board the Fermi Gamma-ray Space Telescope, among which SNRs interacting with nearby MCs appear to be brighter and softer (Liu et al., 2015; Acero et al., 2016).

G51.21+0.11 was classified as an SNR candidate by Anderson et al. (2017) using the combination of 1–2 GHz continuum data from The HI, OH, Recombination line survey of the Milky Way (THOR) and 1.4 GHz Very Large Array Galactic Plane Survey continuum data and observed by Driessen et al. (2018) using the LOw Frequency ARray. Following studies identified it as a complex of two separate SNRs: the compact SNR G51.04+0.07 (Supan et al., 2018) and the shell-type SNR G51.26+0.11 (Dokara et al., 2018). Evidence of nonthermal emission from the two distinct regions in this complex was also provided (Supan et al., 2018; Dokara et al., 2018). Figure 1 displays the radio emission from the complex in the 80–300 MHz GaLactic and Extragalactic All-sky Murchison Widefield Array survey (GLEAM; Wayth et al., 2015; Hurley-Walker et al., 2017) data.

The SNR nature of G51.26+0.11 was further confirmed by Dokara et al. (2021) using the data of the GLObal view of STAR formation in the Milky Way survey that was conducted with the Karl G. Jansky VLA. They measured a degree of polarization of $0.06\pm 0.02$ and a radio spectral index of $\sim-0.4$, which are consistent with the case of SNRs. The measured 200 MHz and 1.4 GHz flux densities are $25.8\pm 3.6$ Jy and $12.4\pm 0.6$ Jy, respectively. Recently, Ranasinghe & Leahy (2022) estimated the distance of $6.6\pm 1.7$ kpc to G51.26+0.11 using the H I absorption spectra constructed with the THOR data. At GeV energies, Araya (2021) performed a study of the $\gamma$-ray emission with 12.3 years of Fermi-LAT data and found an extended $\gamma$-ray source with a significance of 21$\sigma$ above 200 MeV in the region of G51.26+0.11. They suggested that the spectral energy distribution (SED) observed could be fitted by hadronic or inverse Compton emission. In the hadronic scenario, the spectral index of the protons is $\sim$ 2.18 and G51.26+0.11 could be the cause of the $\gamma$-rays. The MCs at the local-standard-of-rest (LSR) velocities +43.4 – +52.2 and +53 – +58.8 km s*-1* were suggested to be compatible with the hadronic scenario. Meanwhile, the pulsar wind nebula (PWN) origin could not be ruled out or confirmed.

In this study, we reanalyze the GeV emission from SNR G51.26+0.11 using Fermi-LAT data and investigate the interstellar molecular environment of the SNR with millimeter CO-line observations. We focus on the morphology of the $\gamma$-ray emission and the possible hadronic contribution resulting from the SNR-MC interaction. In Section 2, we describe the $\gamma$-ray and multi-wavelength observations and data reduction. In Section 3, the corresponding analysis and results are presented. The results are discussed in Section 4 and summarized in Section 5.

2 Observations and Data

2.1 Fermi-LAT Observational Data

The Large Area Telescope (LAT) onboard the Fermi Gamma-ray Space Telescope surveys the $\gamma$-ray sky in the 20 MeV to more than 300 GeV energy range. We collect more than 13.2 years of Fermi-LAT Pass 8 data toward a circular region which is $15^{\circ}$ in radius and centered at the coordinates R.A. = $291.3^{\circ}$, Dec = $16.3^{\circ}$ (J2000), from 2008-08-04 15:43:36 (UTC) to 2021-11-07 02:26:50 (UTC). We analyze the data with the standard softwares Fermitools111http://fermi.gsfc.nasa.gov/ssc/data/analysis/software/ version 2.0.8 released on 2021 January 20, and Fermipy222https://fermipy.readthedocs.io/en/stable/ version 1.0.1 released on 2021 March 12. We select ‘SOURCE’ class and ‘FRONT+BACK’ type events (evclass = 128, evtype = 3) with zenith angle $<90^{\circ}$ to eliminate Earth limb events, and restrict the energy range to 0.2–500 GeV. To choose good time intervals, we apply the recommended filter string ‘(DATA_QUAL==1)&&(LAT_CONFIG==1)’ in tool gtmktime. For ‘SOURCE’ events, the instrument response function (IRF) ‘P8R3_SOURCE_V3_v1’ is used. To build the background model, we include the LAT sources listed in the Fermi-LAT Fourth Source Catalog Data Release 2 (4FGL-DR2, Ballet et al. 2020) in a radius of $25^{\circ}$ around the centre of the region of interest (ROI), as well as the Galactic diffuse emission (gll_iem_v07.fits) and isotropic emission (iso_P8R3_SOURCE_V3_v1.txt).

2.2 CO Line Data

We make use of the data from the Milky Way Imaging Scroll Painting (MWISP) survey project. The CO observation was made in the 12CO ($J=1$ – [math]) line (at 115.271 GHz), and the 13CO ($J=1$ – [math]) line (at 110.201 GHz) in 2016 June using the 13.7 m millimeter-wavelength telescope of the Purple Mountain Observatory at Delingha, China. We focus on a $1.5^{\circ}\times 1.5^{\circ}$ area covering G51.26+0.11 centred at $l=51.26^{\circ}$, $b=+0.11^{\circ}$ in the Galactic coordinate system with a grid spacing $\sim 30^{\prime\prime}$. The half-power beam width of the telescope is about 50*′′* for the two lines. The mean RMS noise level of the main beam brightness temperature is $\sim$ 0.52 K for the 12CO ($J=1$ – [math]) at the velocity resolution of 0.158 km s*-1* and $\sim$ 0.25 K for the 13CO ($J=1$ – [math]) lines at the velocity resolution of 0.166 km s*-1*.

2.3 Other Data

As an aid to constraining the distance to SNR G51.26+0.11, we obtain the 1420 MHz radio continuum and H I line emission data from the Very Large Array Galactic Plane Survey (Stil et al., 2006). The continuum image has a spatial resolution of $18^{\prime\prime}$. The synthesized beam for the H I spectral line images is $18^{\prime\prime}$ and the radial velocity resolution is 0.824 km s*-1*.

3 Multi-wavelength Analysis

3.1 Fermi-LAT Gamma-ray Emission

3.1.1 Morphological Analysis

There are two 4FGL-DR2 catalogue sources, 4FGL J1925.4+1616 and 4FGL J1924.3+1628, toward SNR G51.26+0.11. To analyze the morphology of the GeV emission in the ROI, we remove these two sources from the source models. To reduce uncertainties caused by large point-spread function at lower energies, we select photon events at energies above 5 GeV in this region. Then, we perform the binned maximum likelihood analysis and apply the energy dispersion correction. We free the normalization parameters of the Galactic and isotropic diffuse background components and the sources within $10^{\circ}$ from the ROI centre, as well as the spectral parameters of the sources within $3^{\circ}$ from the ROI centre while fitting the models. Finally, we generate the residual test-statistic (TS) map of the $1.5^{\circ}\times 1.5^{\circ}$ region centred at G51.26+0.11 in the energy range of 5–500 GeV. The test statistic is defined as ${\rm TS}=-2\log(\mathcal{L}_{0}/\mathcal{L}_{1})$, in which $\mathcal{L}_{0}$ is the likelihood of the null hypothesis and $\mathcal{L}_{1}$ is the likelihood of the hypothesis being tested. The detection significance $\sigma$ is usually approximated by the square root of the TS (Wilks, 1938). As shown in Figure 2(a), the excess of the TS map shows an elongated shape toward the northwest, suggesting that an extended symmetric Gaussian may not be sufficient to fit well.

We define the following three source models to describe the $\gamma$-ray emission, where the notation PS stands for point sources, Gauss for the Gaussian template, and Disk for the disk template with a uniform model.

1. i.

1Gauss includes only an extended source with a Gaussian model assuming a power-law (PL)-type spectrum; 2. ii.

1Gauss+1PS includes an extended source with a Gaussian model and a point source, which are assumed to have PL-type spectra; 3. iii.

1Disk+1PS includes an extended source with a uniform disk model and a point source, which are assumed to have PL-type spectra.

We place the extended source and the point source mentioned above (called target sources) at the same positions as the catalogued sources 4FGL J1925.4+1616 and 4FGL J1924.3+1628, respectively. For the models including the point source (4FGL J1924.3+1628), we fix the location of the point source to the original coordinates of 4FGL J1924.3+1628. For all three models, we use method extension of Fermipy to obtain the best-fitting extension and position of the extended source (4FGL J1925.4+1616). We also calculate the likelihood ratio between the best-fitting extended model and the point source hypothesis, defined as TSext $\equiv 2\log({\mathcal{L}_{\rm ext}}/{\mathcal{L}_{\rm PS}})$. The $\gamma$-ray source is considered to be significantly extended only if its TSext exceeds 16 (Lande et al., 2012). In all three models, the measurement results indicate that the extended nature (of 4FGL J1925.4+1616) is genuine (TSext=48.4, 32.2, and 29.3, respectively). In the 1Gauss model, the best-fitting position and 68%-containment radius are R.A. = $291.31\pm 0.02^{\circ}$, Dec = $16.30\pm 0.02^{\circ}$ (J2000) and $\sigma={0.24^{+0.04}_{-0.03}}^{\circ}$, respectively, which are consistent with those reported by Araya (2021). For comparison, we replace the spatial parameters of the 1Gauss model with Araya’s results in the following analysis. The best-fitting spatial parameters of the three models are summarized in Table 1.

We add each of the three templates to the background model and repeat the fit with Fermitools. During the fitting process, we only free the normalization and spectral parameters of the target sources as well as the normalization of the Galactic diffuse emission, while all the other parameters are fixed to their best-fitting values from the above analysis. We use the Akaike Information Criterion (AIC; Lande et al. 2012) to select the best-fitting spatial model. AIC is defined as ${\rm AIC}=2k-2\ln\mathcal{L}$, where $k$ is the number of free parameters in the model and $\mathcal{L}$ is the likelihood of the model tested in the analysis. The model which minimizes the AIC is considered to be the best. Compared with the 1Gauss model, the $\Delta$AIC values for the three models are tabulated in Table 1. As can be seen, 1Disk+1PS could be the best-fitting spatial model, in which we consider 4FGL J1925.4+1616 as a uniform-disk source and 4FGL J1924.3+1628 as a point source. With the inclusion of the 1Disk+1PS template in the background model, there is almost no residual $\gamma$-ray emission in the TS map (see Figure 2(b)). Therefore, we replace the template of 4FGL J1925.4+1616 with its best-fitting parameters and consider it as the GeV counterpart of G51.26+0.11, denoted as ‘Src A’. As we include source 4FGL J1924.3+1628 in the source model, the angular radius of ‘Src A’ ($\sim 0.17^{\circ}$) is a bit smaller than Araya’s (2021) result. The fitted TS values for ‘Src A’ and 4FGL J1924.3+1628 in 5–500 GeV are 95.6 and 22.2, respectively, with the significance of 9.8$\sigma$ and 3.7$\sigma$. The corresponding TS maps which exclude ‘Src A’ and 4FGL J1924.3+1628 from the adjusted model are shown in Figures 2(c) and 2(d), respectively.

3.1.2 Spectral Analysis

We select photon events in the energy range of 0.2–500 GeV for the spectral analysis. After adjusting the background model with the 1Disk+1PS template, we free the normalization and spectral parameters of the sources within 5∘ from the ROI centre. The normalization parameters of the Galactic and isotropic diffuse background components are also set free. To study the spectral properties of ‘Src A’, we test four spectral models: a power-law (PL), an exponentially cutoff power-law (PLEC), a log-parabola model (LogP), and a broken power-law (BPL). The formulae and free parameters of these spectra are listed in Table 2. Since the lowest $\Delta$AIC value points to the best-fitting spectral model, as shown in Table 3, a PL spectrum is preferred for ‘Src A’, and no significant curvature is detected.

Using the best-fitting PL spectral parameters for ‘Src A’, the obtained flux in 0.2–500 GeV is $2.86\times 10^{-11}\ {\rm erg\ cm^{-2}\ s^{-1}}$ (with ${\rm\Gamma=2.19\pm 0.04}$), giving a luminosity of $7.3\times 10^{34}\ {\rm erg\ s^{-1}}$ in 1–100 GeV at a distance of 6.2 kpc (see §4.2 below for the distance estimate). Since we include 4FGL J1924.3+1628 in the background model, the flux of ‘Src A’ is a bit lower than that reported in Araya (2021). The flux of 4FGL J1924.3+1628 in 0.2–500 GeV, also assumed with a PL spectrum, is $5.14\times 10^{-12}\ {\rm erg\ cm^{-2}s^{-1}}$ (with ${\rm\Gamma=1.70\pm 0.17}$). The overall TS values of ‘Src A’ and 4FGL J1924.3+1628 are 381.7 and 26.0, corresponding to the significance of 19.5$\sigma$ and 4.2$\sigma$, respectively.

The SEDs of the two sources are extracted from the maximum likelihood analysis in eight logarithmically spaced energy bins (see Figure 3). In each bin, we fit the normalization parameters of the diffuse backgrounds and the sources within 5∘ from the ROI centre. Additionally, we consider the systematic uncertainty caused by the model imperfection of the Galactic diffuse background. Using an alternative Galactic diffuse emission model (gll_iem_v06.fits), we repeat the above analysis to regenerate the SEDs and calculate the flux deviation as systematic errors. If the TS value of the source in a bin is less than 4, the 95% upper limit of its flux is estimated.

3.2 CO Line Emission

We investigate the potential molecular gas which may contribute to the hadronic emission. Figure 4 shows the 12CO ($J=1$ – [math]) and 13CO ($J=1$ – [math]) line intensity channel maps over successive $\sim$ 2 km s*-1* intervals in the +52 – +58 km s*-1* velocity range. Diffuse CO gas in the interval +54 – +56 km s*-1* seems to pervade the extent of ‘Src A’, which is roughly coincident with SNR G51.26+0.11; we denote this part of molecular gas as ‘Cloud A’. An extended, curved structure (marked as‘Cloud B’ in the middle panels of Figure 4) in the northeast appears in contact with, and tangent to, the shell of SNR G51.26+0.11, where the radio emission is brightened, which may be indicative of the collision of the SNR shell with the cloud. We thus generate a grid of CO-line spectra toward a $5^{\prime}\times 5^{\prime}$ region at the northeastern boundary of the SNR, marked with white boxes in the left panels of Figure 4. As shown in the grid presented in Figure 5, the wings of the line profiles of the CO emission at the LSR velocity $v_{\rm LSR}$ around $\sim$ +54 km s*-1* appear to be broadened in a number of pixels, as exemplified by regions ‘E1’ and ‘E2’.

The average 12CO line profile of ‘E1’, spanning from +45 to +65 km s*-1*, seems to be asymmetrically broadened from +57 km s*-1* to +65 km s*-1* (see the left panel of Figure 6). The peak velocity of 12CO emission ($\sim$ +53 km s*-1*) is also slightly shifted compared to that of 13CO ($\sim$ +54 km s*-1*). Unlike the 12CO line profile of ‘E1’, the average profile of ‘E2’, with a velocity span similar to that of ‘E1’, appears to be symmetrically broadened (see the right panel of Figure 6). The fitted parameters for the $\sim$ +54 km s*-1* emission lines in regions ‘E1’ and ‘E2’ are summarized in Table 4. Combined with the intensity channel maps, we note that ‘E1’ is located in the curved dense structure ‘Cloud B’, while ‘E2’ appears to be in the extent of ‘Cloud A’, which overlaps the ‘Src A’ disk. Considering there is no significant 13CO emission (which is optically thin and represents the quiescent molecular gas) at the wings of the $\sim$ +54 km s*-1* lines, the broadened 12CO profiles suggest that both ‘Cloud A’ and ‘Cloud B’ are likely perturbed by the shock of SNR G51.26+0.11, and thus provide kinematic evidence for the SNR-MC interaction (e.g., Jiang et al., 2010).

3.3 H I Absorption

We follow the methods of Tian et al. (2007) to build the 21 cm H I absorption spectrum. Denote the average H I brightness temperatures of the on-source and the background regions at velocity $v$ as $T_{\rm on}(v)$ and $T_{\rm off}(v)$, respectively. One then has $T_{\rm on}(v)=T_{B}(v)(1-e^{-\tau_{t}(v)})+T^{c}_{s}(e^{-\tau_{c}(v)}-1)$ and $T_{\rm off}(v)=T_{B}(v)(1-e^{-\tau_{t}(v)}+T^{c}_{\rm bg}(e^{-\tau_{c}(v)}-1)$, where $T^{c}_{s}$ and $T^{c}_{\rm bg}$ are the average continuum brightness temperatures for the on-source and background regions, respectively, $T_{B}(v)$ refers to the spin temperature of the H I cloud, $\tau_{t}(v)$ is the total optical depth along line-of-sight, and $\tau_{c}(v)$ is the optical depth from the continuum source to the observer. Then the absorption spectrum can be expressed as $e^{-\tau_{c}(v)}=1-({T_{\rm off}(v)-T_{\rm on}(v)})/({T^{c}_{s}-T^{c}_{\rm bg}})$. The boxes, both with a size of $3^{\prime}\times 3^{\prime}$, in Figure 1 show the on-source and the background regions we select. The obtained H I emission and absorption spectra are displayed in Figure 7.

4 Discussion

4.1 Molecular Environment

To estimate the molecular column density and the mass of the molecular gas, local thermodynamic equilibrium for the molecular gas and optically thick conditions for the 12CO ($J=1$ – [math]) line are assumed. We first calculate the excitation temperature $T_{\rm ex}={5.53}/{\ln[1+{5.53}/({T_{\rm peak}(^{12}{\rm CO})+0.819})]}$ K (Nagahama et al., 1998), where $T_{\rm peak}(^{12}{\rm CO})$ is the peak temperature of the maximum 12CO($J=1$ – [math]) emission point in the cloud region. The peak excitation temperature of the dense gas in ‘Cloud B’ is $T_{\rm ex}\sim 24$ K, which is noticeably higher than the typical temperature of $\sim$ 10 K in quiescent interstellar MCs. Then we estimate the column density via $N({\rm H}_{2})=1.49\times 10^{20}\ {W(^{13}{\rm CO})}/[{1-\exp(-5.29/T_{\rm ex})}]$ cm*-2* (Nagahama et al., 1998), where $W$(13CO) is the integrated intensity of the 13CO ($J=1$ – [math]) line. We take the abundance ratio [H2]/[13CO] as $7\times 10^{5}$ (Frerking et al., 1982). After estimating the $N({\rm H}_{2})$ value in each pixel of the corresponding regions, we obtain the average values for the clouds. The relation $M=2.8m_{\rm H}N({\rm H}_{2})A$ is used to obtain the mass of the molecular gas, where $m_{\rm H}$ is the mass of a hydrogen atom and $A$ is the cross-sectional area of the corresponding clouds. The clouds are simply approximated as spheres when calculating the average number densities of hydrogen nuclei ($n_{\rm H}$). The parameters of ‘Cloud A’ and ‘Cloud B’ are summarized in Table 5.

4.2 Kinematic Distance

Considering the association with the molecular clouds seen at +43.4 – +52.2 and +53 – +58.8 km s*-1*, SNR G51.26+0.11 was placed at a distance between 3.7 and 6.6 kpc (Araya, 2021). The distance was also estimated to be $6.6\pm 1.7$ kpc from the H I absorption spectra that were constructed with the THOR data (Ranasinghe & Leahy, 2022). Based on the association with the $\sim+54$ km s*-1* MCs, we here can provide an independent estimate for the kinematic distance of the SNR. Applying the Galactic rotation curve together with $R_{0}$ = 8.34 kpc and $V_{0}$ = 240 km s*-1* (Reid et al., 2014), the systemic LSR velocity of +54 km s*-1* corresponds to a near distance $4.2\pm 0.5$ kpc and a far distance $6.2\pm 0.5$ kpc (Wenger et al., 2018). Additionally, the predicted tangent point distance toward SNR G51.26+0.11 is at $\sim$ 5.2 kpc and the tangent point LSR velocity is $\sim$ +52 km s*-1*. As shown in Figure 7, the highest H I emission velocity appears at $\sim$ +62 km s*-1*. Whether the tangent point LSR velocity is $\sim$ +52 km s*-1* or $\sim$ +62 km s*-1*, continuous H I absorption features are present from +45 km s*-1* to the tangent point. This suggests that SNR G51.26+0.11 is located beyond the tangent point, and indicates that the SNR and the associated MCs are most probably at the far kinematic distance, $6.2\pm 0.5$ kpc. The distance we determine here is within the ranges given by the previous studies.

4.3 Origin of the Gamma-ray Emission

4.3.1 Pulsar Origin

We first discuss the possibility that the $\gamma$-ray emission comes from pulsars. As shown in Figure 2, three pulsars, PSR J1926+1613, PSR J1924+1631, and PSR J1924+1628, are seen to be located within the ROI when we search in the Australia Telescope National Facility (ATNF; Manchester et al. 2005)333https://www.atnf.csiro.au/research/pulsar/psrcat/ pulsar catalogue and the SIMBAD Astronomical Database (Wenger et al., 2000)444http://simbad.u-strasbg.fr/simbad/sim-fcoo. PSR J1926+1613 is projected inside the disk template for ‘Src A’ and near the edge of the SNR shell. Since the distance ($d$) and the spin-down luminosity of PSR J1926+1613 are unknown, we can not determine the contribution from the pulsar. However, the extended GeV morphology and flat SED of ‘Src A’ seem unfavorable to a pure emission from the pulsar or the putative PWN. For 4FGL J1924.3+1628, both PSR J1924+1631 and PSR J1924+1628 are much closer (with an angular distance of $\sim 0.08^{\circ}$ and $\sim 0.13^{\circ}$, respectively) to it than the SNR is. According to the dispersion measures, PSR J1924+1631 and PSR J1924+1628 are predicted to be at distances of 14.0 kpc (Nice et al., 2013) and 10.5 kpc (Parent et al., 2022), respectively. For 4FGL J1924.3+1628, it would have a luminosity of $\sim 5\times 10^{34}(d/10.5\ {\rm kpc})^{2}$ erg s*-1* in the energy range 0.2–500 GeV, which is two orders of magnitude higher than the spin-down luminosities of $2.4\times 10^{32}$ erg s*-1* for PSR J1924+1628 and $5.7\times 10^{32}$ erg s*-1* for PSR J1924+1631. Thus these two pulsars seem incapable of powering the $\gamma$-ray source 4FGL J1924.3+1628.

4.3.2 SNR G51.26+0.11

Noting that ‘Src A’ essentially coincides SNR G51.26+0.11 in the line of sight, we focus on the SNR origin of ‘Src A’. As shown above, there is also a spatial coincidence between the SNR and the +54 km s*-1* MCs at a distance of $\sim$ 6.2 kpc. Moreover, the broadened profiles of the CO emission lines indicate that the MCs, ‘Cloud A’ and ‘Cloud B’, both are likely perturbed by the shock of SNR G51.26+0.11. Thus a plausible explanation for ‘Src A’ is the hadronic interaction between the relativistic protons accelerated by the shock of SNR G51.26+0.11 and the adjacent MCs. To analyze the possible contribution from SNR G51.26+0.11 to the $\gamma$-ray emission of ‘Src A’, we fit the obtained $\gamma$-ray spectrum in both hadronic and leptonic scenarios using the Naima package (version 0.9.1, Zabalza 2015). We assume the accelerated protons and electrons have a power-law distribution in energy with an exponential cutoff $E_{i,\rm cut}$:

[TABLE]

where $i={\rm p,e}$, $E_{i}$ is the particle kinetic energy, and $\alpha_{i}$ is the power-law index. The normalization is determined by the total energy in particles with energies above 1 GeV ($W_{i}$).

In the hadronic scenario, by adopting a distance of 6.2 kpc to the SNR, the proton index $\alpha_{\rm p}\approx 2.25$ and the lower limit of the proton cutoff energy $E_{\rm p,cut}\approx$ 5 TeV are obtained with the GeV data alone (see the left panel of Figure 8). The proton index ($\alpha_{\rm p}$) we obtain is slightly larger than that (2.11${}^{+0.08}_{-0.12}$) reported by Araya (2021) and hence more consistent with the known gamma-ray-bright interacting SNRs (Acero et al., 2016), and a PWN contribution is not necessary to explain the hard spectrum given 1Disk+1PS as the best template. The cutoff energy can not be well constrained and is fixed as 1 PeV in the calculation. We first consider the contribution from the SNR interaction with ‘Cloud A’, which contains all the diffuse molecular gas in the disk template of ‘Src A’. Because ‘Cloud A’ almost overlaps with the SNR along the line of sight, the fraction of the SNR-MC interaction region to the SNR surface is assumed to be $f\sim 0.5$. With the parameters displayed in Table 5, we use $n_{\rm H}\sim 70$ cm*-3* as the number density of the target gas for proton-proton (PP) hadronic interaction. The energy required in protons is $W_{\rm p}\sim 2.0\times 10^{49}$($d$/6.2 kpc)2($n_{\rm H}$/70 cm*-3*)${}^{-1}(f/0.5)^{-1}$ erg, which is feasible in the context of SNR scenario. For comparison, we also estimate the contribution from the SNR interaction with ‘Cloud B’ in the northeast. The interaction region, assumed as a spherical cap, represents a fraction of the SNR surface with an order of $f\sim 0.01$. We set the value of $n_{\rm H}$ to be 430 cm*-3* and obtain $W_{\rm p}\sim 1.6\times 10^{50}$($d$/6.2 kpc)2($n_{\rm H}$/430 cm*-3*)${}^{-1}(f/0.01)^{-1}$ erg, which is also acceptable for the SNR scenario. But it is likely that the molecular gas surrounding the SNR, including ‘Cloud A’ and ‘Cloud B’, all contribute to the hadronic emission 555It cannot be ruled out that ‘Cloud A’ and ‘Cloud B’ are comprised in the same MC complex. Although no significant spectral line-broadening features can be distinguished in the central regions of ‘Cloud A’, this does not imply that no interaction occurs.. As the $\gamma$-ray flux $F\propto\eta E_{\rm SN}fn_{\rm H}$, where $\eta$ is the acceleration efficiency and $E_{\rm SN}$ is the explosion energy of the supernova, the ratio between the $\gamma$-ray emission from the two clouds, i.e., ‘Cloud A’ and ‘Cloud B’, is estimated to be $\sim 8.14$. Thus the central diffuse gas contributes up to 90% of the observed $\gamma$-ray emission, satisfactorily explaining the morphology of ‘Src A’. On the other hand, if the molecular gas is proximate to the SNR, even if it is not directly shocked, the accelerated protons that escape from the shock front may still hit the molecules (e.g., Li & Chen, 2010; Ohira et al., 2011; Li & Chen, 2012), giving rise to hadronic $\gamma$-rays. Despite the total energy requirements would be increased in this case, the canonical supernova explosion energy is sufficient in our estimation.

For the leptonic scenario, we assume the radio and $\gamma$-ray emission are generated by the same electron population. Their interaction with the magnetic field, i.e., the synchrotron (Sync) radiation, accounts for the radio emission and the inverse Compton (IC) scattering of low energy photons accounts for the $\gamma$-ray emission. Incorporating the radio fluxes of G51.26+0.11 taken from Dokara et al. (2021), we fit the broad-band SED using an exponential cutoff broken power-law with a break energy $E_{\rm break}$ and a cutoff energy $E_{\rm e,cut}$. Due to synchrotron losses of the electrons, we set the spectral index below the break to be $\alpha_{\rm e}$ and above the break to be $(\alpha_{\rm e}$+1). Lack of constraint by the current data, the cutoff energy is fixed 10 TeV based on the X-ray observations on the shell-type SNRs (Reynolds & Keohane, 1999). Besides the cosmic microwave background (CMB), the seed photon fields we use include a far-infrared (FIR) field with a temperature of 30 K and an energy density of 0.5 eV cm*-3* and a star light (SL) field with a temperature of 4000 K and an energy density of 1.0 eV cm*-3*, estimated at a distance of 6.2 kpc from the interstellar radiation field model (Shibata et al., 2011). As shown in the right panel of Figure 8, the data can be fitted with parameters $\alpha_{\rm e}\approx 2.1$ and $E_{\rm break}\approx 50$ GeV, as well as magnetic field $B\approx 3.2\mu$G. The total required energy in primary electrons is $W_{\rm e}\sim 3.4\times 10^{49}$($d$/6.2 kpc)2 erg, which would imply that the leptonic scenario is also feasible energetically. However, the cooling time is about $10^{7}$ yr to produce the break energy of 50 GeV in the magnetic field of 3.2 $\mu$G, which seems to be unphysically long for an SNR with an obvious radio shell morphology. Considering the break feature cannot be explained by the cooling losses alone, alternative models, such as age-limited (Ohira & Yamazaki, 2017; Zhang & Liu, 2019) or escape-limited (Ohira et al., 2010; Ohira et al., 2012; H. E. S. S. Collaboration et al., 2018) acceleration, should be included. This is beyond the scope of this study and not investigated here. In addition, the low magnetic field strength obtained is somewhat lower than expected in the diffusive shock acceleration paradigm. As the scenario is also subject to constraints from X-ray emission, relevant observations of SNR G51.26+0.11 are encouraged.

5 Conclusions

We reanalyze the Fermi-LAT $\gamma$-ray emission in the region of SNR G51.26+0.11 and study its interstellar environment based on the CO-line and H I data. Below are the main results of this study:

1. i.

Using the MWISP survey data, we find some MCs at $v_{\rm LSR}$ $\sim$ +54 km s*-1* overlap with SNR G51.26+0.11 along the line of sight and a dense structure tangent to the northeastern shell of the SNR. Given the morphological coincidences between the molecular gas distribution and the SNR and the asymmetric or broad CO line profiles, it is very likely that the MCs are associated with the SNR. 2. ii.

According to the LSR velocity of the associated MCs and the H I absorption features, SNR G51.26+0.11 is found to be located behind the tangent point at a farther distance, $6.2\pm 0.5$ kpc. 3. iii.

Based on 13.2 years of Fermi-LAT observations, we resolve the extended $\gamma$-ray emission detected in the region of SNR G51.26+0.11 into a uniform-disk $\gamma$-ray source (‘Src A’/4FGL J1925.4+1616) with a significance of about 19.5$\sigma$ and a point source (4FGL J1924.3+1628) with a significance of about 4.2$\sigma$ in 0.2–500 GeV. The best-fitting position and 68%-containment radius of ‘Src A’ are R.A. = $291.38\pm 0.02^{\circ}$, Dec = $16.27\pm 0.02^{\circ}$ (J2000) and $\sigma={0.17^{+0.02}_{-0.02}}^{\circ}$. Its emission can be described by a power-law spectrum with a photon index of ${\rm\Gamma=2.19\pm 0.04}$. The corresponding 0.2–500 GeV flux is $2.86\times 10^{-11}\ {\rm erg\ cm^{-2}\ s^{-1}}$, and the luminosity in 1–100 GeV is $7.3\times 10^{34}$ erg s*-1* at a distance of 6.2 kpc. 4. iv.

A natural explanation for ‘Src A’ is the hadronic scenario in which the energetic protons accelerated by SNR G51.26+0.11 interact with the adjacent MCs. The proton spectrum obtained has an index of $\sim$ 2.25. The leptonic scenario, even if unlikely, cannot be conclusively ruled out.

Acknowledgements

This research made use of the data from the Milky Way Imaging Scroll Painting (MWISP) project, which is a multi-line survey in 12CO/13CO/C18O along the northern galactic plane with PMO-13.7m telescope. We are grateful to all the members of the MWISP working group for their support. W.J.Z. thanks Bing Liu, Qiancheng Liu, Tianyu Tu, Ping Zhou and Xin Zhou for helpful discussions. X.Z. and Y.C. acknowledge the support of National Key R&D Program of China under nos. 2018YFA0404204 and 2017YFA0402600 and the NSFC grants under nos. U1931204, 12173018, and 12121003.

Data Availability

The Fermi-LAT and VLA H I data underlying this work are publicly available, and can be downloaded from https://fermi.gsfc.nasa.gov/ssc/data/access/lat/ and https://www.cadc-ccda.hia-iha.nrc-cnrc.gc.ca/en/search/, respectively. The CO data used in this study will be made available by the corresponding authors upon request.