Studying the nature of the unidentified gamma-ray source HESS J1841-055 with the MAGIC telescopes

We investigate the physical nature and origin of the gamma-ray emission from the extended source HESS J1841-055 observed at TeV and GeV energies. We observed HESS J1841-055 at TeV energies for a total effective time of 43 hours with the MAGIC telescopes, in 2012 and 2013. Additionally, we analysed the GeV counterpart making use of about 10 years of Fermi-LAT data. Using both Fermi-LAT and MAGIC, we study both the spectral and energy-dependent morphology of the source for almost four decades of energy. The origin of the gamma-ray emission from this region is investigated using multi-waveband information on sources present in this region, suggested to be associated with this unidentified gamma-ray source. We find that the extended emission at GeV-TeV energies is best described by more than one source model. We also perform the first energy-dependent analysis of the HESS J1841-055 region at GeV-TeV. We find that the emission at lower energies comes from a diffuse or extended component, while the major contribution of gamma rays above 1 TeV arises from the southern part of the source. Moreover, we find that a significant curvature is present in the combined observed spectrum of MAGIC and Fermi-LAT. The first multi-wavelength spectral energy distribution of this unidentified source shows that the emission at GeV-TeV energies can be well explained with both leptonic and hadronic models. For the leptonic scenario, bremsstrahlung is the dominant emission compared to inverse Compton. On the other hand, for the hadronic model, gamma-ray resulting from the decay of neutral pions ($\pi^0$) can explain the observed spectrum. The presence of dense molecular clouds overlapping with HESS J1841-055 makes both bremsstrahlung and $\pi^0$-decay processes the dominant emission mechanisms for the source.


INTRODUCTION
The unidentified gamma-ray source HESS J1841-055 was first discovered at TeV energies in 2007 by the High Energy Stereoscopic System (H.E.S.S.) during the Galactic plane survey (Aharonian et al. 2008). The observed emission was reported as extended with an elliptical extension of 0.41 • and 0.25 • along the semi-major and semi-minor axes, respectively and centered at Right Ascension (RA):18 h 40 m 55 s and declination (Dec): 5 • 33 00 with a position angle 39 • relative to the RA axis. HESS J1841-055 was detected with a statistical significance of 10.7σ and a flux of (12.8±1.3) × 10 −12 cm −2 s −1 between 0.54 and 80 TeV. The spectrum is best described by a powerlaw with a spectral index of 2.41±0.1 stat ±0.2 sys . These results are compatible with the recent results reported by H.E.S.S. collaboration (H. E. S. S. Collaboration et al. 2018). Using the ARGO-YBJ experiment for energies above 0.9 TeV, Bartoli et al. (2013) reported a similar extension as seen be the H.E.S.S collaboration but a 3 times larger flux due to differing background estimation techniques between the experiments. This region was further investigated by the HAWC observatory also at TeV energies. The source 1HWC J1838-060, from the First HAWC Catalog, was detected at 6.1σ post-trial significance. It is located in the middle of HESS J1841-055 and another known TeV source, HESS J1837-069 (Abeysekara et al. 2016). This detection by HAWC was found to be overlapping with the extension of HESS J1841-055, and the differential flux normalization was compatible with the one reported by the H.E.S.S. collaboration. The second HAWC Catalog also revealed a source, 2 HWC J1837-065 which was likely to be associated to HESS J1841-055 (Abeysekara et al. 2017). Its spectral index varies from -2.90±0.04 for a point-like emission to -2.66±0.03 for a 2 • radius.
Corresponding authors: L. Saha: labsaha@ucm.es, A. López-Oramas: alicia.lopez@iac.es This region was further investigated at other wavelengths to search for possible counterparts. Although no confirmed counterparts of the TeV source HESS J1841-055 at lower energies are known, several possible associations have been suggested. The emission from HESS J1841-055 may be due to either a single extended source or several unresolved sources. Sguera et al. (2009), making use of INTEGRAL data, proposed as counterpart the unidentified transient source 3EG J1837-0423, which was likely to be associated to the Supergiant Fast X-ray Transient (SFXT) AX J1841.0-0536. At X-ray energies, observations of this extended region were done with SUZAKU and an X-ray source was discovered (Nobukawa et al. 2015). The detection of two separate extended sources (FGES J1839.4-0554 andFGES J1841.4-0514) was also reported in this region at energies above 10 GeV using data from the Fermi-Large Area Telescope (LAT, Ackermann et al. 2017;Ajello et al. 2017). Some potential sources at different wavelengths suggested to be associated with HESS J1841-055 is discussed later in detail.
In this paper, we study this complex region using dedicated observations with the MAGIC telescopes at TeV energies. We also explore the GeV counterpart making use of 10-year data of Fermi-LAT. We finally model the GeV-TeV emission to unveil the dominant gamma-ray emission mechanisms at work. The potential counterparts at other frequencies are also investigated. The low energy threshold of MAGIC, which allows to overlap with Fermi-LAT in the GeV domain, combined with the MAGIC capabilities of reaching several TeV, make the MAGIC telescopes a suitable instrument to study this region within a broad energy range. The combination of both MAGIC and Fermi-LAT allows spectral studies of this complex region for almost four decades in energy.
The paper is organized as follows: in section 2, we describe the detailed analyses of the MAGIC and Fermi-LAT data. The results are discussed in section 3. Potential counterparts are proposed in section 4. The multiwaveband modelling of the source is explained in section 5. Finally, we summarize and conclude in section 6.

MAGIC
Very-High-Energy (VHE, E> 100 GeV) gamma-ray observations of HESS J1841-055 are performed using the MAGIC telescopes. MAGIC consists of two 17 m diameter Imaging Atmospheric Cherenkov Telescopes (IACTs) located at the Observatory of Roque de los Muchachos (28 • .8 N, 17 • .9 W, 2200 m above the sea level) on the Canary Island of La Palma, Spain. The energy threshold of the MAGIC stereoscopic system is about 50 GeV, and it is able to detect ∼ 0.6% of the Crab Nebula flux above 250 GeV at 5σ significance in 50 hours of observations at small (<30 • ) zenith angles (Aleksić et al. 2012). HESS J1841-055 was observed between April 2012 and August 2013, for a total of about 43 hours, at zenith angles between 5 • and 50 • ,resulting in an energy threshold for this analysis of ∼150 GeV. To estimate the background simultaneously with the source data, the observations are performed in the so-called wobble-mode (Fomin et al. 1994) at two symmetrical positions, with the source located 0 • .55 off-axis from the center of the camera. After quality cuts, which account for hardware problems, unusual rates, and bad atmospheric conditions, ∼34 hours of high-quality, dark-time data are selected for further analysis. The analysis of the MAGIC data is performed using the standard MAGIC Analysis and Reconstruction Software (Mars; Moralejo et al. 2009;Zanin et al. 2013) and standard analysis procedure.
Given the extension of the source and the possibility of contamination from other nearby sources, we study the region using an iterative maximum likelihood method included in the Skyprism package (Vovk et al. 2018). Skyprism has specifically been developed to perform 2D fitting of IACTs data and has been optimised for MAGIC data. This set of tools compute the instrument response function (IRF) and perform a simultaneous maximum likelihood fit of source models of arbitrary morphology to the sky images. With Skyprism it is then possible to analyse MAGIC data of extended sources of arbitrary morphology and multiple, overlapping sources.
We compute the event count map, the background map, and the instrument response functions which includes Point Spread Function (PSF), energy migration matrix, and exposure map. We use the "exclusion map" method for generating the background map excluding a circular region of 0 • .5 around the center of the HESS J1841-055 (RA = 280 • .23, Dec=−5 • .55) and a circular region of 0 • .3 around a bright spot at the southern edge of the camera (RA = 279 • .4, Dec=−6 • .45). A user-defined source model (2D Gaussian) is used to fit the measured event maps for maximizing the log-likelihood estimate. To calculate the individual spectral parameters of the sources obtained from the modelling of the region, we use the maximum log-likelihood method, as defined in Vovk et al. (2018), assuming a powerlaw model for the source at the pivot energy 1 TeV (energy at which the uncertainty in the normalisation is minimum). The observations for this work are performed at low and medium zenith angles (Z<50 • ). Given the very high signal-to-noise ratio (>0.4), the systematic uncertainties can be considered similar to those reported in Aleksić et al. (2012), defined as 12% in the integral flux for stereoscopic observations.

Fermi-LAT
The Large Area Telescope (LAT) onboard the Fermi Gamma Ray Space Telescope allows for the detection of gamma rays from 30 MeV to > 500 GeV with its large effective area and wide field of view (Atwood et al. 2009). In our analysis, we select nearly ten years (i.e., from 2008 September 1 to 2017 May 5) of Pass 8 SOURCE class (P8R3) LAT events in the reconstructed energy of about 10 GeV to 1 TeV within a 15 • region of interest (ROI) around the fourth Fermi-LAT catalog source 4FGL J1840.9-0532e (associated to 3FHL J1840.9-0532e). TeV source HESS J1841-055 is associated with the Fermi-LAT source 4FGL J1840.9-0532e. The Fermi Science Tools (FST) analysis package 1 version v11r5p3 and the P8R3 − SOURCE − V2 IRFs are used for the analysis. We also use the python-based package Fermipy (version 0.17.4 2 ) to facilitate the analysis of data with the FSTs. We select photons of energies greater than 10 GeV with arrival direction within 105 • from local zenith to remove contamination from the Earth's emission. The PASS 8 source class allows for the use of four different event types which are based on the event-by-event quality of reconstructed direction (PSF) and energy. Hence, the data is separated into these event types to optimize selection of events based on the quality of reconstruction of direction of incoming photons and energy. The Galactic diffuse emission is modeled with the standard Fermi-LAT diffuse emission model (gll_iem_v07.fits). The isotropic emission from extragalactic radiation and residual background models (iso_P8R3_SOURCE_V2_PSF[0/1/2/3]_v1.txt) are also used corresponding to four event types.
We first start with a baseline sky model which includes all 4FGL point and extended sources within the ROI listed in the 4FGL catalogue 3 (Abdollahi et al. 2020). The extended source 4FGL J1840.9-0532e (associated to 3FHL J1840.9-0532e) is our source of interest in the ROI which is associated with two sources from the Fermi Galactic Extended Source Catalog (FGES) (Ackermann et al. 2017) and it is included in the model. Initially, we use the baseline model to optimize parameters of the sources by fitting their flux and spectral parameters. After the initial optimization, we remove all sources for which the values of the predicted number of counts in the model, Npred, are less than 2.0 and we free spectral shapes and normalizations for all the sources which lie within 3 • from the center of the ROI. The isotropic diffused background model is fixed to the value obtained after the first optimization of the ROI but the diffuse Galactic model is kept free for all different configurations or models discussed below. Then we use the binned maximum likelihood method to estimate the best-fit model parameters using a 15 • × 15 • square region centered on 4FGL J1840.9-0532e with a spatial bin-size of 0 • .06 and 10 equally spaced energy bins per decade of energy. We then relocate the source of interest using the maximum likelihood method to find the best source position. As the next step, we use an iterative maximum likelihood-based source finding algorithm to identify new point sources within 0.5 deg from the centre of the ROI. The algorithm finds point sources within this ROI with test-statistics 4 , TS > 16. We continue searching for new sources until all the point sources are added to the baseline model. Following this, we remove all the sources with TS < 16 from the ROI and perform the maximum likelihood method for the best-fit model parameters. 1 fermi.gsfc.nasa.gov/ssc/data/analysis/software/ 2 https://fermipy.readthedocs.io/en/latest/ 3 https://fermi.gsfc.nasa.gov/ssc/data/access/lat/8yr_ catalog/gll_psc_v19.fit 4 The Test Statistic (TS) of a source is evaluated using a likelihood ratio test defined as TS = −2log( L 1 L 0 ), where L 0 and L 1 are the likelihoods of the background model without the source (null hypothesis) and the hypothesis being tested (source plus background), respectively. The detection significance is approximately the square root of the TS.

Morphology
In order to study the energy-dependent morphology of the extended source HESS J1841-055, we produce skymaps for different energy ranges using Skyprism. Fig. 1 shows the relative flux skymaps, with 3σ and 5σ contours extracted from the TS map, produced for energies 50 -500 GeV (low energy, LE, map), 500 GeV -1 TeV (medium energy, ME, map) and >1 TeV (high energy, HE, map), respectively. The relative flux is defined as the excess events divided by the background events. To calculate the extension of the source in each of the energy ranges, we consider a radially symmetrical twodimensional (2D) Gaussian shape. The 1σ standard deviation of the Gaussian is considered to be the extension or radius of the source. The radius is kept as a free parameter and is allowed to change by 0 • .01 over a range of 0 • .1 to 0 • .6 during the maximum likelihood fitting. Moreover, we simultaneously keep changing the origin of the Gaussian by changing RA and Dec by 0 • .01 for both of them over a range of 1 • . The best-fit locations along with the extensions of the source for different energy ranges are shown in Table 1.
The extension of the source at these three energy ranges appears to be the same, however, the overall detection significance of the extended emission reduces at higher energies, revealing only a few hotspots in the southern part of the source (see Fig. 1). The fitted extension is the same (within errors) in the whole energy range (see Table 1). MAGIC observations show that the source has an extension compatible with that measured by H.E.S.S. collaboration at TeV energies. It is also evident from the different maps that the extended region shows several bright hotspots with a significance of more than 5σ. Many bright highly-significant spots are detected at LE and ME energies, while they mostly disappear at HE. These hotspots hint the presence of multiple sources in the region. It also indicates that the most significant emission at higher energies is coming from the southern part of the region. As discussed above, the extended source HESS J1841-055 may potentially consist of multiple sources. To check this, we consider three different source models covering the full energy range, i.e., energies from 50 GeV to above 1 TeV. We first consider a singlesource model where the extended source is considered to be a 2D elliptical Gaussian. We leave the peak position, extension along X and Y direction and angle w.r.t. the X direction free while maximizing the likelihood value of the fit. For the second model, we replace the single-source model with two sources which are modelled as 2D circular Gaussian. The peak position and radius (1σ standard deviation) of the two sources are free parameters of the model. Finally, for the third option, we model the entire source region considering three different sources, one with elliptical disk model and the other two with Gaussian models. The results of the maximum likelihood values are given in Table 2. It is found that both two-source model and three-source models are better than a single-source model. The im-  Table 3.

Spectrum
The spectral energy distribution (SED) is calculated in the energy range of 50 GeV to > 1 TeV, using the Skyprism package. We consider the extended 2D Gaussian template with the extension 0.4 • at the position of the HESS J1841-055 and an isotropic background. The assumed spectrum of the source is considered to follow a simple powerlaw (PL) model which is defined as follows: where N o and α are parameters of the model. The best-fit spectral parameters are N o = (9.43 ± 0.29) × 10 −12 TeV −1 cm −2 s −1 , α = 2.57 ± 0.05. The gamma-ray flux above 50 GeV is, F(> 50 GeV) = 2.23 × 10 −10 cm −2 s −1 . The SED measured by MAGIC is plotted in Fig. 2.
Although morphology studies reveal that the emission region can be modeled better with more than one source, we cannot make any robust estimate on the number of distinct sources due to limitations of the software tool. Hence we do not provide high-quality SEDs associated with these sources.

Morphology
For the morphological analysis of the source, photons with energy above 10 GeV up to 1 TeV are considered to reduce the contamination from nearby pulsars within the ROI. With the baseline model, as discussed in Section 2.2, we perform the binned maximum likelihood analysis on 4FGL J1840.9-0532e and find the best-fit model parameters. To estimate the size of 4FGL J1840.9-0532e, we calculate the TS of the extension (TS ext ) parameter, which is the likelihood ratio of the likelihood for being a point-like source (L pt ) to a likelihood for an assumed extension (L ext ), TS ext = 2log(L ext /L pt ). In order to test the extension of the source of interest, a radially symmetric Gaussian is considered and we vary its sigma from 0 • .01 to 1 • .5 in steps of 0 • .1. We also simultaneously leave the location of the center of the source free within 1σ extension of the Gaussian. We find that the source extension is 0 • .64 ± 0 • .11 with the T ext = 264 which corresponds to a significance of about 16σ. We also consider a radial disk model and found that the source extension is Table 3. Parameters of the best-fit single and multi-source models. The ext x and ext y are extension of the models along X and Y direction respectively. For elliptical Gaussian, they are standard deviation, whereas for elliptical disk model they correspond to semi-major and semi-minor axis respectively.
One   (red) whereas Fermi-LAT energy fluxes are obtained for energy above 10 GeV (blue). The combined fit of the Fermi-LAT and MAGIC SEDs is best described by either a BPL model (green shaded region) or a PLE (blue shaded region) model. 0 • .60 ± 0 • .11 with the T ext = 224, which corresponds to a significance of about 15σ. However, the log-likelihood is maximum for the Gaussian model, hence it will be considered as the preferred model for 4FGL J1840.9-0532e. The resulting Fermi-LAT TS map above between 10 GeV and 1 TeV is shown on Fig. 3.

Spectrum
For the spectral study, we consider data within the energy range 10 GeV-1 TeV. We calculate the SED of HESS J1841-055 using the best model obtained for the morphological study as discussed in Sec. 3.2.1. The SED of HESS J1841-055 is shown in Fig. 2, which is obtained by a fit to the data with the PL model. The best-fit PL model parameters are: prefactor, N 0 = (1.71 ± 0.41) × 10 −14 MeV −1 cm −2 s −1 , spectral index, α = 2.30 ± 0.03 and scale, E 0 = 1 GeV, where the uncertainties are statistical only. The total flux is found to be F(> 10 GeV) = (1.2 ± 0.1) × 10 −8 photons cm −2 s −1 .

Joint fit to MAGIC and Fermi-LAT data
We perform a joint likelihood fit to the observed fluxes from MAGIC and Fermi-LAT to find out the spectral behaviour of the source in the GeV-TeV energy range. We perform a χ 2 fit on the Fermi-LAT-MAGIC spectral points. We consider different spectral shapes as a PL, a PL with exponential cutoff (PLE) and a broken powerlaw (BPL) as spectral shapes for the fit. The PL has already been defined in subsection 3.1.2. The PLE spectral shape is defined as: The BPL model is defined as follows: where A, α 1 , α 2 and E break are parameters of the model. In the case of the BPL model, the spectral break is at 37 GeV, while the cutoff energy in the PLE model is located at 1.8 TeV. Both the BPL and PLE models describe the SED better than a simple PL model, implying that a significant curvature is present in the SED. However, both BPL and PLE models show similar fit probability (p-value), making it difficult to favor one of them the most. A combined fit to the SED with both BPL and PLE is shown in Fig. 2. The parameters of the different models, tested with χ 2 /d.o.f., are given in Table 4.

POTENTIAL COUNTERPARTS
Several point-like sources lie in the FoV of the extended gammaray source HESS J1841-055 and are likely to contribute to the VHE emission. In this section we discuss all these sources and their association with the observed emission. We consider some of the brightest emissions from these sources (see Fig. 4) discussed below to constrain the gamma-ray emission mechanisms in Section 5.

G26.6-0.1
The diffuse hard X-ray source G26.6-0.1 was detected by ASCA in a Galactic plane survey (Bamba et al. 2003) which is located in this region as shown in Fig. 4. The observed X-ray spectrum was found to be featureless and can be fitted with a powerlaw function with photon index 1.3. The diffuse X-ray flux was estimated to be 3.5 × 10 −12 erg cm −2 s −1 from a radius of 12 region in the energy range of 0.7-7.0 keV. We consider this diffuse emission to be associated with HESS J1841-055 and assumed a corresponding scaled X-ray flux from a region with radius 0.4 • similar to the extension of our source, for the multi-wavelength (MWL) modelling in Section 5. The distance to G26.6-0.1 is 1.3 kpc (Bamba et al. 2003). Following this, the distance of HESS J1841-055 is assumed to be 2 kpc.

PSR J1838-0537, PSR J1841-0524 and PSR J1838-0549
Several gamma-ray pulsars lie within the HESS J1841-055 region. Pletsch et al. (2012) discovered the gamma-ray pulsar PSR J1838-0537 in a blind search of Fermi-LAT data. It has been proposed as a potential candidate for the VHE source. It is a radio quiet pulsar. Also, no X-ray pulsation is observed from the location of the pulsar. If it is associated with a nebula, the subsequent observation from this region should have provided a detectable level of radio and X-ray fluxes from this region. The spin down power of the pulsar is  estimated to be E = 5.9 × 10 36 erg s −1 . The integral energy flux of HESS J1841-055 estimated by MAGIC over the range 0.1-10 TeV is lγ ∼ 9.13 × 10 −11 erg cm −2 s −1 . The luminosity for a distance of 2 kpc is, L γ = 4πd 2 lγ= 4.37 × 10 34 erg s −1 for isotropic emission. This implies a conversion efficiency η = L γ / E ∼ 0.7% which is consistent with other suggested pulsar/pulsar wind nebulae (PWN) associations (Hessels et al. 2008). Hence, the pulsar's energetic is likely to power a PWN producing part of the TeV emission. The spectral index derived ∼2.4 is relatively soft in comparison to other PWNe detected at GeV energies by Fermi-LAT. Hence, part of the low-energy emission could have a different origin. There are two other known pulsars PSR J1841-0524 and PSR J1838-0549 (Aharonian et al. 2008) which can contribute to the emission of HESS J1841-055. The estimated E/D 2 values are given as 4.4×10 33 erg s −1 kpc −2 and 4.7×10 33 erg s −1 kpc −2 respectively and they can contribute to the observed emission when considered together. However, if taken separately, each would require approximately 200% efficiency to explain the VHE emission (Aharonian et al. 2008). There is no significant radio emission observed from the location of these pulsars. The observed X-ray emission from these sources is lower than that considered for the multi-waveband modelling. Hence, assuming these constrains, we will not consider them separately for the MWL modelling.
Hence, out of the three pulsars located in the region, it is most likely that only PSR J1838-0537 is contributing to the detected gamma-ray emission.

FGES J1839.4-0554 & FGES J1841.4-0514
Recent Fermi-LAT catalog for Galactic extended sources (FGESs) shows that there are two distinct extended sources in this region with energies above 10 GeV (Ackermann et al. 2017  4. It is evident from the figure that the observed GeV emission is overlapping with the extension found at TeV energies. Therefore, they can be considered as potential counterparts for the TeV emission. Although it appears that there are two different sources, the spectral characteristics at energies above 10 GeV are similar indicating that they may have common origin (Ackermann et al. 2017). In our analysis of Fermi-LAT data in this paper, we consider the entire region which includes both these sources to estimate the SED and we use it for the multi-wavelength modelling.

G27.4+0.00 (Kes 73)
One of the sources studied in radio is a shell-type remnant Kes 73 (G027.4+00.0) which is present at the north-east of the extended emission region (see Fig. 4). The small diameter 5 radio shell is characterized by a steep spectral index (α ∼ −0.68, defined by S ∝ ν α ) between 0.5 to 5 GHz and flux density of 3.5 ± 0.5 Jy at 1.4 GHz (Caswell et al. 1982). Radio studies of the remnant also show an incomplete shell structure with no central engine of Kes 73 and with a radio upper limits of 0.45 mJy and 0.60 mJy at 6 cm and 20 cm radio wavelengths respectively (Kriss et al. 1985). This is considered to be unlikely counterparts due to the very small angular size of 5 and its location on the edge of the extended emission.

AX J1840.4-0537 & AX J1841.4-0536
A weak point-like source, 1RXS J184049.1-054336, is located within G26.6-0.1 and it contributes to less than 10% of the diffuse flux. Hence it is reasonable to exclude this weak X-ray flux from our analysis. The other X-ray point sources AX J1840.4-0537 and AX J1841.4-0536 are located outside the G26.6-0.1 region but well within the extended HESS J1841-055. However the fluxes for these two sources were estimated to be 1.4 × 10 −13 erg cm −2 s −1 and 2.1 × 10 −11 erg cm −2 s −1 , respectively, which are below the level of the scaled diffuse X-ray flux from this extended gamma-ray region. Moreover, due to the point-like morphology of these sources with no associated nebula around them, they can not be considered as potential counterparts of HESS J1841-055. However, a fraction of the total emission could be associated with these sources. Hence, although it is challenging to disentangle which sources are contributing to the observed GeV-TeV emission, we consider that the SNR G26.6-0.1, the pulsar PSR J1838-0537 and the extended FGES J1839.4-9554 and FGES J1841.4-0514 sources are the most promising counterparts, due to their energetics, extension and location within the region.

MODELLING OF THE SPECTRAL ENERGY DENSITY
As already discussed above, there are several sources present in this extended region. Some of them are already argued to be po-tential counterparts at lower energies from the aspects of energetics of the system, while others are excluded due to their very pointlike signatures along with the energetics which can not contribute significantly to the overall emission, considering the extent of the emission. In order to investigate if the multi-wavelength data can be explained self-consistently, we consider that the observed emission is associated with HESS J1841-055. Since the radio and X-ray fluxes from the entire region of the extended emission can not be more than that estimated from different observations, we consider those results as upper-limits in these frequencies after multiplying with a scaling factor attributed to the extended region and the emission regions from where corresponding radio and X-ray measurements are performed. For the GeV-TeV modelling, we will consider the Fermi-LAT and MAGIC data sets from this study and the H.E.S.S. data points from Aharonian et al. (2008). We use the numerical code developed by Saha & Bhattacharjee (2015) for the modelling.

Leptonic model
In Section 3.3, we found that the SED has a spectral curvature and can be better explained with either a BPL model or a PLE model. Since the observed gamma-ray spectra carry imprints of the intrinsic particle distribution, a single population of electrons that follows a BPL type of distribution of electrons is assumed to calculate the Inverse Compton (IC) and bremsstrahlung emission processes.
In general, the electron spectrum might be more complicated than assuming a single population of electrons. For example, for the Crab Nebula, two different population of electrons are considered, namely, radio electrons and wind electrons. Radio electrons are less energetic electrons which reside in the nebular volume throughout its age, and they are mostly responsible for the observed radio fluxes. On the other hand, wind electrons are freshly accelerated electrons and they account for the observed fluxes at X-ray and GeV-TeV energies. However, for simplicity, we consider a single population of electrons that is responsible for the observed emission at GeV-TeV energies.
We first consider that the observed gamma-ray radiation at GeV-TeV energies, is resulting from emission from relativistic electrons through IC and non-thermal bremsstrahlung processes. For IC process, we consider that the high energy photons are produced by the up-scattering of photons from the Cosmic Microwave Background (CMB) and from interstellar dust contribution (Mathis et al. 1983). For bremsstrahlung process, we consider ambient matter density of 100 cm −3 following the estimation discussed in Appendix A. Higher or lower values of the ambient matter densities simply scale the contribution of bremsstrahlung spectrum. Fig. 5 shows both IC and bremsstrahlung spectra for the BPL electron distribution for an ambient matter density of 100 cm −3 . The figure shows that the bremsstrahlung spectrum can explain the observed SED at GeV-TeV energies. On the other hand, the IC emission for the target photons of CMB and star lights can not explain the observed SED Bremsstrahlung emission spectrum, estimated for an ambient matter density of 100 cm −3 , is the dominant one. The parameters of the BPL electron distribution are shown in Table 5.
for the same population of electrons. The parameters of the BPL electron distribution are shown in Table 5. It is to be noted that the electron distribution can be adjusted to explain the observed emission by the IC spectrum. However, the bremsstrahlung spectrum will overestimate the observed flux for the same population of electron due to high ambient matter density. Hence bremsstrahlung becomes dominant emission process within leptonic scenario. In order to check the contribution of the synchrotron spectrum for the electron distribution, we calculate the synchrotron spectrum leaving magnetic field as a free parameter. We find that the synchrotron spectrum for a magnetic field of approximately 5 µG does not overestimate the radio and X-ray limits estimated for this study. The synchrotron component only contributes to radio and X-ray energies. The synchrotron spectrum together with IC and bremsstrahlung spectra is shown in Fig. 6.

Hadronic model
We also introduce a hadronic scenario as an additional component which contributes significantly at gamma rays (GeV-TeV). We calculate the gamma-ray spectrum resulting from the decay of neutral pions following Kelner et al. 2006. The gamma-ray spectrum for the relativistic protons for the BPL model as considered for the leptonic model and for an ambient gas density of n 0 100 cm −3 is shown in Fig. 7. The total energy can be calculated as W p = 5.52 × 10 48 × (100.0/n 0 ) erg. The figure displays that the gamma-ray spectrum resulting from the decay of neutral pions can explain the observed GeV-TeV data very well. The parameters of the model are presented in Table 5.

Gamma rays from extended unidentified sources
The analysis of about 34 hrs of good quality MAGIC data confirms that the gamma-ray emission is as extended as claimed by the H.E.S.S. Collaboration and the source is detected with high Table 5. Parameters for physical models for a single zone particle distribution of a BPL model. The parameters are obtained considering two different models: leptonic and hadronic. Parameters with errors are used as free parameters for the fit.

Parameters
Leptonic Hadronic significance for energy above 50 GeV. In addition to that, we investigate the source morphology as a function of energy. The observed results suggest that at low energies the overall region is detected like a diffuse source, with some few regions around the center of the source where the significance is higher than 5σ. This indicates the possibility that several point-like sources are contributing to the extended emission. At medium energies, between 500 GeV -1 TeV, the emission is concentrated along the center, in a North-South line.
In addition to that, the skymap above 1 TeV shows that the extension of the emission is reduced compared to that of low energies and the most significant flux is located at the southern part of the extended region, with only few hot spots over 5σ. The morphological analysis of MAGIC data also shows that the multiple source model is better over a single-source model. This establishes the fact that several sources are contributing to the extended emission. The morphological study of the source using about 10 years of Fermi-LAT data above 10 GeV also shows that the source is extended. However, the spectral shape is different from that of MAGIC. When comparing with the emission at higher (TeV) energies, HESS J1841-055 displays an extension compatible to that measured by MAGIC.
In the case of the Fermi-LAT detection, the spectrum of the source is best described by a powerlaw. To study the spectral behavior within the entire energy range, from GeV to TeV energies, we performed a joint fit on the spectral data points from MAGIC and Fermi-LAT. We found that the combined SED is best described either by a broken powerlaw model with a spectral break at ∼37 GeV or with a powerlaw with exponential cut-off at 1.8 TeV.

Emission mechanisms
Multi-wavelength modelling of the data indicates that the leptonic model can explain the data well. Due to the higher ambient matter density, the bremsstrahlung spectrum dominates over IC spectrum. The radio and X-ray fluxes put a constraint on the magnetic field in the emission volume when they are accounted with a synchrotron emission process. The magnetic field of 5 µG as mentioned in section 5.1 is very close to that of some other known old PWNe (Reynolds et al. 2012;Kargaltsev et al. 2013). Given the high ambient matter density and presence of molecular clouds, a hadronic emission model is also suitable to explain the observed data at GeV-TeV energies. We find that the hadronic model can explain the data  Table 5. very well for a BPL proton distribution and an ambient matter density of 100/cm 3 . Therefore, both leptonic and hadronic model can explain the data well with the parameters shown in Table 5.
In the whole discussion on multi-wavelength modelling of the data, our assumption was that the observed emission is entirely due to a single source. However, we have already seen that the region is populated by different sources which were established through observations at lower energy bands. Some of the sources are already excluded to be considered as potential gamma-ray emitters while energetic is considered. However, some of them could be potentially associated with the observed emission at GeV-TeV energies. Given the angular resolution of the gamma-ray telescopes at present generation, it is not possible to have an unambiguous association with the sources at other wavebands. One possible scenario for the extension of the emission is the interaction of run-away cosmic particles from the source and the gamma-ray visibility is enhanced due to interac- tion with molecular clouds which are covering the extended source very well as can be seen from Fig. 8 and discussed in appendix A. The presence of molecular clouds along the extension of the source also supports the relatively high ambient matter density required for both leptonic and hadronic model.

The nature of HESS J1841-055
The observations at X-ray energies did not show any bright synchrotron nebula around the pulsars present in this region. However, in this scenario the absence of bright synchrotron nebula can be easily explained. If the TeV source is powered by one or several pulsars present in this region, then pulsars are expected to be relic ones. For such PWNe, IC emission efficiency is more pronounced due to lower magnetic fields. In Section 5, we find that the IC contribution for this source is 10% less compared to the bremsstrahlung spectrum. Therefore, it is reasonable to consider that the synchrotron emission could be even more inefficient, which supports the absence of the bright synchrotron nebula around the pulsars. Therefore, if the bright TeV emission is assumed to be associated with a PWN, then the PWN requires to be a relic one where the remnant of the supernova explosion has already disappeared. It is also discussed in Section 4 that energetically the pulsar PSR J1838-0537 is able to account for the observed GeV-TeV energies. The gamma-ray flux at TeV energies, a factor two lower than the Crab nebula flux, is required to have, S 0 = (L 0 /10 37 erg s −1 )(d/1 kpc) −2 ≥ 10 −3 , where L 0 and d are the luminosity and distance of the source respectively (Aharonian 2004). The two known pulsars PSR J1841-0524 and PSR J1838-0549 cannot contribute to the observed GeV-TeV energies since S 0 is less than 10 −3 . However, S 0 is greater than 10 −3 for PSR J1838-0537 making it a potential counterpart of HESS J1841-055. Since PSR J1838-0537 is not a part of any radio or X-ray nebula, it is also possible to consider that it is an isolated pulsar which has already left the remnant. TeV emission is an effective product of the IC mechanism for such isolated pulsars, with the injection of relativistic electrons in the interstellar magnetic field which is about 3µG. In such a scenario, the bright X-ray and radio synchrotron nebula could be absent. However, the extension of such a source is not readily accepted. Nevertheless, the presence of the molecular clouds along the observed GeV-TeV emission can support its extension within this scenario but through bremsstrahlung processes.
The extension of the source is estimated to be 0.4 • when the source is fitted with a 2D Gaussian with a equal spatial width for both the directions. This extension translates into a radius of approximately 14 pc at a source distance of 2 kpc. The effective diffusion radius can be calculate as , R di f f 2 D(E) t, where D(E) is the diffusion coefficient and can be represented as D(E) = D 0 (E/10GeV) δ (Atoyan et al. 1995). The commonly used diffusion coefficient at 10 GeV is of about D 0 ∼ 10 28 cm 2 s −1 (Berezinsky et al. 1990) and assuming that δ = 0.5 (δ is one of the parameters of the diffusion coefficients; for energy independent diffusion coefficient δ=0), we calculate the diffusion time scale of t di f f = 17 kyr. On the other hand, the lifetime of the bremsstrahlung loss, which is independent of energy, is estimated as t brems 4 × 10 4 (n/1 cm −3 ) −1 kyr = 4 × 10 2 kyr for ambient matter density of 100 cm −3 . Therefore, the dominant emission through bremsstrahlung process for the estimated ambient matter density is a viable solution for the observed extension of the source. We then conclude that the observed emission can be potentially associated with a PWN.
The observed emission can also be considered to be associated with SNRs, since it is seen that there are two SNRs present in and around the source. The first one, G27.4+0.00, is located at the edge of the TeV emission and has a relatively small angular size, hence is unlikely that it can account for the observed gamma rays. G26.6-0.1 is considered to be a potential counterpart for the extended emission. However, there are no strong radio and X-ray nebulae associated with the extent of the emission which is the case for a typical SNR scenario. Hence, a strong association can be made provided that the observed emission is considered due to the particles that escaped the SNR shocks and are interacting with the molecular clouds. Following the diffusion timescale as discussed in the preceding paragraph and the age of the SNR G26.6-0.1 as a middle-age SNR (∼10 3 years; Bamba et al. 2003), it can be considered a possible candidate for at least part of the detected GeV-TeV emission.

SUMMARY AND CONCLUSIONS
We report a deep study of the unidentified gamma-ray source HESS J1841-055 at GeV-TeV energies using about 34 hours of MAGIC and 10 years of Fermi-LAT data. We summarize the results below.
• The results of the detailed analysis show that the observed gamma-ray emission from HESS J1841-055 is significantly extended. The estimated extension of the source using MAGIC data is similar to that reported by the H.E.S.S. Collaboration, found to be ∼ 0.4 • assuming a Gaussian distribution.
• There are several bright hot-spots in the extension of the source which appears to be multiple sources which contribute to the observed emission at GeV-TeV energies. The emission at TeV energies moves towards the south with increasing energy, revealing this region as one of the potential main contributors of the TeV extended emission.
• The extended emission is modelled better with a multi-source model compared to a single-source model.
• The spectral curvature of the SED in the energy range from GeV-TeV is significant and it can either be described by a broken powerlaw model with break at 37 GeV or a powerlaw with exponential cutoff at 1.8 TeV.
• The observed SED can be explained well with both a leptonic (bremsstrahlung) and a hadronic model for the density of ambient matter of 100 cm −3 assuming a BPL distribution of electrons and protons, respectively.
Within the present morphological and spectral studies of this extended source using GeV-TeV data and available MWL information on sources present within the region, we conclude that the extended gamma-ray emission seems to be associated with multiple sources in this region. The GeV-TeV emission is compatible with a PWN scenario, although a fraction of the gamma-ray emission can also be explained within a SNR scenario. However, disentangling these sources at TeV energies (either point sources or extended sources) from one another and quantifying their contribution to the observed morphology of the source demands much better angular resolution compared to the present generation of gamma-ray telescopes. Hence, it becomes naturally an interesting source of study for the next generation of IACT telescopes.

APPENDIX A: TARGET GAS DENSITY
Here we evaluate the target gas density required both for leptonic and hadronic models as discussed in Sec. 5.1 and Sec. 5.2. To evaluate the target gas density, we estimate the densities of each gas phase (neutral hydrogen HI and molecular hydrogen H 2 ) and then sum the estimated values to get the total contribution to the gas density. Under the assumption of the optically thin limit, the HI column density is given by (Dickey & Lockman 1990), N(HI) 1.823 × 10 18 ∫ T b (HI; v r )dv r cm −2 , where T b (HI; v r ) is the brightness temperature of the observed 21 cm line at v r . In order not to overestimate the gas density within the source, we need to integrate over some range of v r . We consider v r in the range of 110 -135 km s −1 corresponding to the distance of about 2 kpc. The average HI density is estimated to be N(HI) = 8.65 × 10 20 cm −2 for a radius of 0.4 • centered on the HESS J1841-055 using the database of the HI4π survey (HI4PI Collaboration et al. 2016). We assume that the HI gas is uniformly distributed within the source. The length of the line of sight along the direction of the HESS J1841-055 is = 2r 0 , where r 0 is the radius of the extended emission. The radius of the source is considered to be 0 • .4 which translates to approximately 14 pc for a distance of 2 kpc to the source. The density of the neutral hydrogen gas is n(HI) = N(HI)/ 10 cm −3 . We use observations (see Fig. 8) of the 13 CO(J = 1 − 0) line, which traces molecular clouds, from the Galactic Ring Survey (GRS; Jackson et al. 2006).The CO spectrum over the range of velocities from +110 to +135 km s −1 are integrated to obtain the velocity integrated CO intensity (W CO ). The W CO averaged over the region with a radius of 0 • .4 covering the extended emission is found to be approximately 63 K km s −1 . To estimate the mass of the molecular cloud, we use the standard CO-to-H2 conversion factor of X CO = N(H2)/WCO = 1.8 × 10 20 cm −2 K −1 km −1 s (Dame et al. 2001). We find N(H2) = 4.8 × 10 22 cm −2 . Therefore, the density of the molecular hydrogen gas, n(H 2 ) = N(H 2 )/ 130 cm −3 . The total gas density, hence, is n(HI) + n(H 2 ) 140 cm −3 . However, for simplicity we consider the gas density of 100 cm −3 for the physical modelling of the source.