© The Authors 2025

1 Introduction

One of the main topics in astroparticle physics is the identification of the sources of cosmic rays (CRs). CRs are particles that originate in astrophysical environments and have energies up to ~10²⁰eV. Since the second half of the last century (Ginzburg & Syrovatskii 1964), supernova remnants (SNRs) have been considered the sources of the bulk of Galactic CRs. This hypothesis has indirectly been supported by recent observations of SNRs at MeV, GeV, and TeV gamma rays. In particular, a striking piece of evidence is the so-called pion-decay bump, which can be observed in the energy range 50 MeV–200 MeV. This spectral feature unequivocally identifies gamma rays produced by π⁰ decay, and consequently, relativistic protons (Ackermann et al. 2013).

The SNR paradigm for the origin of CRs also received support from other observational evidence, including the detection of nonthermal X-ray filaments (Vink 2012) as well as the broadening of Hα lines from Balmer shocks (Morlino 2014). However, several aspects of the problem still remain unsolved. The most relevant aspect concerns the maximum energy: It is unclear from a theoretical perspective whether SNRs can explain the CR spectrum up to the knee energy of ~ PeV, as required by observations (e.g. Schure & Bell 2013; Cardillo et al. 2015; Cristofari et al. 2020). From an observational point of view, the Large High Altitude Air Shower Observatory (LHAASO) collaboration has shown that only a few SNRs are associated with gamma-ray emission beyond ~100 TeV (Cao et al. 2024). Unexpectedly, they are all middle-aged SNRs, while younger objects seem not to show the same high-energy emission. Whether this emission results from particles that escaped from the SNRs when they were younger but are still located around the acceleration region, or if this is connected to pulsar wind nebulae (e.g. Breuhaus et al. 2022) remains to be addressed. Moreover, the final CR spectrum that is released into the Galaxy is poorly constrained. It is in general different from the spectrum that is accelerated at the shock. Both aspects, the maximum energy and the released spectrum, are connected to the way in which particles escape from the SNR, which is difficult to describe given the high nonlinearity that is usually introduced in diffusive shock acceleration (DSA) theory (see, e.g. Blasi 2013). One possible way to study the escape process is to examine middle-aged SNRs and try to detect gamma-ray emission from escaping particles.

In this regard, the case of the middle-aged SNR W44 (~20 kyr) is of great interest. It is embedded in a dense molecular cloud (MC), as observed, for example, by Reach et al. (2005) in the millimeter and near-infrared bands and by Park et al. (2013) in the HI 21 cm line. It was discovered at gamma-ray energies by Abdo et al. (2010) and Giuliani et al. (2011) with the Fermi-Large Area Telescope (LAT) and the Astrorivelatore Gamma ad Immagini LEggero (AGILE) Gamma-Ray Imaging Detector (GRID) instruments, respectively. Ackermann et al. (2013) confirmed the finding that the low-energy cutoff in the gamma-ray spectrum is consistent with emission from decaying π⁰ produced in collisions of accelerated protons. In addition, Uchiyama et al. (2012) showed the presence of gamma-ray emission toward the northwest (NW) and southeast (SE) regions of W44, which has been interpreted as a sign of escaping CRs illuminating the MC complex. More recently, Peron et al. (2020) repeated the Fermi-LAT analysis of the W44 region with a threefold increased exposure compared to the previous study. They found the NW source to be significantly larger than previously thought. This might be caused by a large underlying background component, however, that is not properly included in the background model.

Since W44 is a very bright source in the GeV band, its morphological features affect the features of the surrounding sources. Consequently, it is necessary to describe the morphology of the remnant well.

Although the origin of gamma rays from W44 has been established to be hadronic, the underlying acceleration mechanism is still debated. In its test-particle limit, the DSA model exhibits shortcomings as it tends to predict spectra that are harder than the observed spectra. Thus, following the original idea by Blandford & Cowie (1982), numerous authors have proposed an alternative scenario in which preexisting Galactic CRs are reaccelerated and compressed by the expansion of the SNR forward shock into a dense environment (Uchiyama et al. 2010; Lee et al. 2015; Tang & Chevalier 2015; Cardillo et al. 2016). These models usually require a large fraction of the shock surface (≳50%) to interact with the dense molecular environment (Cardillo et al. 2016). However, it was recently pointed out that some of these models require a circumstellar density that does not fit the evolutionary stage of the SNR (de Oña Wilhelmi et al. 2020; Sushch & Brose 2023).

Another plausible scenario is connected with the escape process. Spectra steeper than E⁻² can result from regular DSA when the time-dependent escape of particles is taken into account together with the fact that the maximum energy of accelerated particles is a decreasing function of time (Celli et al. 2019; Brose et al. 2020; Acciari et al. 2023). Models accounting for particle escape are also better suited to explain the gamma-ray emission from the region outside the SNRs. The escaping flux is tightly connected to the spectrum inside the SNRs, and when the latter is known, the former can therefore be also determined. In contrast, reacceleration models do not make clear predictions of the escaping flux.

The paper is organized as follows. In Sect. 2, we present our analysis of Fermi-LAT data and discuss the morphology (Sect. 2.1) and spectrum (Sect. 2.4). In Sect. 3 we present the novel Major Atmospheric Gamma Imaging Cherenkov (MAGIC) observations and data analysis, and in Sect. 4, we use the escape scenario to model the gamma-ray emission from W44 and the surrounding regions. We summarize our results in Sect. 5.

2 Fermi-LAT observations and data analysis

The Fermi Gamma-ray Space Telescope is a gamma-ray observatory that was launched by NASA on June 11, 2008. It orbits at an altitude of ~565 km with an inclination of ${25}_{\cdot }^{\circ }6$ with respect to the equator. The main instrument on board the satellite is the Large Area Telescope, which is an imaging pair-tracking telescope with a wide field of view that detects photons from 20 MeV to a few TeV (Atwood et al. 2009).

We selected data covering 142 months from a 15° region of interest (RoI) centered on the position of SNR W44 (1: 34.65, b: −0.38) in the energy interval between 100 MeV and 2 TeV.

We performed a binned likelihood analysis in a square region inscribed in the selected RoI. The data were spatially binned with pixels of $0_{\cdot }^{\circ }1$ and divided into energy bins using eight bins/decade. We selected the SOURCE event class, which is the recommended selection for steady point-like or moderately extended sources. In order to increase the Fermi-LAT sensitivity, we divided the dataset according to the event types available in the Pass 8 reprocessing of Fermi-LAT data (Ajello et al. 2021; Atwood et al. 2013; Bruel et al. 2018). In particular, we used the point spread function (PSF) event types, which from PSF0 to PSF3 provide an increasingly better angular reconstruction. We defined components by applying different cuts in event type, energy, and maximum zenith angle z_max. The latter was used to reduce the contamination from the Earth’s limb. A tighter cut was applied at lower energy because the Earth’s limb flux and the LAT PSF increase.

We combined the components using a summed likelihood analysis procedure in which we applied the corresponding set of IRFs, P8R3_SOURCE_V3, to each component. We adopted the following components, as prescribed by Abdollahi et al. (2020):

• PSF2 and PSF3 events with z_max = 90° for the energy range 100–300 MeV;

• PSF1, PSF2 and PSF3 events with z_max = 100° for the energy range 300 MeV–1 GeV;

• PSF0, PSF1, PSF2 and PSF3 events with z_max = 105° for the energy range 1 GeV–2 TeV.

We conducted a morphological analysis of the W44 region by selecting PSF2 and PSF3 events with an energy above 1 GeV and with a maximum zenith angle of 105° in order to exploit the best angular resolution of the instrument¹. On the other hand, we considered the full statistics in the entire energy range 100 MeV–2 TeV to perform a spectral analysis. The correction for the energy dispersion was enabled.

We modeled the RoI considering all the known sources from the 4FGL-DR2 catalog (Ballet et al. 2020) within 20° from the RoI center and the Galactic and isotropic diffuse background models (gll_iem_v07, iso_P8R3_SOURCE_V3_v1²). The normalization parameters of the diffuse components and of sources with a significance higher than 4 σ were fit, while other parameters were kept fixed at the catalog values. Moreover, the correction for energy dispersion was disabled for the isotropic template³. We conducted the analysis using the fermi-tools v2.0.8 (Fermi Science Support Development Team 2019) and the fermipy v1.0.1 (Wood et al. 2017) packages.

2.1 Morphological analysis

We first performed a detailed study of the morphology and spectral shape of the region around the W44 remnant in the energy range 1 GeV–2 TeV. As already stated in the previous sections, we limited the energy range and selected the PSF2+PSF3 event types in order to exploit the better PSF of the instrument.

We first removed all known sources within 1° from the RoI center except for W44. These were

• 4FGL J1856.7+0125c, with flags⁴ 3, 5, 6;

• 4FGL J1857.4+0106, with flags 5, 12;

• 4FGL J1855.8+0150, with flag 3;

• 4FGL J1857.1+0056, with flag 3;

• 4FGL J1857.6+0143, with flags 9, 12;

• 4FGL J1854.7+0153, with flag 3;

• 4FGL J1859.3+0058, (no flags);

• 4FGL J1857.6+0212, with flag 5;

• 4FGL J1858.3+0209, with flag 5.

We then applied an algorithm designed to detect point-like sources within 1° from the center of the RoI. The algorithm used a test source with a power-law (PL) spectrum $$\text{(PL)}\frac{dN}{dE}=N_0{\left(\frac{E}{E_0}\right)}^{-\gamma }$$

with the source position (two parameters) and the normalization and index of the spectrum left free in the fit. The test statistics (TS) was calculated for the test source, where the TS is defined as T S = 2(ln ℒ₁ − ln ℒ₀), with ln ℒ₁ and ln ℒ₀ being the log-likelihood of the models with and without the test source, respectively. Every source found with T S > 25 (equivalent to a significance of approximately 4σ for a fit with four degrees of freedom (Mattox et al. 1996) was tested for evidence of spatial extension and spectral curvature, starting with the most significant source, and the best-fit model for that source was included in subsequent passes of the source-finding algorithm. The procedure was repeated until no new sources were found. The spatial extension test was performed by refitting that source with a disk model in which the radial parameter was left free to vary. A source was considered extended when T S_ext > 25, where T S_ext is defined as twice the log-likelihood difference between a model H1 with an extended disk model and a null hypothesis H0 based on a point-like source: T S_ext = 2 (ln ℒ_disk − ln ℒ_point). The sources were tested for spectral curvature by replacing the PL shape with a curved spectrum. As alternative models, a log-parabola (LP) and a broken power-law (BPL) spectrum⁵ $$\begin{array}{l}\text{(LP)}\frac{dN}{dE}=N_0{\left(\frac{E}{E_b}\right)}^{-(\alpha +\beta \ln \left(E/E_b\right))}\\ \text{(BPL)}\frac{dN}{dE}=N_0\times \{\begin{array}{ll}{\left(\frac{E}{E_b}\right)}^{{\gamma }_1}& \text{ if }E<E_b\\ {\left(\frac{E}{E_b}\right)}^{{\gamma }_2}& \text{ otherwise }\end{array}\end{array}$$

were considered. We evaluated the significance of the curvature using the quantity T S_curv = 2(ln ℒ_LP/BPL − ln ℒ_PL). The quantity ln ℒ_LP/BPL represents the log-likelihood obtained with the LP/BPL model and ln ℒ_PL the one obtained with a PL model. The source was considered curved for T S_curv > 20, and the proper spectral model was chosen as the one providing the highest value of T S_curv.

We point out that the thresholds we used to define an extended or curved source are higher than those adopted in previous works (Ackermann et al. 2018; Abdollahi et al. 2020). This was chosen to avoid an excessive increase in the number of free parameters in the model and to minimize the risk of source confusion.

In order to minimize the influence of close-by extended sources, we also performed an extension and spectral fit for Kes79 and HESS J1857+026, which lie beyond the one-degree circular region.

2.2 W44 morphology

The morphology of the extended SNR W44 has been widely studied in gamma rays for almost 12 years. The first template, an elliptical ring, was derived by Abdo et al. (2010) based on 11 months of data collected by the Fermi-LAT. It is still used as the spatial model for this source in the Fermi catalog (Ballet et al. 2020). Giuliani et al. (2011) confirmed based on AGILE data the quasi-elliptical morphology of the emission, which resembles the 324 MHz radio continuum flux density detected by the Very Large Array (VLA) (Castelletti et al. 2007). Using three years of data of Fermi-LAT, Uchiyama et al. (2012) modeled the spatial distribution of the gamma-ray emission from W44 with a 10 GHz synchrotron radio map. Cardillo et al. (2014) studied the W44 SNR with AGILE data obtained up to June 2012, compared with VLA data and NANTEN2 telescope carbon monoxide (CO) data in the velocity channels 41 km/s and 43 km/s. They highlighted that most of the gamma-ray emission comes from a site of interaction between the SNR and the MC. Finally, Peron et al. (2020) showed in their analysis of almost 10 years of Fermi-LAT data that the W44 SNR morphology in gamma rays can be better described by an elliptical ring than with other spatial shapes they tested (disk, full ellipse, or radio templates derived from Condon et al. 1998 and Egron et al. 2017).

In order to determine the best spatial description of W44 SNR, we repeated the procedure illustrated in the previous section several times and each time adopted a different spatial template for W44. The W44 spectral shape was always described using an LP function, whose parameters were left free in the fitting procedure. We chose the best configuration by evaluating the value of the Akaike information criterion (AIC) (Akaike 1998) for each analysis result. It is defined as AIC = 2k − 2 ln ℒ, where k is the number of free parameters, and ln ℒ is the log-likelihood of the model. The best model is the one for which the AIC is minimized. Since the variation in the number of parameters depends essentially on what happens within 1° from the center of the RoI, we only considered the sources in this restricted region to calculate k.

We considered the following spatial templates (Fig. 1) to describe the morphology of the remnant:

• the 4FGL catalog template, derived from Abdo et al. (2010) after only two years of operations of Fermi-LAT, consisting of an elliptical ring;

• a simple full-ellipse template;

• a radio (1420MHz) template derived from The HI/OH/Recombination line survey of the Milky Way (THOR) (Beuther et al. 2016)⁶.

We derived the dimension and orientation of the full-ellipse template with a dedicated analysis, at first by fitting templates with varying inclination angles with 10° steps in the range [105°, 145°] and semimajor and semiminor axes with $0_{\cdot }^{\circ }0.4$ steps in the intervals [0.3°, 0.46°] and [0.16°, 0.32°]. The centroid of the ellipse was kept fixed at the coordinates taken from the 4FGL catalog template. We used the AIC to determine the best elliptical template. The minimum AIC value was obtained for an inclination angle of 115°, a = $0_{\cdot }^{\circ }38$, b = $0_{\cdot }^{\circ }2$. A finer step of $0_{\cdot }^{\circ }01$ was used for the ellipse semiaxes by fixing the inclination angle to 115°. The log-likelihood was maximized and the AIC minimized for [a = $0_{\cdot }^{\circ }41$ ± $0_{\cdot }^{\circ }01$, b = $0_{\cdot }^{\circ }23$ ± $0_{\cdot }^{\circ }01$]. We also divided the catalog template and the full-ellipse template, similarly to what was done by Ajello et al. (2016), in search for the best division angle, starting from the major axis and considering another 18 angles in the range ±90° around the central one. In Figs. 1b) and 1d) we show the divided templates. The best division angle for the elliptical ring is 115°, and for the full ellipse, it is 125° with respect to the horizontal axis. We fit the two half-ellipses as independent sources (doubling the number of free parameters).

In Table 1 we report the results of the HE (1 GeV–2 TeV) Fermi-LAT analysis obtained for the various models we adopted to describe the RoI. To calculate the number of degrees of freedom (d.o.f.), we considered 3 d.o.f. for each independent source employed to model W44 (since the spatial component was fixed and the spectral model was an LP), 2 d.o.f. for each point-like and 3 d.o.f. for each extended source, 2 d.o.f. for each source with a PL spectral model, and 3 d.o.f. for each source with a LP spectral model. In this energy range, no source was characterized by a BPL spectral model.

The maximum value of log-likelihood and the minimum AIC value were obtained with a model that incorporates the radio template to describe the emission from W44 itself and an additional disk template, hereinafter referred to as the diffuse disk, to describe centrally located diffuse gamma-ray emission in the surroundings of W44 that is otherwise unaccounted for. Four sources were used to describe the region around the remnant in addition to the radio template and the diffuse disk: two moderately extended sources located at opposite edges along the major axis of the SNR, PS J1855.5+0141e and PS J1856.9+0102e, and two additional point-like sources detected by the source-finding algorithm, which correspond to the 4FGL sources 4FGL J1858.3+0209 and 4FGL J1859.3+0058. The extension of PS J1856.9+0102e also includes the locations of the two sources that were previously removed from the RoI, namely 4FGL J1857.4+0106 and 4FGL J1857.1+0056. PS J1855.5+0141e and PS J1856.9+0102e correspond to the previously detected sources SRC-1 and SRC-2 in Uchiyama et al. (2012). We use the convention of the more recent Peron et al. (2020), who referred to them as NW and SE, respectively. As previously discussed, these are spatially modeled as disks. Table 2 reports the best-fit position and radius of the diffuse disk and the NW and SE extended sources.

Fig. 2 shows a deviation probability map⁷ that we obtained with the best-fit model. Following the prescriptions in Abdollahi et al. (2022), we set the PS map parameters to p₀ = 4°, p₁ = 0.9 and p₂ = 0.1°. We also produced another PS map with parameters optimized for extended sources, that is, p₀ = 5°, p₁ = 0.8, and p₂ = 0.5°. No significant deviations (larger than ±2σ) were found in the central part of the RoI in this case either.

Fig. 1 Templates for the spatial modeling of SNR W44. (a-b) Elliptical ring (whole and divided) from Abdo et al. (2010). (c-d) Best full-ellipse (whole and divided). (e) Radio template (1420 MHz) from the THOR survey (Beuther et al. 2016).

Fig. 2 Fermi-LAT deviation probability map, or PS map, of the W44 surroundings with the best model obtained from the likelihood analysis. The white contours represent the W 44 radio template adopted in the analysis, derived from the template reported in Fig. 1e). The red crosses and circles show the new point-like and extended sources we added and fit in this analysis, and the black crosses and circles correspond to sources in the 4FGL-DR2 Fermi-LAT catalog for which the morphology was not changed. The radii of the circles represent the r68 of the extended sources that were modeled with a disk shape.

2.3 Diffuse disk morphology

We also investigated the morphology of the diffuse emission in the surroundings of W44 in detail. In order to find possible associations with gas complexes, we studied ¹²CO and ¹³CO data from the FOREST Unbiased Galactic plane Imaging survey with the Nobeyama 45-m telescope (FUGIN) (Umemoto et al. 2017) in a region within ~1° from the SNR. We chose nine velocity intervals between 30 km/s and 90 km/s from the emission profiles of the gas in correspondence of the large disk and in correspondence of the NW and SE disks. The velocity intervals corresponded to peaks of the gas profile within each of these regions. For each case, we derived a template and repeated the likelihood-fitting procedure for Fermi-LAT data. To do this, we adopted the radio template for the W44.

We refit the other sources in the RoI, including the NW and SE sources, spatially and spectrally in each case. We used the AIC value to decide which was the best configuration. Of the CO templates we derived, the template that provided the best results corresponded to the velocity interval (38.6–49.8) km/s, which is shown in Fig. 3. However, the large disk derived previously (and reported in Table 2) still provided the best description of the RoI, with a ΔAIC=10.6 in favor of the disk hypothesis. Therefore, we also used this latter model as reference spatial model for the MAGIC analysis.

When the CO template is used instead of the large disk to model the large diffuse background, the spatial and spectral parameters of the NW and SE sources did not change significantly. This is compatible with the results in Table 2 within the statistical uncertainties.

Fig. 3 CO template derived from FUGIN survey in the velocity interval (38.6, 49.8) km/s. The template is the sum of the 12CO and 13CO data. For comparison, the white contours and circles represent the W44 radio contours and the extended sources derived from the Fermi-LAT analysis.

2.4 Spectral analysis

For the low-energy spectral analysis (E> 100 MeV), we followed the weighted likelihood approach. This choice was driven by the necessity of accounting for systematic errors, whose main origin is the diffuse background emission. In particular, the Galactic interstellar emission model is affected by systematic uncertainties, which shows imperfections that stem from assumptions made in its construction and uncertainties in the adopted templates. This is a crucial point in particular for the source analysis we performed over energies below a few hundred MeV, as in the case of this work. In this regime, where the PSF is several degrees wide, the source-to-background ratio is low (at the percent level), and the systematic errors on the background model are critical. We carried out the weighted likelihood analysis following the procedure reported in Ballet (2016) and in Appendix B of Abdollahi et al. (2020).

We performed this analysis by fixing the spatial models for the sources of interest to those obtained from the high-energy analysis described above. Since for this spectral analysis we considered a wider energy range, we repeated the curvature fitting procedure on the newly added sources within 1° from the center of the RoI, and we reevaluated the T S_curv for each of them. We also repeated the fit of the spectral parameters of the two extended sources Kes79 and HESS J1857+026. The best spectral model for the diffuse disk is the LP, as shown in Table 3, while BPL is the best model for the NW and SE sources. In Tables 4 and 5 we report the best-fit spectral parameters for the W44 SNR, which has an LP spectral model, and for the three close sources. We evaluated the systematic errors from the uncertainties of the effective area of the instrument. The resulting spectral energy distribution (SED) for W44 is shown in Fig. 4 and the SEDs related to the sources NW and SE are displayed in Fig. 5. We added the systematic errors in quadrature to the errors of the respective source.

In order to probe the presence of a bump in the spectra of the sources NW and SE, which might arise at energies where the spectrum of the W44 SNR starts to decrease, we performed an analysis over a shorter energy range, 300 MeV–30 GeV, and repeated the curvature test, now focused on how well a BPL model can describe the spectra of sources W44, NW, and SE. Fitting a BPL spectral model to the spectrum of the W44 SNR, we obtained a significance T S_curv = 397.1 and a value of E_br = (1.7 ± 0.1) GeV. For source NW, the BPL significance is T S_curv = 8.6 with E_br = (0.47 ± 0.07) GeV, and for source SE, the BPL significance is T S_curv = 22.1 with E_br = (5.4 ± 1.2) GeV.

Fig. 4 SED of SNR W44 (only Fermi-LAT data). The solid line shows the best-fit curve obtained with an LP model, and the dashed lines represent the 1 σ uncertainty.

2.5 Contribution of the Galactic diffuse component to sources NW and SE and to the diffuse disk

We evaluated the contribution of the Galactic diffuse background from the model adopted for the Fermi-LAT analysis (gll_iem_v07) inside the three circular regions identified by sources NW and SE and the diffuse disk. The resulting plots are reported in the Appendix A. In the circular region corresponding to the diffuse disk model (Fig. A.1.3), the Galactic diffuse background dominates all other contributions to the emission, including the diffuse disk source itself. By contrast, the NW and SE sources are larger than the Galactic diffuse background in their corresponding regions. The diffuse disk spectrum is characterized by a cutoff similar to that of the Galactic diffuse model. This suggests that the diffuse disk emission, which corresponds to roughly 32% of the originally fit Galactic diffuse flux, could also be produced by interaction of the Galactic CR sea with gas in the region. However, we do not rule out other possible origins of this emission.

When we assume the same percentage (~32%) as the contamination of the Galactic diffuse background in the NW and SE regions, we obtain that only 10% (20%) of the total emission of NW (SE) can be explained, suggesting that most of the emission has a different origin. A more comprehensive discussion is presented in Sect. 4.

Fig. 5 SEDs of the two sources NW and SE with Fermi-LAT data and MAGIC ULs.

3 MAGIC observations and data analysis

The MAGIC telescopes are a system of two Imaging Atmospheric Cherenkov Telescopes located on the Canary island of La Palma at an altitude of 2200 m above sea level (Aleksić et al. 2016b). MAGIC observed the W44 region between April 2013 and August 2014 for 173.7 h after quality cuts at zenith angles between 25°–45°. The MAGIC observations were centered on the coordinates of the NW source from Uchiyama et al. (2012) using a wobble distance of $0_{\cdot }^{\circ }4$. The low-level analysis was performed with the MAGIC Analysis and Reconstruction Software (MARS) (Zanin 2013; Aleksić et al. 2016a). We used the spatial likelihood analysis package SkyPrism (Vovk et al. 2018) for our subsequent high-level analysis. An estimation of the systematic uncertainties of the MAGIC telescopes can be found in Aleksić et al. (2016a) and Vovk et al. (2018). The uncertainties were directly calculated from the log-likelihood deviation from the obtained best fit.

From the wobble observations, we constructed the background camera exposure model using an exclusion map. We excluded the known sources in our field of view: namely, HESS J1857+026 (Abdalla et al. 2018) using an exclusion zone radius of $0_{\cdot }^{\circ }45$, which is large enough to enclose the two known high-energy sources MAGIC J1857.2+0263 (RA: 18^h57^m${13}_{\cdot }^{\mathrm{\prime }\mathrm{\prime }}0$, Dec: 02°37′31″) and MAGIC J1857.6+0297 (RA: 18^h57^m${35}_{\cdot }^{\mathrm{\prime }\mathrm{\prime }}6$, Dec: 02°58′02″) (Aleksić et al. 2014). We also excluded HESS J1858+020 (Abdalla et al. 2018) with a radius of $0_{\cdot }^{\circ }17$, NW, and SE with their position and extension from the Fermi-LAT analysis above. The centers of the exclusion regions around HESS J1857+026 and HESS J1858+020 were set to the respective maximum in our sky maps. Due to the curved spectra of W44 and the large background components, these sources did not need to be included in the MAGIC energy range. The relative flux map is presented in Fig. 6. No significant detection was obtained with MAGIC. Therefore, we calculated 95% CL upper limits (ULs) for the NW source in each of our energy bins (Fig. 5) following the method described by Rolke et al. (2005). The lowest-energy UL was lower than 3% of the residual background, and we therefore fixed it at this higher level. This is equivalent, for instance, to the analysis of the gravitational wave event GW151226 in Berti (2018). We did not extract ULs for W44 because they would be more than an order of magnitude above the model curve. The SE source is smaller in radius than NW and closer to the edge of the MAGIC field of view. This means that if we were to estimate ULs for this area, they would be even higher than those calculated for the NW region. Therefore, we did not extract ULs for SE to avoid complicating the likelihood fitting unnecessarily.

Fig. 6 Relative flux map of the region observed with the MAGIC telescopes. The modeled sources MAGICJ1857.2+0263, MAGIC J1857.6+0297, and HESS J1858+020 are marked with white diamonds. The center of W44 is marked with a green cross. We also show the location of HESS J1857+026 from Abdalla et al. (2018) with a dark blue triangle. The NW center position is marked with a yellow star, and the extension is indicated by the yellow circle. The 39% and 68% containment contours of the MAGIC PSF are depicted in the bottom left corner.

4 Modeling the gamma-ray emission

In this section, we show that the gamma-ray emission from W44 and its surrounding region can be understood in the framework of particle acceleration and escape from the forward shock.

We note that the NW and SE sources are both located along the major axis of W44, suggesting that the emission from these two sources might be due to accelerated CRs that escape from the W44 in the direction of the local magnetic field and that interact with the interstellar gas. This interpretation is supported by a recent analysis from Liu et al. (2022), who used the velocity gradient technique (VGT) to infer the direction of the large-scale magnetic field⁸. The VGT was applied to formyl ion (HCO⁺), atomic hydrogen (HI), and CO emission lines from MCs that interact with W44, showing that the local magnetic field is oriented along the major axis of the SNR. This result also agrees with radio polarization measurements in the SE region (Goedhart et al. 2024) and with the Plank data at 353 GHz, which measures the dust polarization, at least for regions with high-intensity emission (see discussion in Liu et al. 2022, Sect. 4.2). Hence, the absence of gamma-ray emission from other locations around the SNR could just be the consequence of anisotropic escape of particles in the direction of the local magnetic field, as reported by Peron et al. (2020). A similar scenario, originally proposed by Gabici et al. (2009), was applied to W28 (Nava & Gabici 2013), even though in their case, there was no indication about the direction of the external magnetic field.

Alternatively, the emission from the two sources might result from the interactions of the Galactic CR sea with target gas that is not accounted for in the model for the diffuse emission (Acero et al. 2016). The mass estimated in the diffuse model is affected by an uncertainty of 30% mainly due to the CO-to-H₂ conversion factor (Bolatto et al. 2013). To check this scenario, we estimated the mass of the gas needed to produce a flux compatible with the Fermi-LAT data from NW source assuming that the Galactic CR spectrum is equal to the local spectrum measured by the Alpha Magnetic Spectrometer (AMS)-02 (Aguilar et al. 2015). The corresponding plot is shown in Fig. 7, where we also report MAGIC upper limits. The amount of mass required to fit Fermi-LAT data is ~5.5 × 10⁴ M_⊙. This value corresponds to ~30% of the total mass included in the diffuse model in the NW region integrated along the line of sight. Hence, in principle, the additional mass would be compatible with the mass uncertainty in the diffuse model. The total gamma-ray emission integrated along the line of sight also depends on the mass location, however. The total NW source emission can only be explained when the entire mass is concentrated at the NW source location, resulting in a gamma-ray emission that is three times higher than is predicted from the diffuse model (see Appendix A). In addition, the predicted spectrum does not reproduce the harder Fermi-LAT spectrum at ~10 GeV and slightly overshoots the MAGIC upper limits.

In order to estimate the goodness of this model, we fit the overall flux of the predicted model by keeping the spectral shape fixed. We first only used Fermi-LAT data and then added the MAGIC spectral points. In the latter case, we only used LAT data points up to 100 GeV because the MAGIC upper limits are more constraining for energies above 100 GeV. We calculated the χ² value of the best-fit models in the two cases and report them in Table 6 together with the corresponding p-values. We concluded that the model can be rejected at the ~3.6σ confidence level using Fermi-LAT data alone and at the ~5.8σ confidence level with Fermi-LAT + MAGIC data. For this reason, a scenario in which the gamma-ray emission is due to interactions with the Galactic CR sea is disfavored. In the next section, we focus on the escape scenario of freshly accelerated CRs.

Fig. 7 Gamma-ray emission from source NW compared to the pp emission due to Galactic CRs from a target mass equal to 5.5 × 104 M⊙.

4.1 Model for particle acceleration and escape

We adopted the model developed in Celli et al. (2019) that was applied to γ-Cygni (Acciari et al. 2023) and to the Cygnus Loop (Loru et al. 2021). This model assumes spherical symmetry for acceleration and escape. However, as already noted, we wished to account for particles that only escape along the magnetic flux tube oriented SE-NW. To first approximation, this scenario can be described by our spherical model assuming that all external particles located at distance r from the SNR forward shock are located at the same distance, but in a cylinder with a radius equal to the shock radius. The following equations and assumptions are based on the model described in Celli et al. (2019).

For the SNR properties, we adopted the values that are most commonly used in the literature. The SN explosion kinetic energy was fixed to E_SN = 10⁵¹ erg, the circumstellar density was n₀ = 10 cm⁻³, and the SNR age was t_age = 20 kyr (see the SNRcat9 from Ferrand & Safi-Harb 2012). The distance was estimated in the range 2.1-3.3 kpc; we assumed d = 2.2 kpc following Peron et al. (2020). We then followed Truelove & McKee (1999) to describe the SNR evolution assuming an ejecta power-law index equal to 0. The ejecta mass was fixed to 5 M_⊙ to match the average shock radius of $0_{\cdot }^{\circ }29$, which corresponds to R_sh ≃ 11 pc at a distance of 2.2 kpc.

We accounted for particle acceleration at the forward shock assuming that a fixed fraction of the shock kinetic energy is transferred to nonthermal particles. The instantaneous proton spectrum accelerated at the shock is $$f_{p,0}(p,t)=\frac{3{\xi }_{\mathrm{C}\mathrm{R}}{u}_{\mathrm{s}\mathrm{h}}^2(t){\rho }_0}{4\pi c{\left(m_{p}c\right)}^4\mathrm{\Lambda }(p_{max}(t))}{\left(\frac{p}{m_{p}c}\right)}^{-s}{e}^{-p/p_{max}(t)},$$(1)

where ρ₀ = n₀m_p (m_p is the proton mass), Λ(p_max(t)) is a normalization constant such that the CR pressure at the shock is $P_{\mathrm{C}\mathrm{R}}={\xi }_{\mathrm{C}\mathrm{R}}{\rho }_0{u}_{\mathrm{s}\mathrm{h}}^2,u_{\mathrm{s}\mathrm{h}}$ being the shock velocity. The factor ξ_CR represents the instantaneous acceleration efficiency and was kept constant during the whole evolution of the SNR. The spectral slope s was taken as free parameter and expected to be close to 4. p_max(t) represents the instantaneous maximum momentum (and E_max is the corresponding maximum energy) achieved at time t. Simple considerations on particle confinement suggest that p_max decreases during the Sedov–Taylor (ST) phase (Gabici et al. 2009), and for simplicity, it was assumed to follow a simple power law, $$p_{max}(t)=\{\begin{array}{ll}{p}_{\mathrm{M}}\left(t/t_{\mathrm{S}\mathrm{T}}\right)& \text{ if }t⩽t_{\mathrm{S}\mathrm{T}}\\ p_{\mathrm{M}}{\left(t/t_{\mathrm{S}\mathrm{T}}\right)}^{-\delta }& \text{ if }t>t_{\mathrm{S}\mathrm{T}}\end{array},$$(2)

where t_ST ≃ 420 yr is the beginning of the ST age. The absolute maximum momentum p_M and the slope δ are free parameters that are constrained from the data. The maximum energy was assumed to increase linearly during the ejecta-dominated phase (t < t_ST). However, the final spectrum at the present age (t_age) of the SNR is not very sensitive to the initial phase of the remnant evolution because t_ST ≪ t_age. Hence, the behavior of the maximum energy during the ejecta-dominated phase (i.e. for t < t_ST) cannot be constrained.

Particles with momentum p > p_max(t) escape from the shock and start to diffuse with a diffusion coefficient that we assumed to be proportional to the Galactic one, that is, D_ext = χD_gal, where D_gal(p) = 3.6 × 10²⁸β(p c/GeV)^1/3 cm²s⁻¹.

In summary, the model had seven parameters, ξ_CR, s, p_M, δ, and χ, plus the two masses of the clouds, M_NW and M_SE, that were constrained using the gamma-ray data. The gamma-ray emission from proton-proton interaction was calculated using the parameterization provided by Kafexhiu et al. (2014). The emission has two components: the first one comes from the interior of the SNR, while the second one comes from the exterior and diffuses mainly along the magnetic flux tube. Fig. 8 shows the result of our model compared to data from W44, NW and SE. The emission from W44 is well fit using ξ_CR = 1.3%, s = 4.2, p_M = 100 TeV and δ = 2. As noted by other authors (Ackermann et al. 2013) and as we also found in Sect. 2.4, the spectrum from W44 has a break at a few GeV. In our scenario, this break corresponds to protons with an energy equal to E_max(t_age) = 44 GeV. Above this energy, the spectrum steepens because particles escape from the remnant interior.

It is quite remarkable that the emission from the SE source shows a bump at roughly the same energy at which W44 shows a break, suggesting that the peak is due to escaping particles with p ≳ p_max(t_age). The NW source instead does not show clear evidence of a similar bump, which might have been smoothed due to the broken-shock scenario explained in Sect. 4.2, however. To reproduce the observed flux at a few GeV, we needed to assume cloud masses of M_NW = 9.9 × 10³ M_⊙ and M_SE = 8.7 × 10³ M_⊙ and a diffusion coefficient around the SNR that is suppressed by a factor of 5 compared to the average Galactic diffusion coefficient, namely χ = 0.2. For higher values of D_ext, the gamma-ray flux above ~4 GeV from the W44 interior and the two clouds would be suppressed because particles would escape faster from the region. However, the predicted gamma-ray spectra below ~4 GeV are suppressed because particles with pc ≲ 44 GeV are still confined inside the SNR. The escape model in this basic version is therefore unable to account for the emission from the NW and SE clouds below a few GeV. We tried to include an additional component due to preexisting Galactic CRs (dashed pink lines in Fig. 8) that interact with the cloud masses. In principle, this component is already accounted for in the background model. However, if the total mass along the line of sight is underestimated, a residual component may still be present. By adding this component, the discrepancy with the data is reduced (especially for SE), but is still unsatisfactory.

Before we proceed, we stress that in the model we proposed, the gamma-ray emission is dominated by freshly accelerated particles, while other works invoked the reacceleration of preexisting CRs (Uchiyama et al. 2012; Cardillo et al. 2016). The main difference between these two approaches is connected to the average density of the circumstellar medium in which the shock expands. For a density ≳100 cm⁻³, the shock becomes radiative, and the combination of reacceleration plus compression is enough to account for the gamma-ray emission. Conversely, we assumed a much lower density, n₀ ≃ 10 cm⁻³, and relied on the results obtained by Peron et al. (2020) through the analysis of the CO emission. It is known that the circumstellar medium around W44 contains some denser clumps that interact with the shock, and these easily enter the radiative phase (Cosentino et al. 2018). It is currently unclear which the filling factors of these regions are, however.

Fig. 8 Gamma-ray emission from sources NW (top panel) and SE (bottom) compared to our model. The dashed gray line shows the emission from the high-energy particles escaping from the shock alone. The dashed pink line also includes an additional component due to Galactic CRs interacting with the cloud material, and the dashed orange line shows the sum of the dashed gray line with the contribution from the low-energy particles escaping from the broken-shock region. For comparison, the black points and dashed lines refer to the W44 emission.

4.2 The broken-shock scenario

From the previous section, it is clear that the escape model fails to explain the lowest energies from the external sources, especially for the NW source, where the detected emission is higher by one order of magnitude at least than the model prediction. For this reason, we propose a novel scenario in which some number of particles with an energy lower than E_max(t_age) is allowed to escape from the system. If the shock were a continuous surface, the escape of low-energy particles would not be possible, but in the case of middle-aged SNRs that expand in highly inhomogeneous media, the shock surface could be broken where the upstream medium is denser, resulting in the escape of a small fraction f_br of the total accelerated particles. As a consequence, even low-energy particles might fill the region outside of the shock up to a distance equal to their diffusion length, $l_{\text{diff }}=\sqrt{2D_{\text{ext}}t}$. If the emitting region has a size comparable to or smaller than l_diff, we may expect gamma-ray emission also at an energy lower than ~0.1 E_br. We refer to this model as the broken-shock scenario.

The precise modeling of this scenario is beyond the scope of the present work because it depends on quantities that are not known with sufficient accuracy. In particular, we need to know the density distribution of the circumstellar medium to describe the time evolution of the shock propagating in such a medium.

We limit our discussion here instead to evaluate whether the broken-shock scenario is compatible with the model presented in the previous section in terms of the number of escaped particles, diffusion length, and maximum energy. In particular, concerning the latter, we verified that the fragmentation of the shock surface does not reduce the maximum energy that can be achieved in the system. In other words, the shock surface should be continuous for a typical size l_sh larger than the upstream diffusion length, D(p)/u_sh, of particles with a momentum equal to p_max = 44 GeV/c, as estimated in the previous section. We assumed that the diffusion in the shock region is due to magnetic turbulence that is produced through the resonant streaming instability, such that the diffusion coefficient can be written as D = r_Lc/(3ℱ), where r_L(p) is the Larmor radius, and ℱ is the power of the magnetic turbulence for a logarithmic bandwidth. In quasi-linear theory, this latter quantity is given by (see, e.g. Blasi 2013) $$\mathcal{F}=\frac{\pi }{2}\frac{{\xi }_{\mathrm{C}\mathrm{R}}}{\ln \left(p_{max}/m_{p}c\right)}\frac{u_{\mathrm{s}\mathrm{h}}}{v_{\mathrm{A}}},$$(3)

where v_A = B₀/(4πn₀m_pc)^1/2 is the Alfvén speed upstream of the shock, and B₀ is the large-scale magnetic field. Hence, the condition l_sh > D(p)/u_sh reads $$l_{\mathrm{s}\mathrm{h}}>\frac{2}{3\pi }\frac{\ln \left(p_{max}/m_{p}c\right)}{{\xi }_{\mathrm{C}\mathrm{R}}}\frac{r_L\left(p_{max}\right)c{v}_A}{u_{\mathrm{s}\mathrm{h}}^2}\simeq 0.01\mathrm{p}\mathrm{c},$$(4)

where the last equality was obtained using n₀ = 10 cm⁻³. Interestingly enough, this result does not depend on the value of B₀. The lower limit for the shock size is very low and should not be a limiting condition. Detailed radio maps suggest that l_sh is even much larger (Cotton et al. 2024; Goedhart et al. 2024).

We made the simplifying assumption that in addition to the escape of high-energy particles with E > E_max(t_age), a fraction f_br of low-energy particles (i.e., E < E_max(t_age)) also escapes from the SNR interior and propagates in the circumstellar medium according to the diffusion coefficient D_ext estimated in the previous section. The value of f_br was then determined from the comparison with the gamma-ray data. The dashed orange lines in Fig. 8 show the total gamma-ray emission resulting from the broken-shock scenario assuming a value f_br = 0.2 and 0.1 for NW and SE, respectively. Because a relatively small fraction of 10–20% allows us to increase the expected flux at low energies to levels similar to the measured ones, this option can be considered a possible explanation for the observed flux. The estimated value of f_br is low enough to be considered a weak perturbation to the total number of particles, such that the gamma-ray emission from the W44 interior is not strongly affected. To recover the same gamma-ray flux as in the case of unbroken shock, it is enough to increase the acceleration efficiency by 10–20%. In this scenario, the required masses for the two clouds are similar to the previous estimate, M_NW = 5.4 × 10³ M_⊙ and M_SE = 8.7 × 10³ M_⊙.

As we did for the Galactic CR model (see Sect. 4), we calculated the χ²/N value using the best-fit model and including Fermi-LAT up to 100 GeV + MAGIC data for NW source and Fermi-LAT data up to 100 GeV for SE source. We found values of 33.6/28 and 23.1/22 for NW and SE, respectively, which confirms the good agreement between data and model.

Finally, we verified that the time needed for these particles to fill the two external regions has to be much shorter than the age of W44. This time can be estimated as t_fill ≃ (2R_SRC)²/(2D_ext(p)). For 1 GeV particles, it ranges between 4300 yr and 2400 yr for NW and SE sources, respectively. This corresponds to at most 20% of the age of W44 and is consistent with the idea that the shock only starts to break in the very last phase of the SNR life.

5 Summary

We reported the results of Fermi-LAT and MAGIC observations of the middle-aged SNR W44 and its surrounding environment. Through the Fermi-LAT high-energy analysis, we studied the morphology of the remnant in detail, finding that it resembles an ellipsoid. This confirms the results of previous studies. In particular, a radio template derived from the THOR survey provided the best description of the morphology of the SNR. Moreover, we confirmed the detection of two sources located at the NW and SE edges of the remnant shell. The extensions of the sources NW and SE are $$r_{68}=0_{\cdot }^{\circ }16\pm 0_{\cdot }^{\circ }01$$ and $r_{68}=0_{\cdot }^{\circ }12\pm 0_{\cdot }^{\circ }01$, respectively. The smaller extension obtained for the NW source with respect to previous studies is probably due to the presence of a large diffuse source, which we modeled as a large central disk. We also derived the best spectral models, starting from 100 MeV, for W44 SNR (described with an LP) and for the nearby sources (an LP for the large diffuse disk and a BPL for sources NW and SE).

The observations with the MAGIC telescopes focused on the northwestern region of the W44. We carried out an exclusion region analysis considering as the position and extension of the NW source those derived from the Fermi-LAT high-energy analysis. Since no significant detection was obtained, we derived 95% CL ULs in the energy range above 100 GeV.

We note that no sources coincident with W44, NW, or SE were identified in the first LHAASO catalog of gamma-ray sources (Cao et al. 2024) or in the third catalog of very high-energy gamma-ray sources of the High-Altitude Water Cherenkov gamma-ray Observatory (HAWC) (Albert et al. 2020).

The combined Fermi-LAT and MAGIC spectra of the NW source allowed us to reject a model in which this emission is entirely due to diffuse Galactic CRs. We developed a model to describe the W44, NW, and SE source emissions in the hypothesis of CRs that are freshly accelerated at the remnant shock and escape to the NW and SE clouds. This interpretation was based on the assumption that NW and SE sources are located at the same distance as W44. However, to our knowledge, there are no other known sources along the same line of sight that might cause the observed NW and SE emissions. The spectra of the W44 SNR and the two close-by sources were described with a time-dependent model based on the assumption that CRs are accelerated at the SNR forward shock with a constant efficiency of ~1% up to a maximum energy that decreases in time as E_max ∝ t⁻². At the current age of W 44, the maximum energy is 44 GeV. At any given time, particles with energies >E_max(t) escape anisotropically from the forward shock, following the magnetic field lines along which the NW and SE sources are located. The emission from these sources was then interpreted as due to hadronic interactions of escaping protons with the ambient gas, whose total mass has to be ≈10⁴ M_⊙ for both sources. This model also requires that the diffusion coefficient around W44 has to be suppressed by a factor of about 5 with respect to the average Galactic diffusion coefficient. We note that the NW and SE fluxes can be reasonably well described assuming the same CR acceleration efficiency as is required to fit the gamma-ray flux of W44. Finally, we note that the anisotropic CR escape was proposed in the past for another SNR, namely W28. However, for W44, it is strongly supported by the recently measured direction of the large-scale magnetic field.

Gamma-ray data in the GeV range also indicate the presence of low-energy particles, down to ~1 GeV, in NW and SE sources. At first glance, this seems to contradict our escape model, in which only particles with E ≳ 44 GeV can escape. However, we proposed a novel approach, the broken-shock scenario’, to account for the escape of low-energy particles, in which a small fraction of the shock surface is destroyed (possibly by dense clumps), allowing the escape of particles at all energies. Compared to the total number of accelerated particles, the number of particles that escape from the broken-hock surface has to be ≈10% and 20% for SE and NWsources, respectively. This fraction is small enough for the acceleration process to not be significantly affected. This model proved to agree reasonably well with Fermi-LAT spectral points and MAGIC ULs.

Acknowledgements

Contributions of the authors: R. Di Tria: Fermi-LAT analysis, paper drafting; L. Di Venere: Project coordination, Fermi-LAT analysis, paper drafting; D. Green: Project coordination, CO data analysis, paper drafting; A. Hahn: MAGIC analysis, paper drafting; G. Morlino: theoretical modeling and interpretation, paper drafting; M. Strzys: MAGIC analysis cross-check, paper drafting; E. Bissaldi: paper drafting; S. R. Gozzini, A. López-Oramas: PIs of MAGIC observation campaigns; the rest of the authors have contributed in one or several of the following ways: design, construction, maintenance and operation of the instrument(s) used to acquire the data; preparation and/or evaluation of the observation proposals; data acquisition, processing, calibration and/or reduction; production of analysis tools and/or related Monte Carlo simulations; overall discussions about the contents of the draft, as well as related refinements in the descriptions. The Fermi-LAT Collaboration acknowledges generous ongoing support from a number of agencies and institutes that have supported both the development and the operation of the LAT as well as scientific data analysis. These include the National Aeronautics and Space Administration and the Department of Energy in the United States, the Commissariat ‘a l’Energie Atomique and the Centre National de la Recherche Scientifique /Institut National de Physique Nucleaire et de Physique des Particules in France, the Agenzia Spaziale Italiana and the Istituto Nazionale di Fisica Nucleare in Italy, the Ministry of Education, Culture, Sports, Science and Technology (MEXT), High Energy Accelerator Research Organization (KEK) and Japan Aerospace Exploration Agency (JAXA) in Japan, and the K. A. Wallenberg Foundation, the Swedish Research Council and the Swedish National Space Board in Sweden. Additional support for science analysis during the operations phase is gratefully acknowledged from the Istituto Nazionale di Astrofisica in Italy and the Centre National d’Etudes Spatiales in France. This work performed in part under DOE Contract DE-AC02-76SF00515. We would like to thank the Instituto de Astrofísica de Canarias for the excellent working conditions at the Observatorio del Roque de los Muchachos in La Palma. The financial support of the German BMBF, MPG and HGF; the Italian INFN and INAF; the Swiss National Fund SNF; the grants PID2019-104114RB-C31, PID2019-104114RB-C32, PID2019-104114RB-C33, PID2019-105510GB-C31, PID2019-107847RB-C41, PID2019-107847RB-C42, PID2019-107847RB-C44, PID2019-107988GB-C22, PID2022-136828NB-C41, PID2022-137810NB-C22, PID2022-138172NB-C41, PID2022-138172NB-C42, PID2022-138172NB-C43, PID2022-139117NB-C41, PID2022-139117NB-C42, PID2022-139117NB-C43, PID2022-139117NB-C44 funded by the Spanish MCIN/AEI/ 10.13039/501100011033 and “ERDF A way of making Europe”; the Indian Department of Atomic Energy; the Japanese ICRR, the University of Tokyo, JSPS, and MEXT; the Bulgarian Ministry of Education and Science, National RI Roadmap Project DO1-400/18.12.2020 and the Academy of Finland grant nr. 320045 is gratefully acknowledged. This work was also been supported by Centros de Excelencia “Severo Ochoa” y Unidades “María de Maeztu” program of the Spanish MCIN/AEI/ 10.13039/501100011033 (CEX2019-000920-S, CEX2019-000918-M, CEX2021-001131-S) and by the CERCA institution and grants 2021SGR00426 and 2021SGR00773 of the Generalitat de Catalunya; by the Croatian Science Foundation (HrZZ) Project IP-2022-10-4595 and the University of Rijeka Project uniri-prirod-18-48; by the Deutsche Forschungsgemeinschaft (SFB1491) and by the Lamarr-Institute for Machine Learning and Artificial Intelligence; by the Polish Ministry Of Education and Science grant No. 2021/WK/08; and by the Brazilian MCTIC, CNPq and FAPERJ. We acknowledge the contributions of J. Krause, V. Stamatescu, and P. Colin during the original observation campaigns.

Appendix A Galactic diffuse background contamination

Here we report the E²dN/dE of the Galactic diffuse emission evaluated from the model adopted for the Fermi-LAT analysis (Acero et al. 2016) as described in Sect. 2.5. We considered three disks centered on the coordinates and having the same extensions of the NW (Fig.A.1 top), SE (Fig.A.1 middle) and diffuse disk (Fig.A.1 bottom) sources. The SEDs of the three sources of interest are also shown in the same figures. Also in this case the errors are given by the quadrature sum of the statistical and systematic errors, the latter evaluated taking into account the uncertainties on the effective area.

Fig. A.1 Galactic diffuse background E2dN/dE (dashed orange line) evaluated in correspondence of the sources NW (top panel), SE (middle panel) and Diffuse Disk (bottom panel).

References

All Tables

All Figures

	Fig. 1 Templates for the spatial modeling of SNR W44. (a-b) Elliptical ring (whole and divided) from Abdo et al. (2010). (c-d) Best full-ellipse (whole and divided). (e) Radio template (1420 MHz) from the THOR survey (Beuther et al. 2016).
In the text

Fig. 2 Fermi-LAT deviation probability map, or PS map, of the W44 surroundings with the best model obtained from the likelihood analysis. The white contours represent the W 44 radio template adopted in the analysis, derived from the template reported in Fig. 1e). The red crosses and circles show the new point-like and extended sources we added and fit in this analysis, and the black crosses and circles correspond to sources in the 4FGL-DR2 Fermi-LAT catalog for which the morphology was not changed. The radii of the circles represent the r68 of the extended sources that were modeled with a disk shape.

In the text

	Fig. 3 CO template derived from FUGIN survey in the velocity interval (38.6, 49.8) km/s. The template is the sum of the 12CO and 13CO data. For comparison, the white contours and circles represent the W44 radio contours and the extended sources derived from the Fermi-LAT analysis.
In the text

	Fig. 4 SED of SNR W44 (only Fermi-LAT data). The solid line shows the best-fit curve obtained with an LP model, and the dashed lines represent the 1 σ uncertainty.
In the text

	Fig. 5 SEDs of the two sources NW and SE with Fermi-LAT data and MAGIC ULs.
In the text

Fig. 6 Relative flux map of the region observed with the MAGIC telescopes. The modeled sources MAGICJ1857.2+0263, MAGIC J1857.6+0297, and HESS J1858+020 are marked with white diamonds. The center of W44 is marked with a green cross. We also show the location of HESS J1857+026 from Abdalla et al. (2018) with a dark blue triangle. The NW center position is marked with a yellow star, and the extension is indicated by the yellow circle. The 39% and 68% containment contours of the MAGIC PSF are depicted in the bottom left corner.

In the text

	Fig. 7 Gamma-ray emission from source NW compared to the pp emission due to Galactic CRs from a target mass equal to 5.5 × 104 M⊙.
In the text

Fig. 8 Gamma-ray emission from sources NW (top panel) and SE (bottom) compared to our model. The dashed gray line shows the emission from the high-energy particles escaping from the shock alone. The dashed pink line also includes an additional component due to Galactic CRs interacting with the cloud material, and the dashed orange line shows the sum of the dashed gray line with the contribution from the low-energy particles escaping from the broken-shock region. For comparison, the black points and dashed lines refer to the W44 emission.

In the text

	Fig. A.1 Galactic diffuse background E2dN/dE (dashed orange line) evaluated in correspondence of the sources NW (top panel), SE (middle panel) and Diffuse Disk (bottom panel).
In the text