© The Authors 2024

1 Introduction

Supernova remnants (SNRs) are believed to be the most probable locations for accelerating Galactic cosmic rays (CRs) up to the energy of the spectral knee (Ginzburg & Syrovatskii 1964; Hillas 2005). According to diffusive shock acceleration (DSA) theory, self-produced magnetic turbulence may scatter CRs as they are accelerated at the shock of SNRs. Additionally, the highest-energy CRs near the shock precursor have a poor ability to effectively create their own turbulence and tend to escape from the SNR. These escaped CRs might interact with surrounding dense materials like molecular clouds (MCs) through pp (proton-proton) collisions, leading to significant brightness in the γ-ray range. The decay of neutral pions produced by inelastic collisions between accelerated protons and thick gases in MCs is frequently used to explain gamma-ray (γ-ray) emissions in the GeV range from bright SNRs. Notably, peculiar characteristics connected to π⁰ decay can be seen in several γ-ray spectra with the Fermi Large Area Telescope (Fermi-LAT), such as W44 (Uchiyama et al. 2012; Peron et al. 2020) or W28 (Aharonian et al. 2008; Li & Chen 2010; Hanabata et al. 2014), which are thought to be the clearest indication of relativistic nuclei acceleration in SNRs. The flux of γ-rays depends on the quantity of nuclear cosmic rays released and the diffusion coefficient within the interstellar medium (ISM; Aharonian & Atoyan 1996; Gabici et al. 2009; Aharonian et al. 2004; Rodriguez Marrero et al. 2008) and analyses of the spectra of γ-rays from different regions of SNRs offer valuable insights into the diffusion process and even constrain the energy dependence of the diffusion coefficient in the interstellar medium (Aharonian & Atoyan 1996; Ohira et al. 2011).

PeVatrons refer to astrophysical sources that can accelerate CRs to the PeV energy range around the knee of the CR energy distribution. Many ultra-high-energy (UHE) γ-ray sources have been recently detected owing to unprecedented sensitivities of γ-ray observatories such as the Large High Altitude Air Shower Observatory (LHAASO; Cao et al. 2021b) and the High Altitude Water Cherenkov (HAWC) observatory (Abeysekara et al. 2023, 2021). Most of these sources are associated with pulsars and their UHE gamma rays are likely produced by UHE electrons and positrons via the inverse Compton scattering of low energy photons in the background (Cao et al. 2024). Finding hadronic UHE γ-ray sources (PeVatrons) plays an essential role in addressing the origin of the most energetic Galactic CRs. Although the origin of PeV CRs remains obscure, several sources, including the Boomerang SNR (Tibet ASγ Collaboration 2021), the Cygnus region (Abeysekara et al. 2021; Lhaaso Collaboration 2024), and LHAASO J2108+5157 (Cao et al. 2021a), have been considered as important sources of PeV CRs. In particular, LHAASO J2108+5157 has been shown to be associated with molecular clouds, indicating a hadronic scenario for its γ-ray emission (de la Fuente et al. 2023a,b). Nevertheless, the exact nature of their UHE γ-ray emission is still a matter of debate.

G150.3+4.5 was first observed as an SNR candidate with faint radio emission during the Canadian Galactic Plane Survey (CGPS) given by Gerbrandt et al. (2014). Subsequent detailed observation in radio band reported in Gao & Han (2014) shows that the radio spectral indices of the southeastern and western arcs differ significantly with α_e = −0.40 ± 0.17 and α_w = −0.69 ± 0.24, respectively $(S_{\nu }\propto {\nu }^{{\alpha }_{\text{i}}})$. Its overall GeV spectrum is however hard. Recent γ-ray analysis by Devin et al. (2020) suggested a leptonic scenario for the γ-ray emission and the soft point source 4FGL J0426.5+5434 might be associated with a powerful pulsar, even though no obvious pulsation has been detected. In the TeV energy band, 507-day observations by HAWC find no γ-ray emission (Abeysekara et al. 2017; Devin et al. 2020), while in LHAASO’s recent results, a very large-scale extended emission detected by WCDA (Water Cherenkov Detector Array) is spatially coincident with the extended radio/GeV emission, however, in the higher energy band, the source size reported by KM2A (Kilometer Squared Array) is much smaller centered around the point like soft γ-ray source 4FGL J0426.5+5434 (Cao et al. 2024).

In the present work, we carry out a detailed analysis of the γ-ray emission toward G150.3+4.5 region using 14 yr Fermi-LAT data in the energy range of 100 MeV–1 TeV in Sect. 2. In Sect. 3, the CO (J = 1−0) observation results of molecular clouds in this region is presented to explore the origin of the γ-ray emission. In Sect. 4, we discuss the potential interpretation of the γ-ray emission based on multi-wavelength observational constraints. Lastly, Sect. 5 offers our conclusions.

2 Fermi-LAT data reduction

Fermi-LAT is a pair-conversion instrument that is sensitive to γ-rays from 20 MeV to several hundreds of GeV (an overall description can be seen in Atwood et al. (2009). In the following analysis, the newest Pass 8 data are collected from August 4, 2008 (mission elapsed time 239557418) to August 4, 2022 (mission elapsed time 681264005) to study the GeV emission around SNR G150.3+4.5 labelled in Data Release 3 of fourth Fermi-LAT source catalog 4FGL-DR3 (Abdollahi et al. 2020, 2022). We also select the data with “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 exclude the zenith angle greater than 90° to avoid the earth limb contamination. To derive a better point-spread function, the energies of photon are cut between 2 GeV and 1 TeV to do further morphological analysis, and events with an energy between 100 MeV and 1 TeV are selected to do a more detailed spectral analysis of each spatial component. Events within a 14° × 14° region of interest (ROI) centered at the position of G150.3+4.5 are considered for the binned maximum likelihood analysis (Mattox et al. 1996) together with the instrument response functions (IRF) “P8R3_SOURCE_V3”. The Galactic and isotropic diffuse background models (IEMs; gll_iem_ v07.fits)/(iso_P8R3_SOURCE_V3_v1.txt) are adopted, and all sources listed in the 4FGL-DR3 catalog are included in the background model except for those within 7° from the center of ROI, whose normalisations and spectral indices are set free, and the normalisations of galactic and extra-galactic diffuse emission are also set free with the software make4FGLxml.py¹. The maximum likelihood test statistic (TS) is used to estimate the significance of γ-ray sources, where TS = 2(ln 𝓛₁ − ln𝓛₀), and 𝓛₁ and 𝓛₀ are likelihood values for the background with and without target source (null hypothesis), respectively.

2.1 Morphological analyses

The GeV emission in the SNR G150.3+4.5 region is described by an extended source 4FGL J0425.6+5522e plus a point-like source 4FGL J0426.5+5434 in 4FGL-DR3 (Abdollahi et al. 2022). The former is suggested to be associated with SNR G150.3+4.5. To avoid contamination from the point source 4FGL J0426.5+5434, we use the gtfindsrc command to find its best fitting location (RA = 66.585° ± 0.017°, Dec = 54.594° ± 0.015° with the 68% error radii r₆₈ = 0.025°). Here, the 68% error radii correspond to the radii of the 68% confidence level contours for point-like sources. This measurement is almost the same as the position previously recorded in Devin et al. (2020); Abdollahi et al. (2022). To search for the energy dependent morphology of the GeV emission towards the SNR G150.3+4.5 region, all TS maps shown in Fig. 1 are generated by only considering background fitting without including G150.3+4.5. The γ-ray emission in different energy bands are shown in Fig. 1; for convenience, we indicate 4FGL J0426.5+5434 as SrcX. Above 2 GeV (top left), there are highly asymmetric diffuse γ-ray emission with several bright excesses mainly distributed in the north region, and some of them have good spatial correspondence with the radio intensity. Above 10 GeV (bottom left), the emission located in the southern region gradually disappears, while some parts in the north are still bright. Thus we additionally generated the TS map in 2–10 GeV band (top right) that shows a more or less uniform brightness profile. Above 100 GeV (bottom right), the emission located in the southern part disappears completely, while the northern part is still significant.

To further determine the best spatial template of the GeV emission in this region, we tested several spatial templates under the assumption of a power-law spectral type, with both the normalisation and spectral index set free. We first added one point source (model 1) in the strongest excess location and used the gtfindsrc command to find the best location, recorded as (RA = 67.7207° ± 0.03°, Dec. = 55.9180° ± 0.05°, r₆₈ = 0°.11). The TS value derived is only 81 with 4 degrees of freedom. Then, we added another point source (model 2) into another brightest excess location and we repeated the gtfindsrc command to search for the best-fit location, recorded as (RA = 66.4013° ± 0.02°, Dec = 56.0637° ± 0.03°, r₆₈ = 0°.08). The TS value derived is 147 with 8 degrees of freedom. Considering that there are still several multiple point-like excesses emerging in the 2–10 GeV energy band, we also tested the multi-point sources hypothesis (model 3). In this case, we gradually added point-like sources and repeated the gtfindsrc command to search the best location under the assumption of a power-law spectral shape, until all excess peaks could be removed cleanly and the TS value reaches the maximum. Seven point-like sources were found in total. The TS value derived is 296, with 28 degrees of freedom. In the next step, we adopted the Fermipy tool (Wood et al. 2017) to quantitatively evaluate the extension and location of the extended source scenario, like a single uniform disk (model 4) and a 2D Gaussian (model 5) template. Here, r₆₈ is the 68% containment radius and ${\text{r}}_{{68}_{\text{Disk}}}=0.82\sigma ,{\text{r}}_{{68}_{\text{Gaussian}}}=1.51\sigma$ as suggested by Lande et al. (2012) and the best-fit results are 1) RA = 66.6151°± 0.02°, Dec = 55.3755°± 0.02°, ${\text{r}}_{{68}_{\text{Disk}}}={1.255}^{\circ }\pm {0.03}^{\circ }$, σ = 1.530° ± 0.04° and 2) RA = 66.4258° ± 0.05°, Dec = 55.1577°± 0.04°, ${\text{r}}_{{68}_{\text{Gaussian}}}={1.377}^{\circ }\pm {0.05}^{\circ }$, σ= 0.912° ± 0.03° similar to the results from Devin et al. (2020). In addition, we produced a series of ring templates (model 6) at the center of the uniform disk template with various inner radius from 0º.1 to 0º.8 in steps of 0º.05. The outer radius is fixed to the uniform disk template. We found that the TS values and likelihood values did not increase significantly, and the highest TS value case is shown in Table 1. To search for the spatial correlation between GeV morphology and radio intensity, we also test Urumqi λ6cm map (model 7) as a template. However, its TS value is still low. To further substantiate the potential spectral variations between the northern and the southern region depicted in Fig. 1, we divide the uniform disk template into two half-disks (model 8) along a western and eastern line in celestial coordinates. Then we further conducted additional tests on two half-disk templates with different slope lines, which are not horizontal but represent +30 and −30 degrees (model 9 and model 10), respectively. However, both of them had lower TS value performances compared with model 8. To further determine whether the extended GeV emission is correlated with gas distributions, we tested a spatial template for H₂, which is produced from the carbon monoxide (CO) composite survey described in Sect. 3. Similarly, we also divided the potential H₂ template along the western and eastern line in celestial coordinates, recorded as model 11 and model 12, corresponding to different gas velocity intervals. It is important to note that, in the last two models, the gas distribution template was only applied to the southern half-disk region since molecular clouds are only distributed in the south (shown as Fig. 5, details are discussed in the following Sect. 3), the northern region is still described as a uniform half-disk template as in model 8.

The TS_ext is defined as TS_ext=2(ln 𝓛_ext − ln 𝓛_ps), where 𝓛_ext/𝓛_ps represents the maximum likelihood value for the extend or point-like template, considering only 1 free parameter added by the extended template, and the extension significance is approximately given by $\sqrt{{\text{TS}}_{\text{ext}}}\sigma$. When TS_ext > 16 (Lande et al. 2012), the extended source hypothesis is valid. In model 4 (disk), TS_ext is calculated as 542, which rejects the point-like source hypothesis at 23.3σ. Model 5 (2D Gaussian) shown a marginal improvement, TS_ext = 563 (correspond to 23.7σ). While the results from model 6/7 (ring and radio intensity map) are also similar, and the TS values are less than previous two templates. The two half-disk template (model 8) shows a better performance both in TS (marginal improvement ~2.8σ) and AIC (Akaike information criterion) value, which is defined as AIC = 2k − 2 ln 𝓛 and describe in Akaike (1974), where k is the number of degrees of freedom of the model and 𝓛 is the likelihood value. Additionally, the gas distribution template in different velocity ranges (model 11 & model 12) have similar improvements ~3.7σ. Model 11 has the lowest value of AIC and the highest TS value, and is selected to carry out the subsequent spectrum analysis. We also analysed the energy spectrum of the disk template (model 4) for comparison. Hereafter, we refer to the northern or southern lobes derived from the H₂ template (model 11) as northlobe and southlobe, respectively, and the residual maps obtained after subtracting either the southern, northern, or both lobes independently are shown in the Fig. 2. The above morphology analysis results are summarised in Table 1. Furthermore, we tested the energy-dependent extension of SNR G150.3+4.5 in each energy band, and the results suggest that the extension size does not shrink visibly. No significant energy-dependent morphology is found. The results are summarised in Table 2.

Fig. 1 TS maps in the vicinity of SNR G150.3+4.5 observed with the Fermi-LAT. All figures are generated with square bins of 0.02°. The energy ranges are indicated above the panels. The green circle represents the best-fit 2D Gaussian template from Devin et al. (2020). The white solid and dashed circles show the R68 and radius of the uniform disk model in this work, which is separated by the white horizontal line into two half-spheres: northlobe and southlobe (models 8 & 11 & 12). The white dashed line corresponds to the +30° and −30° scenarios (models 9 and 10). The unnamed white cross in the upper right corner represents a background point source from 4FGL-DR3 (Abdollahi et al. 2020, 2022) and SrcX represents the point source 4FGL J0426.5+5434, which has been treated as a distinct background point source. The cyan contours are extracted from Urumqi λ6 cm results (Gao & Han 2014), showing the complete shell structure. The red and orange circles are R39(σ) size for the LHAASO KM2A and WCDA sources (Cao et al. 2024), respectively.

Fig. 2 Residual maps above 2 GeV, obtained by subtracting southern/northern/both lobe, all labels match those in Fig. 1.

2.2 Energy spectrum

Given the best-fit spatial template model 11 determined above, we performed γ-ray spectral analyses for each components using larger amount of events from 100 MeV to 1 TeV. Considering prior results shown in Devin et al. (2020), we adopted log-parabola (LogPb; $\text{dN}/\text{dE}\propto {\text{E}}^{-(\alpha +\beta \text{log}(\text{E}/{\text{Eb}}_{}))}$) spectra directly for the disk (model 4) and SrcX. For the northlobe and southlobe (model 11), we tested different spectral types and the spectral indices show a significant difference among both the simple power law (PL; dN/dE ∝ E^−α) and broken power law (BPL²) assumptions, as shown in Table 3. In this table, we calculate the estimated systematic errors due to the Galactic diffuse emission by changing the normalisation of the best-fit Galactic diffuse model artificially by ±6% as suggested in Abdo et al. (2010). The sums of statistical and systematic errors are calculated by $\sigma =\sqrt{{\sigma }_{\text{stat}}^2+{\sigma }_{\text{sys}}^2}$ which are shown as the red error bars in Fig. 3. For the northlobe, compared with a single PL model, the BPL model has a notable improvement of the fitting quality, which could be quantify as TS_curve, defined as TS_curve=2(ln 𝓛_BPL – ln𝓛_PL) (Abdollahi et al. 2020). The obtained value of 24 corresponds to a significance level of ~4.9 σ with only one additional free parameter. Thus we conclude that there is an energy spectral break at ~ 146 GeV in the northlobe spectrum. For the southlobe, there is almost no variation in the TS value between PL and BPL cases, thus the southlobe can be described by a PL. Furthermore, to derive the spectral energy distributions (SEDs) of both sources, we divided the events in the 100 MeV–1 TeV energy range into twelve logarithmically equal intervals and repeated the same likelihood fitting analysis for each interval. The normalisations of all sources are left free, and the spectral indices and energy break are fixed to their best-fit value. For bins with TS values less than 5.0, we give upper limits calculated with 95% confidence level using a Bayesian method (Helene 1983). The SED of SrcX exhibits a break around 1 GeV although the flux is insignificant above 10 GeV. Thus the γ-ray excess from SrcX is mainly in the lower GeV band implying a distinct spectral component for the extended TeV emission detected by LHAASO-KM2A. To further test the flux stability for SrcX, we analysed its light curve (LC) in the following part.

Fig. 3 SED and TS values (shaded areas) of the southlobe/SrcX (panel left), northlobe (panel middle), and disk (panel right) as measured by Fermi-LAT in the energy range from 0.1 to 1000 GeV. Upper limit are calculated at 95% confidence level using a Bayesian method for points with TS values less than 5. The statistical errors are shown in black/blue error bars, while the sum of statistical and systematic errors calculated by $\sigma =\sqrt{{\sigma }_{\text{stat}}^2+{\sigma }_{\text{sys}}^2}$ are marked by red error bars. The disk and SrcX can be described with a log-parabola spectrum. The northlobe is best-fitted with a broken power-law spectrum. The southlobe is compatible with a single power law and the black dash lines show the 1 σ statistic error uncertainty.

2.3 Flux variability of SrcX

Considering the complex environment around SrcX, some potential unidentified counterpart that might contributed to γ-ray excess in its 95% error circle through the SIMBAD³ database, such as Galaxy LEDA 2471506 (Paturel et al. 2003), Nebula DSH J0426.0+5433 (Kronberger et al. 2006) and variable star ATO J066.5105+54.5272 detected by Gaia (Gaia Collaboration 2020), which leads its puzzling γ-ray excess, thus it is essential to explore its flux variability. The 14 yr data set, with energies between 0.1–10 GeV, was divided into seven time bins, respectively. For the timing analysis, all spectral indices of the sources included in the model were fixed to the best-fit values and only normalisations are allowed to vary. Then we repeated the likelihood analysis in each time bin. For the time bin with a TS value smaller than 9, an upper limit with 95% confidence level was calculated. The derived light curve are shown in Fig. 4. We also calculated the variability index as defined in Nolan et al. (2012), which under the assumption that with a value in the null hypothesis where the source flux is constant across the full time period and the value under the alternate hypothesis where the flux in each bin is optimised as ${\text{TS}}_{\mathrm{var}}={\displaystyle {\sum }_{i=1}^{N}2}\times (\text{Log}({𝓛}_{\text{i}}({ℱ}_{\text{i}}))-\text{Log}({𝓛}_{\text{i}}({ℱ}_{\text{mean}})))$, where 𝓛_i is the likelihood corresponding to bin i, F_i is the best-fit flux for bin i, and 𝓕_mean is the best-fit flux for the full time assuming a constant flux. For this LC, the source of TS_var ≥ 16.81 was considered to be a variable source with a 99% confidence level, given the TS_var is anticipated to be distributed as ${\chi }_{N-1}^2({\text{TS}}_{\text{var}})$ in the null case, we get TS_var = 6.8 for the flux variability using above energy band, which is very marginal, thus we suggest no significant variability for the γ-ray excess from SrcX. Nevertheless, we notice that the TS values in last two time bins has a insignificant decline. Considering that the data points still locate in 1 σ uncertainty region, it is possible to conclude that the photon flux is stable and the evidence is too weak to draw a conclusion about flux variability, which could be interpreted as statistical fluctuations.

Fig. 4 Light curves of SrcX from 0.1 to 10 GeV. The mean flux and its 1 σ statistic uncertainty are shown as the orange horizontal solid and dashed lines, respectively. The grey histogram represents the TS value for each time bin.

3 CO observations

To track the distribution of H₂ surrounding the SNR G150.3+4.5 region, we used the CO composite survey (Dame et al. 2001) and the standard assumption of a linear relation between the velocity-integrated brightness temperature of CO (J = 1–0) line at 115 GHz (2.6mm) and the molecular hydrogen column density. The HI spectrum obtained from the Leiden/Argentine/Bonn (LAB) survey by Kalberla et al. (2005) reveals the presence of two velocity components within the range of [−9.8, +2.4] km s⁻¹. This observation suggests that clouds in this study may have been disturbed by shocks, as discussed by Castelletti et al. (2013) and Su et al. (2017), which is further confirmed by recent PMO 13.7 m millimeter telescope (Feng et al. 2024), indicating that the signature of SNR associated with surrounding MCs is present in this region. The more precise velocity range ([−9.8, +2.4] km s⁻¹) was considered as whole velocity range for analysis, including the velocity ranges [−9.8, −3.8] km s⁻¹ discussed by Devin et al. (2020) for comparison and depicted in the left/right panel of Fig. 5. The whole velocity interval can be translated into a kinematic distance of $d={0.55}_{-0.55}^{+0.65}\text{ kpc}$ (Devin et al. 2020). This finding is consistent with a kinematic distance of d = 0.8 kpc and an age of 13 kyr reported by Zeng et al. (2023). Given that the SNR is located inside the distance error range of MCs and their potential interactions, here we assume that the molecular clouds are associated with SNR G150.3+4.5 and we adopt the distance d = 0.8 kpc in further calculations. Then θ= 1°.53 corresponds to a physical radius of 21.4 pc. In order to evaluate the column density of H₂ in this region, a conversion factor X_CO = 2 × 10²⁰cm⁻² K⁻¹ km⁻¹ s is used (Bolatto et al. 2013; Dame et al. 2001), and $N_{{\text{H}}_2}$ can be derived as $N_{{\text{H}}_2}=X_{\text{CO}}\times W_{\text{CO}}$· Thus, the mass of the molecular complex can be calculated from W_CO as $$M=\mu m_{\text{H}}{D}^2\text{Δ}{\text{Ω}}_{\text{px}}{X}_{\text{CO}}{\displaystyle \sum _{\text{px}}{W}_{\text{CO}}}\propto N_{{\text{H}}_2},$$(1)

where µ is equal to 2.8 if a relative helium abundance of 25% is assumed, m_H is the mass of the H nucleon, $N_{\text{H}}=2N_{{\text{H}}_2}$ represents the column density of the hydrogen atom in each pixel, and ∆Ω_px corresponds to the solid angle subtended for each pixel in the map (square binning of 0.125° per side). The term Σ_pxW_CO takes into account the binning in velocity of the data cube and is obtained by summing the map content for the pixels in the desired sky region and desired velocity range and scaled by the bin size in velocity. Then, we obtained the mass of the molecular cloud in the whole velocity range through its relation with $N_{{\text{H}}_2}$, assuming a spherical geometry for the gas distribution, the volume is given by ${\text{V}}_{\text{Total}}=\frac{4\pi }{3}{\text{R}}^3$ and V_{Northlobe/southlobe} = 0.5 × V_Total (here R = d × θ), while the mass of gas is calculated to be ${\text{M}}_{\text{Southlobe}}=3.3\times {10}^4{\text{d}}_{0.8}^2{\text{\hspace{0.17em}M}}_{\odot }$, and the average H₂ cubic density in this region is about ${\text{n}}_{\text{Southlobe}}{\text{=66d}}_{0.8}^{-1} {\text{cm}}^{-3}$. For the northlobe, the mass of gas is estimated to be about ${\text{M}}_{\text{Northlobe}}\text{=}4.7\times {10}^3{\text{d}}_{0.8}^2{\text{\hspace{0.17em}M}}_{\odot }$, the average H₂ cubic density is about ${\text{n}}_{\text{Northlobe}}{\text{=9d}}_{0.8}^{-1}{\text{\hspace{0.17em}cm}}^{-3}$. For the whole-disk region, the total mass of gas within 1.53° is estimated to be ${\text{M}}_{\text{Total}}=3.8\times {10}^4{\text{d}}_{0.8}^2{\text{\hspace{0.17em}M}}_{\odot }$, average density is approximately ${\text{n}}_{\text{Total}} \text{=} {\text{37d}}_{0.8}^{-1} {\text{cm}}^{-3}$, which is roughly consistent with the measurements from Zeng et al. (2023) and Feng et al. (2024); however, it differs noticeably from the 10⁻³ order estimated from X-ray measurements with ROSAT in Devin et al. (2020). We suggest that the actual ambient density of about several tens cm⁻³ derived from gas detector measurements will be more accurate. For the SrcX region, we estimate that the mass of gas within 0. 1° around the central location of SrcX is calculated to be ${\text{M}}_{\text{SrcX}}\text{=2}.5\times {10}^2{\text{d}}_{0.8}^2 {\text{M}}_{\odot }$, corresponding to ${\text{n}}_{\text{SrcX}}{=89\text{d}}_{0.8}^{-1}{\text{\hspace{0.17em}cm}}^{-3}$, which indicates that there are some higher-density area in the southern region, but due to the resolution limit (0.125°) of CO composite survey data, we are currently unable to obtain the significantly fiber structure of gas distribution, and more high-resolution telescopes may be needed to observe this area.

As shown in Fig. 5, there is good spatial correspondence between the gas distribution and the location of γ-ray emission in southern region, which might be the main contributor to southlobe. While in the northern region, there is almost no distribution of molecular cloud. To further investigate the origins of γ-ray emission in each part, we model the γ-ray fluxes in the next section.

Fig. 5 Integrated CO (J = 1–0) emission intensity (Kkm s−1) toward SNR G150.3+4.5 at 115 GHz using CO composite survey (Dame et al. 2001). Panel left: integrated on the whole velocity range [−9.8, +2.4] km s−1, and overlaid with the radio shell structure extracted from Urumqi λ6cm results (Gao & Han 2014). Panel right: integrated on [−9.8, −3.8] km s−1 while overlaid with the GeV emission contour line generated from the top right panel in Fig. 1. Other labels are same as above.

4 A two-zone model

4.1 Hadronic emission dominating in the southern region

The above Fermi-LAT data analysis shows that the GeV γ-ray emission in the southern region consists of two separate sources: the extended southlobe and the point-like SrcX. Among them, southlobe has a flat spectrum with an index of ~2.0. Considering that there is a good spatial correspondence between the southern γ-ray emission and molecular cloud distribution shown in Fig. 5, the southlobe could originate from interaction of accelerated CRs with MCs. A number of arguments suggest that SrcX can be interpreted as molecular clouds illuminated by escaping CRs with their projection located within the shell of the SNR (Li et al. 2023a). First, for SrcX, there is no detected X-ray emission and no pulse signal (Devin et al. 2020). Moreover, it is in spatial coincidence with molecular clouds and the high density environment calculated in Sect. 3. Finally, the fact that its GeV spectrum index is much softer than others region from the SNR itself and that its flux is so stable (as shown in Sect. 2.3) further support this scenario. The parameters adopted for this scenario are summarised in Table 4.

To fit the γ-ray spectra of the southlobe and SrcX. We assume an instantaneous injection of protons into an uniform emission zone at T = 13 kyr ago. The injected proton spectrum is a broken power law: $$Q(E)=Q_0\frac{{(E/E_{p,\text{br}})}^{-{\gamma }_1}}{1+{(E/E_{p,\text{br}})}^{{\gamma }_2-{\gamma }_1}}.$$(2)

Here, γ₁ and γ₂ are the spectral indices below and above the break energy E_p,br, respectively. A part of the injected protons will be trapped inside while others can escape via diffusion. Such a broken power-law spectrum is often seen in γ-ray SNRs (Zeng et al. 2017, 2019, 2021; He et al. 2022; Li et al. 2023a). Given the flat γ-ray spectrum in GeV energy band for southlobe, we have γ₁ = 2.0 and γ₂ = γ₁ + 1, the break energy of protons E_p,br is fixed at 10 TeV.

Assuming that the shock-cloud interaction site is far away from the shock front and CRs is released at the time of the supernova explosion, the distribution of the escaped protons follows (Thoudam & Hörandel 2012; Liu et al. 2020; Li et al. 2023b) $$N_p(E,r_{\text{s}},T)=\frac{Q(E)}{{[4\pi D(E)T]}^{\frac{3}{2}}}\exp \left[\frac{-r_{\text{s}}^2}{4D(E)T}\right],$$(3)

where the diffusion coefficient is assumed to be uniform and taken to be D(E) = χD₀(E/E₀)^δ for E > E₀, where D₀ = 1 × 10²⁸ cm² s⁻¹ at E₀ = 10 GeV and δ = 1. With a distance between the Earth and SNR G150.3+4.5 d = 0.8 kpc, the γ-ray SNR has a radius of R = θ × d = 21.4 pc, the same as in Sect. 3. Considering the physical distance between SrcX and SNR central cannot be constrained by observations for the projection effect, we assume r_s as the distance between the injection site (SNR centre) and the emitting molecular clouds in the subsequent model building part to fit the SED of SrcX. For an injected source spectrum given by Q(E) ∝ E^−Γ and D(E) ∝ E^δ, Eq. (3) shows that N_p(E) will follow $N_p(E)\propto E^{-\left(\text{Γ}+\frac{3}{2}\delta \right)}$ at high energies, while displaying a lower energy spectral cut-off at E_b where $\sqrt{4D(E_b)T}\simeq r_{\text{s}}$. The bottom right panel of Fig. 6 shows the dependence of the SED on the SNR age and r_s. The ~21 pc (magnetic solid line) scenario is perfectly compatible with the SED trend, while under the larger r_s assumption, the peak values of the fit lines are difficult to match with the observed data points (even with a higher injection energy). Thus we propose that there is such a constraint on the approximate value of r_s.

The SED trend is perfectly compatible with the 21 pc (magnetic solid line) scenario, but the peak values of the fit lines under the bigger r_s assumption are difficult to match the observed data (even with a greater inject energy). For this reason, we propose that there is a constraint regarding the approximate value of Rs.

The total energy of injected protons is assumed to be W_inj = η E_SN, the ratio η of the kinetic energy converted into accelerated protons with a standard value of 0.1, the E_SN is the kinetic energy of the SNR with a typical value of 10⁵¹erg. The correction factor χ of the diffusion coefficient is set to be free parameter. The corresponding γ-ray fluxes are calculated using the naima package (Zabalza 2015) with $$\frac{\text{d}{N}_{\gamma }}{\text{d}{E}_{\gamma }}=\frac{n_i{v}_{i}c}{4\pi d^2}{\displaystyle \int \frac{\text{d}{\sigma }_{pp}}{\text{d}{E}_{\gamma }}}\left(E_{\gamma },E\right)N_{\text{p}}(E,r)\text{d}E,$$(4)

where the differential proton-proton inelastic cross section for γ-ray production, dσ_pp/dE_γ, is adopted from Kafexhiu et al. (2014). Also, ni/vi represents for the total gas density/volume in different emission zones (i = southlobe, X) and are the same as those given in Sect. 3. For the southern region, the resulting γ-ray flux with the parameter χ = 0.3 could perfectly explain the observational data for T = 13 kyr as indicated by the magenta line in the top left panel of Fig. 6. And the total energy of escaped protons above 1 GeV in the X emission zone is calculated to be W_escaped,X = 2.21 × 10⁴⁷(n_X/89cm⁻³)⁻¹ erg, while the total energy of trapped protons above 1 GeV in the southlobe emission zone is calculated to be W_tapped = 2.35 × 10⁴⁷(n_S/66 cm⁻³)⁻¹ erg. We note that the diffusion coefficient is one order of magnitude lower than the standard Galactic value, and a smaller diffusion coefficient (e.g. χ = 0.1) would leads to a lower value of W_escaped,X to explain the observed flux. However, the soft GeV spectrum requires a value of δ as high as 1, corresponding to Bohm-like diffusion, which is different from typical turbulence models such as Kolmogorov and Kraichnan found in others SNRs. The presence of Bohm diffusion around this SNR could be due to CR-driven instabilities from CRs leaving that source and modifying the magnetic field in the ISM around the SNR (Bell 2004). It is known from numerical simulations (Reville & Bell 2013) that the magnetic field in the upstream of SNR shocks leads to Bohm diffusion for CRs. It is therefore plausible that CRs could undergo Bohm-like diffusion around their SNR. A similar diffusion coefficient was also found in HESS J1912+101 (Li et al. 2023a).

Fig. 6 Modelling of the multi-wavelength spectra of the SrcX southlobe (top-left), northlobe (top-right), and disk (bottom-left). The blue and purple butterfly shadow regions come from recent LHAASO results (Cao et al. 2024). The orange data points come from 4FGL J0426.5+5434 as measured by Devin et al. (2020). Top left: Solid, dashed, and dash-dotted magenta lines show the spectra produced by escaping ions with different values of χ indicated, and the black solid line shows the hadronic model for trapped protons with a broken power-law distribution. Top right: Leptonic models, where the power law (ECPL) and broken power law (ECBPL) are presented by red dashed and grey solid lines, while the orange and green dashed lines show a special ECPL and ECBPL case where the cut-off energy is fixed to 25 TeV. The radio fluxes and the X-ray upper-limit come from Zeng et al. (2021) and Devin et al. (2020), respectively. Bottom left: Hybrid model for the total γ-ray flux shown as black solid, red dashed, orange dashed lines added by different leptonic spectral types. Blue data points are derived from the best-fit results measured by Devin et al. (2020) for comparison. Bottom right: Dependence of the SED on the SNR age and rs. The solid, dotted, dashed, dash-dotted lines correspond to rs = 21.4/35/45/55 pc, respectively.

4.2 Leptonic emission dominating in the northern region

On the inner side of the shock front, the γ -ray emission is assumed to be mainly from the shock-accelerated electrons, which are also producing the synchrotron emission in the spectral range from radio to X-rays, mainly owing to the consistency with the hard GeV γ-ray spectrum similar to typical PWNe (pulsar wind nebulae) such as HESS J1825-137 (Principe et al. 2020), HESS J1640–165 (Mares et al. 2021), and so on. The radio spectral index −0.69 ≤ α ≤ −0.40 (Gao & Han 2014) in the SNR G150 region is much softer than that of a typical PWN (−0.3 ≲ α ≲ 0, Gaensler & Slane 2006) and is more compatible with a SNR scenario (α ≲ −0.40, Green 2017). Moreover, there is no pulsation (Barr et al. 2013), nor X-ray emission detected (Devin et al. 2020), and the radio shell structure shown in Fig. 1 is very symmetric. Therefore, we suggest that the γ-ray emission in the northlobe also comes from SNR G150 and is dominated by the inverse-Compton scattering process. In this case, we performed a multi-wavelength modelling using the naima package (Zabalza 2015) with a simple one-zone leptonic model for the northlobe with a Markov chain Monte Carlo algorithm to give the numerical constraint. The best-fit parameters are summarised in Table 5. On the other hand, we note the potential irregular and elongated structure in the northern region, indicating the possibility of (filamentary) diffusion caused by local magnetic turbulence (Giacinti et al. 2012); thus, the high-energy electrons escape scenario cannot be ruled out. Combined with the age of SNR G150.3+4.5 in the ranges 10– 100 kyr, there is an accordance with the potential pulsar halo evolution stage 2 that is depicted in Giacinti et al. (2020), which means that under this assumption, the potential pulsar halo/PWN has not yet fully escaped from the SNR region. Thus, the γ-ray emission might be the result of a SNR plus PWN contribution (Zeng et al. 2023). Such a scenario cannot be totally ruled out and more high-resolution observations are needed to reveal the nature in this region.

Some evidence was shown that the radio continuum spectral indices for the northern and southern parts of G150.3+4.5 might be different (Gao & Han 2014). We therefore treat the two parts separately. According to the procedure introduced in Gao et al. (2011), the large-scale Galactic diffuse emission in the observation image is filtered out by using the technique of ‘background filtering’ (Sofue & Reich 1979). For a further adjustment, a twisted hyper-plane is constructed based on the pixel values in the four image corners without any obvious structures and subtracted to remove the local-background Galactic diffuse emission, which ensures that the radio continuum emission seen in the target area all comes from G150.3+4.5. Finally, we estimate the flux densities for the two parts which are separated by the white line as marked in Fig. 1. Firstly, the above procedure is applied to the Effelsberg 21 cm image (Reich et al. 1997) as shown in the Fig. 1 of Gao & Han (2014). The results are S_{21 cm} = 5.1 Jy and S_{21 cm} = 6.3 Jy for the northern and southern parts of G150.3+4.5, respectively. The flux density at 21 cm for the northern part (~5.1 Jy) could be regarded as a lower limit, since some of the faint emission seen at 6 cm in this part is not detected at 21 cm. Then the same routine is further applied to Urumqi 6 cm data, and we obtain S_{6 cm} = 3.3 Jy, and S_{6 cm} = 3.5 Jy accordingly. Considering the uncertainty, i.e. 10% in the measurements (Gao & Han 2014), the radio flux densities given by Urumqi 6 cm in both parts are almost equal. Thus, we adopted 0.5 times radio flux for the north lobe model building calculation shown in Fig. 6 (top-right panel).

In the leptonic model, the inverse Compton scattering (ICS) or bremsstrahlung processes of relativistic electrons is considered, especially for the ICS process, the radiation field includes the CMB (cosmic microwave background) and IR component from interstellar dust and gas with T = 30 K and u = 1 eV cm^-3 (Porter et al. 2006, 2008). The distance is also adopted as 0.8 kpc. The spectra of the electrons was assumed to be a power law with an exponential cut-off (ECPL) taking the following form: $$\frac{\text{d}{N}_e}{\text{d}E}\propto {\left(\frac{E}{E_0}\right)}^{-{\alpha }_{e1}}\exp \left(-\frac{E}{E_{\text{e,cut}}}\right)$$(5)

and a broken power law (ECBPL) with an exponential cut-off that follows (Zeng et al. 2019): $$\frac{\text{d}{N}_e}{\text{d}E}\propto \exp \left(-\frac{E}{E_{\text{e,cut}}}\right)\left\{\begin{array}{ll}{\left(\frac{E}{E_0}\right)}^{-{\alpha }_{e1}}\hfill & ;E<E_{\text{e,break}}\hfill \\ {\left(\frac{E_{\text{e,break}}}{E_0}\right)}^{{\alpha }_{e2}-{\alpha }_{e1}}{\left(\frac{E}{E_0}\right)}^{-{\alpha }_{e2}}\hfill & E\ge E_{\text{e,break}}\hfill \end{array}\right\}$$(6)

where α_e and E_e,cut are the spectral index and the cut-off energy, E_e,break represents the break energy, and α_e2 = α_e1 + 1.0. Then we set the synchrotron cooling timescale being equal to the age of SNR G150.3+4.5 to reduce the number of free parameters, at the cut-off energy of E_e,cut = 1.25 × 10⁷ ${\text{t}}_{\text{age,yr}}^{-1}{\text{\hspace{0.17em}B}}_{\text{μG}}^{-2}$. We combined it with the radio flux data points and the upper limit in X-ray to constrain the total energy of electrons above 1 GeV and magnetic field strength, which is calculated as W_e = 1.73 × 10⁴⁷ erg and B ≃ 7.0 µG in broken power law case, the corresponding E_e,cut is calculated as 19.6 TeV. For the power law case (red dashed line), we derived similar results, namely: We = 1.75 × 10⁴⁷ erg and B ≃ 7.6 µG, while E_e,cut is 16.5 TeV. For the PL (BPL) cases with a fixed cutoff energy equal to 25 TeV, W_e = 2.02 (1.74) ×10⁴⁷ erg and the magnetic strength is calculated as B ≃ 6.2 µG. Based on the above best-fit parameters, we derived a similar magnetic strength value compared with Devin et al. (2020). However, due to the ambient gas densities being calculated from different surveys, the estimated age of the SNR varies, and we suggest that using gas detector measurements will provide a more accurate density estimate than X-ray measurements. Since the best-fit results for freshly accelerated electrons and protons is measured in each lobe, the electron-to-proton ratio can be calculated as K_ep = 6 × 10⁻³ at 1 GeV in the whole SNR region. This is consistent with the expectation in the range of 10⁻⁵−10⁻² (Cristofari et al. 2013, 2021), while the typical value of K_ep = 10⁻² was derived from the CRs measured at Earth (Zhang & Chen 2016). Also, a lower limit of 10⁻³ was obtained according to the radio observations of the SNRs in nearby galaxies (Katz & Waxman 2008).

Then we added the γ-ray contributions from the southlobe and the northlobe together to derive the total γ-ray flux in the SNR G150.3+4.5 region, which is shown in the bottom left panel of Fig. 6. The obtained total γ-ray flux can well explain the apparent discontinuity at the spectral boundary between WCDA and KM2A given by LHAASO (Cao et al. 2024). Notably, the red dashed line, with a lower cut-off energy value, is perfectly consistent with both KM2A and WCDA measurements. In contrast, the black solid line and the orange and green dashed lines, which have higher cut-off energies, do not fit well with the WCDA/KM2A results. This means that the γ-ray emission in the lower TeV energy band detected by WCDA jointly comes from the hadronic and leptonic components; as the energy increases, the leptonic component quickly disappears. In the higher TeV energy band, only the hadronic component remains in the high density region. This leads to the softer spectrum detected by KM2A, as well as the smaller extended source, which can smoothly connect with the SED from the southlobe. This evidence suggest that SrcX is not the GeV counterpart of the KM2A source, even though they have similar positions. Correspondingly, the KM2A source should be associated with southlobe. We also noticed that the disappearance of the leptonic components is limited by multi-wavelength data. The leptonic model can match with LHAASO observations only when the cut-off energy is below the KM2A energy band (< 25 TeV). Therefore, this conclusion is to some extent constrained by observations. Conducting further theoretical research about the mechanism of electron cooling in this region is worthwhile.

5 Discussion and conclusions

We analysed the MeV-GeV γ-ray emission in the vicinity of SNR G150.3+4.5, using 14 yr of Fermi-LAT data, and we found that the northern half-sphere has a spectral break near 146 GeV with a significance of TS ~24 (4.9 σ). The photon indices below and above the energy break are 1.54 ± 0.04_stat ± 0.07_syst and 2.28 ± 0.08_stat ± 0.12_syst, respectively. The southern half-sphere could be described well by a single power law with an index of 1.95 ± 0.07_stat ± 0.09_syst even though emission above 100 GeV is insignificant. Molecular clouds are only found in the southern half-sphere, implying soft GeV spectra in high density regions, a feature that has been observed in RX J1713.7-3946 (Sano et al. 2013; Tanaka et al. 2020), SNR Puppis A (Xin et al. 2017), and (more recently) HESS J1912+101 (Li et al. 2023a). In light of the discovery of spectral evolution of high energy particles in SNRs (Zeng et al. 2017, 2019), these results suggest that shocks evolving in a high density environment produce high energy particle distribution similar to those in old SNRs.

The northern half-sphere has a soft radio spectrum (Gao & Han 2014) and the multi-wavelength SED can be fitted with the leptonic scenario for the γ-ray emission, under a ECPL/ECBPL spectral assumption, the total energy of electrons above 1 GeV is calculated as 1.75/1.73 × 10⁴⁷ erg, respectively; both of them reveal a high energy cut-off around several tens of TeV, which is consistent with measurements towards the whole SNR region (Zeng et al. 2023). The bright Fermi source SrcX (4FGL J0426.5+5434) has a soft spectrum and is spatially coincident with a puzzling source discovered by LHAASO-KM2A and a bright radio arc with a hard spectrum in the southern half-sphere (Gerbrandt et al. 2014; Gao & Han 2014). It can be attributed to molecular clouds illuminated by CRs injected at the birth of the SNR, and its SED can be fitted by a model with the typical SNR kinetic energy E_SN = 10⁵¹erg and energy convert fraction η = 0.1, while the best-fit diffusion coefficient is about one order of magnitude lower than the standard Galactic value (χ = 0.3). We found that the multi-wavelength SED of the southern half-sphere can be fitted with the hadronic scenario for the γ-ray emission with the emitting ions having a broken power-law spectrum, reminiscence spectra of middle age SNRs. The γ-ray emission detected by the LHAASO-KM2A is then produced by CRs interacting with ions in the background. Although the soft spectrum of SrcX may be attributed to CRs escaping from the SNR and interacting with a dense MCs, it may also be associated with a soft spectrum of CRs accelerated as the shock of SNR G150.3+4.5 encounters a dense MCs and slows down (Li et al. 2023a). Further observations are needed to reveal of the nature of SrcX.

Overall, we found that the GeV spectral indices differ significantly between the north and south sides of the SNR G150.3 + 4.5. Combining the spatial correspondence between the distribution of MCs and the γ-ray emission in southern region, we suggest that the softer spectral index in the high-density southern region originates from hadronic emission processes. The puzzling point-like source SrcX can be explained by an escape model. In the northern region, the harder spectrum is dominated by the leptonic inverse Compton scattering process. Combining multi-wavelength data, we constrained the fitting parameters for each source in this region. The fitting results show that the hybrid model with both leptonic and hadronic contributions to the γ-rays is more reasonable than a purely leptonic model. We suggest that the hybrid γ-ray emission is mainly due to evolution of the SNR in an inhomogeneous environment. Shocks sweeping through the low-density northern region lead to a harder spectral index dominated by the leptonic emission, while the high-density southern region has a softer spectral index dominated by the hadronic component. In higher energy bands probed by LHAASO-KM2A, only high-density regions bombarded by high energy CRs can be observed.

Acknowledgements

We thank Shaoqiang Xi for the help on the molecular cloud data analysis approach. We thank Xuyang Gao for the help on the radio continuum data based on the Urumqi and Effelsberg observations. We also would like to thank Houdun Zeng, Yang Su, Xi Liu and P.P.Delia for invaluable discussions. This work is supported by the National Natural Science Foundation of China under the grants No. 12393853, 12375103, 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).

Appendix A Possibility of other potential theoretical models for multi-wavelength spectra

In contrast to the preceding discussion, we calculated the hadronic contribution in the northern low gas density region (n_Northlobe $=9{\text{d}}_{0.8}^{-1}{\text{\hspace{0.17em}cm}}^{-3})$; also, the contribution of lepton components in the southern region. In this scenario, both the northern and southern regions are dominated by the hybrid scenario. The noteworthy point is that due to the existence of dense molecular clouds interacting with SNRs in the southern region, the magnetic field strength in this area could be significantly amplified around the clump (Zirakashvili & Aharonian 2010; Inoue et al. 2012) because of both the field compression and the magnetohydrodynamical (MHD) instabilities that progress in the shock-clump interaction (Celli et al. 2019). For the sake of simplicity and clarity in the model, here we adopt the electron spectrum with the format of ECPL in Table 5, doubling the magnetic field strength (~ 15.24 µG), then we set the synchrotron cooling timescale equaling to the age of SNR G150.3+4.5 (E_e,cut = 1.25 × 10⁷ ${\text{t}}_{\text{age,yr}}^{-1}{\text{\hspace{0.17em}B}}_{\text{μG}}^{-2}\text{\hspace{0.17em}TeV})$, then the cut-off energy in electron spectrum should be reduced by four times(~ 4.1 TeV). The total energy of electrons above 1 GeV is calculated as W_e = 2.66 × 10⁴⁶ erg, which is nearly one order of magnitude lower than the protons energy in the southern region or the electrons energy in the northern region. This indicates that the leptonic contribution in the southlobe can be nearly neglected (Fig. A1 panel left), similar to the hadronic contribution in the northlobe (Fig. A1 panel middle). In this scenario, the total contribution added from northlobe and southlobe could also be consistent with γ-ray observation results (Fig. A1 panel right), suggesting the real magnetic field strength in southlobe should be constrained around two times that of the northlobe. Since there is an absence of X-ray data points (only the upper limit detected), the synchrotron emission cannot be well constrained. Therefore, more high-resolution X-ray surveys ought to be conducted in this region.

Fig. A1 Modelling of a hybrid scenario in southlobe (left), northlobe (middle), and disk (right). In particular, the leptonic contribution in the northlobe is in agreement with the red dashed line shown in the top-right panel of Fig. 6.

Appendix B Chance coincidence with a background γ-ray source

Although the above arguments support SrcX origins from the escaped CRs illuminated surrounding molecular clouds, we still cannot totally rule out the possibility that it is a background point-like source just located in the LOS of SNR G150 further spatially coincidence with the surrounding molecular cloud clump. Here we using a Poisson distribution, the chance probability of observing a background point-like source located inside SNR G150.3+4.5 and spatially coincidence with the nearest molecular cloud clump, which can be calculated as $P_{ch}=1-\exp \left[-\pi \left(R_0^2+4{\sigma }_{\gamma }^2+R_c^2\right){\displaystyle \sum \left(>F_{th}\right)}\right]$ (Yuan et al. 2018; Xi et al. 2020), where ∑(> F_th) is the surface density of sources with fluxes higher than F_th, σ_γ is the 68% position uncertainties of SrcX, and R₀ is the angular distance between the location of SrcX and the target SNR, similarly, R_C is the angular distance between the location of SrcX and surrounding MCs clump. Since the density distribution of the 4FGL sources in the box defined by Galactic latitude |b| < 5° is uniform with respect to the angle, we estimate a number density of ∑(> F_th) = 129.5 sr⁻¹ (corresponding to 0.04 degree⁻²) above the flux F_th = 1.47 × 10⁻¹¹ erg cm⁻² s⁻¹ (0.1 − 100 GeV), using the power-law fitting of the cumulative numbers of sources as a function of the threshold fluxes for the 4FGL sources, the chance coincidence probability is estimated to be 0.0018 for R₀ = 0° (SrcX located inside target SNR and its γ-ray emission region), σ_γ = 0.025° and R_C = 0.11° (brightest pixel projection distance between the location of SrcX and molecular cloud clump). This results indicate that SrcX, as a background point source, located precisely within the interior of SNR G150 and spatially corresponds to the position of a molecular cloud clump with very low probability.

Fig. B1 4FGL source number distribution in galactic latitude (b) shown in the top-left panel. Bottom left panel: 4FGL source number distribution in galactic longitude for |b| < 5°. Right panel: Cumulative numbers of sources and the threshold fluxes for LAT 4FGL sources (histograms) with the power-law fitting results of their high-flux trends (red dashed line). The grey dotted line represents for the 95% upper limit obtained from flux analysis in Sect. 2.3.

References

All Tables

All Figures

Fig. 1 TS maps in the vicinity of SNR G150.3+4.5 observed with the Fermi-LAT. All figures are generated with square bins of 0.02°. The energy ranges are indicated above the panels. The green circle represents the best-fit 2D Gaussian template from Devin et al. (2020). The white solid and dashed circles show the R68 and radius of the uniform disk model in this work, which is separated by the white horizontal line into two half-spheres: northlobe and southlobe (models 8 & 11 & 12). The white dashed line corresponds to the +30° and −30° scenarios (models 9 and 10). The unnamed white cross in the upper right corner represents a background point source from 4FGL-DR3 (Abdollahi et al. 2020, 2022) and SrcX represents the point source 4FGL J0426.5+5434, which has been treated as a distinct background point source. The cyan contours are extracted from Urumqi λ6 cm results (Gao & Han 2014), showing the complete shell structure. The red and orange circles are R39(σ) size for the LHAASO KM2A and WCDA sources (Cao et al. 2024), respectively.

In the text

	Fig. 2 Residual maps above 2 GeV, obtained by subtracting southern/northern/both lobe, all labels match those in Fig. 1.
In the text

Fig. 3 SED and TS values (shaded areas) of the southlobe/SrcX (panel left), northlobe (panel middle), and disk (panel right) as measured by Fermi-LAT in the energy range from 0.1 to 1000 GeV. Upper limit are calculated at 95% confidence level using a Bayesian method for points with TS values less than 5. The statistical errors are shown in black/blue error bars, while the sum of statistical and systematic errors calculated by $\sigma =\sqrt{{\sigma }_{\text{stat}}^2+{\sigma }_{\text{sys}}^2}$ are marked by red error bars. The disk and SrcX can be described with a log-parabola spectrum. The northlobe is best-fitted with a broken power-law spectrum. The southlobe is compatible with a single power law and the black dash lines show the 1 σ statistic error uncertainty.

In the text

	Fig. 4 Light curves of SrcX from 0.1 to 10 GeV. The mean flux and its 1 σ statistic uncertainty are shown as the orange horizontal solid and dashed lines, respectively. The grey histogram represents the TS value for each time bin.
In the text

Fig. 5 Integrated CO (J = 1–0) emission intensity (Kkm s−1) toward SNR G150.3+4.5 at 115 GHz using CO composite survey (Dame et al. 2001). Panel left: integrated on the whole velocity range [−9.8, +2.4] km s−1, and overlaid with the radio shell structure extracted from Urumqi λ6cm results (Gao & Han 2014). Panel right: integrated on [−9.8, −3.8] km s−1 while overlaid with the GeV emission contour line generated from the top right panel in Fig. 1. Other labels are same as above.

In the text

Fig. 6 Modelling of the multi-wavelength spectra of the SrcX southlobe (top-left), northlobe (top-right), and disk (bottom-left). The blue and purple butterfly shadow regions come from recent LHAASO results (Cao et al. 2024). The orange data points come from 4FGL J0426.5+5434 as measured by Devin et al. (2020). Top left: Solid, dashed, and dash-dotted magenta lines show the spectra produced by escaping ions with different values of χ indicated, and the black solid line shows the hadronic model for trapped protons with a broken power-law distribution. Top right: Leptonic models, where the power law (ECPL) and broken power law (ECBPL) are presented by red dashed and grey solid lines, while the orange and green dashed lines show a special ECPL and ECBPL case where the cut-off energy is fixed to 25 TeV. The radio fluxes and the X-ray upper-limit come from Zeng et al. (2021) and Devin et al. (2020), respectively. Bottom left: Hybrid model for the total γ-ray flux shown as black solid, red dashed, orange dashed lines added by different leptonic spectral types. Blue data points are derived from the best-fit results measured by Devin et al. (2020) for comparison. Bottom right: Dependence of the SED on the SNR age and rs. The solid, dotted, dashed, dash-dotted lines correspond to rs = 21.4/35/45/55 pc, respectively.

In the text

	Fig. A1 Modelling of a hybrid scenario in southlobe (left), northlobe (middle), and disk (right). In particular, the leptonic contribution in the northlobe is in agreement with the red dashed line shown in the top-right panel of Fig. 6.
In the text

Fig. B1 4FGL source number distribution in galactic latitude (b) shown in the top-left panel. Bottom left panel: 4FGL source number distribution in galactic longitude for |b| < 5°. Right panel: Cumulative numbers of sources and the threshold fluxes for LAT 4FGL sources (histograms) with the power-law fitting results of their high-flux trends (red dashed line). The grey dotted line represents for the 95% upper limit obtained from flux analysis in Sect. 2.3.

In the text