© The Authors 2023

1. Introduction

Cosmic rays (CRs) with energies up to the knee (∼1 PeV) are believed to be produced in hadronic PeVatrons, cosmic accelerators located in our Galaxy (for a review see e.g., Gabici et al. 2019). Despite substantial observational efforts in the last decade, the origin of the highest-energy galactic CRs remains unknown, mainly because of the difficulty in reconstructing the direction of their origin, as they are subject to deflection in the Galactic magnetic field. When accelerated protons interact with ambient matter, they emit γ rays through π⁰ decay. Similarly, electrons and positrons produce γ rays via inverse Compton (IC) scattering on low-energy photon fields via bremsstrahlung when colliding with atomic nuclei of ambient matter, or via synchrotron radiation when interacting with magnetic fields. Studying very high-energy (VHE, 0.1 < E < 100 TeV) and ultrahigh-energy (UHE, E > 0.1 PeV) cosmic γ-rays and disentangling the different origins of the radiation are therefore essential in order to search for cosmic PeVatrons (CTA Consortium 2019, and in prep.).

Diffusive shock acceleration (DSA; Bell 1978) taking place in supernova remnants (SNRs) and pulsar wind nebulae (PWNe) has been proposed as a possible mechanism to accelerate CRs (Bell 2013). In several SNRs, a characteristic spectral feature known as a “pion-decay bump” has been detected, providing evidence that proton acceleration takes place in these sources (Ackermann et al. 2013; Jogler & Funk 2016; Hess & Abdalla 2018a; Ambrogi et al. 2019; Abdollahi et al. 2022). However, none of the γ-ray spectra of the sources firmly identified as SNRs extend beyond 100 TeV (Aharonian et al. 2019; Zeng et al. 2019), which suggests that these sources are probably not capable of proton acceleration up to PeV energies.

The search for PeVatrons continues, and in the last few years several new candidates showing γ-ray emission above 100 TeV have been discovered by the Tibet Air Shower (AS) collaboration (Amenomori et al. 2019; Tibet ASgamma Collaboration 2021) and the High Altitude Water Cherenkov (HAWC)

observatory (Abeysekara et al. 2020). Recently, the Large High Altitude Air Shower Observatory (LHAASO) collaboration, exploiting the unprecedented sensitivity of the LHAASO-KM2A instrument in the UHE range, reported the discovery of 12 UHE γ-ray sources reaching energies up to 1.4 PeV (Cao et al. 2021b). Among these, there is only one unidentified source with an as-of-yet undetected TeV counterpart, namely LHAASO J2108+5157.

LHAASO J2108+5157 is the first γ-ray source directly discovered in the UHE band, and was detected with a post-trial significance of 6.4σ above 100 TeV (Cao et al. 2021a). The position of the source is RA = 317.22° ±0.07° , Dec = 51.95° ±0.05°. The source is reported to be point-like with a 95% confidence level upper limit on its extension of 0.26° with a two-dimensional symmetrical Gaussian shape assumption. The spectrum of LHAASO J2108+5157 above 25 TeV can be described by a single power law (PL) with a photon index of 2.83 ± 0.18. There is no VHE or X-ray counterpart to the source, but Cao et al. (2021a) identified a close high-energy (HE) point source, 4FGL J2108.0+5155 (Abdollahi et al. 2020), at an angular distance of 0.13°. A dedicated analysis suggested that the HE source might be spatially extended (4FGL J2108.0+5155e) with extension of 0.48°. Its spectrum can be described with a single PL with a photon index of 2.3 between 1 GeV and 1 TeV. However, the physical connection between the spectral energy distributions (SEDs) of the HE and UHE sources is not particularly clear because of the very different spectral indices. Cao et al. (2021a) found the source to be coincident with the position of a molecular cloud [MML2017]4607 (Miville-Deschênes et al. 2017), which would support the hypothesis that the emission has a hadronic origin, if CR protons collide with the ambient gas. The authors suggested that the extended emission of 4FGL J2108.0+5155e could be related to an old SNR, while the point-like UHE emission could be due to interaction of the escaping CRs from the SNR with the molecular cloud. Alternatively, the authors proposed that the relativistic CRs could be accelerated in one of the nearby open stellar clusters, but confirmation of these hypotheses is complicated because of the unknown distance of the source. A lepto-hadronic emission scenario was also proposed by Kar & Gupta (2022), whereby shock-accelerated electrons and protons were injected in the molecular cloud several thousand years ago during an explosion.

TeV halos in the vicinity of pulsars were recently established as a class of extended VHE sources (Linden et al. 2017; Abeysekara et al. 2017; López-Coto et al. 2022a), featuring bright TeV emission and a hard spectrum (Sudoh et al. 2019). Γ-ray emission in such sources can be produced in IC scattering of ambient photons by VHE electrons and positrons accelerated by the pulsar-wind termination shock (Sudoh et al. 2019). Even though IC γ-ray emission beyond 100 TeV is suppressed because of the Klein-Nishina effect, it has been shown that IC can still dominate UHE emission in radiation-dominated environments (e.g., Vannoni et al. 2009; Breuhaus et al. 2021). This mechanism was used to explain UHE γ-ray emission of extended sources detected by HAWC, and three LHAASO sources associated with pulsars (Breuhaus et al. 2021, 2022). In addition, the study of Albert et al. (2021) suggests that UHE γ-ray emission may be a generic feature in the vicinity of pulsars with high spin-down powers Ė > 10³⁶ erg s⁻¹. The same limit on Ė was also derived from first principles, showing that only very energetic pulsars can power PeV γ-ray emission (de Oña Wilhelmi et al. 2022).

According to the ATNF database¹ (Manchester et al. 2005) there is no detected pulsar within a 1° radius around LHAASO J2108+5157. This does not a priory exclude the PWN/TeV halo scenario, as the pulsar might remain undetected if its beam is not pointing towards us. The spectral analysis of Cao et al. (2021a) showed that a PWN scenario can also explain the observed UHE emission of LHAASO J2108+5157. With a lack of other observational data, and especially the missing VHE counterpart, the nature of the source remains unknown.

In this paper, we present the results of a dedicated observation of the source region with the first Large-Sized Telescope (LST-1) and XMM-Newton², and results of a dedicated analysis of Fermi-LAT data, providing strong constraints on LHAASO J2108+5157 γ-ray and X-ray emission and the physical nature of the source. We also use all these available data to carry out a detailed modeling of the multiwavelength emission of the source, including a discussion of possible emission scenarios.

The paper is structured as follows: in Sect. 2 we describe our detailed analysis of LST-1 and Fermi-LAT data, followed by a dedicated analysis of the XMM-Newton data, and an analysis of ¹²CO(1–0) emission lines in the direction of the source that were collected in a composite survey by Dame et al. (2001). In Sect. 3, we present and discuss the results of multiwavelength spectral modeling of the source. Finally, a summary and conclusions can be found in Sect. 4.

2. Observations and data analysis

2.1. LST-1

LHAASO J2108+5157 was observed with LST-1 (López-Coto et al. 2021, LST Collaboration, in prep.) for 91 hours over 49 nights from June to September 2021. The data were taken in Wobble mode – which allows the background to be evaluated from the same observations (Fomin et al. 1994) – using four positions centered at RA = 317.15°, Dec = 51.95°, that is, at coordinates that lie between the LHAASO source and the possible Fermi-LAT counterpart. The offset of each Wobble with respect to the center of the field of view (FoV) was 0.5°, instead of the standard 0.4°, in order to decrease the number of excess events leaking into the background regions in case the source was found to be extended. The observations presented were taken up to > 55° zenith angle, and with the Moon below horizon. We applied quality cuts based on the stability of the trigger rate, the atmospheric transmission (using MAGIC LIDAR measurements Fruck 2022), and the rate of CR events, resulting in 49.3 h of good-quality data (dead-time-corrected) used for the analysis.

The data calibration and shower reconstruction was carried out using the standard pipeline implemented in lstchain v0.9 (Lopez-Coto et al. 2022b). To separate pixels containing Cherenkov light emitted by the atmospheric shower particles from background noise, we carried out image cleaning, taking into account the pixel-wise noise fluctuations (LST Collaboration, in prep.). This method exploits background level information provided by dedicated interleaved pedestal events (containing only noise) acquired during observations at a rate of 100 Hz to reduce the effect of increased pixel noise due to stars in the FoV for example (LST Collaboration, in prep.). The cleaned shower images were parameterized with Hillas parameters (Hillas 1985). We then applied a random forest (RF) algorithm trained on Monte Carlo (MC) simulated images of γ and proton events in order to separate the γ rays from the hadrons (resulting in an estimation of the so-called Gammaness parameter for each event)³ and to reconstruct the energy and arrival direction of each event, using the Hillas parameters as inputs to the RF. As shower development for a primary particle or photon of given energy depends on the specific orientation of each shower relative to the local magnetic field, and the amount of Cherenkov light collected with the telescope depends on the zenith angle of the shower, we trained the RFs on MC events simulated along a declination line that follows the approximate path of the source in zenith and azimuth during the night (LST Collaboration, in prep.). While a constant low night-sky background (NSB) level was assumed in MC simulations, the real observations are performed in a wide range of NSB conditions. To reach the best possible performance in this particular analysis, we tuned the NSB levels in the training sample of MC-simulated images to that of real data, adopting the so-called “noise padding” method (for further details see LST Collaboration, in prep.).

Global selection cuts on Gammaness⁴ and the squared angular distance between the reconstructed event direction and the source (θ²) – assuming a point-like nature for the source – were optimized on 36 high-quality runs of Crab Nebula observations taken in 2021, applying the same selection criteria as in the LHAASO J2108+5157 source analysis. Crab Nebula is the brightest persistent γ-ray source in the sky, which makes it an ideal target for LST-1 calibration and validation of new data-reconstruction methods (e.g., López-Coto et al. 2021; Alispach et al. 2022; Juryšek et al. 2021). The Crab Nebula detection significance was evaluated on a grid of Gammaness ∈(0.5, 0.98) with a step of 0.02 and θ² ∈(0.01, 0.1) deg² with a step of 0.002 deg², resulting in the best global selection cuts of Gammaness > 0.84 and θ² < 0.04 deg² used in the spectral analysis. In order to reach the best possible performance for a potential source detection, we also optimized the Gammaness cut on Crab detection significance in individual energy bins (five bins per decade) keeping the θ² cut fixed.

Figure 1 shows a θ² plot for three OFF regions reflected with respect to the center of the FOV after selection cuts optimized on Crab detection significance. There is no significant source detection in either of the four energy bins, but at the highest energies (3–100 TeV), we see a hint of a source with a detection significance of 3.67σ (Li & Ma 1983) and a signal-to-noise ratio (S/N) of 46% under the point-like source assumption. This signal is insufficient to claim a detection of the source in the TeV energy range. However, if the excess above 3 TeV proves to be significant in future observations, then no excess seen at lower energies may suggest a hard spectral index of the potential VHE source.

Fig. 1. ON (blue) and OFF (orange) counts detected by the LST-1 telescope after selection cuts in 49.3 hours of effective observation time in four blindly selected energy bins. Number of excess events in the first two θ2 bins for the highest energies is 45 ± 13 with a Li and Ma detection significance of 3.67σ.

We performed a 1D spectral analysis using the Gammapy (Deil et al. 2017) package, adopting the source coordinates reported by Cao et al. (2021a). As LST-1 is in the commissioning phase and the analysis software and methods are still under intensive development, we only performed a preliminary 2D analysis using an acceptance model taken from real LST data, which was then used to correct for radial dependencies in the background models in Gammapy. The results of the 2D analysis can be found in Appendix A, but a detailed analysis of the source morphology is left for future studies.

The 1D spectral analysis was performed in the energy range of 100 GeV to 100 TeV using the ‘reflected regions background’ method to estimate the number of background events in the signal region. We assumed the source spectrum to follow a single PL defined as dN/dE = N₀(E/E₀)^−Γ, with fixed reference energy E₀ = 1 TeV. Having no morphological information and no sign of a potential source extension in the LST-1 data, the best-fit spectral parameters were found under the assumption of a point-like source. The resulting spectral parameters are listed in Table 1. Despite the large statistical uncertainties, the resulting photon indices suggest a relatively hard spectrum in the TeV range. We estimated the systematic uncertainty due to the energy resolution on the fitted photon index Γ to be ∼2%; this is negligible compared to the statistical uncertainty.

We also performed a joint likelihood fit of the LST-1 data and LHAASO flux points to find the spectrum of the source in the multi-GeV to multi-TeV range. Provided the hard TeV and soft multi-TeV spectrum, for the spectral shape, we considered a PL with an exponential cutoff (ECPL) defined as dN/dE = N₀(E/E₀)^−Γexp(−E/E_cutoff). The best-fit parameters are listed in Table 1.

As a second step, we performed a maximum-likelihood estimation of the source flux in six logarithmically spaced energy bins between 100 GeV and 100 TeV, using ECPL spectral parameters fitted in the previous step. We did not reach a significant source detection in the TeV range, and therefore calculated 95% confidence level differential flux upper limits (ULs) in individual energy bins, which are shown in Fig. 2. In the first energy bin (0.1–0.316 TeV), the telescope effective area drops below 10%, which we used as a safe threshold in the analysis, and therefore the Gammapy flux point estimator only calculated the flux in five energy bins. Table 2 summarizes the LST-1-measured flux ULs and the corresponding TS in each energy bin.

Fig. 2. Spectral energy distribution of the LHAASO J2108+5157 source observed with LST-1. The green confidence band represents the best-fitting PL spectral model of LST-1 data and its statistical uncertainties. The blue confidence band shows a joint likelihood fit of the LST-1 data and LHAASO flux points with an ECPL spectral model. The ECPL spectral model was used to estimate the 95% confidence level ULs on the differential fluxes shown in all energy bins.

Although no significant point-like source was detected at energies above 300 GeV (2.2σ for PL spectral model), the LST-1 ULs provide strong constraints on the source emission, which is further discussed in Sect. 3. The resulting integral 95% UL on the flux of the source is F(E > 300 GeV) < 5.0 × 10⁻¹³ ph cm⁻² s⁻¹.

2.2. Fermi-LAT

We performed a dedicated binned analysis of the region around LHAASO J2108+5157 with the publicly available Fermi-LAT science tools⁵. We followed the standard recommendations from the Fermi-LAT team to analyze data from 1 GeV up to 500 GeV within 11° of the LHAASO source in order to take into account the broad instrument point spread function (PSF). We specifically selected data with evclass=128 and evtype=3 and a zenith angle < 90°. The selected data cover the time interval from 4 August 2008 until 31 January 2022, and we filtered out bad time intervals using (DATA_QUAL>0)&&(LAT_CONFIG==1). The data were reprocessed with the most recent P8R3_SOURCE_V3 instrument response function (IRF), and we used gll_iem_v07⁶ for the most updated Galactic diffuse emission model and iso_P8R3_SOURCE_V3_v1 for the isotropic emission. We used a spatial binning of 0.1°/pixel and eight logarithmically uniform energy bins per decade. As opposed to the LHAASO Collaboration, who used the 10-year 4FGL-DR2 catalog (Ballet et al. 2020; Abdollahi et al. 2020), we used the more recent 12-year 4FGL-DR3 catalog (Fermi-LAT Collaboration 2022) in our analysis to create the source model, which includes more sources. More specifically, the LHAASO team manually added two new γ-ray sources in their source model to describe the region of interest (ROI) found with the “find source” tool of fermipy (Wood et al. 2017). One of the two sources visible in their plot (see bottom left square in the left of Fig. 6 of Cao et al. 2021a) is now confirmed in the new DR3 catalog as an interstellar gas clump, but with slightly shifted coordinates (4FGL J2115.2+5219c). The other source does not appear in the new Fermi-LAT catalog, and is also not detected in our analysis. As recommended, we included sources in our model up to 5° beyond our ROI, keeping their parameters frozen in the likelihood fit.

Any source in the catalog with a nonzero flags parameter is affected by systematic errors. Sources with these indicators should be used with great care. These correspond to significant excesses of photons, but such excesses can result from residual extended emission or confused source pile-up. In the 4FGL-DR3 catalog, three sources within less than 1.5° of LHAASO J2108+5157 are flagged with “c”, which stands for “interstellar gas clump”. More than 50% of the neighboring sources within less than 4° of LHAASO J2108+5157 present at least one nonzero flag (for more details see Fermi-LAT Collaboration 2022), which make the low-energy Fermi-LAT analysis nontrivial.

Looking at all sources present in the 4FGL-DR3 catalog and their 95% positional errors, there are no counterparts overlapping with the uncertainty position provided by LHAASO (see Fig. 3). The closest source is 4FGL J2108.0+5155, which lies 0.13° away. Above 1 GeV, the source spectrum is well fitted by a log parabola with a normalization of (9.8 ± 0.9)×10⁻¹³ ph cm⁻² s⁻¹ MeV⁻¹, α = 2.5 ± 0.2 and β = 0.37 ± 0.18, assuming the same E_b value as that of the catalog (i.e., E_b = 1580.67 MeV). The other three 4FGL sources visible in Fig. 3 are fainter and present a softer log-parabolic spectrum with a turnover at lower energies, with the only exception being J2109.5+5238c, whose spectrum is a PL with a photon index of 2.6, which locally overtakes the flux of 4FGL J2108.0+5155 above a few tens of GeV.

Fig. 3. Fermi-LAT TS map in Galactic coordinate above 2 GeV, which shows the sources present in the 4FGL-DR3 catalog with their 95% positional errors (magenta and red ellipses). The small green rectangle indicates the position of the LHAASO source with statistical uncertainty on RA and Dec derived from a two-dimensional Gaussian model, while the smaller green circle represents 95% position uncertainty of 0.14° reported by Cao et al. (2021a). The larger green circle indicates the 95% UL on the source extension (0.26°). The white cross highlights the position of a new potential hard source, whereas the yellow contour indicates the FoV of the previously discussed XMM observation.

The spectrum of the closest source to LHAASO J2108+5157, namely 4FGL J2108.0+5155, presents a steep decrease above a few GeVs, which is not compatible with the UHE LHAASO points. Therefore, its physical relation to the UHE source is challenging (see the discussion in the following sections). By rerunning the analysis, extending the low-energy threshold to 500 MeV and to 300 MeV, and properly increasing and adapting the selected ROI, the fitted spectra that we obtain present some scatter at low energy, which is due to the large instrument PSF. Although it depends on how much freedom we allow in the fit to the neighboring sources and to the Galactic diffuse emission, in all cases the trend converges toward a unique and consistent behavior above a few GeVs.

In order to verify the goodness of the used source model at high energies, we constructed a 15° ×15° TS map centered on the LHAASO source, removing the source 4FGL J2108.0+5155 from the model. We computed the TS map above different threshold energies, from 1 GeV to 10 GeV, and we used a PL spectrum for the putative source, assuming different γ indices (from −1.5 to −3). Some of these TS maps are reported in Fig. 4. Each TS map has been smoothed with a Gaussian with a standard deviation equal to 68% of the Fermi-LAT containment angles at each different threshold energy. From this analysis, we can clearly see that, assuming a very soft photon index, above 1 GeV the peak of the TS map coincides with the position of 4FGL J2108.0+5155, whereas using a harder photon index moves the peak toward the southeast (in Galactic coordinates). This trend becomes even more evident when we move towards higher energies. Already above 2 GeV, the excess of the TS maps assumes an elongated shape toward the southeast, and can no longer be considered as point-like, nor can it be reproduced by an extended symmetric Gaussian. These TS maps (Fig. 4) confirm the very soft spectral behavior of 4FGL J2108.0+5155, whose flux steeply drops above a few GeVs, and suggest the presence of two different sources with clearly distinct spectra, located at two different positions separated by ∼0.4°. One of these sources is 4FGL J2108.0+5155, which is already included in the 4FGL-DR3 catalog, whereas the other is a new hard source (hereafter HS), approximately located at l = 92.35° and b = 2.56°, not included in the catalog. Such sources are difficult to distinguish from one another at low energies because of the relatively large PSF of the Fermi-LAT instrument, and it is not trivial to spatially disentangle them. On the contrary, they are clearly distinguishable above a few GeVs, where the PSF becomes smaller than the two source separations⁷. The existence of two distinct peaks is also evident in the nonsmoothed TS maps. Assuming a flat spectrum, the excess at the position of HS dominates over that of 4FGL J2108.0+5155 above ∼4 GeV. If instead we assume a harder spectrum, a similar transition occurs at even lower energies. It is important to mention that the new HS source does not spatially correlate with the local structure of the diffuse Galactic emission model.

Fig. 4. Fermi-LAT TS maps computed assuming various photon indexes for the putative source, and above different threshold energies. Each TS map is smoothed with a Gaussian whose sigma is equal to 68% of the Fermi-LAT PSF containment radius measured at each corresponding threshold energy. This size is reported with a cyan segment in the bottom left corner of each plot. Black contours are overplotted with a linear scale to better localize the position of the TS peaks. Green, red, magenta, and white elements are the same as those used and described in Fig. 3. The small white ticks on both axes are in units of 0.1°. Each subplot has been renormalized to its own maximum value to make its color-scale and isocontours comparable.

Adding the new HS source in the original source model and rerunning the likelihood fit analysis provides slightly different results for the spectral shape of 4FGL J2108.0+5155, which is now fitted with a log parabola with a normalization of (9.9 ± 0.9)×10⁻¹³ ph cm⁻² s⁻¹MeV⁻¹, α = 2.7 ± 0.2 and β = 0.32 ± 0.16, assuming the same fixed value for E_b. The new HS source is detected with a significance of ∼4σ, and its spectrum can be fitted with a PL with a normalization of (1.5 ± 0.9)×10¹³ ph cm⁻² s⁻¹MeV⁻¹ and a photon index of Γ = 1.9 ± 0.2, using an energy scale of E₀ = 1 GeV. If we fix the photon index, the normalization accuracy of HS improves to (1.5 ± 0.5)×10⁻¹³ ph cm⁻² s⁻¹MeV⁻¹. Due to the HS small flux at low energies, its inclusion in the model does not significantly affect the spectral results of the neighboring sources, in particular at low energies. Using a different model to represent the HS source, such as a log parabola or ECPL does not improve the likelihood fitting, and so the simple PL is preferred, which presents fewer degrees of freedom. The angular separation of this HS from the LHAASO J2108+5157 source is 0.27°, which is larger than the 95% upper limit of the extension provided in Cao et al. (2021a), and is therefore unlikely to be its counterpart.

The SED points of J2108.0+5155 and HS shown in Fig. 5 were computed by running a separate independent likelihood analysis in each smaller energy band, replacing the source of interest with a simple PL spectrum. The normalization of this spectrum was let free to vary in the fit, whereas its photon index was fixed to the local slope (α) of the log parabola in the case of J2108.0+5155, and to the previous obtained photon index Γ in the case of the HS source. The error bar represents 1σ statistical error. The confidence band represents the 1σ error obtained from the covariance matrix of the fit.

Fig. 5. Fermi-LAT SED for J2108.0+5155 (blue) and HS (red) analyzed with the energy threshold of 300 MeV (see the text for details). Fluxes in energy bins with TS > 10 are drawn as flux points, and for lower TS 95% confidence level ULs are shown.

The discrepancy between our flux and that provided by Cao et al. (2021a) can arise from the several differences present between the two analyses, which we highlight in this article. In particular, we used a more recent IRF, a more recent source catalog, and a more recent isotropic diffuse emission component. Furthermore, Cao et al. (2021a) provided the integral flux value assuming a symmetric Gaussian extended source with a radius of 0.48°, and our TS map results suggest this is not a correct assumption (see Fig. 4).

2.3. XMM-Newton

The field surrounding LHAASO J2108+5157 was observed by XMM-Newton on June 11, 2021, for a total of 13.6 ks. The observation was centered on RA(J2000) = 317.0170°, Dec(J2000) = +51.9275°. We reduced the data from the European Photon Imaging Camera (EPIC)⁸ of XMM-Newton using XMMSAS v19.1 and the X-COP data-analysis pipeline (Eckert et al. 2017; Ghirardini et al. 2019). After screening the data and creating calibrated event files using the standard chains, we used the XMMSAS tasks pn-filter and mos-filter to filter out time periods affected by strong soft proton flares. After excising the flaring time periods, the clean exposure time is 4.7 ks (MOS1), 4.9 ks (MOS2), and 3.0 ks (pn). From the clean event files, we extracted images in the soft (0.5–2 keV) and hard (2–7 keV) bands, and used the eexpmap task to create effective exposure maps accounting for vignetting, bad pixels, and chip gaps. To estimate the nonX-ray CR-induced background (NXB), we made use of the unexposed corners of the detectors to rescale the filter-wheel-closed event files available in the calibration database. We then reprojected the filter-wheel-closed data to match the attitude file of our observation and extracted model-particle background images from the rescaled filter-wheel-closed event files. Finally, we summed up the images of the three EPIC instruments as well as the NXB maps. We also created a total EPIC exposure map by summing up the exposure maps of the three instruments, weighted by their respective effective area.

We ran the XMMSAS task ewavelet to detect X-ray sources in the field. The resulting catalog contains 11 sources in the soft band and 2 in the hard band. The brightest source in the field is associated with the eclipsing binary star V* V1061 Cyg (RX J2107.3+5202). We extracted the spectrum of the source and fitted it in XSPEC. We find that the eclipsing binary exhibits a very soft spectrum, which appears better represented by a thermal model with kT ∼ 0.8 keV than with a PL. The total model flux of the source is 3.7 × 10⁻¹³ erg cm⁻² s⁻¹ (0.5–2 keV) and 2.0 × 10⁻¹⁴ erg cm⁻² s⁻¹ (2–7 keV). Assuming a thermal nature of the spectrum, this is an unlikely counterpart for the LHAASO source. The other sources in the field are substantially fainter than the fluxes reported here.

If the source were associated with leptonic emission from a SNR or PWN, we would expect the corresponding X-ray source to be extended. We therefore attempted to place upper limits on the possible X-ray flux of an extended source with varying radius and centered on the LHAASO source position. To this end, we used the public code pyproffit (Eckert et al. 2020) to extract the brightness profile of unresolved emission across the field. The NXB-subtracted, vignetting-corrected profiles are consistent with a flat brightness distribution, which is compatible with the expected sky background emission made of a combination of the local hot bubble, the Galactic halo, and the cosmic X-ray background. If the source is smaller than the XMM-Newton FoV (30′ in diameter), we can use the measured count rates in the outer regions of the field to estimate the local sky background emission and set the maximum allowed number of source counts on top of the background. For aperture photometry performed in the Poisson regime, the S/N is given by

$$\begin{array}{c}\hfill S/N=\frac{N_{\mathrm{src}}}{\sqrt{N_{\mathrm{src}}+N_{\mathrm{bkg}}}},\end{array}$$(1)

where N_src is the number of (background-subtracted) source counts and N_bkg is the total number of background counts (sky background + NXB). A 2σ upper limit can therefore be set by computing the number of source counts for which S/N = 2, that is

$$\begin{array}{c}\hfill N_{2\sigma }=2+2\sqrt{1+N_{\mathrm{bkg}}}.\end{array}$$(2)

For a uniform, circular source, the corresponding upper limit depends on the assumed source radius because the background expectation is proportional to the area. We computed such upper limits for three possible apertures: a radius of 0.5 arcmin (point-like source) and extended sources with radii of 3 and 6 arcmin, respectively, which was the maximum possible radius considering the relatively small FoV of the EPIC camera and that the XMM-Newton observation is not centered at the position of the UHE source. The estimated 2σ upper limits on the allowed number of counts and count rates are given in Table 3.

To convert the upper limits on the count rates into fluxes, a model spectrum needs to be assumed; it is also especially important that we understand whether or not the source is affected by Galactic extinction. In the case where the source is located at the opposite side of the Galaxy, the measured X-ray fluxes will be attenuated by the Galactic N_H, which is substantial along the Galactic plane. The Galactic HI column density at the position of the LHAASO source is 1.05 × 10²² cm⁻² (HI4PI Collaboration 2016), which implies that the majority of the low-energy photons will be absorbed along the way, such that the actual emitted flux can be substantially higher. Conversely, if the source is local, the source spectrum will be mostly unabsorbed. Here, we provide upper limits for the two extreme cases of a completely absorbed and a completely unabsorbed source to bracket all possible scenarios. We simulated XMM-Newton spectra assuming that the source spectrum is a PL with a photon index of 2.0 and computed the conversion between count rate and flux in the two energy bands and for the absorbed and unabsorbed scenarios alike. The conversion factors were then used to determine the corresponding flux upper limits. Table 3 provides a list of upper limits for a wide range of scenarios. The allowed values are in the range 10⁻¹⁵ − 10⁻¹³ erg cm⁻² s⁻¹ depending on the assumed source size and spectrum.

In order to estimate a flux upper limit in the full region of a possible extension of the UHE source from the data available (95% extension upper limits has a radius of 16 arcmin), we scaled the 6 arcmin upper limits by the ratio of the two radii. However, we note that there is no surface brightness gradient across the XMM-Newton FoV, which would be expected for a very extended source centered at the UHE source coordinates. This means that any low-surface-brightness extended source would need to be not only very extended but also uniform across the XMM-Newton FoV, which is not very likely.

2.4. Molecular clouds

LHAASO J2108+5157 is located in the direction of relatively dense molecular clouds. Cao et al. (2021a) searched the database of molecular clouds in the Galactic plane (Miville-Deschênes et al. 2017) and found its position close to the direction of two molecular clouds [MML2017]2870 and [MML2017]4607 with distances of 1.43 kpc and 3.28 kpc, and masses of 3.5 × 10⁴ M_⊙ and 8.5 × 10³ M_⊙, respectively.

To estimate the total gas density in the direction of LHAASO J2108+5157, we considered a contribution of H₂ molecular density and of neutral hydrogen HI density. We used ¹²CO(1–0) line-emission observations collected during the composite survey by Dame et al. (2001) to estimate the density and distance of the molecular clouds. We note that, given the Galactic longitude of the source, namely l > 90°, the kinematic distance of each molecular cloud has only a single solution (see e.g., Roman-Duval et al. 2009). In order to get a mean radial-velocity spectrum of brightness temperature T_B(v) in the UHE source region, we averaged T_B(v) in all pixels within 0.26°, which is 95% ULs on the UHE source extension reported by Cao et al. (2021a). Consistently with Cao et al. (2021a), we identified three peaks in the T_B(H₂, v) spectrum, which can be well described by three Gaussians at v₁ ≈ −11.8 km s⁻¹, v₂ ≈ −2.7 km s⁻¹, and v₃ ≈ 8.4 km s⁻¹. Two of these correspond to the centroid velocities of the molecular clouds, of namely, −1.1 km s⁻¹ and −13.7 km s⁻¹, identified by Miville-Deschênes et al. (2017), but the origin of the last peak remains unknown, and we further consider only the clouds with negative radial velocities.

For determination of the kinematic distances d to the molecular clouds, we adopted a PL rotation curve (Russeil et al. 2017) with a galactocentric radius of R₀ = 8.34 kpc and with the same orbital velocity as the Sun, V₀ = 240 km s⁻¹, leading to d₁ ≈ 3.1 kpc and d₂ ≈ 2.0 kpc. The angular radius of the UHE emission can then be converted to a physical radius of the source, which gives r₁ = 7.1 pc and r₂ = 4.5 pc, depending on its distance, if spatially coincident with the molecular clouds.

In order to obtain the average densities of the molecular clouds in the source region, we integrated the individual Gaussian contributions W_CO over the velocity ranges v₁ ∈ ( − 23, −1) km s⁻¹ and v₂ ∈ ( − 9, 4) km s⁻¹, given by 3σ ranges of the two Gaussian contributions. W_CO can then be converted into the column density as N(H₂) = X_COW_CO/N_px, where X_CO = 2 × 10²⁰ cm⁻² s K⁻¹ km⁻¹ (Bolatto et al. 2013), and N_px is the number of pixels on the sky within the source extension. Assuming a spherically symmetrical UHE source emission region, we estimate the number density of molecular hydrogen to be n₁(H₂) = N₁(H₂)/2r₁ = 51 cm⁻³ and n₂(H₂) = 170 cm⁻³, for the distant and the close molecular cloud, respectively. A sky map of the velocity-integrated ¹²CO intensity is shown in Fig. 6 for both molecular clouds.

Fig. 6. Velocity-integrated 12CO intensity (WCO) of two molecular clouds spatially coincident with the direction of LHAASO J2108+5157. Left: Integrated velocity of the first Gaussian component peaking at v1 ≈ −11.8 km s−1, with corresponding distance of d1 ≈ 3.1 kpc. Right: Integral of the second Gaussian component at v2 ≈ −2.7 km s−1 and d1 ≈ 2.0 kpc. The white contour represents 1420 MHz continuum emission from the Canadian Galactic Plane Survey (Taylor et al. 2003). The position of LHAASO J2108+5157 is marked with a black cross, and 95% UL on its extension (0.26°) is indicated with a black circle (Cao et al. 2021a). Bilinear interpolation is used to smooth out the contributions from individual pixels.

To estimate the density of neutral hydrogen HI, we used the brightness-temperature velocity spectrum T_B(HI, v) in the direction of the source averaged over its angular extension from the HI4π survey (HI4PI Collaboration 2016). The column density of HI can be estimated under the optically thin limit assumption as N(HI) = 1.823∫T_B(HI, v) dv cm⁻² (Dickey & Lockman 1990). Assuming the same distance for both gas phases, we integrated over the same velocity ranges as for the H₂ molecular clouds. Following the same geometrical assumptions as for H₂, we estimate the number density of neutral hydrogen to be n₁(HI) = 64 cm⁻³ and n₂(HI) = 70 cm⁻³. The total number density of the gas in the source region can therefore be estimated for both clouds as n₁ = n₁(HI)+n₁(H₂) = 115 cm⁻³ and n₂ = 240 cm⁻³.

3. Spectral modeling and discussion

3.1. Leptonic scenario of emission

As most of the VHE sources in our Galaxy have been identified as PWNe (Hess & Abdalla 2018c), we first examine this possible scenario of emission. Here, we provide quantitative estimates of and limits on the leptonic scenario of emission derived from the data, and compare them with the physical properties of known PWNe/TeV halos.

We used the naima⁹ package (Zabalza 2015) to derive a parent electron distribution reproducing IC dominated (Khangulyan et al. 2014) VHE to UHE emission of LHAASO J2108+5157. For the sake of simplicity, we assume a single electron spectrum in the form of an ECPL f(E)∼E^−αexp(−E/E_c). As the source was not detected in the X-ray range, we only considered IC emission using the LST-1 and LHAASO flux points (Cao et al. 2021a). The target photon field for IC is expected to consist of cosmic microwave background (CMB) and far-infrared (FIR) radiation of the dust, with temperatures of 2.83 K and 20 K, and energy densities of 0.26 eV cm⁻³ and 0.3 eV cm⁻³, respectively, which are typical values at the galactocentric radius of the Sun (Hinton & Hofmann 2009; Popescu et al. 2017).

The electron distribution that best describes the observations is shown in Fig. 7; it has a cutoff at $E_c={100}_{-30}^{+70}\mathrm{TeV}$ and a spectral index of α = 1.5 ± 0.4. All parameters of the model are summarized in Table 4. If there is magnetic field, the same population of electrons must emit synchrotron radiation (Aharonian et al. 2010), which together with the XMM-Newton upper limits on the X-ray emission allow us to put constraints on the strength of the magnetic field in the PWN. Figure 7 shows 95% XMM-Newton upper limits on the absorbed emission in a circular region with r = 6′, which is a reasonable upper limit on the angular size of a Galactic PWN emitting X-rays at a distance d ≤ 10 kpc (Bamba et al. 2010). The comparison with the 95% CL of synchrotron emission constrains the magnetic field to B ≲ 1.2 μG, which is lower than a typical value of Galactic magnetic field B_G ≈ 3 μG. We note that for unabsorbed X-ray upper limits, which are relevant if the source is relatively close, the constraints on the magnetic field would be even stronger, namely B ≲ 0.5 μG. Given its Galactic latitude of b ≈ 3°, the source is close to the Galactic plane if it is not too distant from the Sun, and one should not expect a background magnetic field strength significantly below the typical level; therefore the absorbed case is favored. The possibility of greater extension of the undetected PWN – which would potentially lead to more relaxed constraints on its magnetic field – cannot be excluded. However, we note that even the approximate absorbed X-ray flux ULs scaled on the full UHE source extension lead to a relatively low B ≲ 1.9 μG compared to the average Galactic magnetic field (also shown in Fig. 7 for reference).

Fig. 7. Multiwavelength SED of LHAASO J2108+5157 showing a leptonic scenario of emission. Observations with different instruments are represented by data points of different colors: XMM-Newtonr = 6′ (blue), r = 16′ (green), Fermi-LAT (red), LST-1 (purple), and LHAASO-KM2A (yellow). The black solid line represents the best-fitting IC-dominated emission of LST-1 and LHAASO data. The corresponding synchrotron radiation of the same population of electrons is represented with dashed and dash-dotted lines for B = 1.2 μG and B = 1.9 μG, respectively. The dotted line represents a phenomenological model of a tentative pulsar: the best-fit PL with a subexponential cutoff on the Fermi-LAT data.

Such a weak magnetic field, needed to suppress the synchrotron emission of LHAASO J2108+5157, is on the lower end of the typical range seen for B_PWN, which is, 1 − 100 μG (Martin et al. 2014; Zhu et al. 2018). However, we note that a relatively weak magnetic field is needed to explain a leptonic UHE emission, which is only possible in radiation-dominated environments (Vannoni et al. 2009; Breuhaus et al. 2021, 2022). MILAGRO (Abdo et al. 2009) and HAWC (Abeysekara et al. 2017) detected an extended 2° TeV emission surrounding the pulsar Geminga, leading to the recent establishment of a new class of TeV-halo sources (Linden et al. 2017; Sudoh et al. 2019). Resulting from propagation of relativistic electrons that already left the PWN in the interstellar medium (ISM), magnetic field in the TeV halos can be expected to follow the level of the magnetic field in the ISM. However, Liu et al. (2019) obtained an upper limit on the magnetic field in the halo of Geminga of B < 1 μG, and therefore the TeV halo scenario for LHAASO J2108+5157 is also feasible. The TeV nebula surrounding Geminga has a large angular extension, but this pulsar is also relatively close (d = 250 pc Faherty et al. 2007). In the Geminga-like scenario, the lower limit on the distance of LHAASO J2108+5157 is approximately 2 kpc in order not to violate the source-extension UL of 0.26° provided by Cao et al. (2021a).

Inverse-Compton-dominated radiation of a single electron population cannot explain the soft GeV emission of 4FGL J2108.0+5155, which is spatially coincident with LHAASO J2108+5157. There are 117 γ-ray pulsars identified in the Fermi-LAT data showing similar spectral properties to 4FGL J2108.0+5155 (Abdo et al. 2013). We therefore put forward the hypothesis that the GeV emission is the signature of a γ-ray pulsar. Saz Parkinson et al. (2016) applied machine learning methods to classify sources in the Third Fermi-LAT catalog in two major classes: AGNs and pulsars. 3FGL J2108.1+5202, which is the Third Fermi-LAT general catalog (Acero et al. 2015) counterpart of 4FGL J2108.0+5155, was classified consistently with logistic regression (LR) and RF classifiers as a pulsar, which provides support for our hypothesis. However, we note that the resulting LR and FR probabilities are relatively low, that is, only about 30%, and therefore we cannot exclude the possibility of misclassification of the source, and an extragalactic origin of the HE emission cannot be excluded (for further details see Saz Parkinson et al. 2016).

Gamma-ray pulsars are characterized by soft spectra, with the flux steeply falling above a few GeV (e.g., MAGIC Collaboration 2020). In the Fermi-LAT energy band, the typical differential spectrum can be described with a PL with a subexponential cutoff dN/dE = N₀(E/E₀)^−Γexp(−(E/E_cutoff)^b), where E₀ is the energy scale, Γ the photon index, E_cutoff the cutoff energy, and b the cut-off strength (Leung et al. 2014; MAGIC Collaboration 2020). In order to reduce the degeneracy of the model parameters, considering that there are only three significant Fermi-LAT flux points, we fixed b = 0.7, which is the cut-off strength of the PL with a subexponential cutoff model of the Geminga pulsar SED in the GeV band (MAGIC Collaboration 2020). The best fit of the Fermi-LAT data consistent with XMM-Newton ULs shown in Fig. 7 has $\mathrm{\Gamma }=1.5_{-0.2}^{+0.1}$ and E_cutoff = 0.9 ± 0.2 GeV. Despite the large uncertainty, the photon index is consistent with that of γ-ray pulsars with a spin-down power of Ė = 10³⁴ − 10³⁷ erg s⁻¹ (Abdo et al. 2013). The γ-ray luminosity of 4FGL J2108.0+5155 of between 1 and 100 GeV is L_{1 − 100 GeV} = 2 × 10³³(d/1 kpc)² erg s⁻¹. One should note that, assuming a Galactic origin for the source, the Galactic longitude of the source, of namely l = 92.2°, implies an UL on the source distance of d ≲ 8 kpc, because of the geometry between the position of the Sun and the edge of the Galaxy at the Galactic longitude of the source. Such UL on the source distance leads to an UL on the source luminosity of L_{1 − 100 GeV} ≲ 1.3 × 10³⁵ erg s⁻¹, which is consistent with the population of γ-ray pulsars for any possible distance to the source.

The allowed spin-down power range of the tentative pulsar also implies limits on the luminosity in the TeV band for the nebula powered by the pulsar of L_{1 − 10 TeV} ≈ 10³¹ − 10³⁵ erg s⁻¹ when compared with the sample of known TeV PWNe (see Fig. 7 in Hess & Abdalla 2018b). The TeV luminosity of LHAASO J2108+5157 is L_{1 − 10 TeV} ≈ 6 × 10³²(d/1 kpc)² erg s⁻¹, which makes the pulsar scenario valid for any possible distance of the source within the Galaxy. In the Geminga-like scenario, assuming the source extension of the order of 10 pc, the minimum possible distance of 2 kpc implies a TeV PWN luminosity of L_{1 − 10 TeV} > 2 × 10³³ erg s⁻¹ and therefore a pulsar spin-down power of Ė > 10³⁵ erg s⁻¹, which is consistent with the estimates derived from the SED of the Fermi-LAT counterpart.

The total energy in electrons dominating the TeV emission is E(> 1 GeV)≈1 × 10⁴⁵(d/1 kpc)² erg, which can be compared with the total energy released by the pulsar during relativistic electron cooling time, given by

$$\begin{array}{cc}\hfill t_{\mathrm{cool},\mathrm{yr}}& ={(1/t_{\mathrm{IC},\mathrm{yr}}+1/t_{\mathrm{syn},\mathrm{yr}})}^{-1},\hfill \end{array}$$(3)

$$\begin{array}{cc}\hfill t_{\mathrm{IC},\mathrm{yr}}& \approx 3.1\times {10}^5{(1+40E_{\mathrm{e},\mathrm{TeV}}{(kT)}_{\mathrm{eV}})}^{1.5}{u}_{\mathrm{rad},\mathrm{eV}{\mathrm{cm}}^{-3}}^{-1}{E}_{\mathrm{e},\mathrm{TeV}}^{-1},\hfill \end{array}$$(4)

$$\begin{array}{cc}\hfill t_{\mathrm{syn},\mathrm{yr}}& \approx 1.3\times {10}^7{E}_{\mathrm{e},\mathrm{TeV}}^{-1}{B}_{\mu \mathrm{G}}^{-2},\hfill \end{array}$$(5)

where E_e is electron energy, and u_rad and T are the energy density and temperature of the radiation field, respectively (Moderski et al. 2005; Hinton & Hofmann 2009). For the cutoff-energy electrons E_c = 100 TeV, CMB radiation field, and B = 1.9 μG, we get t_cool ≈ 20 kyr. The total energy released by the pulsar would be in the range of E ≈ 6 × 10⁴⁶ − 6 × 10⁴⁸ erg, for the spin-down power in the range of Ė ≈ 10³⁵ − 10³⁷ erg s⁻¹. Provided that the distance to the source is d ≤ 5 kpc, only a small fraction of the total energy released by the pulsar during its lifetime needs to be invested in acceleration of the HE electrons, even for Ė on the lower end of the allowed spin-down powers. For Ė ≳ 10³⁶ erg s⁻¹, the pulsar could power the IC emission for any possible distance of the source in the Galaxy.

We also note that there is an unidentified point-like radio source, NVSS 210803+515255 (Condon et al. 1998) or WENSS B2106.4+5140 (de Bruyn et al. 2000), well within the 95% error ellipse of the source 4FGL J2108.0+5155. We did not use the radio flux points to further constrain the HE emission of the tentative pulsar as, in general, it does not follow PL from radio to sub-GeV.

3.2. Hadronic scenario of emission

The absence of an X-ray counterpart supports a hadronic emission mechanism, where UHE γ rays are produced through inelastic collisions of hadronic relativistic CR with thermal protons followed by π⁰ decay (see Kafexhiu et al. 2014). There is an open stellar cluster, Kronberger 80, and an open cluster candidate, Kronberger 82 (Kronberger et al. 2006), in the close vicinity of LHAASO J2108+5157. These open clusters may potentially act as sources of accelerated PeV protons, as young massive stars seem to contribute significantly to the Galactic CR (Aharonian et al. 2019). Interaction of the PeV protons with ambient gas leads to γ-ray emission via π⁰ decay, and therefore may be responsible for the observed emission of LHAASO J2108+5157, which is spatially coincident with two relatively dense molecular clouds. However, the 95% confidence interval for the distance to Kronberger 80 derived from Gaia DR2 data is 7.9–13.7 kpc (Cantat-Gaudin & Anders 2020), which puts the mutual distance of the cluster and the more distant cloud (d₁ ≈ 3.1 kpc) at more than 4.8 kpc, excluding the possibility of protons accelerated in this stellar cluster interacting with both molecular clouds. We note that the same can be concluded for the previously reported distance of Kronberger 80, which is 4.98 kpc (Kharchenko et al. 2016). On the other hand, the distance to Kronberger 82 remains unknown and one cannot exclude that this open cluster is close to one of those molecular clouds. Another source of CR protons are SNRs, which have been shown to accelerate protons emitting γ rays of multi-TeV energies (Ackermann et al. 2013), which diffuse away from their acceleration sites and interact with ISM and molecular clouds, as observed in several cases (i.e., W28; Aharonian et al. 2008 or IC443 Albert et al. 2007).

Assuming a π⁰ decay-dominated origin of the UHE emission, we used the jetset v1.2.2¹⁰ (Tramacere et al. 2009, 2011; Tramacere 2020) package and its implemented frequentist fitting routines (based on iminuit) to fit the LST-1 and LHAASO. The hadronic proton-proton (pp) model implemented in jetset is based on the parametrization of Kelner et al. (2006), and takes into account the γ-ray emission from π⁰ decay, radiation (synchrotron, IC, and bremsstrahlung) from the secondaries (evolved to equilibrium) of charged pions, and neutrinos. This means that, for each step of the minimization, secondaries pairs temporarily evolve to equilibrium state, taking into account synchrotron radiation and escape timescales. In Fig. 8, we show for reference the π⁰ decay emission model with α = 2. and γ_min = 1 proposed by Cao et al. (2021a), which we find to be ruled out by the Fermi-LAT and LST-1 ULs, suggesting a very hard spectrum or a high value of the low-energy cut-off in the proton distribution responsible for the pp-induced γ-ray emission. Based on this observational evidence, we frame our model in the phenomenological scenario proposed by Celli et al. (2019), that is, we assume that the hard spectrum (dictated by the LST-1 upper limits) is generated by protons escaping a shock around a middle-aged SNR, and illuminating the molecular cloud. Celli et al. (2019) showed that the escaped proton distribution can be very narrow, with a low-energy cut-off of up to ≈10⁵ GeV. In order to test such a scenario, we model the proton distribution with a simple PL, we leave the value of γ_min free during the fit, and we freeze the value of the spectral index to α = 2.75, which is similar to that of the galactic CR in this energy range (Acero et al. 2016a). In this simplified scenario, γ_min is mimicking the energy break investigated in Celli et al. (2019), and α represents the index of the protons escaping from the acceleration region. This index is expected to soften with respect to the index of the confined accelerated proton (α_acc; Celli et al. 2019). Table 5 shows the best-fit model parameters shown in Fig. 8 for both molecular clouds in the direction of the source, where the distances, total number densities of the target protons, and sizes of the emitting regions were fixed on the values reported in Sect. 2.4, In particular, we notice that the value of γ_min and α fit well within the parameter space investigated in Celli et al. (2019) for the case of α_acc ∈ [2, 2.3], and are in agreement with the expectations from standard DSA theory (e.g., Bell 1978). The total required energy of all relativistic protons interacting with molecular cloud is E_T, 1 = 7.5 × 10⁴⁶ erg and E_T, 2 = 1.5 × 10⁴⁶ erg, assuming the interaction of protons with the more distant and closer molecular cloud, respectively. This is well below the energy content of CR protons interacting with molecular clouds in the vicinity of W28 and IC 443, which is 1%−10% of the total energy of a typical SN, which is E_SN ≈ 10⁵¹ erg (Ackermann et al. 2013; Cui et al. 2018).

Fig. 8. Multiwavelength SED of LHAASO J2108+5157 with hadronic models of emission. Observations with different instruments are represented by datapoints of different colors: XMM-Newtonr = 16′ (blue), Fermi-LAT (green), LST-1 (red), and LHAASO-KM2A (purple). The best-fitting hadronic model of VHE-UHE emission (solid line) with fixed spectral index of the proton PL distribution α = 2.75 has γmin = 1.6 × 105 for both clouds. Dashed line represents the total neutrino flux for both clouds. Black and red dash-dotted lines represent the synchrotron emission of secondary particles for cloud 1 and cloud 2, respectively. The gray dash-dotted line represents the π0 decay emission model with α = 2 and γmin = 1 shown for reference (Cao et al. 2021a).

The total neutrino flux resulting from π^+/− decay is comparable with the γ-ray flux in the TeV range (see Fig. 8), which makes this source an interesting candidate for a follow-up analysis of data from a neutrino detector in this region. However, we note that the sensitivity of current neutrino detectors is about an order of magnitude lower than the predicted neutrino flux, and only future instruments will have the potential to definitively confirm or reject the hadronic emission hypothesis (e.g., Grant et al. 2019).

The HE γ-ray emission cannot be explained in a single-component hadronic scenario together with VHE-UHE emission. Cao et al. (2021a) suggested that the spectrum of the extended source 4FGL J2108.0+5155e may be associated with an old SNR, which usually features a soft spectrum above 1 GeV (Acero et al. 2016b; Li et al. 2021). However, our dedicated analysis of Fermi-LAT data shows that the Gaussian extended-source assumption is not correct. Fitting the SED of 4FGL J2108.0+5155 in the Fermi-LAT energy band above 1 GeV with a single PL provides a photon index of Γ = 3.2 ± 0.2, which in turn tends to be too soft compared to the observations of old SNRs interacting with dense molecular clouds (see Yuan et al. 2012). One might also consider a significant contribution of HE emission from the sea of CRs. The energy density of CRs at Galactocentric radii > 8 kpc is u_CR(E > 1 GeV)≈0.5 eV cm⁻³ (Yang et al. 2016), which is lower than observed in the Solar System. Considering the angular resolution of Fermi-LAT at 1 GeV, which is about 0.9°, the UL on the radius of 4FGL J2108.0+5155 to still appear as a point-like source can be written as R < dtan(0.9° ). This results in rather weak limits on the proton energy density u_p, 1(E > 1 GeV) > 0.14 eV cm⁻³ and u_p, 2(E > 1 GeV) > 0.10 eV cm⁻³, for the more distant and closer molecular cloud, respectively, and a hadronic origin for the HE emission therefore cannot be excluded.

4. Summary and conclusions

In this work, we present a multiwavelength study of the unidentified UHE γ-ray source LHAASO J2108+5157, which has not yet been found to be associated with any PWN, SNR, or pulsar. Dedicated observations of the source with LST-1 – yielding 49 hours of good-quality data – resulted in a hint (2.2σ) of hard-spectrum emission at energies between 300 GeV and 100 TeV, which can be described with a single PL with a photon index of Γ = 1.62 ± 0.23. Our data analysis with selection cuts optimized for source detection show a possible excess (3.7σ) at energies E > 3 TeV. Although a confirmed detection of the VHE emission would require deeper observations, the LST-1 data provide important constraints on the source emission in the TeV range.

The VHE-UHE γ-ray emission can be well described with IC-dominated emission of relativistic electrons with a spectral index of α = 1.5 ± 0.4 and a cut-off energy of $E_{\mathrm{c}}={100}_{-30}^{+70}$ TeV, favouring the PWN scenario. However, there is no sign of any X-ray source in 13.6 ks of dedicated observation with XMM-Newton, which puts strong constraints on the magnetic field in the emission region, B ≲ 1.2 − 1.9 μG, depending on the angular extension of the X-ray-emitting region. Such a weak magnetic field is on the lower end of a typical magnetic field in PWNe, and also compatible with a magnetic field in the TeV halo around Geminga. A detailed morphological study of the region with a high-resolution instrument, such as the future completed CTA observatory, or a deeper X-ray observation would shed more light on the nature of the source and help to distinguish between the PWN and TeV-halo hypotheses.

The lack of detection of a pulsar associated with the UHE source presents another challenge for the PWN/TeV-halo scenario. Our dedicated analysis of the 12 years of Fermi-LAT data allowed us to precisely determine the spectral properties of the HE source 4FGL J2108.0+5155, which is spatially consistent with LHAASO J2108+5157. The spectral analysis shows that the HE emission is compatible with the spectral properties of the population of known γ-ray pulsars, which are characterized with soft spectra. Limits on the total energy released by the tentative pulsar derived from the luminosity of the source and relativistic electron cooling time significantly exceed the total energy in relativistic electrons, and such an object would therefore be able to power the VHE-UHE emission.

The presence of two molecular clouds in the direction of LHAASO J2108+5157 supports a hadronic scenario of emission, with the UHE γ rays being produced by the interaction of relativistic protons accelerated in a close stellar cluster or SNR with one of the two molecular clouds. In such a scenario, X-ray ULs on the synchrotron radiation of secondaries would allow significantly higher levels of magnetic field in the source compared to the leptonic scenario of emission. We performed a detailed spectral modeling under the assumption of π⁰ decay-dominated emission based on the phenomenological scenario proposed by Celli et al. (2019), that is, assuming that the hard spectrum in the LST spectral window is due to protons escaping a shock around a middle-aged SNR and illuminating the molecular cloud. The best-fit value of γ_min ≈ 1.5 × 10⁵ agrees with one of the scenarios propose by Celli et al. (2019), which predicts a narrow distribution of escaped protons, with a low-energy cut-off up to ≈10⁵ GeV. The proposed model reproduces the observed broad-band SED reasonably well, but a followup analysis is still needed to get a better constraint on the spectrum of the escaped proton illuminating the molecular cloud, and will be presented in a future publication. For example, a narrower proton distribution and/or a softer spectrum in the escaped protons can still account for the observed SED, with implications as to the age of the possible SNR. Nonetheless, the current analysis shows the potential of the LST in providing robust observational constraints useful for testing theoretical frameworks, and reducing the degeneracy in the parameter space.

The HE γ-ray emission in the hadronic scenario cannot be explained by a single-component model together with VHE-UHE emission. While the old SNR scenario seems to be unlikely due to the very soft spectral index of the emission, the lower limit on the relativistic proton energy density is compatible with an interaction between the sea of CRs and the molecular clouds.

The total neutrino flux predicted in our model is comparable with the γ-ray flux in the TeV range, which makes LHAASO J2108+5157 an interesting candidate for future neutrino experiments of sufficient sensitivity to either confirm or reject the hadronic scenario of emission.

Acknowledgments

We gratefully acknowledge financial support from the following agencies and organisations: Conselho Nacional de Desenvolvimento Científico e Tecnológico (CNPq), Fundação de Amparo à Pesquisa do Estado do Rio de Janeiro (FAPERJ), Fundação de Amparo à Pesquisa do Estado de São Paulo (FAPESP), Fundação de Apoio à Ciência, Tecnologia e Inovação do Paraná – Fundação Araucária, Ministry of Science, Technology, Innovations and Communications (MCTIC), Brasil; Ministry of Education and Science, National RI Roadmap Project DO1-153/28.08.2018, Bulgaria; Croatian Science Foundation, Rudjer Boskovic Institute, University of Osijek, University of Rijeka, University of Split, Faculty of Electrical Engineering, Mechanical Engineering and Naval Architecture, University of Zagreb, Faculty of Electrical Engineering and Computing, Croatia; Ministry of Education, Youth and Sports, MEYS LM2015046, LM2018105, LTT17006, EU/MEYS CZ.02.1.01/0.0/0.0/16_013/0001403, CZ.02.1.01/0.0/0.0/18_046/0016007 and CZ.02.1.01/0.0/0.0/16_019/0000754, Czech Republic; CNRS-IN2P3, France; Max Planck Society, German Bundesministerium für Bildung und Forschung (Verbundforschung/ErUM), Deutsche Forschungsgemeinschaft (SFBs 876 and 1491), Germany; Istituto Nazionale di Astrofisica (INAF), Istituto Nazionale di Fisica Nucleare (INFN), Italian Ministry for University and Research (MUR); ICRR, University of Tokyo, JSPS, MEXT, Japan; JST SPRING - JPMJSP2108; Narodowe Centrum Nauki, grant number 2019/34/E/ST9/00224, Poland; The Spanish groups acknowledge the Spanish Ministry of Science and Innovation and the Spanish Research State Agency (AEI) through the government budget lines PGE2021/28.06.000X.411.01, PGE2022/28.06.000X.411.01 and PGE2022/28.06.000X.711.04, and grants PGC2018-095512-B-I00, PID2019-104114RB-C31, PID2019-107847RB-C44, PID2019-104114RB-C32, PID2019-105510GB-C31, PID2019-104114RB-C33, PID2019-107847RB-C41, PID2019-107847RB-C43, PID2019-107988GB-C22; the “Centro de Excelencia Severo Ochoa” program through grants no. SEV-2017-0709, CEX2019-000920-S; the “Unidad de Excelencia María de Maeztu” program through grants no. CEX2019-000918-M, CEX2020-001058-M; the “Ramón y Cajal” program through grant RYC-2017-22665; the “Juan de la Cierva-Incorporación” program through grant no. IJC2019-040315-I. They also acknowledge the “Programa Operativo” FEDER 2014-2020, Consejería de Economía y Conocimiento de la Junta de Andalucía (Ref. 1257737), PAIDI 2020 (Ref. P18-FR-1580) and Universidad de Jaén; “Programa Operativo de Crecimiento Inteligente” FEDER 2014-2020 (Ref. ESFRI-2017-IAC-12), Ministerio de Ciencia e Innovación, 15% co-financed by Consejería de Economía, Industria, Comercio y Conocimiento del Gobierno de Canarias; the “CERCA” program of the Generalitat de Catalunya; and the European Union’s “Horizon 2020” GA:824064 and NextGenerationEU; We acknowledge the Ramon y Cajal program through grant RYC-2020-028639-I; State Secretariat for Education, Research and Innovation (SERI) and Swiss National Science Foundation (SNSF), Switzerland; the research leading to these results has received funding from the European Union’s Seventh Framework Programme (FP7/2007-2013) under grant agreements No 262053 and No 317446. This project is receiving funding from the European Union’s Horizon 2020 research and innovation programs under agreement No 676134. This research has made use of the observations obtained with XMM-Newton, an ESA science mission with instruments and contributions directly funded by ESA Member States and NASA. This research has made use of the SIMBAD database, operated at CDS, Strasbourg, France. Author Contribution: M. Balbo: Fermi-LAT data analysis, paper drafting and edition; D. Eckert: XMM-Newton data analysis, paper drafting and edition; J. Juryšek: project leadership, LST data analysis, paper drafting and edition, theoretical modelling, theoretical interpretation; A. Tramacere: theoretical modelling, theoretical interpretation; G. Pirola: LST analysis cross-check; R. Walter: project initiator & leadership, theoretical interpretation, XMM ToO. 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; discussion and approval of the contents of the draft.

References

Appendix A: Two-dimensional LST-1 analysis

Here, we present the results of a first attempt to build significance and excess sky maps. At present, LST-1 does not have a background model that could be used for background-count predictions. Nevertheless, a tool¹¹ has been developed for the creation of an acceptance model from real data that can be used for radial corrections in Gammapy background models. In our case, the background method adopted for the 2D analysis is the ring-background technique: in this method the OFF region is defined as a ring around a trial source position, and is used to provide a background estimate (Berge et al. 2007). As opposed to the reflected-region-background model, in which we use the same offset between the ON and OFF regions and the pointing direction, with this method the detector acceptance cannot be assumed to be constant within the ring, because it covers areas with different offsets from the pointing position (Berge et al. 2007).

As mentioned, the radial acceptance model (which assumes the instrument response to be symmetrical under rotations around the pointing direction) is extracted directly from the data set under analysis on a run-by-run basis and is then projected onto the sky in order to evaluate the background for the whole data set (de Naurois 2021). We excluded the three putative γ-ray sources (LHAASO J2108+5157, 4FGL J2108.0+5155 and HS) in order to avoid contamination of the acceptance (de Naurois 2021).

As expected, we noticed that the background acceptance rapidly decreases with energy; on the other hand, it does not significantly change with the offset from the center of the FoV at the high energies we are looking at, because showers from high-energy gamma rays have a larger Cherenkov photon density on the ground than the ones produced by low-energy γ-rays, and the telescope is triggered by showers with larger impact parameters and hence large angular offsets (Berge et al. 2007).

Finally, in order to take into account the acceptance dependence on the zenith angle variation across the FoV, we divided the sample into four zenith bins and stacked them together according to the average run zenith angle.

In Figure A.1, the significance and excess maps show an excess of signal above 3 TeV up to significance values above 4σ in the proximity of LHAASO J2018+5157. The Gaussian fit over the significance distribution for the off bins (see Figure A.1) gives an average value consistent with 0, but the distribution is slightly asymmetric, meaning that the background modeling for our dataset is not perfect and could be improved.

Fig. A.1. Statistical significance (left) and excess (center) maps in a region of 2° ×2° around LHAASO J2108+5157 in the range of energy between 3 TeV and 100 TeV. The map was produced using the Gammapy tool for the ring-background model: the OFF region is defined as a 0.3° radius ring, with a width of 0.3°, around a 0.2° radius circular ON region. The correlation radius for the building of the maps was fixed at 0.2° according to the value of the PSF (68% containment) in this range of energy, averaged for the different zenith bins used for the production of the IRFs. The same exclusion region was used as that adopted to build the acceptance model. Right: Significance distribution for the off bins (purple) and for all bins (pink). The black line represents a Gaussian fit of the background, providing the mean value μ = −0.03, and standard deviation σ = 1.15.

An LST-1 background modeling program is still under development, but it is already possible to use it to extract the radial acceptance from real data and use it for radial corrections on the ring-background tool implemented in gammapy. Although the obtained background model is not found to be perfect for our dataset, we find significance values in the skymap that are reasonably consistent with those achieved through the 1D signal extraction above 3 TeV under the point-like-source assumption (see Figure 1).

All Tables

All Figures

	Fig. 1. ON (blue) and OFF (orange) counts detected by the LST-1 telescope after selection cuts in 49.3 hours of effective observation time in four blindly selected energy bins. Number of excess events in the first two θ2 bins for the highest energies is 45 ± 13 with a Li and Ma detection significance of 3.67σ.
In the text

Fig. 2. Spectral energy distribution of the LHAASO J2108+5157 source observed with LST-1. The green confidence band represents the best-fitting PL spectral model of LST-1 data and its statistical uncertainties. The blue confidence band shows a joint likelihood fit of the LST-1 data and LHAASO flux points with an ECPL spectral model. The ECPL spectral model was used to estimate the 95% confidence level ULs on the differential fluxes shown in all energy bins.

In the text

Fig. 3. Fermi-LAT TS map in Galactic coordinate above 2 GeV, which shows the sources present in the 4FGL-DR3 catalog with their 95% positional errors (magenta and red ellipses). The small green rectangle indicates the position of the LHAASO source with statistical uncertainty on RA and Dec derived from a two-dimensional Gaussian model, while the smaller green circle represents 95% position uncertainty of 0.14° reported by Cao et al. (2021a). The larger green circle indicates the 95% UL on the source extension (0.26°). The white cross highlights the position of a new potential hard source, whereas the yellow contour indicates the FoV of the previously discussed XMM observation.

In the text

Fig. 4. Fermi-LAT TS maps computed assuming various photon indexes for the putative source, and above different threshold energies. Each TS map is smoothed with a Gaussian whose sigma is equal to 68% of the Fermi-LAT PSF containment radius measured at each corresponding threshold energy. This size is reported with a cyan segment in the bottom left corner of each plot. Black contours are overplotted with a linear scale to better localize the position of the TS peaks. Green, red, magenta, and white elements are the same as those used and described in Fig. 3. The small white ticks on both axes are in units of 0.1°. Each subplot has been renormalized to its own maximum value to make its color-scale and isocontours comparable.

In the text

	Fig. 5. Fermi-LAT SED for J2108.0+5155 (blue) and HS (red) analyzed with the energy threshold of 300 MeV (see the text for details). Fluxes in energy bins with TS > 10 are drawn as flux points, and for lower TS 95% confidence level ULs are shown.
In the text

Fig. 6. Velocity-integrated 12CO intensity (WCO) of two molecular clouds spatially coincident with the direction of LHAASO J2108+5157. Left: Integrated velocity of the first Gaussian component peaking at v1 ≈ −11.8 km s−1, with corresponding distance of d1 ≈ 3.1 kpc. Right: Integral of the second Gaussian component at v2 ≈ −2.7 km s−1 and d1 ≈ 2.0 kpc. The white contour represents 1420 MHz continuum emission from the Canadian Galactic Plane Survey (Taylor et al. 2003). The position of LHAASO J2108+5157 is marked with a black cross, and 95% UL on its extension (0.26°) is indicated with a black circle (Cao et al. 2021a). Bilinear interpolation is used to smooth out the contributions from individual pixels.

In the text

Fig. 7. Multiwavelength SED of LHAASO J2108+5157 showing a leptonic scenario of emission. Observations with different instruments are represented by data points of different colors: XMM-Newtonr = 6′ (blue), r = 16′ (green), Fermi-LAT (red), LST-1 (purple), and LHAASO-KM2A (yellow). The black solid line represents the best-fitting IC-dominated emission of LST-1 and LHAASO data. The corresponding synchrotron radiation of the same population of electrons is represented with dashed and dash-dotted lines for B = 1.2 μG and B = 1.9 μG, respectively. The dotted line represents a phenomenological model of a tentative pulsar: the best-fit PL with a subexponential cutoff on the Fermi-LAT data.

In the text

Fig. 8. Multiwavelength SED of LHAASO J2108+5157 with hadronic models of emission. Observations with different instruments are represented by datapoints of different colors: XMM-Newtonr = 16′ (blue), Fermi-LAT (green), LST-1 (red), and LHAASO-KM2A (purple). The best-fitting hadronic model of VHE-UHE emission (solid line) with fixed spectral index of the proton PL distribution α = 2.75 has γmin = 1.6 × 105 for both clouds. Dashed line represents the total neutrino flux for both clouds. Black and red dash-dotted lines represent the synchrotron emission of secondary particles for cloud 1 and cloud 2, respectively. The gray dash-dotted line represents the π0 decay emission model with α = 2 and γmin = 1 shown for reference (Cao et al. 2021a).

In the text

Fig. A.1. Statistical significance (left) and excess (center) maps in a region of 2° ×2° around LHAASO J2108+5157 in the range of energy between 3 TeV and 100 TeV. The map was produced using the Gammapy tool for the ring-background model: the OFF region is defined as a 0.3° radius ring, with a width of 0.3°, around a 0.2° radius circular ON region. The correlation radius for the building of the maps was fixed at 0.2° according to the value of the PSF (68% containment) in this range of energy, averaged for the different zenith bins used for the production of the IRFs. The same exclusion region was used as that adopted to build the acceptance model. Right: Significance distribution for the off bins (purple) and for all bins (pink). The black line represents a Gaussian fit of the background, providing the mean value μ = −0.03, and standard deviation σ = 1.15.

In the text