© The Authors 2025

1 Introduction

Supernova remnants (SNRs) are widely recognized as the leading candidates for accelerating Galactic cosmic rays (CRs) up to energies below the so-called knee in the CR spectrum, at about the Peta-electronvolt (PeV) scale (Hillas 2005). In these extreme environments, high-energy CRs give rise to γ-ray emission through different radiative processes. In the leptonic scenario, γ-rays are produced via inverse-Compton scattering or bremsstrahlung of relativistic electrons. Conversely, the hadronic scenario involves inelastic collisions between accelerated protons and ambient matter that result in the production of neutral pions that in turn decay into γ-rays. A particularly distinctive signature of hadronic CR acceleration is the detection of γ-ray emission from dense molecular clouds (MCs) that are illuminated by CRs that have escaped nearby SNRs. Observational support for this scenario has been provided by the Fermi Large Area Telescope (Fermi-LAT), which has reported γ-ray emission that is consistent with CR-MC interactions in several well-known SNRs, including IC 443 (Ackermann et al. 2013; Abdo et al. 2010), W44 (Uchiyama et al. 2012; Peron et al. 2020), W28 (Aharonian et al. 2008; Li & Chen 2010; Hanabata et al. 2014), W51C (Abdo et al. 2009b), G150.3+4.5 (Li et al. 2024), and G15.4+0.1 (Li et al. 2023b). These remnants exhibit prominent GeV γ-ray emission that is generally attributed to hadronic interactions between shock-accelerated protons and nearby dense gas. In particular, the γ-ray spectra of W44 and W51C reveal a characteristic spectral bump that is associated with neutral pion decay. This provides some of the most compelling evidence to date that relativistic protons are accelerated in SNRs (Giuliani et al. 2011; Jogler & Funk 2016).

The SNR γ-Cygni (G78.2+2.1) is a puzzling γ-ray source located in the complex Cygnus region. The multiwavelength analyses toward this well-studied region have never ceased, and the origins of its γ-ray emission is still not well understood. Ladouceur & Pineault (2008) revealed that the spectral index of the radio flux −α_v varies between ~−0.8 and ≲−0.4 across the whole SNR region. The softest index is found in the bright southeastern part, and the northwest part is harder (~0.55) (Ladouceur & Pineault 2008; MAGIC Collaboration 2023). Near the center lies the energetic GeV-bright pulsar PSR J2021 + 4026, whose pulse profile has been provend to vary (Wang et al. 2024; Razzano et al. 2023). The unknown distances make the γ-ray emission in this region also more puzzling. In the very high-energy (VHE) energy band, VERITAS¹ observations report an extended source of ~0.23° in the northwest part (Aliu et al. 2013) that was characterized using an ellipse template (0.29° × 0.19°) with a power-law (PL; dN/dE ∝ E^−α) index of 2.79 ± 0.39_stat ± 0.20_sys (Abeysekara et al. 2018). The location of the VERITAS source is roughly consistent with part of the MAGIC² source (MAGIC Collaboration 2023), and both of them are spatially inconsistent with 2HWC J2020+403 measured by HAWC³ (Abeysekara et al. 2017b). Furthermore, MAGIC Collaboration (2023) found that there is a bright arclike TeV source located west of the SNR that is not detected in the GeV band. This makes the origins of the γ-ray emission in this region even more puzzling.

Pulsar halos are a new type of γ-ray source that was identified in recent years. As pulsar wind nebulae (PWNe) age, high-energy electrons and positrons may escape from the nebula and diffuse in the surrounding interstellar medium (ISM) after ~100 kyr (Giacinti et al. 2020). The electrons and positrons that diffuse out of the PWN produce γ-ray emission through inverseCompton scattering of ambient photon fields (López-Coto et al. 2022). The recent detection of extended γ-ray emission around Geminga and Monogem (Abeysekara et al. 2017b,a; Di Mauro et al. 2020) and the slow-diffusion region around them has sparked a widespread discussion about TeV halos, while the nature of the emission is still uncertain. It was also noted that pulsar halos may exhibit an asymmetric morphology when the coherence length is large enough (López-Coto & Giacinti 2018), which indicates potentially anisotropic diffusion inside the pulsar halo. Limited by the resolution and systematic uncertainties, however, an asymmetric morphology of these pulsar halos has not yet been confirmed by experiments (Abeysekara et al. 2017a; H.E.S.S. Collaboration 2023).

We conduct a comprehensive analysis of the GeV γ-ray emission toward the γ-Cygni area based on 15 years of data from the Fermi-LAT instrument. The photon events were selected in an energy range from 100 MeV to 1 TeV, and the results are presented in Sect. 2. In Sect. 3 we present the results of the gas observation of ¹²CO(J = 1–0) that we obtained from the Milky Way imaging scroll painting (MWISP). Sect. 4 discusses potential origins of the γ-ray emission by comparing multiwavelength observations with theoretical expectations. Last, Sect. 5 provides our conclusions.

2 Fermi-LAT data reduction

In the subsequent analysis, we adopt the standard LAT analysis software Fermitools v2.2.0 together with the Fermipy v1.1.6 (Wood et al. 2017) to quantitatively calculate the extension and position of extended sources. The most recent Fermi-LAT (Atwood et al. 2009) Pass 8 data were collected from August 4, 2008 (mission elapsed time 239557418) to August 4, 2023 (mission elapsed time 712800005) to study the GeV emission around γ-Cygni. We selected the data with a source event class P8R3_SOURCE (evclass=128) and event type FRONT + BACK (evtype=3), with the standard data-quality selection criteria (DATA_QUAL > 0)&&(LAT_CONFIG == 1), and we excluded the zenith angle over 90° to avoid contamination from the Earth limb. To derive a better point-spread function and reduce contamination from the pulsar, we limited the photon energies to between 20 GeV and 1 TeV for the further morphological analysis. Additionally, events with energies between 100 MeV and 1 TeV were selected for a more detailed spectral analysis. All photon events within a region of interest (ROI) of 20° × 20° centered at the position of γ-Cygni were modeled using the Fermi-LAT 4FGL catalog data release 4 (4FGL-DR4; Abdollahi et al. 2020; Ballet et al. 2023) in a binned maximum likelihood analysis (Mattox et al. 1996). The instrument response functions (IRF) P8R3_SOURCE_V3⁴, along with the Galactic and isotropic diffuse background models (IEM, gll_iem_v07.fits) and (iso_P8R3_SOURCE_V3_v1.txt) were adopted. All sources listed in the 4FGL-DR4 catalog were included in the background model, all sources within 7° from the center of ROI and two diffuse backgrounds above were set free, ass generated by the script make4FGLxml.py⁵. We adopted the likelihood test statistic (TS) to calculate the significance of the γ-ray sources, which is defined as TS = 2(ln ℒ₁ − ln ℒ₀), where ℒ₁ and ℒ₀ are maximum likelihood values for the model with and without the target source. Furthermore, the TS_ext is defined as TS_ext = 2(ln ℒ_ext − ln ℒ_ps), where ℒ_ext and ℒ_ps represent the maximum likelihood values for the extended and point-like templates, respectively. This calculation only considered one additional free parameter that was introduced by the extended template, and the significance is approximately given by $\sqrt{{\mathrm{T}\mathrm{S}}_{\mathrm{e}\mathrm{x}\mathrm{t}}}$ in units of σ.

Fig. 1 Folded pulse profile and two-dimensional phaseogram in 32 phase bins derived for PSR J2021+4026. The grayscale represents the weights of photons in each bin. White indicates a higher weight, and black indicates a lower weight. The dotted and dashed lines mark the minimum and maximum phases of the off-pulse phase interval, respectively.

2.1 Timing analysis

We performed a timing analysis of the 0.1–1000 GeV LAT data of the PSR J2021+4026 to separate its γ-ray emission from other sources in the region. We first included the LAT photons within an aperture radius of 3° and weighted them by their probability of originating from the pulsar (using gtsrcprob in the Fermitools). The LAT photons with low weights(<0.001) were excluded from the analysis. We assigned pulse phases to the LAT photons using the Fermi plugin of TEMPO2 (Hobbs et al. 2006; Edwards et al. 2006), according to the known ephemeris given in the LAT third pulsar catalog (Smith et al. 2023). This ephemeris covers the time period of MJD 54689–58175, and we therefore only included the LAT data during this time period. The weighted pulse profile and the two-dimensional phaseogram are plotted in Figure 1. We defined the phase 0.28–0.47 as the off-pulse phase range and the remaining phase as the on-pulse phase range.

Then, we assigned pulse phases to all of the LAT events in the ROI. With LAT events during different phase ranges, we generated counts maps in different energy ranges to verify whether we were able to avoid the contamination from PSR J2021+4026. The resulting count maps are depicted in the top row of Figure 2. The top left panel of Figure 2 shows that most >10 GeV photon events are concentrated around PSR J2021+4026 within the full-phase (which includes the on- and off-pulse) range. Almost no γ-ray excess is detected when >20 GeV photon events are selected, however, as shown in the top right panel of Figure 2. Additionally, we calculated the TS maps in the vicinity of PSR J2021 + 4026 for both the on-pulse (bottom left panel of Figure 2) and off-pulse phase using >10 GeV photon events. The significant excess around the pulsar is still present even when we only included events during the off-pulse phase range (shown in the bottom right panel of Figure 2). These results indicate that it is not enough to separate the pulsar contamination in the >10 GeV energy band, where the TS value of the pulsar still exceeds 500, while when we adopted >20 GeV data events, the TS value was only ~4, in contrast to ~249 000 measured with >100 MeV data events. Therefore, we adopted >20 GeV full-phase range data events to conduct a quantitative analysis of the spatial components within the region.

Fig. 2 Top row: Fermi-LAT count maps above 10 GeV and 20 GeV within the full-phase range. Bottom row: Fermi-LAT TS maps above 10 GeV during the on- and off-pulse phases, respectively. The black diamond indicates the position of the gamma-ray pulsar PSR J2021+4026 (Abdo et al. 2009a).

2.2 Morphological analysis

In the beginning, we attempted to verify the γ-ray emission position offset with energy and the energy-dependent morphology. Combined with the above results, we only selected photon events above 20 GeV to generate TS maps, and all TS maps were generated by only considering background fitting, but did not include the extended source γ-Cygni. The γ-ray contribution from PSR J2021 +4026 was removed, which was modeled by a point-like source whose best-fit location in the full-phase range given by gtfindsrc was marked as RA = 305.384° ± 0.009°, Dec = 40.441° ± 0.012°. This is consistent with the position given in the 4FGL-DR4 source list, and the extended γ-Cygni source is described as a uniform disk with 68% error radius R₆₈ = 0.517°. Since Fraija & Araya (2016) and MAGIC Collaboration (2023) argued that there are two different spectrum components inside the γ-Cygni region, we directly started with two-component templates to fit the γ-ray emission. First, we used a point-like source plus a uniform disk and 2D Gaussian template as model 2 and model 3 (listed in Table 1). Their TS values were calculated to be 839 and 846, respectively, and both improved significantly compared to the single-disk template given in the 4FGL-DR4 catalog (model 1).

Then, we further tested the two extended Gaussian templates as model 4, where the location of the smaller Gaussian was perfectly consistent with TeV observation results given by the VERITAS Collaboration (Aliu et al. 2013). The best-fit results for these two Gaussian sources were recorded as RA = 305.277°, Dec. = 40.431°, r₆₈ = 0.622° and R.A. = 305.113°, Dec = 40.832°, r₆₈ = 0.151°, and they are shown as the solid yellow and green circles in Figure 3 (hereafter, the 0.622° Gaussian source that spatially corresponds to the whole SNR region is called SrcFermiThe 0.151° Gaussian source spatially corresponds to the VERITAS region and is called Src_VERITAS). In this scenario, the TS_ext for Src_VERITAS was calculated to be approximately 26, which rejects the point-like source hypothesis at the 5.1σ level. According to Lande et al. (2012), when TS_ext > 16, the extended-source hypothesis is valid. These results suggest that the γ-ray emission from the γ-Cygni SNR can be described by two extended components.

To search for the correlation between GeV and TeV emission, we further tested other templates that also included two extended components, but directly adopted the location and extension parameters from the VERITAS measurement results. They are labeled model 5 and model 6, and they are represented by the best-fit 0.622° Gaussian template plus a 0.23° circle given by Aliu et al. (2013) or by the 0.29° × 0.19° ellipse given by Abeysekara et al. (2018). These sources are shown as a dashed cyan circle and a solid cyan ellipse in Figure 3, respectively (the uniform ellipse template was generated by the pyFits⁶). Then, we tested the 0.622° Gaussian template plus two extended components presented by MAGIC (MAGIC Collaboration 2023), labeled model 7, which is shown as red circles in Figure 3. The TS values obtained by models 5, 6, and 7 are lower than the two Gaussian templates (model 4). Last, to test whether we missed other potential sources, we added an extended source to model 4 to fit three Gaussian components simultaneously. The extension and location were all set as free parameters, but the TS value did not increase (model 8). The results are summarized in Table 1. In this table, the Akaike information criterion (AIC) value is defined as AIC = 2k − 2 ln ℒ and was described by Akaike (1974), where k is the number of degrees of freedom of the model, and ℒ is the likelihood value. The model with the lowest AIC value is preferred. We therefore adopted model 4 for the following spectral analysis.

Fig. 3 Fermi-LAT TS maps in the vicinity of γ–Cygni. The energy range for each panel is shown above them. The solid yellow and green circles show the extension size of the best-fit two Gaussian template (R68). The red circles show the two extended TeV sources measured by MAGIC (MAGIC Collaboration 2023), called MAGIC-1 and MAGIC-arc, respectively. The cyan circle and ellipse show the VERITAS results (Aliu et al. 2013; Abeysekara et al. 2018). The orange circle shows the center of the HAWC source with the 1 σ error radius (Abeysekara et al. 2017b). The black diamond shows the position of the gamma-ray pulsar PSR J2021+4026 (Abdo et al. 2009a) as in Figure 2. The white contours come from the CGPS 1420 MHz radio observation of the SNR shell structure (Ladouceur & Pineault 2008).

Fig. 4 γ-ray SED obtained from the events in 0.1–1000 GeV shown as black data points. The gray solid line shows the best-fit PL spectra for each component. The dashed lines show the 68% confidence bands. The arrows indicate the 95% upper limits when the TS value is lower than 4 in a given bin. The blue data points are the VERITAS results (Abeysekara et al. 2018). The red butterfly was extracted from Abeysekara et al. (2017b). The gray histogram denotes the TS value for each bin.

2.3 Energy spectrum

In Section 2.2, model 4 performed best and was adopted for the further spectral analyses. After using events from 100 MeV to 1 TeV, we obtained different spectral types, such as a simple PL and broken power law (BPL⁷), and the results are then summarized in Table 2. For Src_Fermi, the BPL model does not provide a notable improvement of the fitting quality compared with a PL assumption. This could be quantified as TS_curve, defined as TS_curve = 2(ln ℒ_BPL − ln ℒ_PL) (Abdollahi et al. 2020), and the derived value of 4.1 corresponds to a significance level of ~2.0 σ. We therefore suggest that there is no energy break in the Src_Fermi spectrum. For Src_VERITAS, its TS_curve = 3.6 is even lower than 2 σ, and its calculated TS value shows barely any difference between a PL and a BPL assumption. We therefore suggest that its spectrum can also be described by a single PL. Then, we obtained their spectral energy distributions (SEDs) by separating the events in the 100 MeV–1 TeV energy range into 12 logarithmically equal intervals and repeated the likelihood evaluation analysis for each interval. In this part, the normalization values of all sources were left free and the other parameters were fixed. For bins with a TS values lower than 5.0, we calculated the upper limits using a Bayesian method (Helene 1983) at a 95% confidence level. The derived SEDs are shown in Fig. 4, together with the global best-fit spectra in the energy range of 100 MeV–1 TeV.

3 Gas observations

We used the data from the Milky Way Imaging Scroll Painting (MWISP⁸) project with the high-resolution carbon monoxide (CO) survey along the northern Galactic plane with the 13.7 m telescope of the Purple Mountain Observatory (Su et al. 2019). Furthermore, Ladouceur & Pineault (2008) argued that the HI emission at velocities between [−19.0, −11.0] km s⁻¹ and [3.0, 10.0] km s⁻¹ that was measured by the Canadian Galactic Plane Survey (CGPS) appears to circumscribe the radio continuum emission, which is considered to be a signature of association with γ-Cygni, and the average HI spectrum (Fig. 17 in Ladouceur & Pineault 2008) shows that the cloud might be disturbed by the shock front in the velocity range [−16.0, +16.0] km s⁻¹ (see also Figure 37, 38 and 39 in Ladouceur & Pineault 2008). In the direction of γ-Cygni, the standard velocity-distance relation becomes unreliable as a result of complex gas kinematics and significant noncircular motions. As discussed by Reid et al. (2009, 2014), the Galactic bar, spiral arm shocks, and streaming motions near star-forming regions introduce strong deviations from the standard Galactic rotation curve. In particular, the line of sight toward the Cygnus X region (l ~ 75° − 85°) intersects multiple overlapping molecular cloud complexes (Schneider et al. 2006, 2007; Gottschalk et al. 2012), which makes it difficult to associate the radial velocity with a unique distance. Moreover, large-scale motions such as superbubble expansions and spiral-arm streaming further distort the local velocity field (Reid et al. 2014). Therefore, instead of relying on a kinematic distance model, we considered the HI studies by Ladouceur & Pineault (2008), which provided a more accurate representation of the gas structures that are likely associated with γ-Cygni in this complex region. Based on this, we integrated the radial velocity range [+9.0, +12.0] km s⁻¹ and [−16.0, −11.0] km s⁻¹ to calculate the molecular hydrogen column mass.

We assumed that the molecular cloud and γ-Cygni are associated and then adopted the distance d = 1.5 kpc for further calculations (which were derived from the assumption that γ-Cygni is associated with the nearby star-forming region G077.901+01.769; (Solin et al. 2012); while in the Cygnus region direction, we cannot simply obtain the actual distance through the velocity). Under the assumption of θ = 0.622°, the obtained physical radius is ~16.2 pc. To determine the column density of H₂ in this region, we employed a conversion factor X_CO = 2 × 10²⁰ cm⁻² K⁻¹ km⁻¹ s (Bolatto et al. 2013; Dame et al. 2001). Using this factor, the column density N_H2 was calculated as N_H2 = X_CO × W_CO. Consequently, the mass of the molecular complex can be derived from the integrated intensity (W_CO), $$M=\mu m_{\mathrm{H}}{D}^2\mathrm{\Delta }{\mathrm{\Omega }}_{\mathrm{p}\mathrm{x}}{X}_{\mathrm{C}\mathrm{O}}\sum _{\mathrm{p}\mathrm{x}}{W}_{\mathrm{C}\mathrm{O}}\propto N_{{\mathrm{H}}_2}.$$(1)

In this formula, μ was set to 2.8, reflecting a relative helium abundance of 25%. The mass of a hydrogen nucleon is denoted as m_H. The solid angle subtended by each pixel is given by ΔΩ_px. The term ∑_pxW_CO accounts for the velocity binning of the data cube, which is calculated by summing the map content for the pixels within the target sky region and the specified velocity range, and then scaled by the velocity bin size.

As shown in the left panel of Figure 5, in the velocity range [+9.0, +12.0] km s⁻¹, only the VERITAS source region spatially corresponds to the gas distribution and the location of the γ-ray emission. The green circle indicates the radius of Src_VERITAS. Since part of the molecular cloud (MC) is located outside of the projection size of Src_VERITAS characterized in Sect. 2.2, the magenta circle was used to calculate the total MC mass in this velocity range, which includes the MC mass that might interact with the escaped CRs. In the right panel of Figure 5, within the velocity range of [−16.0, −11.0] km s⁻¹, the brightest part of the SNR shell in the southeast coincides well spatially with the gas distribution. The red circle denotes the total MC mass in this velocity range, and the orange circle indicates the 1 σ uncertainty radius of the HAWC source. Details about their correlation can be found in Sect. 4. Combined with the HI spectrum, which also supports the association between the SNR and the gas in this velocity range (MAGIC Collaboration 2023; Ladouceur & Pineault 2008), we calculated the density and mass of each cloud to derive a more precise normalization factor for the subsequent hadronic model.

By using the estimate made for N_H2, we calculated in the velocity ranges [+9.0, +12.0] km s⁻¹ the mass of the molecular cloud within the Src_VERITAS region (green circle) to be about ${\text{M}}_{\text{veritas\hspace{0.17em}}}=102{\text{d}}_{1.5}^2{\text{M}}_{\odot }$, and the gas mass within the magnetic circle region was calculated to be ${\mathrm{M}}_{\text{CloudV\hspace{0.17em}}}=171{\mathrm{d}}_{1.5}^2{\mathrm{M}}_{\odot }$. Assuming a spherical geometry of the gas distribution, we estimated the volume to be ${\mathrm{V}}_{\text{veritas\hspace{0.17em}}}=\frac{4\pi }{3}{\mathrm{R}}^3$, where R = d × θ, and the average H₂ cubic density in this region is about ${\mathrm{n}}_{\text{veritas\hspace{0.17em}}}=16{\mathrm{d}}_{1.5}^{-1}{\mathrm{c}\mathrm{m}}^{-3}$. Correspondingly, in the velocity ranges [−16.0, −11.0] km s⁻¹, the gas mass within the spherical orange and red regions was calculated to be ${\text{M}}_{\text{hawc\hspace{0.17em}}}=63{\text{d}}_{1.5}^2{\text{M}}_{\odot }$ and ${\text{M}}_{\text{CloudH}}=542{\text{d}}_{1.5}^2{\text{M}}_{\odot }$, respectively. Since the >100 GeV TS map in Figure 3 shows a decent overlap between the γ-ray peak and the location of the HAWC source, with the excess remaining above 3 σ significance and extending beyond 200 GeV, we assumed that the gas within the two source regions (VERITAS and HAWC) corresponds to the trapped component. The gas masses we adopted to calculate the γ-ray fluxes within the trapped emission zones are therefore M_veritas and M_hawc. Correspondingly, the gas located outside the source region but within the region that is affected by escaped CRs (CloudV and CloudH shown in Figure 5) corresponds to the escaped component, whose gas mass we adopted to calculate the γ-ray flux within the escaped emission zone. Therefore, we subtracted the gas mass in the source region from that in the cloud region to estimate the mass of the escaped component. This yielded approximate values of ${\mathrm{M}}_{\text{escapeH\hspace{0.17em}}}={\mathrm{M}}_{\text{CloudH\hspace{0.17em}}}-{\mathrm{M}}_{\text{hawc}}=479{\mathrm{d}}_{1.5}^2{\mathrm{M}}_{\odot }$ and ${\mathrm{M}}_{\text{escapeV\hspace{0.17em}}}={\mathrm{M}}_{\text{CloudV\hspace{0.17em}}}-{\mathrm{M}}_{\text{veritas\hspace{0.17em}}}=69{\mathrm{d}}_{1.5}^2{\mathrm{M}}_{\odot }$, which represents the gas mass associated with the escaped components in different emission zones. Considering the radius of each spherical region, we calculated their densities to be about several dozen cm⁻³. To be more precise, they fluctuate between approximately $4{\mathrm{d}}_{1.5}^{-1}{\mathrm{c}\mathrm{m}}^{-3}$ and $29{\mathrm{d}}_{1.5}^{-1}{\mathrm{c}\mathrm{m}}^{-3}$.

Fig. 5 Integrated 12CO (J = 1–0) emission intensity (K km s−1) toward γ-Cygni region in the velocity ranges [+9.0, +12.0] km s−1 (left panel) and [−16.0, −11.0] km s−1 (right panel) using MWISP data. The solid circles are the same as in Figure 3. The dashed magenta and red circles correspond to different regions that are affected by escaped CRs, called CloudV and CloudH, respectively.

Fig. 6 Modeling the γ-ray spectrum in the hadronic scenario. Left panel: total contribution (solid black line) from the sum of the dashed red and orange lines, corresponding to the escaped and trapped ions, respectively. The red butterfly was extracted from Abeysekara et al. (2017b). Middle panel: total γ-ray contribution calculated with different rs. Right panel: predicted γ-ray flux in the VERITAS region, obtained by adding the escaped and trapped components. The blue data points and butterfly were taken from Abeysekara et al. (2018).

4 Discussion of the possible origins of the γ-ray emission

4.1 SNR-MC interaction

Since the MC distribution shown in Section 3 spatially corresponds well with both the γ-ray emission measured by FermiLAT and part of the SNR shell, we suggest that the northwestern and southeastern shock fronts interact with nearby dense gas, similar to what has been observed in RX J1713.7-3946 (Sano et al. 2013; Tanaka et al. 2020). These SNR–MC interactions that occur in lower-density environments tend to result in relatively harder spectral indices in the GeV energy band (Yuan et al. 2012). This scenario implies that only parts of the SNR shell interact with molecular clouds of varying densities, leading to inhomogeneous γ-ray emission. The soft spectral component (the escaped contribution shown in Figure 6) observed in the lower GeV band can be interpreted as escaped cosmic rays (CRs) illuminating nearby MCs, even if their projected positions appear to be located within the SNR shell (Li et al. 2023a). In contrast, the flatter spectrum component (the trapped contribution shown in Figure 6) is likely produced by CRs that are still trapped within the shock region and interact with the local dense gas. Moreover, based on the gas-mass calculations presented in Section 3, the gas-mass ratios ϵ in the two different cloud regions shown in Figure 5 are ϵ_V = M_veritas/M_escapeV = 102/69 ≈ 1.48 and ϵ_H = M_hawc/M_escapeH = 63/479 ≈ 0.13, respectively. This gas-mass ratio ϵ differs by more than an order of magnitude, indicating that the γ-ray emission from these two emission regions is mainly contributed by trapped ions and escaped ions. This difference also naturally accounts for the discrepancy in Src_Fermi and Src_VERITAS source detections and their different spectrum characteristics.

We considered a scenario in which protons are injected instantaneously into a uniform emission zone (corresponding to CloudV and CloudH mentioned in Section 3) approximately 7000 years ago (same as the SNR age argued by Zeng et al. 2019), and the spectrum of the injected protons can be characterized as a broken power-law spectrum, $$Q(E)=Q_0\frac{{\left(E/E_{p,\text{break\hspace{0.17em}}}\right)}^{-{\gamma }_1}}{1+{\left(E/E_{p,\text{break\hspace{0.17em}}}\right)}^{{\gamma }_2-{\gamma }_1}}.$$(2)

E_p,break represents the break energy, where γ₁ and γ₂ are the spectral indices below and above the energy break, respectively. Because of the flat γ-ray spectrum ~2.0 in the GeV energy band for the VERITAS region, we set γ₁ = 2.0 for simplicity and γ₂ = γ₁ + 1.0 (Zeng et al. 2019). E_p,break was fixed to 5 TeV, and we then calculated the distribution of the escaped protons in the emission region following the method outlined by Liu et al. (2020), $$N_p(E,r_{\mathrm{s}},T)=\frac{Q(E)}{[4\pi D(E)T{]}^{\frac{3}{2}}}\exp \left[\frac{-r_{\mathrm{s}}^2}{4D(E)T}\right].$$(3)

The diffusion coefficient was assumed to be uniform, following the form D(E) = χD₀(E/E₀)^δ for E > E₀, where D₀ = 1 × 10²⁸ cm² s⁻¹ at E₀ = 10 GeV and δ = 1/2, which is consistent with Kraichnan turbulence (Kraichnan 1965; Blasi 2013). As a result of projection effects, the actual distance between the gas complex and the supernova remnant (SNR) remains uncertain. Furthermore, in the direction of the Cygnus region, the relation between kinematic distance and gas velocity can no longer be derived by the method described by Reid et al. (2009, 2014). We therefore adopted r_s as a free parameter representing the distance between the injection site and the illuminated molecular clouds. With the injected source spectrum defined by Q(E) ∝ E^−Γ and D(E) ∝ E^δ, Equation (3) reveals that N_p(E) shows a lower-energy spectral cutoff at E_p,break when $\sqrt{4D\left(E_b\right)T}\simeq r_{\mathrm{s}}$. Correspondingly, N_p(E) follows $N_p(E)\propto E^{-(\mathrm{\Gamma }+\frac{3}{2}\delta )}$. The total energy of injected protons is denoted W_inj = η E_SN, where η is the efficiency of converting kinetic energy into accelerated protons, which is typically set to 0.1. The kinetic energy of the SNR, E_SN is generally taken to be 10⁵¹ erg (Blasi 2013). Additionally, the correction factor χ for the diffusion coefficient was set to 1.0, which corresponds to the standard value of the Galactic diffusion coefficient (Blasi 2013). The corresponding γ-ray fluxes produced in escaped and trapped emission zones were calculated with the naima package (Zabalza 2015).

For the total γ-Cygni region, the resulting γ-ray flux can be well described by adding the escaped ions together with the trapped ion contribution, shown as dashed red and orange lines in the left panel of Figure 6. The total energy of escaped and trapped protons above 1 GeV was calculated to be W_escapedH = 6.32 × 10⁴⁸(M_escapeH/479M_⊙)⁻¹ erg and W_trappedH = 8.98 × 10⁴⁹(M_hawc/63M_⊙)⁻¹ erg, respectively. As mentioned above, assuming a distance of 1.5 kpc, we calculated the source radius to be r_s = 16.2 pc. The precise value for r_s cannot be constrained because of the projection effect, however. We therefore added four scenarios in which r_s = 10/20/30/40 pc to show the dependence between r_s and the γ-ray flux. We plot the results in the middle panel in Figure 6. For the VERITAS region, we suggest that its GeV emission with a flat spectrum is mainly dominated by the ongoing shock-cloud interaction (Liu et al. 2022), and the escaped component contributes as well. Because the ratio value (ϵ_V) of the gas mass between the trapped emission zone and the escaped emission zone is relatively high, however, the escaped component can only contribute little. Correspondingly, the low ratio value (ϵ_H) of the gas mass between the trapped emission zone and the escaped emission zone would lead to a higher escaped γ-ray flux component in the lower-energy band, similar to the left panel in Figure 6. The total energy of each part can be calculated as W_escapedV = 4.77 × 10⁴⁸(M_escapeV/69 M_⊙)⁻¹ erg and W_trappedV = 4.31 × 10⁴⁹(M_veritas/102 M_⊙)⁻¹ erg. We also calculated the theoretical prediction for the trapped particles in the higher-energy band within model C proposed by Bell (2015). In this model, the higher-energy particles are trapped in a relatively small region compared with the lower-energy particles. More specifically, Figure 5 of Bell (2015) predicts that the CRs with energies ≳100 TeV (corresponding to the ~10 TeV photons in HAWC energy band) are trapped in a small region inside the SNR, whose radius is roughly one-fifth of the SNR shock radius. The shock radius is ~0.5 degree, as shown in Figure 5, and therefore, this model predicts that the high-energy CRs are confined in a region with a radius of ~0.1 degree. This is consistent with the measurement from HAWC (Abeysekara et al. 2017b) that 2HWC J2020+403 is a point-like source, with a 1 σ uncertainty radius of 0.11 degree. These high-energy CRs correspond to the CRs that were accelerated early during the expansion phase of the SNR when its shock speed was higher and that were advected inside the SNR. They remain trapped in tangled magnetic fields from CR-driven instabilities.

Fig. 7 Energy density of TeV sources calculated using two different approaches. The top panels are calculated as ${\epsilon }_{\mathrm{e}}=\dot{E}{\tau }_c/V$, and the bottom panels are calculated following the approach in Section 3.2 of Giacinti et al. (2020). The newly identified pulsar halo associated with PSR J0622+3749 (Aharonian et al. 2021), indicated by a star, is shown along with the Geminga halo and Monogem halo (Abeysekara et al. 2017a). The circles representing several γ-ray PWNe are based on data from Giacinti et al. (2020), and γ-Cygni under different distance assumptions is shown as squares.

4.2 Pulsar wind nebula and pulsar halo

The detected GeV-bright pulsar might mean that the γ-ray emission in this region might be powered by high-energy electrons and positrons generated by a pulsar wind nebula (PWN). Like MSH 15–52 (H.E.S.S. Collaboration 2018) and HESS J1825–137 (Principe et al. 2020), for example, typical γ-ray PWNe detected by Fermi-LAT are always driven by energetic pulsars with higher $\dot{E}$ between 10³⁶ ~ 10³⁹ erg s⁻¹ (Acero et al. 2013). The value of 10³⁵ erg s⁻¹ is lower by at least one order of magnitude, which corresponds to the PSR J2021 + 4026 associated with the γ-Cygni. The measured photon spectrum indices inside the γ-Cygni region are also much softer than typical PWNe with ~1.6 photon index in the GeV energy band, for example, MSH 15–52 (H.E.S.S. Collaboration 2018) and HESS J1640-465 (Mares et al. 2021). These arguments disfavor the scenario according to which the γ-ray emission originates from a PWN.

We determined whether the TeV emission might be interpreted as inverse-Compton scattering of ambient photon fields by relativistic electrons and positrons that have escaped from a potential PWN and diffuse into the ISM (López-Coto et al. 2022), leading to a pulsar halo around the pulsar. To do this, we adopted the two different approaches described by Giacinti et al. (2020) to estimate the energy density ε_e from the pulsar properties and the γ-ray luminosity. The region ε_e ≪ ε_ISM (area shaded in gray in Figure 7) corresponds to when the energy density in the relativistic electrons and positrons is negligible compared with that of the ISM. This corresponds to the region in which pulsar halos may form. In contrast, ε_e ≳ ε_ISM corresponds to the case when the emitting electrons are contained in the region that is dynamically dominated by the pulsar (Giacinti et al. 2020).

To present a more intuitive comparison, we also included the identified pulsar halo LHAASO⁹ J0622+3749 (Aharonian et al. 2021) in our figure, which is located near Geminga and PSR B0656+14 (Abeysekara et al. 2017a). Its parameters are summarized in Table 3. For the γ-Cygni region, we adopted the TeV extension size of θ = 0.23° in our calculations. Since the distance to the pulsar is unknown, we considered 1.0 kpc and 2.0 kpc as the minimum and maximum distances, assuming PSR J2021+4026 is associated with the SNR. The corresponding radii for each scenario are 4.014 pc and 8.028 pc, respectively, and they are shown with different squares in the Figure 7. The lower panel of Figure 7 shows that it is hard to definitively attribute the TeV emission to a potential pulsar halo. We note, however, that the green square, representing the larger distance (2.0 kpc), is situated near the boundary of shadow region, but not within it. Considering that the pulsar age is only 77 kyr, we suggest that this source might be in an intermediate stage of its transition from a PWN to a pulsar halo. This implies that the halo has not yet fully formed, causing its $\dot{E}$ and ε_e to fall between those of a typical PWN and a mature pulsar halo. Additionally, we observe that the calculated ε_e(L_γ) is approximately 0.17 eV/cm³ below the 2.0 kpc assumption. This suggests that even when the actual distance is much greater (e.g., 4.0 kpc), the energy density would not decrease sufficiently to fall entirely within the shaded area. These results contrast sharply with those for the three identified pulsar halos shown with stars in Figure 7, which are all located in the shaded areas in the bottom panels. Recent numerical simulations suggested, on the other hand, that the morphology of a pulsar halo can exhibit highly asymmetric and distorted shapes, especially when the coherence length of the interstellar turbulence exceeds 10 pc (López-Coto & Giacinti 2018). Furthermore, the simulations by Bao et al. (2025) showed an offset between the brightest TeV excess and the pulsar position. This implies that in extreme cases, the actual extension size of the halo may be difficult to measure because of the filamentary structure that is caused by the asymmetric diffusion of the cosmic rays (Giacinti et al. 2012, 2013). Therefore, further investigations of this region using high-resolution detectors at other wavelengths are essential.

5 Conclusions

We analyzed the γ-ray emission in the vicinity of γ-Cygni using 15 years of Fermi-LAT data and confirmed that the GeV counterpart of the VERITAS source is quite different from other parts inside the SNR. Because the molecular cloud clump in this region spatially coincides with part of the GeV emission and the radio shell, we suggest that the flatter component might be attributed to trapped ions, while the low-energy soft component may be due to escaped cosmic rays that interact with nearby molecular clouds. The GeV spectrum of γ-Cygni therefore consists of a mixture of these two components. The different gas density ratios between the escaped emission zone and the trapped zone result in different γ-ray flux characteristics in different regions. This scenario matches the prediction from Bell (2015) well. In this case, high-energy ions are trapped within a relatively small region around the center of the SNR, where they are confined by an enhanced magnetic field and are unable to escape. The calculated spatial extension of the trapped region agrees well with the observed size of the HAWC source within 1 σ uncertainty. When combined with gas observation results, the resulting γ-ray flux can explain the measured multiwavelength data reasonably well. On the other hand, considering the high spin-down luminosity of the central γ-ray pulsar PSR J2021+4026, we estimated its energy density ε_e using two approaches, as shown by Giacinti et al. (2020). The calculated energy density shown as the squares in Figure 7 lies well within the range of known PWNe, which is located near the boundary of known pulsar halos. This suggests that a pulsar halo origin cannot be entirely ruled out. The spatial correlation between dense gas, the radio SNR shell, and the γ-ray emission provides stronger support for the SNR scenario. Taking into account recent numerical simulation work (López-Coto & Giacinti 2018; Bao et al. 2025), we suggest that it remains plausible that the SNR overlaps with a twisted-morphology pulsar halo.

Acknowledgements

This research made use of the data from the Milky Way Imaging Scroll Painting (MWISP) project, which is a multiline survey in 12CO/13CO/C18O along the northern galactic plane with PMO-13.7 m telescope. We are grateful to all the members of the MWISP working group, particularly the staff members at PMO-13.7 m telescope, for their long-term support. MWISP was sponsored by National Key R&D Program of China with grants 2023 YFA1608000 & 2017YFA0402701 and by CAS Key Research Program of Frontier Sciences with grant QYZDJ-SSW-SLH047. We also would like to thank P.P. Delia, Yang Su for invaluable discussions. This work is supported by the National Natural Science Foundation of China under the grants Nos. 12393853, U1931204, 12103040, 12147208, and 12350610239, the Natural Science Foundation for Young Scholars of Sichuan Province, China (No. 2022NSFSC1808), and the Fundamental Research Funds for the Central Universities (No. 2682022ZTPY013).

References

All Tables

All Figures

	Fig. 1 Folded pulse profile and two-dimensional phaseogram in 32 phase bins derived for PSR J2021+4026. The grayscale represents the weights of photons in each bin. White indicates a higher weight, and black indicates a lower weight. The dotted and dashed lines mark the minimum and maximum phases of the off-pulse phase interval, respectively.
In the text

	Fig. 2 Top row: Fermi-LAT count maps above 10 GeV and 20 GeV within the full-phase range. Bottom row: Fermi-LAT TS maps above 10 GeV during the on- and off-pulse phases, respectively. The black diamond indicates the position of the gamma-ray pulsar PSR J2021+4026 (Abdo et al. 2009a).
In the text

Fig. 3 Fermi-LAT TS maps in the vicinity of γ–Cygni. The energy range for each panel is shown above them. The solid yellow and green circles show the extension size of the best-fit two Gaussian template (R68). The red circles show the two extended TeV sources measured by MAGIC (MAGIC Collaboration 2023), called MAGIC-1 and MAGIC-arc, respectively. The cyan circle and ellipse show the VERITAS results (Aliu et al. 2013; Abeysekara et al. 2018). The orange circle shows the center of the HAWC source with the 1 σ error radius (Abeysekara et al. 2017b). The black diamond shows the position of the gamma-ray pulsar PSR J2021+4026 (Abdo et al. 2009a) as in Figure 2. The white contours come from the CGPS 1420 MHz radio observation of the SNR shell structure (Ladouceur & Pineault 2008).

In the text

Fig. 4 γ-ray SED obtained from the events in 0.1–1000 GeV shown as black data points. The gray solid line shows the best-fit PL spectra for each component. The dashed lines show the 68% confidence bands. The arrows indicate the 95% upper limits when the TS value is lower than 4 in a given bin. The blue data points are the VERITAS results (Abeysekara et al. 2018). The red butterfly was extracted from Abeysekara et al. (2017b). The gray histogram denotes the TS value for each bin.

In the text

	Fig. 5 Integrated 12CO (J = 1–0) emission intensity (K km s−1) toward γ-Cygni region in the velocity ranges [+9.0, +12.0] km s−1 (left panel) and [−16.0, −11.0] km s−1 (right panel) using MWISP data. The solid circles are the same as in Figure 3. The dashed magenta and red circles correspond to different regions that are affected by escaped CRs, called CloudV and CloudH, respectively.
In the text

Fig. 6 Modeling the γ-ray spectrum in the hadronic scenario. Left panel: total contribution (solid black line) from the sum of the dashed red and orange lines, corresponding to the escaped and trapped ions, respectively. The red butterfly was extracted from Abeysekara et al. (2017b). Middle panel: total γ-ray contribution calculated with different rs. Right panel: predicted γ-ray flux in the VERITAS region, obtained by adding the escaped and trapped components. The blue data points and butterfly were taken from Abeysekara et al. (2018).

In the text

Fig. 7 Energy density of TeV sources calculated using two different approaches. The top panels are calculated as ${\epsilon }_{\mathrm{e}}=\dot{E}{\tau }_c/V$, and the bottom panels are calculated following the approach in Section 3.2 of Giacinti et al. (2020). The newly identified pulsar halo associated with PSR J0622+3749 (Aharonian et al. 2021), indicated by a star, is shown along with the Geminga halo and Monogem halo (Abeysekara et al. 2017a). The circles representing several γ-ray PWNe are based on data from Giacinti et al. (2020), and γ-Cygni under different distance assumptions is shown as squares.

In the text