11institutetext: Erlangen Centre for Astroparticle Physics, Friedrich-
Alexander-Universität Erlangen-Nürnberg, Nikolaus-Fiebiger-Straße 2, 91058
Erlangen, Germany
11email<EMAIL_ADDRESS>
# LHAASO J2108+5157 as a Molecular Cloud Illuminated by a Supernova Remnant

A.M.W. Mitchell https://orcid.org/0000-0003-3631-564811

(Received –; accepted –)

###### Abstract

Context. The search for Galactic PeVatrons - astrophysical accelerators of
cosmic rays to PeV energies - has entered a new phase in recent years with the
discovery of the first Ultra-High-Energy (UHE, $E>100$ TeV) gamma-ray sources
by the HAWC and LHAASO experiments. Establishing whether the emission is
leptonic or hadronic in nature, however, requires multiwavelength data and
modelling studies. Among the currently known UHE sources, LHAASO J2108+5157 is
an enigmatic source without clear association to a plausible accelerator, yet
spatially coincident with molecular clouds.

Aims. We investigate the scenario of a molecular cloud illuminated by cosmic
rays accelerated in a nearby supernova remnant (SNR) as an explanation for
LHAASO J2108+5157. We aim to constrain the required properties of the SNR as
well as which of the clouds identified in the vicinity is the most likely
association.

Methods. We use a model for cosmic ray acceleration in SNRs, their transport
through the interstellar medium and subsequent interaction with molecular
material, to predict the corresponding gamma-ray emission. The parameter space
of SNR properties is explored to find the most plausible parameter combination
that can account for the gamma-ray spectrum of LHAASO J2108+5157.

Results. In the case that a SNR is illuminating the cloud, we find that it
must be young ($<10$ kyr) and located within $40-60$ pc of the cloud. A SN
scenario with a low Sedov time is preferred, with a maximum proton energy of 3
PeV assumed. No SNRs matching these properties are currently known, although
an as yet undetected SNR remains feasible. The galactic CR sea is insufficient
to solely account for the observed flux, such that a PeVatron accelerator must
be present in the vicinity.

###### Key Words.:

Gamma ray: ISM – ISM: supernova remnants – ISM: clouds – (ISM:) cosmic rays

##  1 Introduction

Cosmic Rays (CRs) are energetic particles originating from astrophysical
accelerators and continuously arriving at Earth. The all-particle CR spectrum
exhibits a spectral softening at $\sim$ 1-3 PeV, known as the ‘knee’,
generally understood to indicate the start of the transition from galactic to
extragalactic accelerators being responsible for the bulk of CRs (Ginzburg &
Syrovatskii, 1964; Hillas, 2006; Apel et al., 2013; Parizot, 2014).
Astrophysical sources capable of accelerating particles to PeV energies are
known colloquially as ‘PeVatrons’. Gamma rays produced as a consequence of
particle interactions at the source typically have energies a factor $\sim
1/10$ that of the parent particle population (Blumenthal & Gould, 1970).
Hence, the detection of gamma rays with $E>100$ TeV indicates the presence of
particles with PeV energies, corresponding to the CR ‘knee’.

Definitive evidence for the presence of PeVatrons in our galaxy has, however,
proven elusive. Although diffusive shock acceleration of CRs at supernova
remnants (SNRs) can account for the energy budget of CRs in our galaxy, their
gamma-ray spectra cut-off at energies below 100 TeV (Bell, 1978; Lagage &
Cesarsky, 1983; Bell, 2013; H.E.S.S. Collaboration et al., 2018). This
indicates that active acceleration of particles to PeV energies is not
occurring at these SNRs, although the detection of the characteristic ‘pion-
bump’ signature of neutral pion decay in several SNRs indicates that the
emission is hadronic in origin (Bell et al., 2013; Ackermann et al., 2013;
Jogler & Funk, 2016).

Indications for PeVatron activity from the galactic centre were found (HESS
Collaboration et al., 2016), yet only in recent years have experimental
facilities been capable of measuring gamma-rays with energies $>100$ TeV.
Water Cherenkov and particle detector based facilities in particular, such as
HAWC (Abeysekara et al., 2017), LHAASO (Aharonian et al., 2021) and Tibet-
AS$\gamma$ (Amenomori et al., 2009) have contributed significantly to this
advance.

Until 2023, the ultra-high-energy (UHE, $>100$ TeV) gamma-ray sky comprised a
mere $\sim 15$ sources, with the Crab nebula one of the first identified
(Amenomori et al., 2019; LHAASO Collaboration et al., 2021). A further 31 UHE
sources were announced in the first LHAASO catalogue (Cao et al., 2023).
Twelve UHE sources were reported by LHAASO in 2021 (Cao et al., 2021c), the
majority of which are spatially coincident with known Very-High-Energy (VHE,
$\gtrsim 100$ TeV) gamma-ray sources. In particular, several associations with
energetic pulsar wind nebulae from which the emission is understood to be
dominantly leptonic in origin. Despite Klein-Nishina suppression of inverse
Compton scattering at the highest energies, this suppression is relaxed in the
case of high radiation field environments, and a leptonic scenario remains a
viable interpretation for the UHE sources associated with known energetic
pulsars (Vannoni et al., 2009; Breuhaus et al., 2021, 2022).

There is, however, one source reported in Cao et al. (2021c), for which the
gamma-ray emission was first discovered at UHE and without any known
counterpart accelerators, such as pulsars or supernova remnants (SNRs). LHAASO
J2108+5157 is an enigmatic source, spatially coincident with molecular clouds
yet with the accelerator mechanism remaining unidentified (Cao et al., 2021b).
In the wake of the LHAASO discovery, follow-up observations were conducted by
several facilities, including in the radio and X-ray bands as well as by
gamma-ray experiments. The Fermi-LAT source 4FGL J2108.0+5155 is spatially
coincident with the UHE emission, but due to the differing spectral properties
a physical association remains unclear (Abdollahi et al., 2020). A re-analysis
of the Fermi-LAT data found a potential spatial extension of $0.48^{\circ}$
angular size of the source, designated 4FGL J2108.0+5155e (Cao et al., 2021b).

A $3.7\sigma$ signal of gamma-ray emission was measured at $E>3$ TeV by the
Large Sized Telescope (LST-1), a prototype telescope for the forthcoming
Cherenkov Telescope Array (CTA) (Acharya et al., 2013). They derive upper
limits in the energy range 0.32 TeV to 100 TeV that considerably constrain
model scenarios for the origin of the emission (Abe et al., 2023).

The HAWC observatory recently reported a $\sim 7\,\sigma$ detection in $\sim
2400$ days of data (Kumar et al., 2023). However, observations and analysis by
the VERITAS IACT array did not result in a detection, with constraining upper
limits being reported, consistent with those from the LST-1 (Kumar et al.,
2023; Abe et al., 2023).

Although there is little observational evidence for SNRs currently acting as
PeVatrons, it remains feasible that SNRs act as PeVatrons only for a
comparatively short period during their lifetimes, such that the rate of
currently detectable SNR PeVatrons is low (Cristofari et al., 2020). Particle
escape from the shock region occurs in an energy-dependent manner, such that
the most energetic particles will also be the first to leave the shock region
(Celli et al., 2019a). Evidence for PeV particles may therefore be found not
at the location of the accelerator, but rather from subsequent interactions of
these particles with target material in the ambient medium, such as nearby
molecular clouds (Gabici & Aharonian, 2007; Morlino et al., 2009; Inoue et
al., 2012; Celli et al., 2019b). This scenario has been proposed as a possible
explanation for the UHE emission from LHAASO J2108+5157 (Cao et al., 2021b;
Kar & Gupta, 2022; De Sarkar, 2023; Abe et al., 2023).

In contrast to previous models for LHAASO J2108+5157, in this study we scan
the parameter space in two free variables, namely SNR age and the distance
between the cloud and the SNR, to determine the range of plausible values for
the required properties of the responsible SNR. We investigate the influence
of uncertainties in the cloud properties on the resulting gamma-ray flux for
the best-matched models. The corresponding expected neutrino flux is
estimated, and the plausibility of the best-matched models is discussed.

##  2 Method

We adopt the model of Mitchell et al. (2021, 2023), based on Aharonian &
Atoyan (1996) and Kelner et al. (2006), to investigate the scenario of a SNR
illuminating molecular clouds as a possible explanation for LHAASO J2108+5157.
Whilst there are several clouds identified in the vicinity, we focus on clouds
that are spatially coincident with the $\gamma$-ray emission and located
closest to the best-fit centroid of LHAASO J2108+5157 at
($l=92.2148^{\circ}$,$b=2.9359^{\circ}$). Cloud 4607 from the Miville-
Deschênes et al. (2017) catalogue based on data from the 12CO survey of Dame
et al. (2001) has been considered in previous models of the region (Kar &
Gupta, 2022; De Sarkar, 2023), whilst recently a newly identified cloud has
been detected in the region (de la Fuente et al., 2023). Adopting the
convention of prior works, we henceforth refer to these two clouds as
MML[2017]4607 and FKT[2022] respectively. Table 1 summarises the key physical
properties of the clouds relevant for this study.

Table 1: Properties of the molecular clouds considered in this study. $d$ is the distance to the cloud and $n$ is the total number density of hydrogen gas.  Cloud | MML[2017]4607 | FKT[2022]
---|---|---
($l$ (deg), $b$ (deg)) | (92.272, 2.775) | (92.4, 3.2)
$d$ (kpc) | 3.28 | 1.7 $\pm$ 0.6
$n$ (cm-3) | 30 | 37 $\pm$ 14
diameter (deg) | 0.5 | 1.1 $\pm$ 0.2

For convenience, we summarise here the key features of the model from Mitchell
et al. (2021, 2023) adopted for this work.

  * –

Protons are accelerated impulsively with a power-law spectrum of slope
$\alpha$.

  * –

The particle probability density function $f(E,r^{\prime},t^{\prime})$ is
taken from equation (3) of Aharonian & Atoyan (1996), and is a function of the
particle energy $E$, distance travelled from the SNR $r^{\prime}$ and time
since escape from the SNR $t^{\prime}$.

  * –

The SNR radius, $R_{\rm SNR}$, expands with time (t) adiabatically during the
Sedov-Taylor phase as $R_{\rm SNR}\propto t^{2/5}$ (Truelove & McKee, 1999).

  * –

Particles escape from the SNR at a time $t_{\rm esc}$ in a momentum-dependent
manner, following $t_{\rm esc}\propto(p/p_{M})^{1/\beta}$ where $p_{M}$ is the
maximum particle energy reached, assumed to be 3 PeV/c at the Sedov time,
$t_{\rm sed}$ (Celli et al., 2019a).

  * –

Particles are either transported diffusively through the ISM to reach the
cloud or are injected directly into the cloud if the SNR is sufficiently
expanded.

  * –

Diffusion within the intervening ISM is assumed to be slow with respect to the
Galactic average value due to the local accelerator activity (Gabici et al.,
2007).

  * –

Within the cloud, diffusion is suppressed with respect to the ISM by a factor
$\chi$ that relates to local turbulence.

For our default scenario, the Sedov time is assumed to commence at $1.6$ kyr
corresponding to the case of type II supernovae. In the case of type IA
supernovae, the Sedov time commences at a mere 234 yr, which is considered as
an alternative scenario (Celli et al., 2019a).

Details of the SNR forward shock interaction with the cloud, in the case that
the SNR is in close proximity and sufficiently evolved, are neglected. The
diffusion coefficient $D(E)$ is considered to have a power-law dependence on
the energy $E$ as:

$D(E)=\chi D_{0}\left(\frac{E/\mathrm{GeV}}{B(n)/3\mathrm{\mu
G}}\right)^{\delta}\leavevmode\nobreak\ ,$ (1)

where $\delta$ and $\chi$ relate to the properties of the magnetic field
turbulence in the region, and $B(n)$ describes the dependence of the magnetic
field strength $B$ on the cloud density $n$ (see Mitchell et al. (2021)).
Values adopted for $\delta$, $\chi$ and the diffusion coefficient
normalisation $D_{0}$ at 1 GeV are listed in table 2. Within the ISM, $\chi$
is taken to be 1, whilst a value of 0.1 is adopted to account for suppressed
diffusion within the clouds. From the above ingredients, the particle spectrum
as a function of energy $E$, age $t$ and distance from the accelerator $r$ is
obtained, $f(E,r^{\prime},t^{\prime})$.

Table 2: Assumed parameters of the model. The value adopted for $D_{0}$ is taken from Gabici et al. (2007).  Description | Parameter | Value
---|---|---
Proton power law spectral index | $\alpha$ | 2.0
Index characterising energy | $\delta$ | 0.5
-dependence of diffusion |  |
Index characterising momentum | $\beta$ | 2.5
-dependence of particle escape |  |
Diffusion suppression factor | $\chi$ | 0.1 or 1
due to turbulence |  |
Diffusion coefficient | $D_{0}$ | $3\times 10^{26}\frac{\rm cm^{2}}{\rm s}$
normalisation at 1 GeV |  |
ISM particle density | $n_{\rm ISM}$ | 1 cm-3
Maximum particle momentum | $p_{M}$ | 3 PeV/c
Sedov time (type II SNR) | $t_{\rm sed}$ | 1.6 kyr
Sedov time (type Ia SNR) | $t_{\rm sed}$ | 234 yr

Experimental measurements are, however, bound to neutral messengers such as
$\gamma$-rays and neutrinos as the signatures for the presence of energetic
hadronic particles. For comparison to data, the proton spectrum can then be
converted into a gamma-ray emissivity
$\Phi_{\gamma}(E_{\gamma},r^{\prime},t^{\prime})$ (in ph cm-3 s-1 TeV-1) by
using the expressions from Kelner et al. (2006):

$\Phi_{\gamma}(E_{\gamma},r^{\prime},t^{\prime})=cn\int_{E_{\gamma}}^{\infty}\sigma_{\mathrm{inel}}(E)f(E,r^{\prime},t^{\prime})F_{\gamma}\left(\frac{E_{\gamma}}{E},E\right)\frac{dE}{E}\leavevmode\nobreak\
,$ (2)

for which we adopt the parameterisation of the inelastic cross-section for
proton-proton interactions $\sigma_{\mathrm{inel}}(E)$ from Kafexhiu et al.
(2014), noting that due to high uncertainties below $\sim\,100$ GeV, we take
this as an energy threshold and restrict our model predictions to energies
$>100$ GeV only.

Lastly, we obtain the $\gamma$-ray flux $F(E_{\gamma},t)$ at a distance $d$
away from the cloud (i.e. at Earth) taking into account the volume of the
molecular cloud traversed by particles $V_{c}$ via:

$F(E_{\gamma},t)=\Phi_{\gamma}(E_{\gamma},t)V_{c}/(4\pi
d^{2})\leavevmode\nobreak\ .$ (3)

The diffusive galactic CR flux permeates the entire Galaxy, and as such will
also contribute to the total particle flux interacting with the molecular
clouds. To take this contribution into account, we include the proton flux as
measured by the Alpha Magnetic Spectrometer on the International Space Station
(Aguilar et al., 2015). This flux is added to the particle flux arriving at
the cloud, $f$ in equation (2), enabling the relative contributions of
accelerator and the diffuse CR sea to be evaluated.

Kelner et al. (2006) also provide expressions for the neutrino production via
charged pion and muon decay via $F_{\nu}\left(\frac{E_{\nu}}{E},E\right)$ for
the total production of electron and muon neutrinos from the same proton
interactions. By analogy with equations (2) and (3) the corresponding total
neutrino flux can be obtained.

In the next section, we use this model to generate predictions for the gamma
ray flux arising from a hypothetical SNR illuminating the molecular clouds
identified in the vicinity of the LHAASO J2108+5157. Additionally, we consider
the contribution from the galactic CR sea, to establish whether it is
sufficient to account for the observed gamma-ray emission without requiring a
nearby accelerator. The model is compared to measurements from LHAASO and
HAWC, and upper limits from the LST-1 and VERITAS (Cao et al., 2021b; Abe et
al., 2023; Kumar et al., 2023).

Figure 1: Dependence of gamma-ray flux on properties of a type II supernova.
Left: Variation in supernova age for a fixed distance of 24 pc. Right:
Variation with distance between the SNR and the cloud for a fixed age of 4
kyr. Solid lines correspond to MML[2017]4607 and dashed lines correspond to
FKT[2022].

Figure 2: Dependence of gamma-ray flux on properties of a type Ia supernova.
Left: Variation in supernova age for a fixed distance of 10 pc. Right:
Variation with distance between the SNR and the cloud for a fixed age of 1
kyr. Solid lines correspond to MML[2017]4607 and dashed lines correspond to
FKT[2022].

##  3 Results

###  3.1 Scan over SNR parameter space

As the properties of the molecular clouds are known (table 1), we vary the
properties of a hypothetical SNR to investigate the required values to account
for the $\gamma$-ray flux of LHAASO J2108+5157. We assume that the SNR is
located at the same distance from Earth as the cloud. The SNR age is varied in
ten logarithmically spaced steps between 1 kyr and 500 kyr, for a fixed
separation distance between the SNR and cloud of 24 pc. Similarly, the
separation distance is independently varied in ten logarithmically spaced
steps between 10 pc and 500 pc for a fixed SNR age of 4 kyr. For type Ia
supernovae, the fixed reference values were reduced to 10 pc and 1 kyr. These
curves are shown in figures 1 and 2 for type II and type Ia supernovae
respectively.

In the case of type II supernova remnants shown in figure 1, the predicted
flux is comparable to the data for cloud FKT[2022], yet the flux predicted for
MML[2017]4607 consistently falls below the measured flux. Younger ages are
preferred, with the flux at energies $\lesssim 1$ TeV becoming over predicted
between $\sim$7 kyr and 30 kyr for FKT[2022]. The key features of the model
are that the highest energy particles escape the SNR at earlier times and are
first to arrive at the cloud. The spectral energy distribution hence rises at
the highest energies at earlier times (and for shorter distances). The
particle distribution then cools as a function of age, with the peak shifting
towards lower energies.

In the case of type Ia supernova remnants shown in figure 2, a separation
distance larger than 24 pc is required for FKT[2022] to avoid over predicting
the flux in the $<10$ TeV range. MKT[2022]4607 is better able to account for
the gamma-ray flux under the type Ia scenario, yet only for an optimum
combination of low distance and young age. As $t_{\rm sed}$ is lower for the
type Ia scenario, the spectral energy distribution is more highly populated at
an earlier stage.

###  3.2 Contribution from the galactic CR sea

As described above, the contribution from diffusive galactic CRs is included
in the model, assuming that the particle flux is comparable to that measured
at Earth (Aguilar et al., 2015). From the parameter scan, we find that the
contribution from the nearby SNR dominates over that from galactic CR sea in
most cases. Indeed, the contribution from diffusive galactic CRs only exceeds
that from the SNR if either the cloud-SNR distance is $\gtrsim 200$ pc (for
young $\lesssim 10$ kyr SNRs), or if the SNR is old, $\gtrsim 400$ kyr (for
nearby $\lesssim 50$ pc SNRs).

In order to test whether the diffuse galactic CR sea could be solely
responsible for the measured gamma-ray flux, the normalisation of the galactic
flux contribution was varied, in the absence of considering any hypothetical
SNR. To match the observed emission at TeV energies using the molecular clouds
considered, the normalisation must be of order $\sim 10^{3}$ higher than that
measured at Earth. This enhancement is unlikely to be achieved without the
presence of an accelerator nearby.

Next, we consider all possible combinations of SNR age and separation distance
within the aforementioned ranges. A chisquare evaluation of the model curve to
the LHAASO data points only is used to establish which model curves provide
the closest match to the data. Due to the large number of free parameters
entering into the model, we do not perform a minimisation, as there will be
multiple local minima in the parameter space able to account for the data.
Rather, we aim provide a plausible range of allowed values for the specific
case of this model, with assumed fixed parameters as in table 2.

Figure 3: Model curves corresponding to the parameter combinations that best
match the data as listed in table 3.3.1. Left: type II and Right: type Ia
supernova remnant. Solid lines and shaded uncertainty band correspond to
MML[2017]4607. Dashed lines and hatched region correspond to FKT[2022].

###  3.3 Best-matched models

####  3.3.1 Clouds MML[2017]4607 and FKT[2022]

For each cloud, model curves corresponding to the two best matching
combinations of SNR age and separation distance are shown in figure 3. Model
curves for MML[2017]4607 were consistently below the data points for
$\chi=0.1$ within the cloud, as seen in figures 1 and 2. To obtain parameter
values within comparable agreement to the data as for FKT[2022], we neglected
the suppressed diffusion within the cloud for MML[2017]4607 (and for this
section only) by setting $\chi=1$. This corresponds to the most optimistic
case in which CRs can freely penetrate the cloud, although we note that
$\delta$ was kept fixed to $0.5$ and we did not investigate the effect of
altering the energy dependence of the diffusion coefficient in equation (1).

The best matching combinations are summarised in table 3.3.1. In general, the
SNR age was found to have a stronger influence on the curve shape and hence
quality of the match to LHAASO data than the separation distance. FKT[2022]
yielded more parameter combinations with a lower $\chi^{2}$ than
MML[2017]4607, where model curves for the same age yet for smaller distances
were essentially consistent. This is supported by figures 1 and 2 \- for a
fixed age, provided the distance is sufficiently small that CRs have had time
to traverse the cloud, the gamma-ray flux remains constant with decreasing
distance. (Equivalently, the gamma-ray flux drops with increasing distance.)
Overall, the type Ia scenario (i.e. a lower $t_{\rm sed}$) is preferred.

One might ask whether or not a finer-resolution of values covering the
reasonable parameter space would lead to a model that better matched the data.
Whilst this may be the case, we first consider the effect of propagating the
uncertainties in the measured properties of the clouds (table 1) through the
model. An upper bound to the flux is obtained by adopting the 1 $\sigma$
deviation $d-\sigma_{d}$ and $n+\sigma_{n}$, whilst a lower bound is similarly
obtained from the model evaluated with $d+\sigma_{d}$ and $n-\sigma_{n}$,
where we intrinsically assume that the uncertainties are Gaussian distributed.
Increasing $n$ will increase the target material and hence flux as per
equation (2), whilst increasing $d$ will decrease the flux as per equation
(3). For FKT[2022] uncertainties are reported in de la Fuente et al. (2023),
whilst for MML[2017]4607 uncertainties are not provided in the case that near
and far estimates agree (Miville-Deschênes et al., 2017). We therefore adopt a
20% uncertainty in $d$ and $n$ for MML[2017]4607 as a rough estimate, given
that the true uncertainty and subsequent variation in the model is unknown.
Resulting uncertainty bands corresponding to the parameter combinations
reported in table 3.3.1 are shown in figure 3.

Combinations of SNR age, $t$ and separation distance, $\Delta d$ for the model
curves that best match the LHAASO data, listed in ranked order. These curves
are shown in figure 3.  Cloud $t$ (kyr) $\Delta d$ (pc) SN type $\chi^{2}$
MML[2017]4607 1 37 Ia 5.1 FKT[2022] 4 37 * Ia 6.7 FKT[2022] 4 57 Ia 9.2
FKT[2022] 4 57 II 15.5 FKT[2022] 8 24 ** II 17.0 MML[2017]4607 4 24 ** II 24.4
MML[2017]4607 2 37 II 25.0 MML[2017]4607 1 24 Ia 28.2

  * *

Model curves for the same SNR age yet with smaller distances provided a
comparable fit to the LHAASO data, but severely overestimated the LST-1 upper
limits and are hence not shown.

  * **

Model curves for the same SNR age yet with a distance of 10 pc and 15 pc were
comparable to the 24 pc distance quoted.

Figure 3 clearly illustrates that the uncertainty introduced to the model from
experimental measurements (or the adopted 20% uncertainty) on the cloud
properties leads to variation in predicted flux comparable to that seen by
varying the input age and distance of the parameter scan. Therefore, a more
finely-grained exploration of the SNR parameter space is not well-motivated.

####  3.3.2 Corresponding neutrino flux

For two of the best matching models from table 3.3.1, we show the
corresponding total neutrino flux in figure 4. For MML[2017]4607 this is for
$t=1$ kyr and $\Delta d=37$ pc, whilst for FKT[2022] we show $t=4$ kyr and
$\Delta d=57$ pc, both for the SN Ia case. Although $\Delta d=37$ pc yielded a
lower $\chi^{2}$ for FKT[2022] with respect to the LHAASO data, this curve is
disfavoured as it exceeds the upper limits provided by LST-1 (upper dashed
curve in figure 3). These curves essentially scale with the $\gamma$-ray flux,
yet still lie at least an order of magnitude in flux below the sensitivity of
current neutrino experiments suited for the detection of astrophysical
neutrinos, such as IceCube (Aartsen et al., 2019).

Figure 4:  Gamma-ray and neutrino fluxes for two of the best-matching
scenarios for a type Ia SN from table 3.3.1. Solid lines correspond to the
expected gamma-ray flux and dashed lines to the neutrino flux.

##  4 Discussion

LHAASO J2108+5157 is an intriguing UHE gamma-ray source with no known
counterparts yet spatially coincident with molecular clouds. In this study, we
investigate a scenario whereby the molecular cloud is illuminated by energetic
protons accelerated at a SNR in the vicinity. By scanning the parameter space
of SNR age and separation distance between the hypothetical SNR and the cloud,
we obtain model predictions that can be compared to data, thereby constraining
the most likely SNR properties. Consistently, we find that a comparatively
young ($<10$ kyr) and nearby ($d\lesssim 40-60$ pc) SNR is required.

There are currently no known SNRs matching this description. From the SNR
catalogue (Ferrand & Safi-Harb, 2012), the two closest SNRs are G094.0+01.0
and G093.7-00.2, at angular distances of more than $2.4^{\circ}$ from LHAASO
J2108+5157. At the 3.28 kpc distance of MML[2017]4607 this corresponds to 140
pc and 190 pc separation from the cloud respectively, whilst at the 1.7 kpc
distance of FKT[2022] the SNRs are situated 80 pc and 110 pc away from the
cloud. Additionally, G094.0+1.0 has an estimated age of 25 kyr, far older than
the SNR ages preferred by our model. We conclude that neither SNR is
associated to LHAASO J2108+5157.

Nevertheless, it remains plausible that there are further, as yet undiscovered
SNRs located in the region. Recent results from the EMU/POSSUM survey,
performed using the Australian Square Kilometer Array Pathfinder (ASKAP)
observed a region of the galactic plane containing 7 known SNRs, and found 21
candidates, of which 13 were new discoveries Ball et al. (2023). This supports
the notion that radio surveys to date may not be sufficiently sensitive to
detect all SNRs within a given region.

Several molecular clouds have been identified in the region, two from Miville-
Deschênes et al. (2017) based on Dame et al. (2001) (MML[2017]4607 and
MML[2017]2870) and most recently a new cloud FKT[2022] reported by de la
Fuente et al. (2023). Model parameters were explored for the two clouds
spatially coincident with LHAASO J2108+5157, namely MML[2017]4607 and
FKT[2022].

Both type II and type Ia supernova explosion scenarios were considered, where
the main difference is in the assumed time for transition to the Sedov-Taylor
phase ($t_{\rm sed}$). Although a better match could be achieved under the
type Ia scenario, we consider this unlikely. Type Ia supernoave occur in older
systems where at least one member of a binary system has sufficiently evolved
to become a white dwarf, generally corresponding to environments not rich in
molecular material. Type II supernovae, however, occur in younger environments
where an abundance of molecular material can be expected, similar to that
observed in the vicinity of LHAASO J2108+5157. Hence, we rather interpret
these results as indicating that an earlier transition into the Sedov-Taylor
phase is preferred, which may reflect (e.g.) properties of the ambient medium
rather than the nature of the progenitor (Truelove & McKee, 1999).

In all model curves, the highest energy data point at $\sim$500 TeV could not
be well matched with a maximum energy of the proton spectrum of 1 PeV.
Therefore, throughout this study we assumed a maximum energy at the Sedov time
of 3 PeV.

For MML[2017]4607 to account for the data, we neglected an additional
suppression within the cloud due to turbulence compared to the ISM (i.e.
$\chi=1$). With $\chi=0.1$ within the cloud, MML[2017]4607 consistently under
predicted the data in our model (figures 1 and 2). Our model assumed locally
suppressed diffusion compared to the Galactic average also in the intervening
medium between the SNR and the cloud, a reasonable assumption for regions of
active particle acceleration (Gabici et al., 2007; D’Angelo et al., 2018;
Inoue, 2019).

Suppressed diffusion and young SNR age as preferred model parameters is in
agreement with the 4.5 kyr age obtained by Kar & Gupta (2022), although De
Sarkar (2023) suggests an older SNR age of 44 kyr, obtained with a different
spectral index for the particle population. A young SNR may still be a
comparatively weak producer of synchrotron emission, or could be of small size
and remain embedded within (or obscured by) molecular clouds in the region.
Given the angular size of the molecular clouds, a young SNR could be
completely hidden behind the clouds along the line of sight. Using the
relation $R_{\rm SNR}\propto t^{2/5}$ for evolution in the Sedov-Taylor phase,
an SNR younger than 12 kyr for MML[2017]2870 and 19 kyr for FKT[2022] would be
small enough to be obscured by the cloud. This is consistent with the
preferred $<10$ kyr SNR age.

Nonetheless, other scenarios for the origin of LHAASO J2108+5157 remain
plausible. Young stellar clusters have been hypothesised as suitable Galactic
PeVatrons, with particle acceleration occurring at the termination shock of
the collective wind (Aharonian et al., 2019; Morlino et al., 2021). There are
two known young stellar clusters nearby to LHAASO J2108+5157: although the
distance to Kronberger 80 is known to be at least 4.8 kpc or larger
(Kharchenko et al., 2016; Cantat-Gaudin & Anders, 2020), disfavouring an
association with the molecular clouds in the region, and the distance to
Kronberger 82 remains unknown (Kronberger et al., 2006). As such, a stellar
cluster is a potential alternative accelerator, also capable of illuminating
molecular clouds with CRs, but not well motivated in this region.

Given the spatial correlation of LHAASO J2108+5157 with molecular clouds a
leptonic scenario for the emission seems unlikely, nevertheless it has been
demonstrated that powerful pulsar wind nebulae are capable of accelerating
leptons to beyond 1 PeV and can account for UHE gamma rays, especially in high
radiation field environments (Vannoni et al., 2009; Breuhaus et al., 2022).
The lack of a pulsar counterpart, or of X-ray synchrotron emission that would
indicate the presence of a pulsar wind nebula also in cases where the pulsed
emission is mis-aligned, disfavours such a scenario.

With the advent of current generation detectors such as LHAASO sensitive to
UHE gamma rays, we may expect other enigmatic sources to emerge, corresponding
to clouds illuminated by unknown accelerators. Other unidentified gamma-ray
sources for which no known counterpart has been identified to date, such as
LHAASO J0341+5258, may have a similar origin (Cao et al., 2021a). The first
LHAASO catalogue reported no fewer than seven further new sources that seem to
be “dark” in nature, without any known counterparts (Cao et al., 2023).
Undoubtedly, further follow-up studies, both in terms of observation and
interpretation, are necessary to determine the origin of these enigmatic
gamma-ray sources.

##  5 Conclusion

LHAASO J2108+5157 is a dark UHE gamma-ray source spatially coincident with two
molecular clouds. We find that the gamma-ray emission can be accounted in
terms of molecular cloud illumination by CRs from a nearby ($\lesssim 40-60$
pc) young ($<10$ kyr) SNR. Although no SNR is currently known matching these
criteria, such an SNR could be obscured by other material along the line of
sight, or simply lie below the detection threshold of previous surveys (Ball
et al., 2023). Interactions of the diffuse galactic CR sea with the molecular
clouds is found to be insufficient to explain the observed gamma-ray flux.

As the exposure of current survey instruments increases, and with the advent
of future facilities such as CTA and SWGO, we can anticipate further such
discoveries, potentially unveiling a population of UHE sources tracing the
presence of PeV particles (Acharya et al., 2013; Albert et al., 2019;
Casanova, 2022; Cao et al., 2023). The key to identifying PeVatrons may lie
not in emission from the accelerators themselves, but rather from evidence of
energetic particles that have escaped the source region.

###### Acknowledgements.

The author is grateful to G. Rowell & C. van Eldik for fruitful discussions
and especially to A. Specovius for reading the manuscript. This work is
supported by the _Deutsche Forschungsgemeinschaft, DFG_ project number
452934793.

## References

  * Aartsen et al. (2019) Aartsen, M. G., Ackermann, M., Adams, J., et al. 2019, ApJ, 886, 12
  * Abdollahi et al. (2020) Abdollahi, S., Acero, F., Ackermann, M., et al. 2020, ApJS, 247, 33
  * Abe et al. (2023) Abe, S., Aguasca-Cabot, A., Agudo, I., et al. 2023, A&A, 673, A75
  * Abeysekara et al. (2017) Abeysekara, A. U., Albert, A., Alfaro, R., et al. 2017, ApJ, 843, 39
  * Acharya et al. (2013) Acharya, B. S., Actis, M., Aghajani, T., et al. 2013, Astroparticle Physics, 43, 3
  * Ackermann et al. (2013) Ackermann, M., Ajello, M., Allafort, A., et al. 2013, Science, 339, 807
  * Aguilar et al. (2015) Aguilar, M., Aisa, D., Alpat, B., et al. 2015, Phys. Rev. Lett., 114, 171103
  * Aharonian et al. (2021) Aharonian, F., An, Q., Axikegu, et al. 2021, Chinese Physics C, 45, 025002
  * Aharonian et al. (2019) Aharonian, F., Yang, R., & de Oña Wilhelmi, E. 2019, Nature Astronomy, 3, 561
  * Aharonian & Atoyan (1996) Aharonian, F. A. & Atoyan, A. M. 1996, A&A, 309, 917
  * Albert et al. (2019) Albert, A., Alfaro, R., Ashkar, H., et al. 2019, arXiv e-prints, arXiv:1902.08429
  * Amenomori et al. (2019) Amenomori, M., Bao, Y. W., Bi, X. J., et al. 2019, Phys. Rev. Lett., 123, 051101
  * Amenomori et al. (2009) Amenomori, M., Bi, X. J., Chen, D., et al. 2009, ApJ, 692, 61
  * Apel et al. (2013) Apel, W. D., Arteaga-Velàzquez, J. C., Bekk, K., et al. 2013, Phys. Rev. D, 87, 081101
  * Ball et al. (2023) Ball, B. D., Kothes, R., Rosolowsky, E., et al. 2023, MNRAS, 524, 1396
  * Bell (1978) Bell, A. R. 1978, MNRAS, 182, 147
  * Bell (2013) Bell, A. R. 2013, Astroparticle Physics, 43, 56
  * Bell et al. (2013) Bell, A. R., Schure, K. M., Reville, B., & Giacinti, G. 2013, MNRAS, 431, 415
  * Blumenthal & Gould (1970) Blumenthal, G. R. & Gould, R. J. 1970, Reviews of Modern Physics, 42, 237
  * Breuhaus et al. (2021) Breuhaus, M., Hahn, J., Romoli, C., et al. 2021, ApJ, 908, L49
  * Breuhaus et al. (2022) Breuhaus, M., Reville, B., & Hinton, J. A. 2022, A&A, 660, A8
  * Cantat-Gaudin & Anders (2020) Cantat-Gaudin, T. & Anders, F. 2020, A&A, 633, A99
  * Cao et al. (2021a) Cao, Z., Aharonian, F., An, Q., et al. 2021a, ApJ, 917, L4
  * Cao et al. (2021b) Cao, Z., Aharonian, F., An, Q., et al. 2021b, ApJ, 919, L22
  * Cao et al. (2023) Cao, Z., Aharonian, F., An, Q., et al. 2023, arXiv e-prints, arXiv:2305.17030
  * Cao et al. (2021c) Cao, Z., Aharonian, F. A., An, Q., et al. 2021c, Nature, 594, 33
  * Casanova (2022) Casanova, S. 2022, Universe, 8, 505
  * Celli et al. (2019a) Celli, S., Morlino, G., Gabici, S., & Aharonian, F. A. 2019a, MNRAS, 490, 4317
  * Celli et al. (2019b) Celli, S., Morlino, G., Gabici, S., & Aharonian, F. A. 2019b, MNRAS, 487, 3199
  * Cristofari et al. (2020) Cristofari, P., Blasi, P., & Amato, E. 2020, Astroparticle Physics, 123, 102492
  * Dame et al. (2001) Dame, T. M., Hartmann, D., & Thaddeus, P. 2001, ApJ, 547, 792
  * D’Angelo et al. (2018) D’Angelo, M., Morlino, G., Amato, E., & Blasi, P. 2018, MNRAS, 474, 1944
  * de la Fuente et al. (2023) de la Fuente, E., Toledano-Juarez, I., Kawata, K., et al. 2023, PASJ, 75, 546
  * De Sarkar (2023) De Sarkar, A. 2023, MNRAS, 521, L5
  * Ferrand & Safi-Harb (2012) Ferrand, G. & Safi-Harb, S. 2012, Advances in Space Research, 49, 1313
  * Gabici & Aharonian (2007) Gabici, S. & Aharonian, F. A. 2007, ApJ, 665, L131
  * Gabici et al. (2007) Gabici, S., Aharonian, F. A., & Blasi, P. 2007, Ap&SS, 309, 365
  * Ginzburg & Syrovatskii (1964) Ginzburg, V. L. & Syrovatskii, S. I. 1964, The Origin of Cosmic Rays
  * H.E.S.S. Collaboration et al. (2018) H.E.S.S. Collaboration, Abdalla, H., Abramowski, A., et al. 2018, A&A, 612, A3
  * HESS Collaboration et al. (2016) HESS Collaboration, Abramowski, A., Aharonian, F., et al. 2016, Nature, 531, 476
  * Hillas (2006) Hillas, A. M. 2006, arXiv astro-ph, 0607109
  * Inoue (2019) Inoue, T. 2019, ApJ, 872, 46
  * Inoue et al. (2012) Inoue, T., Yamazaki, R., Inutsuka, S.-i., & Fukui, Y. 2012, ApJ, 744, 71
  * Jogler & Funk (2016) Jogler, T. & Funk, S. 2016, ApJ, 816, 100
  * Kafexhiu et al. (2014) Kafexhiu, E., Aharonian, F., Taylor, A. M., & Vila, G. S. 2014, Phys. Rev. D, 90, 123014
  * Kar & Gupta (2022) Kar, A. & Gupta, N. 2022, ApJ, 926, 110
  * Kelner et al. (2006) Kelner, S. R., Aharonian, F. A., & Bugayov, V. V. 2006, Phys. Rev. D, 74, 034018
  * Kharchenko et al. (2016) Kharchenko, N. V., Piskunov, A. E., Schilbach, E., Röser, S., & Scholz, R. D. 2016, A&A, 585, A101
  * Kronberger et al. (2006) Kronberger, M., Teutsch, P., Alessi, B., et al. 2006, A&A, 447, 921
  * Kumar et al. (2023) Kumar, S., Martin, M., & Wang, X. 2023, in 38th International Cosmic Ray Conference (ICRC2023), 914
  * Lagage & Cesarsky (1983) Lagage, P. O. & Cesarsky, C. J. 1983, A&A, 125, 249
  * LHAASO Collaboration et al. (2021) LHAASO Collaboration, Cao, Z., Aharonian, F., et al. 2021, Science, 373, 425
  * Mitchell et al. (2021) Mitchell, A. M. W., Rowell, G. P., Celli, S., & Einecke, S. 2021, MNRAS, 503, 3522
  * Mitchell et al. (2023) Mitchell, A. M. W., Rowell, G. P., Celli, S., & Einecke, S. 2023, MNRAS, 520, 300
  * Miville-Deschênes et al. (2017) Miville-Deschênes, M.-A., Murray, N., & Lee, E. J. 2017, ApJ, 834, 57
  * Morlino et al. (2009) Morlino, G., Amato, E., & Blasi, P. 2009, MNRAS, 392, 240
  * Morlino et al. (2021) Morlino, G., Blasi, P., Peretti, E., & Cristofari, P. 2021, MNRAS, 504, 6096
  * Parizot (2014) Parizot, E. 2014, Nuclear Physics B Proceedings Supplements, 256, 197
  * Truelove & McKee (1999) Truelove, J. K. & McKee, C. F. 1999, ApJS, 120, 299
  * Vannoni et al. (2009) Vannoni, G., Gabici, S., & Aharonian, F. A. 2009, A&A, 497, 17