LAT observations and data reduction

The LAT instrument aboard the Fermi spacecraft that collects photons in the energy range of 20 MeV to 1 TeV was employed for this study. We used approximately 16 years of Fermi-LAT Pass 8 Release 3 (P8R3) data69,70, analyzed with the latest instrument response functions (P8R3_SOURCE_V3), including both front- and back-converting events. The data span mission elapsed times from 239557417 to 750035897, corresponding to coordinated universal time from 4 August 2008 (15:43:36.0) to 7 October 2024 (23:18:12.0). All standard cleaning processes aimed at avoiding data contamination by transient events were followed. Bad time intervals were identified and excluded from the analysis. Earth limb contamination was effectively treated, by excluding data sets that exceed a zenith angle of 90∘ and 105∘, respectively, above 100 MeV and 1 GeV.

Data processing, including event selection and exposure calculations, was performed using the fermipy analysis package (v1.3.1)71, a Python interface to the fermitools software suite (v2.2.0; https://fermi.gsfc.nasa.gov/ssc/data/analysis/software/) designed primarily for binned analyses of Fermi-LAT data. The analysis employed up-to-date background models, including the Galactic diffuse emission (gll_iem_v07.fits) and the isotropic component accounting for residual instrumental and extragalactic backgrounds (iso_P8R3_SOURCE_V3_v1.txt).

The models were retrieved from the Fermi Science Support Center (FSSC; https://fermi.gsfc.nasa.gov/ssc/data/access/lat/BackgroundModels.html). All 4FGL-DR4 sources embedded within a region of 15∘ radius from the center of the region of interest (ROI) were included in the model. Following the above conventions, we performed a binned analysis over a 10∘ × 10∘ region centered at the G189.6+3.3 SNR (ROI). The spatial analysis was performed  > 1 GeV, and up to 800 GeV, while the spectral analysis was conducted in the broader 100 MeV-800 GeV energy range.

Morphological analysiseROSITA data

Due to the lack of a complete radio shell, the motivation behind this project arose from the newly acquired eROSITA data72,73. eROSITA is the main instrument aboard the Russian-German Spektrum Roentgen Gamma (SRG) observatory74. It collects photons in the 0.2-10 keV energy range and it consists of seven parallel aligned mirrors of 1∘ field of view (FOV) each (TM1-7). Reference 12 reported the first complete view of the G189.6+3.3 SNR with eROSITA. In this section, we describe the eROSITA imaging analysis used to construct spatial templates for the subsequent Fermi-LAT morphological analysis.

Aiming at a direct spatial comparison of the remnant in X-rays and gamma-rays and subsequently at examining the implications of such a spatial correlation in the two distinct energy bands, we employed publicly available data from the first eROSITA all-sky survey (eRASS1) in the c020 processing version. We performed standard data reduction, cleaning, and analysis processes employing the corresponding tools of eSASSusers_201009 version of eSASS (eROSITA Standard Analysis Software)75. We report on the eROSITA mosaic intensity map displayed on the left panel of Fig. 2 in the 0.2-5.0 keV energy band at and around the remnant’s position. The latter energy range was optimized based on spectral analysis results; these results, however, are beyond the scope of this paper. All point sources have been properly masked to improve visibility of diffuse X-ray emission from G189.6+3.3 and IC 443 SNRs. Extended diffuse X-ray emission from G189.6+3.3, which is completely thermal and has a distinct shell-type morphology that overlaps with the IC 443 to the west, is detected in eRASS1.

LAT data

We performed a binned spatial analysis, of 0.02∘ bin size,  > 1 GeV (since at lower energies the instrument PSF size exceeds 1∘), handling events types PSF0, PSF1, PSF2, and PSF3 separately in the SummedLikelihood framework. We simultaneously fitted the spectral parameters of all sources within 3∘ from the ROI center, along with the normalizations of the Galactic and isotropic background components, as well as the index of the Galactic diffuse component.

Three 4FGL-DR4 sources fall within the well-defined X-ray emitting shell extension detected with eROSITA, namely J0618.9+2240c, J0620.1+2246, and J0620.9+2201, as shown on the left panel of Fig. 2. An additional point source, namely 4FGL J0622.7+2248, lies just outside of the remnant’s X-ray emitting shell. Within a few degrees of the remnant, there are no other nearby Fermi sources. All four sources do not have known counterparts that would indicate potential association with other astrophysical objects. In the fitting process, GeV gamma-ray emission from all four 4FGL-DR4 sources was included to reproduce the emission from the G189.6+3.3 SNR.

To spatially inspect the GeV emission from the G189.6+3.3 SNR and to identify any additional sources within the ROI, we computed a test statistic (TS) map excluding from the model the four aforementioned Fermi point sources. As such, the significance of a source is examined with a standard E−2 spectrum at the position of each pixel against the null hypothesis (TS = 2(ln\({{{{\mathcal{L}}}}}_{1}\)-ln\({{{{\mathcal{L}}}}}_{0}\)), with \({{{{\mathcal{L}}}}}_{1}\) tested and \({{{{\mathcal{L}}}}}_{0}\) null hypothesis likelihoods). The obtained TS map displays gamma-ray emission from the site of the IC 443 remnant, which is most likely due to mismodeling of this prominent gamma-ray emitting SNR. At present, the gamma-ray emission from IC 443 is described using two catalog components: the extended source 4FGL J0617.2+2234e (associated with IC 443) and 4FGL J0616.5+2235 (a second gamma-ray source located within the angular extent of the remnant and also likely associated with SNR). An improved spatial template would be necessary to account for the apparent residuals, but such a study falls beyond the scope of this work.

Instead, to account for the apparent residuals and prevent mismodeling of IC 443 from affecting the analysis of G189.6+3.3 two additional point sources (RAJ2000: 94.46∘, DecJ2000: 22.37∘ and RAJ2000: 94.11∘, DecJ2000: 22.44∘) were added to the model. PS maps were generated and examined before and after the addition of the two point sources to ensure that no significant negative residuals emerged from their locations. A PS map aims to assess data with model agreement and is sensitive to both positive and negative deviations, allowing us to confirm the validity of the model by computing the p-value, which is more reliable than residual count maps76. After validating the background model using the PS maps, we recomputed the TS map above 1 GeV. The obtained TS map  > 1 GeV, shown on the left panel of Fig. 3, reveals a gamma-ray source to the east of IC 443.

However, the PSF size at the lower energy cut exceeds the source size. To confirm the extended nature of the detected GeV source and report on a clear statement of a separate gamma-ray source, which is not associated with IC 443, we restricted the spatial analysis  > 5 GeV (despite the lower statistics) to ensure an improved PSF size. The TS map  > 5 GeV, shown on the middle panel of Fig. 3, shows extended GeV emission, detected  > 10σ significance level to the northeast and at 10σ significance level to the southeast using the spatial template 13 – best fit spatial template as obtained from the morphological analysis performed in this work (refer to Table 1). This is perfectly coincident with the eastern part of the remnant’s well-defined X-ray emitting shell, as depicted in the composite image shown in Fig. 1. Furthermore, the clear separation of the GeV emission to the southeastern region of G189.6+3.3, which coincides with the enhanced X-ray emission at the eastern boundary of the remnant’s shell (Fig 1), strongly suggests the presence of a separate gamma-ray source rather than gamma-ray emission from IC 443 protons diffusing in the surrounding areas. This finding is further corroborated by the nature of the emission from the southeastern region, which appears to be leptonically induced, as discussed in subsections ‘Spectral analysis of the LAT data’ and ‘Multiwavelength SED modeling’.

Table 1 Gamma-ray spatial model comparison

Since the gamma-ray emission conventionally associated with IC 443 is predominantly interpreted as hadronically induced26, the presence of a distinct leptonic component further supports the interpretation that the detected emission originates from a separate source associated with G189.6+3.3 rather than from particles escaping IC 443. A similar conclusion is also suggested by the morphology of the northern GeV component. The enhanced GeV emission is not broadly distributed around IC 443 or simply correlated with external gas illuminated by runaway particles, but is instead spatially confined to the northern boundary of G189.6+3.3, where it coincides with the Hα filament, UV emission, optical forbidden-line diagnostics of radiative shocks, and dense material associated with S249. This close correspondence with local shock tracers favors gamma-ray production at the G189.6+3.3/S249 shock-cloud interface rather than interaction of unrelated surrounding material with particles escaping from IC 443.

To further strengthen our argument regarding a separate gamma-ray source and provide a direct comparison with13, we extended the spatial analysis to even higher energies and present the remnant’s TS map above 10 GeV – where the substantially improved Fermi-LAT PSF ( ≲ 0. 1∘ containment radius) allows a much clearer spatial separation from IC 443, as illustrated on the right panel of Fig. 3. In addition, we construct a substantially larger 5. 7∘ × 5. 7∘Fermi-LAT PS map  > 10 GeV, as shown on the right panel of Supplementary Fig. 3, without modeling the gamma-ray emission from the IC 443 region. This map provides compelling evidence for the presence of an extended source to the east and southeast of IC 443, thereby excluding the possibility of IC 443 mismodeling that could result in residual emission in its vicinity, which might be erroneously interpreted as gamma-ray emission from G189.6+3.3.

Overall, Fig. 3 and Supplementary Fig. 3 (right panel) provide a firm detection of gamma-ray emission from the G189.6+3.3 SNR that appears to be in spatial coincidence with the eastern half of the well-defined X-ray shell detected with eROSITA. A particularly notable feature of this gamma-ray emission morphology is the pronounced enhancement coincident with the Hα filament (shown as a cyan contour in Fig. 3), which appears clearly separated from IC 443 and more prominent than the gamma-ray emission from the remainder of G189.6+3.3. This feature will be discussed further in a subsequent section.

Having established the presence of extended gamma-ray emission spatially associated with G189.6+3.3, we next investigated its morphology in greater detail by first removing the 4FGL sources within the region and iteratively introducing alternative spatial source models. We explored several single-component spatial templates, including a 2D Gaussian, a radial disk, the eROSITA X-ray image (right panel of Fig. 2), and a flat (on-source region: pixel values set to 1, off-source region: pixel values set to zero) polygonal region that encompasses the X-ray emission detected with eROSITA (highlighted as a green contour on the right panel of Fig. 3 and in Supplementary Fig. 3). For the 2D-Gaussian and the radial disk templates (single source geometrical models), we initially set the added source center coordinates to (RAJ2000: 94.98∘, DecJ2000: 22.15∘; moderately shifted to the east compared to the center of the X-ray shell since we only focused on the analysis of gamma-ray emission from the eastern half of the remnant, where there is no contamination from IC 443) and performed both a position and extension fitting and an extension-only fitting by computing a likelihood ratio test with respect to the no-extension (point-source) hypothesis. We note that fitting both the position and extension of the 2D-Gaussian and radial disk templates results in the center being shifted towards the northern brighter blob and the extension, sigma or radius, being significantly reduced (localized to the region of the northern bright blob). Consequently, a portion of the gamma-ray emission from the southern region of the remnant remains inadequately modeled. The resulting positions and extensions for the 2D-Gaussian and radial-disk models are summarized in Table 1. The Akaike Information Criterion (AIC)77 is used to compare the goodness of the fit between the different tested models as follows, ΔAIC=\({{{{\rm{AIC}}}}}_{4{{{\rm{point}}}}{{{\rm{sources}}}}}-AI{C}_{i}=2(\Delta d.o.f.-\Delta ln{{{\mathcal{L}}}})\). According to both the fit quality (refer to Table 1) and the residual count maps, no single-component model can simultaneously reproduce the southern extension of the gamma-ray emission and the strong concentration of emission observed at the northern boundary of the remnant.

Consequently, since single-component models fail to adequately reproduce the observed gamma-ray morphology, and motivated by the apparent flux variation between the northern and southern regions below 10 GeV as well as by the fact that only the northern part overlaps with the S249 HII region, we tested several two-component, hereafter named G189N and G189S, source models (models 7–13 in Tabl. 1) to determine the morphology that best reproduces the observed emission. The implications of such a treatment in the characterization of the particle population(s) responsible for the gamma-ray radiation (hadronic and/or leptonic) are discussed in subsection ‘Spectral analysis of the LAT data’. The two spatial templates are motivated by distinct physical and geometric considerations. The southern component is defined by the large-scale geometry of the SNR shell as traced by the eROSITA X-ray emission, whereas the northern component is motivated by the presence of a localized Hα filament tracing the interaction with dense gas associated with the S249 HII region. For this reason, the Hα filament itself was the first spatial template tested for the northern gamma-ray emission. Notably, the Hα template used in some of those two-component templates was constructed similarly to the flat eROSITA polygonal template by setting pixel values to 1 on filament regions and zero elsewhere.

In this framework, the terms “northern” and “southern” are used as descriptive labels indicating the relative spatial location of the dominant emission regions rather than to imply strictly separate extraction areas. The two components partially overlap (see Fig. 3), and this overlap is explicitly accounted for in the likelihood analysis. Specifically, the three Fermi-LAT point sources located toward the northeastern part of the remnant are replaced by a single extended source represented by a conspicuous optical filament detected at the northern part of the remnant. A second extended component, modeled using either the eROSITA X-ray emission map and/or a radial disk template, is assigned to the southeastern part of the remnant, replacing the fourth Fermi point source, which is molecular gas-free. It is noteworthy that the observed emission from the southern part of the remnant is detected at approximately 10σ,   > 1 GeV. The variations in the TS values obtained for the gamma-ray emission detected at the southern region of the remnant, as reported in the second column of Table 1, are contingent upon the model applied to the northern region with the addition of the second component (the Hα filament, a second disk, and a 2D-Gaussian).

From the morphological analysis, we conclude that despite the apparent spatial correlation of the X-ray shell with the diffuse GeV emission, the X-ray eROSITA template alone, shown on the right panel of Fig. 2, after properly masking IC443’s X-ray emission, was not obtained as the best fit spatial model to the GeV gamma-ray data. This is primarily because the gamma-ray emission appears to extend further to the north of the X-ray shell, the X-ray emission is purely thermal (i.e., the regions of enhanced emission for X-rays and gamma-rays differ), and the gamma-ray emission appears strongly enhanced at a specific position at the north edge of the remnant. Among the above tested models, the best ΔAIC value (significantly improved compared to four point sources or single source models that account for emission from the entire remnant) is obtained while adopting a spatial template that consists of two distinct but overlying components: the polygonal region obtained from the eROSITA X-ray image of the remnant (primary contributions from its southern part) and a 2D-Gaussian at its northern boundary, spatially coincident with the Hα filament. A 5. 7∘ × 5. 7∘ PS map,  > 1 GeV, was generated and examined after running the analysis with our best fit spatial model to ensure that no significant positive or negative residuals emerged from the nearby locations (refer to the left panel of Supplementary Fig. 3). The latter model demonstrates that the gamma-ray emission is spatially aligned with the eROSITA X-ray shell and peaks at the location of the filament (center of the north bright blob). This result is evident of particle acceleration within the filament and subsequent diffusion in the surrounding medium, as further discussed in subsection ‘Spectral analysis of the LAT data’ and subsection ‘Multiwavelength SED modeling’.

To further visually demonstrate that the gamma-ray emission from G189.6+3.3 is more accurately described by extended source components rather than by the currently adopted four-point-source representation, we investigated the residual emission obtained under different source-modeling configurations. Specifically, we constructed TS maps  > 1 GeV by removing one, two, three, and then all four Fermi sources within the region of interest: 4FGL J0618.9+2240c detected with TS=192.8, 4FGL J0620.1+2246 detected with TS=124.3, 4FGL J0620.9+2201 detected with TS=30.1, and 4FGL J0622.7+2248 detected with TS=36.9. This procedure reveals significant residual emission from the southeastern part of the shell, as shown in the left panel of Supplementary Fig. 1. Even when both 4FGL J0618.9+2240c and 4FGL J0620.1+2246 sources are modeled (the two Fermi sources with the highest significance among the four), diffuse gamma-ray emission is detected from the remnant’s eastern half  > 10σ level. In fact, even after modeling all four Fermi sources, strong residuals in the form of diffuse emission from the eastern and southeastern half of the shell persist ( > 10σ), as shown on the left panel of Supplementary Fig. 1. In contrast, no comparably significant residual structures remain when adopting the best-fit extended spatial model discussed above (right panel of Supplementary Fig. 1). This direct comparison further validates the significant improvement in goodness of fit (given in Table 1) when utilizing extended models compared to the currently adopted four-point-source model.

Spectral analysis of the LAT data

A two-component source model needs to be employed to account for the entire gamma-ray emission from the remnant’s location. To this end, the G189N (originating from the long Hα filament but best represented by a 2D-Gaussian) and the G189S (represented by the polygonal region obtained from the eROSITA X-ray image – illustrated as a green contour on the right panel of Fig. 3 and in Supplementary Fig. 3) spatial templates were employed to investigate the nature of the gamma-ray emission from different parts of the remnant (spatial model 13), as discussed in subsection ‘LAT data’. The spectral analysis is performed in the broader 100 MeV to 800 GeV energy range.

Two additional years of data compared to 4FGL-DR4 have been employed for the purposes of this study. To account for potential additional faint sources not included in the 4FGL-DR4 catalog and a potential change of the morphology of the gamma-ray emission from the remnant, we construct and inspect the TS map from the remnant’s location  > 100 MeV. The shape of the gamma-ray emission from the remnant’s location remains intact – nearly identical to the one obtained  > 1 GeV. In addition, no new sources within 2∘ from the analysis ROI were detected that could directly contaminate our spectral data, thus no additional entries were included in our model prior to running the spectral analysis process.

As a next step, we tested both a simple power-law and a LogParabola spectral model for each component, using the curvature test statistic (TScurv) implemented in fermipy to quantitatively discriminate between them. The TScurv is defined as twice the difference in the log-likelihood between a curved spectral model (e.g. LogParabola) and a simple power-law model. We find that a powerlaw emerges as the best-fit spectral model for G189S. Conversely, a LogParabola describes best the energy distribution of the gamma-ray photons from G189N. The corresponding spectral parameters of the north and south components are summarized in Table 2.

Table 2 Gamma-ray spectral parameters

We then divided the 100 MeV to 800 GeV energy range into eight logarithmically spaced bins, with the highest-energy bin covering a narrower interval of 320–800 GeV. A likelihood fitting analysis was performed to derive the gamma-ray SED by computing the photon flux in each bin separately for the G189N and G189S components, with the aid of a model that contains both the G189N and G189S components (model 13) instead of the entire source. In the above process the spectral index, and α,  β spectral parameters of the source components (where applicable) were kept fixed to the best fit values. The normalizations of all 4FGL-DR4 sources within 3∘, as well as the normalizations of the background components, were let to vary. We set a TS threshold of 4 and provide 95% upper limits for intervals that exhibit lower values. The derived spectral points for both source components are illustrated in Fig. 4.

Regarding the systematics, we considered three distinct types of systematic uncertainties and added them in quadrature. These encompass uncertainties in the Galactic diffuse background, which primarily affect the lower energy regime; uncertainties in the effective area; and uncertainties in the optimal spatial model. To determine the former two types of uncertainties, we adhere to the methodology delineated in ref. 78. The latter source of systematic uncertainties is typically evaluated by extracting the spectrum for the second-best fit model, which in this case comprises a radial disk to the south and a 2D-Gaussian to the north (model 12). However, to avoid underestimating the systematic error, given that the 2D Gaussian template used to model G189N in model 12 is nearly identical, we opted for model 11 instead, which is not significantly inferior to model 13. The resultant SED, incorporating both statistical and systematic uncertainties, is presented in red in Fig. 4.

Estimates on the remnants’ intrinsic propertiesDistance estimates

G189.6+3.3: Typical luminosities of SNRs at 1.4 GHz have values in the range of 5 × 1014 < L1.4GHz < 1017 W/Hz79. To measure a lower limit for the distance, it is assumed that G189.6+3.3 has L1.4GHz = 5 × 1014 W/Hz, yielding a \({d}_{\min }=1.4\) kpc. A maximum distance estimate can also be provided considering that SNRs typically do not exceed a diameter of 100 pc80. Thus, given the angular size of G189.6+3.3 (approximately 1. 6∘ in diameter based on its visible eROSITA X-ray shell), we estimate a maximum distance of \({d}_{\max }=3.5\) kpc. We note that a slightly larger angular extent of about 1. 8∘ is suggested by eROSITA spectral evidence indicating that G189.6+3.3 and IC 443 completely overlap12.

Additionally, the absorption column density values derived from the remnant’s eROSITA X-ray analysis (X-ray spectral fit)12 and the up-to-date 3D extinction maps obtained from the combination of the 2MASS and GAIA data21 were utilized to provide an additional independent distance consistency check. Reference 81 employed a large sample of Chandra observations of Galactic SNRs to derive a statistical relation between the X-ray absorption column density (NH) and the mean color excess (EB-V)/extinction (Aν), as follows:

$${N}_{{{{\rm{H}}}}}/{E}_{{{{\rm{B-V}}}}}=8.9\times 1{0}^{21}\,{{{{\rm{cm}}}}}^{-2}\cdot {{{{\rm{mag}}}}}^{-1}\\ {N}_{{{{\rm{H}}}}}/{A}_{{{{\rm{\nu }}}}}=2.87(\pm 0.12)\times 1{0}^{21}{{{{\rm{cm}}}}}^{-2}\cdot {{{{\rm{mag}}}}}^{-1}$$

(1)

Considering the X-ray absorption column density measured for G189.6+3.3, NH = 4 × 1021 cm−2, which is found to be rather uniform across the remnant12, we derive an optical extinction of \({A}_{\nu }=1.3{9}_{-0.05}^{+0.06}\). The quoted uncertainty accounts for the range of column densities, NH = (3 − 6) × 1021cm−2, obtained from the spatially resolved X-ray spectral analysis of different sub-regions and spectral models across the remnant. Inputting such a result in the 3D extinction maps of21, implemented via the EXPLORE G-Tomo tool, yields a distance uncertainty of \({{{\rm{d}}}}=1.{8}_{-0.04}^{+0.06}\) kpc (gray shaded region on the right panel of Fig. 8). This measurement was derived by considering the cumulative extinction estimate towards the center of the remnant.

However, because G189.6+3.3 is an extended source, it is important to assess how sensitive this estimate is to the specific line of sight adopted. The EXPLORE G-Tomo tool is a visualization and query interface that extracts line-of-sight extinction profiles and two-dimensional slices from up-to-date three-dimensional dust extinction maps. Because G-Tomo is intrinsically a pencil-beam estimator (i.e., it applies to a single set of sky coordinates), its application to an extended source requires some care. We therefore queried the extinction-distance relation at multiple positions across the angular extent of G189.6+3.3 in order to assess how strongly the inferred distance varies across the remnant and to quantify the associated systematic uncertainty. Specifically, conducting the distance estimation by considering the extinction observed towards the eastern boundary of the SNR results in a distance estimate closer to 1.9 kpc. The large majority of positions sampled across the remnant yield distance estimates clustered between 1.8 and 1.9 kpc, indicating that this interval can be considered a confident distance estimate for G189.6+3.3. By incorporating the uncertainty in the cumulative extinction and the full range of positional variation across the SNR area, a broader distance uncertainty range of 1.7 to 2.3 kpc can be obtained.

IC 443:

The distance to IC 443 based on kinematic measurements of optical Hα emission is estimated to be approximately 1.9 kpc82, while independent extinction-based measurements yield a consistent distance of 1.8 ± 0.05 kpc83. A general 1.5–2.0 kpc distance for IC 443 is broadly used in the literature (e.g.,84, and references therein).

Claims that the S249 HII region lies several hundred parsecs in front of IC 443 are based on kinematic distance estimates inferred from small differences in radial velocity82. However, near the Galactic anti-center the velocity-distance relation is shallow, and kinematic distances become intrinsically uncertain and highly sensitive to modest non-circular motions and projection effects. As such, in this direction small differences in radial velocity—such as those reported for S249 and IC 443 in ref. 82—do not map cleanly to distance. In this regime, assigning a unique systemic velocity to S249 is model-dependent, particularly along a line of sight that contains multiple kinematic components.

In contrast, multiple observational tracers point to a physical association between IC 443 and S249, including the morphological correspondence between the northeastern rim/optical filaments of IC 443 and the S249 HII region37, and enhanced radio/continuum brightness and morphological distortion of IC 443 toward the northern (S249-facing) side46. In addition, gamma-ray emission from the IC 443 environment is widely interpreted as arising from particle acceleration followed by interactions with dense gas in its surroundings44, consistent with the presence of a dense, structured medium in the IC 443-S249 complex45. When combined with the detected gamma-ray emission from G189.6+3.3 reported here, these observations strongly support a common distance for IC 443, S249, and G189.6+3.3.

Age estimates

G189.6+3.3: There is no single precise, widely-accepted age for G189.6+3.3 in the literature; however, the available estimates consistently classify it as a middle-aged or evolved SNR, older than IC 443. The initial age of G189.6+3.3 was estimated to be around 105 years based on its X-ray appearance and the very soft emission detected in ROSAT data9. This estimate was, however, derived from empirical considerations and is subject to large uncertainties. Assuming instead that the remnant is in the Sedov–Taylor phase, the same work reports a younger age of 3 × 104 yr. Subsequent X-ray observations with Suzaku revealed recombining plasma along the eastern edge of the remnant, providing clear evidence that G189.6+3.3 is in a relatively advanced evolutionary stage11. Such recombining plasma signatures are commonly observed in middle-aged SNRs, including IC 44385. Taken together, these observational constraints favor an age of order a few  × 104 yr for G189.6+3.3. However, it is important to note that, depending on the ambient density and thermal history, recombining plasma features may also persist in older remnants; and therefore an age of order 105 yr remains a reasonable estimate for this SNR.

In this work, we employed different methodologies aiming at providing updated age estimates. From radio data, one can estimate an upper limit for the age of the remnant when knowing its linear radio extension (R=22.6 pc, assuming a distance of 1.8 kpc as computed above) and applying86 model, which provides an age estimate of: \(t={({R}_{pc}/21.9)}^{3.23}\times 1{0}^{5}=1.1\times 1{0}^{5}\) yr. We additionally applied the python calculator for SNR evolution as provided in ref. 87 to perform an updated age estimation. For the obtained absorption column density of the entire remnant, NH = 4 × 1021cm−2 and a distance of 1.8 kpc, we derive an average ISM number density of n0 = 0.72 cm−3. This value represents the local density at the location of the SNR, under the assumption that the line-of-sight average density is representative of the ambient medium. Considering as inputs the derived local ISM density as calculated above, a typical explosion energy on the order of 1051 erg and maintaining the remaining parameters at default input values, we derived an age of 4.5 × 104 yrs. This value is an approximation since the density at the region of G189.6+3.3 can be significantly impacted by the local medium at the remnant’s distance. As such, we also adopt 0.5 cm−3 and 0.15 cm−3 as the ISM ambient density, as computed in ref. 11 at the distance of the remnant, which yields ages of 3.2 × 104 yr and 1.8 × 104 yr, respectively.

Finally, we confirmed the aforementioned results by performing a series of computations based on the Sedov-Taylor self-similar solution as outlined in ref. 88:

$$R=0.314\cdot {\left(\frac{{E}_{51}}{{n}_{0}}\right)}^{\frac{1}{5}}\cdot {t}_{{{{\rm{yr}}}}}^{\frac{2}{5}}\,{{{\rm{pc}}}},$$

(2)

where R is the remnant’s linear radius (shock radius), E51 the kinetic energy of the explosion in units of 1051 erg, n0 the ambient density, and tyr the remnant’s age. Assuming that the SNR is still in the Sedov phase, we obtained an age estimate of 3.1 × 104 yr for n0 = 0.5 cm−3 and 1.7 × 104 yr for n0 = 0.15 cm−3.

We emphasize that both the Sedov-Taylor relations and the SNR evolution calculator87 provide global, average estimates of the time since explosion. These approaches are not intended to describe the detailed dynamical state of individual regions within the remnant. In an inhomogeneous environment, shocks propagating into dense material can locally decelerate and enter the radiative regime significantly earlier than shocks expanding into lower-density regions, while still originating from the same explosion at the same epoch. Consequently, the presence of a locally radiative shock at the interaction region at the northern boundary of the SNR does not contradict the use of a Sedov-like framework to estimate the global explosion age of the remnant.

For the present analysis, the quantity of primary relevance is the explosion epoch, as this constrains the possible time delay between the supernova events associated with G189.6+3.3 and IC 443. We therefore adopt a 20–110 kyr age estimate for G189.6+3.3. Using this age range, one can compute the shock speed as follows:

$${\upsilon }_{s}=\frac{2}{5}\times \frac{R}{t},$$

(3)

obtaining υs = 285 km/s for n0 = 0.5 cm−3 and υs = 520 km/s for n0 = 0.15 cm−3. Both sets of (n0, vs) estimates, together with the assumed distance, were employed to perform the multiwavelength SED modeling, as outlined in subsection ‘Multiwavelength SED modeling’.

IC 443:

The age of IC 443 has been calculated to range between 3 − 30 kyr38,47,54,89,90,91. On the lower limit, ref. 89 proposed an age of 3 kyr based on the shock velocity suggested by X-ray measurements of the remnant. However, the measured plasma temperature does not probe the current shock velocity. In fact, shock velocities derived from optical filaments imply a significantly older age of order 2 × 104 yr age. A comparably young age of 4 kyr was also derived by ref. 54 under the assumption that the SNR is expanding within a pre-existing cavity. On the upper limit, ref. 90 proposed an age of 3 × 104 yr under the assumption that the SNR is expanding in a dense medium by using a model for shock propagation in dense molecular material. Additionally,91 obtained a reliable age estimate of 3 × 104 yr for what is believed to be the associated neutron star to IC 443. Finally, ref. 38 reported an age estimate of 2 × 104 yr assuming a break out SNR. The most recent age estimation of IC 443 at the time of this study is derived from the dynamical modeling analysis presented in ref. 47, wherein the authors developed a 3D Hydrodynamical model to elucidate the interaction between IC 443 and its environment. This study concludes that the remnant’s age is approximately 8–9 kyr.

Comparing the aforementioned methods utilized in the various studies for age estimation, it is evident that an age as young as 3–4 kyr can be justified only if the SNR explodes in a large cavity. In the absence of a pre-explosion cavity formation, the 3–4 kyr estimate serves merely as a lower limit on the age of IC 443. Indeed, the morphology of the remnant supports an explosion in a dense medium where ejected material got trapped. Specifically, IC 443 presents as a two-shell asymmetric structure, with the two shells exhibiting different sizes. The larger radio shell expanding towards the southwest demonstrates a smooth decline in its radio continuum intensity radial profile. This indicates that the shock towards the southwest is propagating in a relatively uniform medium after breaking out, whereas the northeastern shell continues to expand in dense medium. The northeastern shell exhibits an inverse radio continuum intensity radial profile that peaks at its northeastern boundary as it continues to interact with dense medium, resulting in a smaller size. This preferred scenario, supported by observational data and consistent with the 3D Hydrodynamical model established to describe the evolution of IC 44347, provides an age estimate of approximately 8–9 kyr, which aligns well with the likely associated neutron star.

In synthesizing literature reports and the estimates presented in this study, an age uncertainty range of 3–30 kyr and 30–100 kyr has been identified for IC 443 and G189.6+3.3, respectively. This indicates that IC 443 represents the more recent explosion within the system, with a discernible time delay of up to approximately 100 kyr between the two supernova explosions.

Multiwavelength SED modeling

In this section, we investigate the physical origin of the gamma-ray emission from G189.6+3.3 through multiwavelength SED modeling. Given the distinct spectral and multiwavelength characteristics of the northern and southern regions, we model the G189N and G189S components separately. The sparseness of the available radio data which are only indicative of the non-thermal nature of the radio emission (radio synchrotron radiation), as reported in ref. 10, together with the absence of a detected non-thermal X-ray component, limits the degree to which a fully constrained multiwavelength SED modeling of the SNR can be achieved. We nevertheless carry out multiwavelength modeling to probe the origin and physical properties of the gamma-ray emission and to strengthen our conclusions.

Consequently, we utilized data from the DRAO synthesis telescope and obtained a flux density of 2.0 Jy at 1.4 GHz for the entire SNR. This value is treated as an upper limit; the radio flux-density extraction is discussed in subsection ‘Proton re-acceleration’. The data quality at 0.4 GHz does not permit a similar estimate. Furthermore, despite the purely thermal nature of the X-ray emission from the remnant’s location, we computed an upper limit on the non-thermal X-ray flux of 1.75 × 10−12 erg/cm2/s in the 2.3–5.0 keV energy range. This estimate assumes an absorbed power-law model with a typical photon index of Γ = 2 and a column density of NH = 4 × 1021 cm−2, as inferred from the X-ray spectral fit to the thermal emission.

We note that the absence of a detected non-thermal X-ray synchrotron component is not unexpected, even in the presence of a relatively strong magnetic field. In an evolved SNR interacting with dense material, synchrotron-emitting electrons are expected to suffer strong radiative losses and rapid spectral steepening, shifting any non-thermal X-ray emission to energies above the range where eROSITA is most sensitive. Indeed, eROSITA is optimized for the soft X-ray band, with a rapidly declining sensitivity and effective area above 3 keV73. Moreover, the dominance of thermal X-ray emission in the eROSITA band implies that any synchrotron component, if present, must lie below the derived upper limits and therefore remains subdominant compared to the thermal X-ray flux.

Consistent with this expectation, and given the age of the remnant, an X-ray synchrotron component might be expected to decrease sharply with a Γ = 3 spectral index. We therefore adopted, for the multiwavelength SED modeling, an alternative spectral model consisting of an absorbed power law with NH = 4 × 1021 cm−2 and Γ = 3, which yields an upper limit on the non-thermal X-ray flux of 1.66 × 10−12 erg/cm2/s in the 2.3–5.0 keV energy range, corresponding to a discrepancy of approximately 5% relative to the flatter spectral assumption. To derive the upper limits for our preferred spectral model, we used the methodology outlined in ref. 92, and specifically in Appendix A.

Contextualizing these measured fluxes in combination with the obtained Fermi-LAT SEDs, we performed a multiwavelength SED modeling using the naima package for non-thermal radiative modeling93, which yielded results consistent with a hadronic-dominated gamma-ray component for the northern part and a leptonic-dominated gamma-ray component for the southern part. We explain below our SED modeling which uses the distance of 1.8 kpc derived in subsection ‘Estimates on the remnants’ intrinsic properties’. For both components, G189S and G189N, we adopted the simplest assumption that GeV and TeV gamma rays are emitted by a population of accelerated protons and electrons distributed in the same region characterized by constant density and magnetic field strength. We described the proton spectrum by a power law with an exponential cut-off (defining the maximum energy reached by the particles) and the electron spectrum by a broken power law with an exponential cut-off. The break in the electron spectrum is due to synchrotron cooling, and occurs at the energy Eb above which the particle spectral index is se,2 = se,1 + 1, with se,1 being the spectral index below the energy break. We emphasize that this break should be regarded as a phenomenological, volume-averaged representation of radiative losses, rather than a sharp physical discontinuity, as spatial and temporal variations in the shock and magnetic field are expected to smear out such features in more realistic multi-zone or hydrodynamical descriptions.

To make the derivation of the principal model parameters explicit, we briefly summarize the key relations used to estimate the characteristic particle energies. Following the procedure described in detail in ref. 88, the electron cooling break and the maximum particle energies are obtained by comparing the relevant acceleration and radiative loss timescales with the characteristic age of the remnant, as commonly adopted in standard one-zone broadband SED modeling of SNRs. In the framework of DSA, the acceleration time can be written as

$${t}_{{{{\rm{acc}}}}}=30.6\,{k}_{0}\left(\frac{E}{1\,{{{\rm{TeV}}}}}\right){\left(\frac{B}{100\mu {{{\rm{G}}}}}\right)}^{-1}{\left(\frac{{v}_{s}}{1000{{{\rm{km}}}}{{{{\rm{s}}}}}^{-1}}\right)}^{-2}{{{\rm{yr}}}},$$

(4)

where E is the particle energy, B is the magnetic-field strength, vs is the shock velocity, and k0 represents the ratio between the particle mean free path and the gyroradius. The synchrotron cooling time for relativistic electrons is given by

$${\tau }_{{{{\rm{sync}}}}}=1.25\times 1{0}^{3}{\left(\frac{E}{1{{{\rm{TeV}}}}}\right)}^{-1}{\left(\frac{B}{100\mu {{{\rm{G}}}}}\right)}^{-2}{{{\rm{yr}}}}.$$

(5)

The cooling break energy of the electron spectrum is obtained by equating the synchrotron cooling time with the age of the remnant,

$${\tau }_{{{{\rm{sync}}}}}({E}_{b})={t}_{{{{\rm{age}}}}}.$$

(6)

The maximum particle energy is determined by equating the acceleration time to either the remnant age or the dominant loss timescale,

$${t}_{{{{\rm{acc}}}}}({E}_{\max })=\min \left({t}_{{{{\rm{age}}}}},\,{\tau }_{{{{\rm{sync}}}}}\right).$$

(7)

For protons we consider the age-limited case, while for electrons the maximum energy is obtained by equating the acceleration and synchrotron loss timescales, tacc(Emax,e) = τsync(Emax,e). These estimates depend on the unknown ratio between the mean free path and the gyroradius (k0). For young SNRs, we expect k0  ≃ 1, while for evolved systems, we expect k0 > 1. Below, we fixed this value at 1, providing therefore an upper limit on the derived maximum energy. In subsection ‘Age estimates’, we derived the age / velocity estimates of 3.1 × 104 yr / υs = 285 km/s for n0 = 0.5 cm−3, and 1.7 × 104 yr / υs = 520 km/s for n0 = 0.15 cm−3. These values are used below to constrain the characteristic particle energies in the SED modeling.

For the G189S component, the constraining radio upper limit implies a low magnetic field of B = 5μG and a hard injection index of 2.0 consistent with DSA. The shock speeds derived above for the remnant imply Emax,p = 4.2 TeV and Emax,p = 7.5 TeV, for densities of 0.5 cm−3 and 0.15 cm−3, respectively. When adopting these two distinct density/shock speed values, the only resulting difference is in the maximum energy (Emax,p/e), which has a negligible impact on the resulting SED. For this reason, we adopt an intermediate value of 5 TeV. The physical parameters used to reproduce the multiwavelength data for the G189S component are summarized in Supplementary Table1 as representative solutions, while the modeled SED is presented in Fig. 5, showing that the dominant radiation at gamma-ray energies is IC scattering of accelerated electrons on ambient photon fields due to the low density of the medium in this region which suppresses proton-proton interaction. Very-high-energy (TeV) measurements from HAWC and LHAASO are shown for qualitative comparison only and are not included in the spectral modeling, owing to current uncertainties in source morphology and association within the IC 443 complex. We note that while such a low magnetic field implies long synchrotron cooling times and could in principle allow for electron escape, efficient escape would be expected to produce a broader or spatially displaced gamma-ray morphology, which is not observed. Within the angular resolution of the Fermi-LAT data, the southern gamma-ray emission remains spatially correlated with the eROSITA X-ray shell, indicating that electron escape does not significantly affect the adopted spatial modeling.

Due to the lack of multiwavelength data on the G189N component, we focus on the constraints imposed by the gamma-ray spectrum derived in this work. The density of the medium within the S249 HII region amounts to n0 = 50 cm−3 to 2000 cm−3 as reported by ref. 16. We used the most conservative value of 50 cm−3, fixed the particle injection index at 2.0 and the electron-proton ratio at 0.01, as done above for the G189S component to reproduce the gamma-ray data. The steepness of the gamma-ray spectrum above a few GeV implies a low energy cut-off value of 100 GeV as clearly demonstrated in Fig. 6 (left panel). Alternatively, the gamma-ray data can be reproduced using a steeper injection spectrum of 2.3 and a higher energy cutoff at 0.6 TeV (middle panel). This low energy cutoff with respect to the southern component can be due to the recent interaction of the supernova remnant with the dense material of the S249 HII region, which also strongly influences the energetics and radiative properties of the system.

Within the framework of these standard DSA scenarios, both spectral solutions require a comparable energy injected into relativistic protons and, given the broad uncertainty in the ambient gas density (n0) discussed above, the latter can be expressed in a scalable form as

$${W}_{p}\simeq 1.9\times 1{0}^{48}\left(\frac{50\,{{{{\rm{cm}}}}}^{-3}}{{n}_{0}}\right){{{\rm{erg}}}}\quad \left({{{\rm{equivalently}}}}\,{W}_{p}\simeq 9.5\times 1{0}^{49}\left(\frac{1\,{{{{\rm{cm}}}}}^{-3}}{{n}_{0}}\right){{{\rm{erg}}}}\right).$$

(8)

Although the gamma-ray spectrum in the northern region can be reproduced using proton spectra compatible with standard diffusive shock acceleration, the dense environment and inferred interaction with the S249 HII region imply strongly decelerated and radiative shock conditions. Under such conditions, fresh particle acceleration to high energies is expected to be inefficient, and the observed low cutoff energy is naturally explained by compression and re-acceleration of pre-existing cosmic rays in the dense post-shock gas. As such, alternative acceleration mechanisms associated with dense gas, and capable of reproducing the observed emission in radiative regions, are explored in subsection ‘Proton re-acceleration’.

The dense environment also strongly constrains alternative leptonic interpretations of the northern gamma-ray emission. Reproducing the observed gamma-ray flux through a Bremsstrahlung-dominated scenario would require an electron-to-proton ratio approaching unity, significantly larger than typically inferred for Galactic cosmic rays and SNRs. Likewise, reproducing the observed gamma-ray emission through an IC-dominated scenario would require unrealistically low ambient densities and unusually high electron-to-proton ratios in order to suppress both the proton–proton interaction and Bremsstrahlung emission. Such conditions are difficult to reconcile with the dense environment of the S249 interaction region and with the strong spatial correlation between the gamma-ray emission, molecular material, and the prominent Hα filament.

Proton re-acceleration

The intensity variation observed between the northern (hadronic) gamma-ray blob, spatially coincident with the Hα filament, and the leptonically induced gamma-ray emission from the remainder of the remnant’s area suggests a scenario of proton re-acceleration followed by contraction within the thin Hα filament. To further investigate this scenario, we performed crushed-cloud modeling of the northern hadronic gamma-ray component, following the framework introduced by ref. 94 and subsequently developed in dedicated numerical studies (e.g.95,96,97), with the aim of assessing whether a re-acceleration scenario can reproduce the observed emission.

In particular, we suggest that the gamma-ray emission from the G189N component is the result of the shock reacceleration of the energetic nuclei that constitute the ubiquitous Galactic cosmic ray background, followed by their adiabatic compression in the radiative layer downstream of the shock. We assumed that the SNR shock is still evolving in the Sedov phase and first encountered the S249 cloud when its radius (age) was about 16 pc (13 kyr) and its velocity approximately 500 km/s. When the SNR shock hits the cloud, a strong and much slower (50 km/s) shock is driven into the dense gas. If the cloud has a density equal to 50 cm−3, the current value for the total mass of the gas shocked so far is  ≳ 300M⊙. After being reaccelerated at the shock and compressed downstream of it, cosmic rays interact with the very dense radiative layer and produce gamma-rays via hadronic interactions with the ambient gas. A comparison between the prediction of the model and data is shown on the right panel of Fig. 6, where we assumed that the maximum energy of particles accelerated at the shock is 150 GeV and that all reaccelerated particles remain confined downstream of the shock. This indicates that a simple crushed cloud model as the one developed in ref. 94 can provide a quite natural explanation of the gamma-ray emission from G189N.

To further support our claim, we searched for UV signatures from the area where the remnant interacts with the molecular material. In principle, UV emission serves as an effective indicator of the radiative shocks that might be present in interstellar clouds with high densities. Unlike non-radiative shocks, radiative shocks are slower and do not efficiently accelerate particles via the DSA mechanism. As such, the detection of prominent UV filaments from the region where the remnant interacts with dense material would further substantiate the model involving the re-acceleration of preexisting ambient cosmic rays.

No Galaxy Evolution Explorer (GALEX) observations of the remnant have been performed; however, we found three Swift UVOT observations of small portions of the northern part of the remnant (ObsID: 00038720001, 00084926003, 00084926004), one of which (00084926003) partially covers the prominent Hα filament. In Supplementary Fig. 2, we present the Swift UVOT level 2 processing data (reduced and corrected as provided by the Barbara A. Mikulski Archive for Space Telescopes [MAST]) recorded with the U filter (central wavelength: 346.5 nm), which corresponds to the near-UV (UVA) energy range. The UV emission closely traces the morphology of the optical filament. When considered together with the strong optical forbidden-line emission and elevated S II to Hα ratios measured in this region16, the UV detection provides independent and complementary evidence that the filament traces radiative, cooling post-shock gas. This multiwavelength consistency strengthens the interpretation that the northern boundary of the remnant, where it interacts with the S249 HII region, is dominated by radiative shocks, thereby supporting a scenario in which the observed gamma-ray emission arises from the re-acceleration of pre-existing cosmic rays rather than an interpretation linked to freshly accelerated particles.

If the gamma-ray emission from G189N is indeed associated with compression and re-acceleration of pre-existing cosmic rays within the radiative shock interacting with the S249 region, then enhanced radio synchrotron emission is also expected from the same interaction site, since the processes responsible for re-accelerating relativistic protons should likewise act on relativistic electrons, albeit potentially with different efficiencies and energy-dependent acceleration rates. In this context, our objective was also to report on the predicted radio synchrotron flux density from the radiative shock and to compare it with observational DRAO data at 1.4 GHz from the filament location. Our ultimate goal is to establish a reference radio synchrotron flux density value from G189N to facilitate future research endeavors. To this end, we have performed an aperture-based radio flux-density estimate for the filament region (G189N) using the CGPS/DRAO 1.42 GHz data. As expected in the re-acceleration scenario discussed above, the filamentary radio structure, confirmed to be non-thermal in nature by ref. 10, is clearly enhanced relative to both its immediate surroundings and the broader SNR emission (Supplementary Fig. 4), despite the inherent background uncertainty and potential confusion between thermal and non-thermal radio contributions from the S249 region, which may however, affect the precise determination of the discrepancy between the predicted and observed radio synchrotron flux density.

Given this localized radio enhancement, it becomes important to estimate the synchrotron flux density specifically associated with the filament region rather than with the remnant as a whole. In instances where an SNR exhibits a highly non-uniform distribution of radio emission, the total flux density may not accurately represent the flux density of particularly brighter regions. To quantify this, we utilize the DRAO-1.42 GHz map (as applied in the computation of radio flux density from the entire SNR) and focus our analysis on the northern bright radio boundary of the SNR (G189N). We select as on-source flux density extraction region the entire wide arc feature (comparable to the size of the G189N gamma-ray component) seen in the DRAO maps encompassing the Hα filament, which is evidently even brighter. We employ two distinct background control regions in the vicinity of the SNR but well outside its extension. One is positioned in a region devoid of thermal emission from S249, while the second is placed within S249 to account for potential thermal emission contributions from the on-source region. Applying the Rayleigh-Jeans law we obtained a radio synchrotron flux density of 2.0 Jy (accounting for thermal contributions from S249 and/or apparent radio contributions that might be unrelated to S249 to the north of G189.6+3.3, as if G189N is embedded at the most enhanced radio emission regions detected across S249 area) and 6.3 Jy (considering no contributions from the S249 region).

We emphasize that extracting the radio synchrotron flux density from the entire remnant, based on the available observational data, presents significant challenges, primarily due to the background uncertainty of a large-area SNR, the overlap with IC 443 to the west, and the presence of a thermal radio component (S249) potentially overlapping with the northern regions of the SNR. It is also evident from the radio surveys (refer to Supplementary Fig. 4) that the radio emission from the full shell is significantly dimmer compared to the northern boundary of the remnant coincident with G189N and is close to, if not at, the background level. Thus, we conservatively adopt the lower end of the estimated flux-density range (2 Jy) when discussing the radio synchrotron emission associated with the remnant, while future radio observations with improved sensitivity, angular resolution, and multi-frequency coverage are expected to enable improved background characterization and a more robust disentanglement of thermal and non-thermal emission, even if the full extent of the remnant is not completely recovered.

Although these flux-density estimates provide a useful reference for comparison with the re-acceleration model, the radio spectral index could in principle serve as an additional diagnostic of whether the filament emission is consistent with a particle re-acceleration scenario. However, the currently available radio data do not permit a sufficiently robust determination. While radio spectral indices can often be robustly derived for bright, isolated filaments where background uncertainties largely cancel, this assumption does not hold for the faint filamentary emission discussed here. In this case, the structure is detected with varying significance across the available radio surveys, each characterized by different angular resolution, spatial filtering, and background properties. As a result, uncertainties in background subtraction do not cancel in spectral-index measurements and instead dominate the inferred slope. In addition, while a radio spectral index can formally be derived from two flux-density measurements, although additional measurements at multiple frequencies are strongly preferable, a robust determination requires sufficiently detected emission in matched-resolution datasets with well-controlled systematics.

In the present case, although the filament is clearly identifiable in the available radio surveys, the emission remains relatively faint and heterogeneous across sufficiently separated frequencies, making the inferred slope highly sensitive to systematic effects. In addition to the factor-of-three uncertainty associated with the CGPS background treatment, any corresponding estimate from the GB6 survey is expected to carry even larger uncertainties because of differences in angular resolution, spatial filtering, sensitivity to diffuse emission, and local background structure. Under these conditions, any derived spectral-index estimate would be dominated by systematic rather than statistical uncertainties and therefore would not provide a sufficiently robust physical diagnostic. We therefore refrain from overinterpreting any spectral-index estimate as a diagnostic for the re-acceleration scenario and instead focus on the consistency between the observed radio flux density (estimated above) and that predicted by the re-acceleration model, which provides a more robust assessment given the limitations of the current radio data.

To estimate the level of synchrotron emission expected under the radiative shock re-acceleration scenario, we adopted the crushed-cloud framework previously applied to evolved interacting SNRs. We applied pressure equilibrium to derive the radiative compression factor, s = (compressed cloud density/cloud density)/rsh, with the shock compression ratio rsh = 4. For an upstream magnetic field value of 10 μG in the cloud, a cloud density of 50 cm−3, and a 50 km/s shock propagating within the cloud, we adopted the methodology outlined in refs. 41,94 to compute the relevant physical parameters. This approach yielded a compressed density of 1662 cm−3 and a compressed magnetic field of 271μG (refer to eq. 3, 4, and 5 in refs. 41). Consequently, the value of s is determined to be 8.31.

Using these compressed physical conditions, we estimated the expected synchrotron-to-π0-decay ratio by comparison with the Cygnus Loop, which represents a similarly evolved interacting SNR. Synchrotron emission scales with the magnetic field strength raised to the power of (1+α), where α represents the electron energy index (here α=0.5). In contrast, the π0 decay scales with ambient density. Given that the Cygnus Loop is a SNR with conditions similar to our case, we compute the ratios of the compressed density and magnetic field of the two remnants. Those are estimated to be 5.67 and 1.11, respectively. Therefore, at this first-order approximation, the ratio of synchrotron to π0 for G189.6+3.3 should be approximately 0.206 of that in the Cygnus Loop. Considering the predicted radio flux density (at 1.4 GHz) for the re-acceleration scenario of the Cygnus Loop41 and the gamma-ray flux density estimates in the 0.1-100 GeV energy range for both the Cygnus Loop41 and G189N, 1.89 ± 0.14 × 10−5 MeV cm−2 s−1 (or 3.03 ± 0.22 × 10−11 erg cm−2s−1) as derived in this work, we predict a synchrotron flux density at 1.4 GHz of 4.5 ± 0.6 Jy.

Consequently, the apparent radio synchrotron enhancement, as depicted in all radio maps of different frequencies in Supplementary Fig. 4, which is evidently spatially coincident with the thin filament, and the broadly consistent results between the observational radio flux-density estimate from the filament region and the predicted radio synchrotron flux density from the crushed-cloud model, further support the re-acceleration scenario.