Energy-dependent Analyses of the Gamma-Ray Emission from HESS J1857+026 with Fermi-LAT

We report the discovery of the energy-dependent morphology of the GeV gamma-ray emission from HESS J1857+026 with more than 13 yr of Fermi Large Area Telescope data. The GeV gamma-ray emission from this region is composed of two extended components. The hard component with an index of 1.74 ± 0.07 in the energy range of 0.5–500 GeV is spatially coincident with HESS J1857+026, and its 68% containment radius varies from ∼0.°44 below 40 GeV to ∼0.°30 above 140 GeV. The hard GeV gamma-ray spectrum and the energy-dependent morphology of HESS J1857+026 make it favor a pulsar wind nebula origin, which is associated with the energetic pulsar PSR J1856+0245. The soft component with an index of 2.70 ± 0.16 and another extended gamma-ray source detected in this region, 4FGL J1857.9+0313e, with an index of 2.55 ± 0.07, are spatially coincident with two molecular clumps in the northeast and southwest of HESS J1857+026, which favors the hadronic process, and the protons could be accelerated by the hypothetical supernova remnant associated with PSR J1856+0245.


INTRODUCTION
More than two hundred of very-high-energy (VHE; > 100 GeV) gamma-ray sources have been detected with the operations of ground-based Cherenkov telescopes, such as H.E.S.S. (Aharonian et al. 2004), MAGIC (Aleksić et al. 2016), HAWC (Abeysekara et al. 2013) and LHAASO (Cao et al. 2023).The gamma-ray emission could be produced by the hadronic interactions in which π 0 decays into two gamma-ray photons subsequently, via the inverse Compton (IC) scattering process or through non-thermal bremsstrahlung radiation from high energy electrons.Most of the VHE gamma-ray sources have been identified to be pulsar wind nebulae (PWNe), supernova remnants (SNRs), X-ray binaries, TeV halos and so on.However, some TeV gamma-ray sources are still unidentified, and the multi-wavelength studies of these sources are crucial for revealing their nature and probing the origin of cosmic rays.
HESS J1857+026 was first discovered as a VHE gamma-ray source with an extension of (0.11±0.08) • ×(0.08±0.03)• (Aharonian et al. 2008).In H. E. S. S. Collaboration et al. (2018), the morphology of HESS J1857+026 was described by a two dimensional (2D) Gaussian component with an approximate size of 0.26 • ± 0.06 • .Hessels et al. (2008) discovered an energetic pulsar, PSR J1856+0245, in the direction of HESS J1857+026, which makes HESS J1857+026 to be a potential PWN candidate.The period and spin-down luminosity of PSR J1856+0245 are P = 81 ms and Ė = 4.6 × 10 36 erg s −1 , respectively, with the characteristic age of τ c = 21 kyr.The distance of PSR J1856+0245 was first estimated to be ∼ 9 kpc, derived by the dispersion measure (DM) with electron density model of Cordes & Lazio (2002).And, an updated distance of ∼ 6.3 kpc is given by ATNF Pulsar Catalog 1 according to the electron density model of Yao et al. (2017).MAGIC carried out follow-up observations and presented energy-dependent morphology of this region (MAGIC Collaboration et al. 2014).In the energy range of 300 GeV-1 TeV, the morphology observed by MAGIC is compatible with H.E.S.S. observation.However, two separate gamma-ray sources, named MAGIC J1857.2+0263(hereafter MAG1) and MAGIC J1857.6+0297(hereafter MAG2), were detected with the data above 1 TeV.MAG1 is an extended source with the intrinsic extension of (0.17 ± 0.03 stat ± 0.02 sys ) • × (0.06 ± 0.03 stat ± 0.02 sys ) • , while MAG2 is compatible with 1 https://www.atnf.csiro.au/research/pulsar/psrcat/arXiv:2312.11189v1[astro-ph.HE] 18 Dec 2023 a point source.MAGIC Collaboration et al. (2014) interpreted MAG1 as a PWN powered by PSR J1856+0245, and MAG2 may be associated with a molecular cloud complex containing an Hii region located at ∼ 3.7 kpc and a possible gas cavity.Another pulsar, PSR J1857+0300, was discovered in the direction of MAG2 (Lyne et al. 2017).The characteristic age and spin-down luminosity of PSR J1857+0300 are τ c ∼ 4.6 × 10 6 yr and Ė ∼ 2.3 × 10 32 erg s −1 , with the distance of ∼ 6.7 kpc.Meanwhile, an elliptical superbubble, was detected with neutral gas observation (Petriella et al. 2021), which is also spatially coincident with HESS J1857+026.The kinematical distance of the superbubble is about 5.5 kpc and is close to the DM distance of PSR J1856+0245 (∼ 6.3 kpc).Petriella et al. (2021) concluded that the TeV emission of HESS J1857+026 originates from the superbubble, and PSR J1856+0245 is located inside the superbubble.In addition, they found five molecular components in the velocity interval between 78 and 90 km s −1 with 13 CO(J=1-0) observations, which are probably associated with the superbubble.And they favor a single gamma-ray source scenario instead of the superposition of two gamma-ray sources.Petriella et al. (2021) also performed radio observations at 1.5 GHz and 6.0 GHz with VLA.Nevertheless, no significant radio emission was detected in this region.
In the energy range of 1-25 TeV, the Water Cherenkov Detector Array (WCDA) of Large High Altitude Air Shower Observatory (LHAASO) detected an extended gamma-ray source, 1LHAASO J1857+0245, which is spatially consistent with HESS J1857+026 (Cao et al. 2023).The spatial template of 1LHAASO J1857+0245 is described by a 2D Gaussian with σ = 0.24 • ± 0.04 • , and its gamma-ray spectrum in 1-25 TeV is modeled by a power law with an index of 2.93 ± 0.07.
The GeV gamma-ray emission from HESS J1857+026 was first detected using a point source hypothesis by Fermi-LAT (Atwood et al. 2009), while no gamma-ray pulsation from PSR J1856+0245 was observed (Rousseau et al. 2012).Since no obvious X-ray emission was detected, only an upper limit of F 1−10 keV = 2 × 10 −12 erg cm −2 s −1 was obtained (Rousseau et al. 2012).Considering the uncertain origin and the complexity of this region, the detailed analysis with more GeV observational data will be helpful to investigate the origin of the gamma-ray emission.
Taking advantage of more than 13 years of Fermi-LAT data, we performed an energy-dependent analyses of the region around HESS J1857+026 and discussed the nature of this source.The paper is organized as follows.In Section 2, we describe the data analysis routines and present our results.In Section 3, we discuss the radiation mechanisms and the nature of HESS J1857+026.And the conclusion for this work is presented in Section 4.

FERMI-LAT DATA AND RESULTS
The Fermi-LAT Pass 8 data we analyzed are collected from August 4, 2008 to March 4, 2022 with energies from 500 MeV to 500 GeV.The region of interest (ROI) is a 10 • × 10 • square centered at the TeV gamma-ray centroid of HESS J1857+026 (R.A. = 284.296• , Dec. = 2.667 • ; H. E. S. S. Collaboration et al. 2018).To reduce the contamination from the Earth Limb, events with zenith angle greater than 90 • are eliminated.The standard analysis software of Fermitools 2 is used with the instrumental response function (IRF) of P8R3 SOURCE V3.The models we used to describe the Galactic and isotropic diffuse emissions include gll iem v07.fits, iso P8R3 SOURCE V3 v1.txt and iso P8R3 SOURCE V3 PSF3 v13 .In addition, all sources in the incremental version of the fourth Fermi-LAT source catalog (4FGL-DR3; Abdollahi et al. 2020Abdollahi et al. , 2022) ) within a radius of 15 • centered at HESS J1857+026 are included in the model.And the binned maximum likelihood analysis method with gtlike is applied.During the fitting procedure, the spectral parameters of all sources located in the ROI are left free, together with the normalizations of the two diffuse backgrounds.

Average Spatial Extension
In the 4FGL-DR3 catalog, the gamma-ray emission of HESS J1857+026 is described by an uniform disk (named as 4FGL J1857.7+0246e)centered at R.A. = 284.449• , Dec. = 2.774 • with the 68% containment radius of r 68 = 0.50 • , which is given by the analysis of Fermi-LAT extended Galactic sources (FGES; Ackermann et al. 2017).Meanwhile, there is a point source named as 4FGL J1857.9+0313clocated in the north of the disk, which has no identified counterpart (Abdollahi et al. 2020).And in the south of the disk, two point sources named as 4FGL J1857.6+0212 and 4FGL J1858.3+0209 are identified to be associated with SNR G35.6-0.4 and HESS J1858+020 (Cui et al. 2021;Zhang et al. 2022).To obtain the spatial template of HESS J1857+026, we performed the spatial extension analysis in the energy ranges of 1-3 GeV (low energy) and 10-500 GeV(high energy), respectively.And in the low energy band, only "PSF3" type (evclass=128 & evtype=32) data with better angular resolution are selected to reduce the contamination from nearby sources, while the data with "SOURCE" type (evclass=128 & evtype=3) are used in the high energy range.After subtracting the background sources included in the model (except 4FGL J1857.9+0313c in left panel and 4FGL J1857.7+0246e in the middle and right panel), we created three 2.5 • × 2.5 • TS maps centered at HESS J1857+026 with the different energy ranges, as shown in Figure 1.
As shown in the left panel of Figure 1, the position of 4FGL J1857.9+0313cgiven by 4FGL-DR3 is not coincident with the the gamma-ray peak of this region.In addition, there are discrepancies between the morphologies of the GeV emission around HESS J1857+026 and the spatial templates given by FGES both in the low and high energy bands, as shown in the middle and right panels of Figure 1.Therefore, we refined the spatial templates of HESS J1857+026 of the low (hereafter referred to be "Src A") and high ("Src T") energies, and the morphology of 4FGL J1857.9+0313c was also re-analyzed.The uniform disk, 2D Gaussian, as well as the H.E.S.S. image model are tested.Meanwhile, a point source model with the best-fit coordinate calculated by gtfindsrc is also applied for Src A and 4FGL J1857.9+0313c to test the spatial extension of it.And the MAGIC image above 1 TeV is also adopted as the spatial template of Src T to explore whether or not the gamma-ray emission above 10 GeV is a superposition of two sources, which is similar to MAG1 and MAG2 (MAGIC Collaboration et al. 2014).The centroids and extensions of the 2D Gaussian and uniform disk are fitted by Fermipy4 , a PYTHON package that automates analyses with the Fermi Science Tools (Wood et al. 2017).And the results of the spatial analysis with the different energy ranges are listed in Table 1.From Table 1, we can see that both the uniform disk and 2D Gaussian we analyzed can describe the GeV emission of Src A and 4FGL J1857.9+0313c(hereafter renamed as 4FGL J1857.9+0313e).And for Src T, a 2D Gaussian model is the best-fit template instead of the MAGIC image, indicating that the GeV emission from Src T is in favor of one-source scenario instead of two separate gamma-ray sources.In the following analyses, the uniform disk is adopted as the spatial template of Src A and 4FGL J1857.9+0313e, while Src T is described by the 2D Gaussian model.As shown in the middle panel of Figure 1, the centroid of the gamma-ray emission from Src A is far from the position of PSR J1856+0245/PSR 1857+0300 and is on the edge of the TeV emission of HESS J1857+026, which suggests that there is no spatial coincidence between Src A and the pulsars or TeV emission.For Src T, the extension of the best-fit 2D-Gaussian template is much smaller than the result given by FGES, which could be attributable to the improvement of the Galactic diffuse background or the newly detected gamma-ray sources.Considering the comparable extension and the spatial coincidence between Src T and H.E.S.S. observation (H.E. S. S. Collaboration et al. 2018), we suggest that the GeV emission of Src T has the same origin of the TeV emission.In addition, we found that the spectral index of Src A in the energy range of 1-3 GeV is Γ = 2.50 ± 0.34, while the spectral index of Src T for the data above 10 GeV is Γ = 1.84 ± 0.10.Considering the distinct GeV spectra and morphologies, the origins of Src A and Src T are probably different.And both Src A and Src T are considered in the following analysis.⋆ r68 is the 68% containment radius, where r68 = 1.51σ for 2D Gaussian model and r68 = 0.82σ for uniform disk model (Lande et al. 2012).† Calculated with respect to the spatial model used in 4FGL-DR3.

Energy-Dependent Extension Analysis of Src T
To further explore the energy-dependent behavior of Src T, we performed the extension analyses in the energy ranges of 10-40 GeV, 40-140 GeV and 140-500 GeV, respectively.For Src T, the 2D-Gaussian template is adopted, while the centroid and extension in each energy band are refitted with Fermipy.The results with the different energy ranges are listed in Table 2, and the corresponding TS maps are presented in Figure 2.
The energy-dependent analysis show that the extension of Src T in the energy range of 10-40 GeV is larger than that of higher energy bands, with r 68 varying from 0.44 • below 40 GeV to ∼ 0.30 • above 140 GeV.Similar phenomenons are also observed in the typical PWN HESS J1825-137 at TeV and GeV energies (Aharonian et al. 2006; H. E. S. S. Collaboration et al. 2019;Principe et al. 2020).In addition, it should be noted that the gamma-ray emission regions in the energy range of 40-140 GeV and 140 GeV-500 GeV seem to be different, as shown in the middle and right panels of Figure 2.However, such phenomena could be attributed to the limited statistics of gamma-ray photons, and more observational data will be helpful to explore the energy-dependent behavior of Src T.

Spectral Analysis
The spectral analysis was performed in the energy range of 500 MeV -500 GeV.And similar to the spatial analysis, the events of "PSF3" type with a better angular resolution were selected for the data below 3 GeV.During the analysis process, the summed likelihood analysis method was adopted, and both Src A and Src T were included in the model with their spatial templates given in Sect.2.1.1.The spectra of Src A, Src T and 4FGL J1857.9+0313e are adopted to be the power law models.And we also tested the GeV spectral curvature for each source by adopting the log-parabola spectrum, while no significant improvement was obtained compared with the power law models.The global fitting gives a hard spectrum for Src T with the photon index of Γ = 1.74 ± 0.07, which is comparable to the result in Rousseau et al. (2012).And the integral photon flux is estimated to be (5.04 ± 1.08) × 10 −9 ph cm −2 s −1 .While the spectrum of Src A is very soft, and the photon index is fitted to be Γ = 2.73 ± 0.10, with the integral photon flux of (1.98 ± 0.19) × 10 −8 ph cm −2 s −1 .And the best-fit photon index of 4FGL J1857.9+0313e is Γ = 2.55 ± 0.07 with the integral photon flux of (1.41 ± 0.13) × 10 −8 ph cm −2 s −1 .
To study the spectral energy distributions (SEDs) of these three sources, the data are divided into nine logarithmically space energy bins.And the summed likelihood analysis is repeated in each energy bin, with only the normalizations of sources located in ROI and the diffuse backgrounds in the model left free, while the spectral parameters are fixed to be the global-fitting values.For the energy bin with TS value lower than 4.0, an upper limit with 95% confidence level is calculated.And the SEDs are shown in Figure 3.The GeV spectrum of Src T could connect with the TeV SED of HESS J1857+026 smoothly, which suggests the same physical origin.2014).And the cyan butterfly indicates the global power law spectrum of 1LHAASO J1857+0245 detected by LHAASO-WCDA (Cao et al. 2023).Right: The GeV SED of 4FGL J1857.9+0313e.The gray solid and dashed lines show the best-fit power law spectrum and its 1σ statistic error.And the hadronic models with Ep,cut = 3 PeV and Ep,cut = 1 TeV for 4FGL J1857.9+0313e are shown as the red solid and dashed lines, respectively.The spatial and spectral data analyses above reveal that the diffuse GeV gamma-ray emission around HESS J1857+026 could be distinguished into two separate extended components: one of them, namely Src A, has no spatial and spectral coincidence with HESS J1857+026.And another one, namely Src T, shows both spatial and spectral consistency with the TeV measurement of HESS J1857+026, which supports Src T as the GeV counterpart of HESS J1857+026.
Although MAGIC observation revealed two sources in the region of HESS J1857+026, MAG1 and MAG2, with the events above 1 TeV (MAGIC Collaboration et al. 2014), there is no evidence show that Src T is composed by two gamma-ray sources limited by the event statistics and the PSF of Fermi-LAT.Therefore, we suggest that the GeV emission of HESS J1857+026 originates from a single gamma-ray source.The right panel of Figure 1 shows that the centroid of GeV gamma-ray emission is consistent with PSR J1856+0245 associated with MAG1, not PSR J1857+0300 in the direction of MAG2, which supports the same origin for HESS J1857+026 and MAG1.And for the possible origin of the gamma-ray emission from HESS J1857+026, MAGIC Collaboration et al. ( 2014) suggested a PWN powered by PSR J1856+0245 with Ė = 4.6 × 10 36 erg s −1 , and the spin-down luminosity of PSR J1857+0300 with Ė ∼ 2.5 × 10 32 erg s −1 is too low to power a gamma-ray PWN (Abdo et al. 2013;Acero et al. 2013).While Petriella et al. (2021) revealed the existence of a superbubble with the analysis of atomic gas in this region, and suggested the superbubble origin for the TeV emission from HESS J1857+026.However, the hard GeV gamma-ray spectrum of HESS J1857+026 makes it to be different from the typical superbubble with an index of ∼ 2.2, e.g.Cygnus Cocoon (Ackermann et al. 2011;Aharonian et al. 2019), but is similar to the typical PWNe, e.g.MSH 15-52 (Abdo et al. 2010a), HESS J1825-137 (Grondin et al. 2011).In addition, the GeV gamma-ray emission of HESS J1857+026 also shows the energy-dependent morphology, with the emission radius varying from 0.44 • below 40 GeV to ∼ 0.30 • above 140 GeV.And the centroid of the GeV emission moves towards PSR J1856+0245 with increasing energies.Such characteristics are also detected in the PWN HESS J1825-137 and HESS J1303-631 (Aharonian et al. 2006;H. E. S. S. Collaboration et al. 2019H. E. S. S. Collaboration et al. , 2012;;Principe et al. 2020).All these evidences above support that the gamma-ray emission of HESS J1857+026 could originate from a PWN associated with PSR J1856+0245.
For PWNe, the emission from radio to X-rays are normally produced by the synchrotron emission, whereas the gamma-ray emission are explained by the IC scattering process (leptonic process).Here, a simple one-zone leptonic model is applied for HESS J1857+026.The IC background photon fields include the cosmic microwave background (CMB) and infrared (IR) photon from dust with the temperature of T ∼ 30 K and density of 1 eV cm −3 et al. 2012).And the distance of HESS J1857+026 is adopted to be 6.3 kpc, derived from the dispersion measurement of PSR J1856+0245 (Yao et al. 2017).Considering the absence of radio detection for the PWN, we used the sensitivity of VLA in Petriella et al. (2021) to calculate the upper limits of radio flux by assuming that the spatial extension of the radio PWN is same to the that of the gamma-ray emission with a radius of 0.4 • .And the radio upper limits are estimated to be 9.2×10 −14 erg cm −2 s −1 at 1.5 GHz and 1.5×10 −13 erg cm −2 s −1 at 5 GHz, respectively.The electron spectrum is assumed to be a broken power law with an exponential cutoff in the form of dNe dE ∝ (E/E br ) −α 1 1+(E/E br ) α 2 −α 1 exp − E Ee,cut (Grondin et al. 2011;Xin et al. 2018), where E br , E e,cut , α 1 and α 2 are the break energy, cutoff energy, and indices of electron spectrum, respectively.The model fitting is performed using the naima package (Zabalza 2015).As shown in Figure 4, the gamma-ray spectrum of HESS J1857+026 can be reproduced with α 1 ∼ 2.2, α 2 ∼ 3.2, a break energy with E br ∼ 3.5 TeV, and a cutoff energy with E e,cut ∼ 70 TeV, respectively.The total energy of electrons above 1 GeV, W e , is calculated to be ∼ 1.1 × 10 49 (d/6.3kpc) 2 erg.The cooling timescale at break energy is estimated to be t cool ≈ 56(E/3.5 TeV) −1 [(U ph + U B )/1.66 eV cm −3 ] −1 kyr, here U ph = 1.26 eV/cm 3 and U B = B 2 /8π.And such value is about two times larger than the characteristic age of PSR J1856+0245, which suggests that the break structure could be an intrinsic characteristic of injected electronic spectrum, not be produced by the cooling effect (Gaensler & Slane 2006).With the electron spectrum, the radio and X-ray upper limits constrain the magnetic field strength to be lower than ∼4µG, which is consistent with the typical values for gamma-ray PWNe (e.g.Grondin et al. 2011Grondin et al. , 2013)).It should be noted that the X-ray upper limit was calculated with a radius of 0.1 • by Rousseau et al. (2012).And with the same radius of 0.4 • for the calculation of radio upper limits, the X-ray flux would be much larger, which would have no affect to the SED fitting.
Along with the evolution of the PWN into the interstellar medium (ISM), the energetic electrons could escape and transport dominated by diffusion to form a detectable halo around the pulsar, which is defined as a pulsar halo.Such halos were first detected around Geminga and PSR B0656+14 with the TeV gamma-ray emission Abeysekara et al. (2017).And Di Mauro et al. (2019) also claimed to detect the corresponding GeV gamma-ray emission around Geminga.Based on the definition of an electron halo in Giacinti et al. (2020), namely that of the overdensity of relativistic electrons around pulsar compared with ISM, we estimate the electronic energy density of HESS J1857+026, ε e , and compared it with the typical value of ISM, ε ISM = 0.1 eV cm −3 .The gamma-ray emission region of HESS J1857+026 is assumed to be a sphere with a radius of ∼ 0.40 • , which corresponds to a physical radius of ∼ 44 pc for the distance of 6.3 kpc.Using the total energy of electrons we derived, the electronic energy density of ε e ∼ 0.60 eV cm −3 is obtained, which is much larger than that of ISM.Therefore, we suggest that the relativistic electrons are still contained in a region energetically and dynamically dominated by the pulsar.2021) studies the molecular clouds in the range of 79 -90 km s −1 and derived a kinetic distance of ∼5.5 kpc, which is compatible, within the errors, with the DM distance of PSR J1856+0245 (∼6.3 kpc).By adopting the value of the conversion factor of X CO = 2 × 10 20 cm −2 (K km s −1 ) −1 (Bolatto et al. 2013), we estimate the total mass contents of clumpA and clumpB with the different distances, respectively.And the total masses of clumpA and clumpB are estimated in the region of 0.35 • and 0.28 • sky integration radii, corresponding to the 68% containment radii of the extended gamma-ray emission of SrcA and 4FGL J1857.9+0313e.respectively.For the velocity range of 50 -65 km s −1 with distance of 3.7 kpc, the total mass contents of clumpA and clumpB are calculated to be ∼ 2.1 × 10 5 d 2 3.7 M ⊙ and ∼ 1.7 × 10 5 d 2 3.7 M ⊙ , corresponding to the average gas number densities of n gas,A = 175 cm −3 and n gas,B = 280 cm −3 by assuming a spherical geometry of the gas distribution.For the velocity range of 79 -90 km s −1 with the compatible distance of PSR J1856+0245 with 6.3 kpc, the total mass content of clumpA and clumpB are calculated to be ∼ 1.5 × 10 5 d 2 6.3 M ⊙ and ∼ 1.6 × 10 5 d 2 6.3 M ⊙ .And the corresponding average gas number densities are about n gas,A = 27 cm −3 and n gas,B = 50 cm −3 , respectively.The spatial coincidence between the extended gamma-ray emission and the molecular gas suggest the hadronic origin for SrcA and 4FGL J1857.9+0313e.
In the hadronic scenario, the proton spectrum is assumed to be a single power law with an exponential cutoff in the form of And the cutoff energy of protons cannot be well constrained and was first adopted to be the energy of the cosmic ray knee with E p,cut = 3 PeV.The hadronic model for Src A is shown as the gray solid line in Figure 4, and the corresponding proton spectrum should be much soft with γ ∼ 2.8.The total energy of protons above 1 GeV is estimated to be W p ∼ 1.6 × 10 50 (n gas,A /27 cm −3 ) −2 (d/6.3 kpc) 2 erg or W p ∼ 1.6 × 10 48 (n gas,A /175 cm −3 ) −2 (d/3.7 kpc) 2 erg.For 4FGL J1857.9+0313e, the spectral index of protons with γ ∼ 2.7 and total energy of W p ∼ 5.8×10 49 (n gas,B /50 cm −3 ) −2 (d/6.3 kpc) 2 erg or W p ∼ 6.4×10 47 (n gas,B /280 cm −3 ) −2 (d/3.7 kpc) 2 are needed to explain the GeV gamma-ray emission, which is shown as the red solid line in the right panel of Figure 3.In addition, considering the fact that the cut-off energy of protons in the middle-aged SNRs is usually much lower than PeV (Gelfand et al. 2013;Supán et al. 2023), we decreased the cutoff energies of protons to estimate the different hadronic models.And the allowed minimum value of cutoff energy is about 1 TeV for SrcA and 4FGL J1857.9+0313e, with the total energies of protons of ∼ 1.5 × 10 50 (n gas,A /27 cm −3 ) −2 (d/6.3 kpc) 2 erg and ∼ 5.7 × 10 49 (n gas,B /50 cm −3 ) −2 (d/6.3 kpc) 2 erg, respectively.
For the molecular cloud in 50 -65 km s −1 with the distance of 3.7 kpc, there is no candidates as the origin of high energy protons.While for the molecular clouds in 79 -90 km s −1 , the compatible distance with PSR J1856+0245 suggests that the hypothetical SNR associated with PSR J1856+0245 could provide enough power assuming that ∼10% of the supernova kinetic energy of ∼10 51 erg is transferred to the energy of particles.Moreover, the soft GeV gammaray spectra of SrcA and 4FGL J1857.9+0313e are also similar to that of the old SNRs interacting with molecular clouds (e.g.IC 443 & W44 Abdo et al. 2010b;Ackermann et al. 2013).The absence of the detection of associated SNR in this region suggests that it may have already dissipated into the ambient gas.And further observations, especially in the radio and X-ray bands, will be crucial to reveal the physical origin of the gamma-ray emission in this region.

CONCLUSION
Using more than 13 years of Fermi-LAT observations, we studied the GeV gamma-ray emission in the direction of HESS J1857+026, and found that the GeV emission around HESS J1857+026 is composed of two extended gamma-ray sources: Src A and Src T. Src T is spatially coincident with HESS J1857+026 and its hard GeV gamma-ray spectrum could connect with the TeV SED of HESS J1857+026 smoothly, indicating that SrcT could be the GeV counterpart of HESS J1857+026.In addition, we performed the energy-dependent analyses of the GeV gamma-ray emission from HESS J1857+026, and its extension decreases towards higher energies.The energy-dependent morphology and the hard GeV gamma-ray spectrum of HESS J1857+026 make it favor the PWN origin.And one-zone leptonic model with a broken power law electronic spectrum can well describe the multi-wavelength data of HESS J1857+026.SrcA and another extended gamma-ray source, 4FGL J1857.9+0313e with soft GeV gamma-ray spectra, have no identified counterparts in other wavelength.However, two molecular clumps in the northeast and southwest of HESS J1857+026 are spatially coincidence with the GeV gamma-ray emission from SrcA and 4FGL J1857.9+0313e, which suggest the hadronic process for their gamma-ray emission.With the single power law model for protons, the GeV spectra of SrcA and 4FGL J1857.9+0313e could also be explained with the soft proton spectra.And these high energy protons could be produced by the hypothetical SNR associated with PSR J1856+0245, which may have already dissipated into the ambient gas.
HESS J1857+026 is one of the peculiar gamma-ray sources, which shows the energy-dependent morphology.More detailed observations by LHAASO (Cao et al. 2019) andCTA (Cherenkov Telescope Array Consortium et al. 2019) would be helpful to explore the particle transport mechanisms, and the future radio and X-ray observations are crucial to investigate the origin of the gamma-ray emission in this region.

Figure 1 .
Figure 1.2.5 • × 2.5 • TS maps centered on HESS J1857+026.Each pixel is 0.02 • and the TS maps are smoothed with a Gaussian kernel of σ = 0.04 • .The green contours indicate the 12 CO(J=1-0) emission integrated between 78 and 90 km s −1 by FUGIN (Umemoto et al. 2017).The magenta contours are the H.E.S.S. observation of HESS J1857+026 (H.E. S. S. Collaboration et al. 2018).The cyan dashed circle is the r68 of uniform disk of HESS J1857+026 given by FGES(Ackermann et al. 2017).The red solid circle shows the r68 of 2D Gaussian for 1LHAASO J1857+0245 detected by LHAASO-WCDA(Cao et al. 2023).And the blue and white crosses are the position of PSR J1856+0245 and PSR 1857+0300, respectively.Left: TS map for the data below 3 GeV with the diffuse backgrounds and 4FGL-DR3 sources subtracted except for 4FGL J1857.9+0313e.The black plus shows the position of 4FGL J1857.9+0313c as a point source given by 4FGL-DR3, and the white circle sign the r68 of best-fit uniform disk we analyzed here.Middle: TS map for the data below 3 GeV with the diffuse backgrounds and 4FGL-DR3 sources subtracted except for Src A. The r68 of the uniform disk for SrcA is described by the white circle.Right: TS map with data above 10 GeV.The white circle marks the r68 of the best-fit 2D Gaussian model for SrcT.

Figure 3 .
Figure 3. Left: The GeV SEDs of SrcA and SrcT are marked by the brown and black circles, respectively.The red and purple circles indicate the H.E.S.S. observations by Aharonian et al. (2008) and H. E. S. S. Collaboration et al. (2018), respectively.And the blue circles show the MAGIC observation in MAGIC Collaboration et al. (2014).And the cyan butterfly indicates the global power law spectrum of 1LHAASO J1857+0245 detected by LHAASO-WCDA(Cao et al. 2023).Right: The GeV SED of 4FGL J1857.9+0313e.The gray solid and dashed lines show the best-fit power law spectrum and its 1σ statistic error.And the hadronic models with Ep,cut = 3 PeV and Ep,cut = 1 TeV for 4FGL J1857.9+0313e are shown as the red solid and dashed lines, respectively.

Figure 4 .
Figure 4.The broadband SED of HESS J1857+026 with the leptonic model and SrcA with hadronic model.The radio and X-ray flux upper limits are shown as the green and pink dots (Rousseau et al. 2012; Petriella et al. 2021).The blue dashed line shows the synchrotron component for HESS J1857+026.The red and green dotted lines represent IC scattering of different seed photons.And the black solid curve is the sum of different leptonic radiation components for HESS J1857+026.The gray solid and dashed lines indicate the hadronic models with Ep,cut = 3 PeV and Ep,cut = 1 TeV for SrcA.

Figure 5 .
Figure 5. 12 CO(J=1-0) intensity maps (in the unit of K km s −1 ) in the velocity ranges of 50-65 km s −1 (left) and 78-90 km s −1 (right).The green, white and magenta circles show the GeV gamma-ray extensions of SrcA, HESS J1857+026, and 4FGL J1857.9+0313e,respectively.The position of PSR J1856+0245 is marked by the yellow cross.And the cyan contours indicate the MAGIC gamma-ray flux map above 1 TeV (MAGIC Collaboration et al. 2014).

Table 1 .
Spatial Properties for the GeV Emission in the Direction of HESS J1857+026