Using Interstellar Clouds to Search for Galactic PeVatrons: Gamma-ray Signatures from Supernova Remnants

Interstellar clouds can act as target material for hadronic cosmic rays; gamma rays subsequently produced through inelastic proton-proton collisions and spatially associated with such clouds can provide a key indicator of efficient particle acceleration. However, even in the case that particle acceleration proceeds up to PeV energies, the system of accelerator and nearby target material must fulfil a specific set of conditions in order to produce a detectable gamma-ray flux. In this study, we rigorously characterise the necessary properties of both cloud and accelerator. By using available Supernova Remnant (SNR) and interstellar cloud catalogues, we produce a ranked shortlist of the most promising target systems, those for which a detectable gamma-ray flux is predicted, in the case that particles are accelerated to PeV energies in a nearby SNR. We discuss detection prospects for future facilities including CTA, LHAASO and SWGO; and compare our predictions with known gamma-ray sources. The four interstellar clouds with the brightest predicted fluxes>100 TeV identified by this model are located at (l,b) = (333.46,-0.31), (16.97,0.53), (110.43,1.89) and (336.73,-0.98). These clouds are consistently bright under a range of model scenarios, including variation in the diffusion coefficient and particle spectrum. On average, a detectable gamma-ray flux is more likely for more massive clouds; systems with lower separation distance between the SNR and cloud; and for slightly older SNRs.


INTRODUCTION
Understanding the origin of Galactic cosmic rays (CRs) continues to be a very active research area, with many open questions (Blasi 2013;Grenier et al. 2015). The search for evidence of Galactic PeVatrons -astrophysical accelerators capable of achieving the CR knee at ∼ 1 PeV (10 15 eV) -is a major focus of several present day astrophysical experiments. SNRs have long been thought to provide the bulk of the CR flux at these energies, however the gamma-ray spectra observed from SNR to date typically cut-off at energies of a few tens of TeV. Cosmic rays at PeV energies, however, would produce ∼100 TeV gamma rays. These CRs are sufficiently energetic to escape the SNR shock and travel through the interstellar medium (ISM).
Interstellar clouds provide suitable targets for CR interactions generating pions (among other particles), that may subsequently produce a detectable gamma-ray flux through the decay of neutral pions. Similarly, a detectable neutrino flux may be produced through the decay of charged pions. Both gamma rays and neutrinos are neutral messengers, providing a direct indication of their origin. The arrival direction of charged cosmic rays typically does not indicate the direction of origin, due to deflection by interstellar magnetic fields. In ★ E-mail: alison.mw.mitchell@fau.de this paper, we focus on the gamma-ray signatures left by SNRs on nearby clouds, a scenario that was first investigated in the search for PeVatrons by .
The amount of gamma-ray flux produced grows both with the mass of the cloud, corresponding to the amount of target material available for interactions, and with the CR over-density. A local over-density of CRs may be present due to CR production from a source in the vicinity of the cloud, such as an SNR; or due to a strong suppression of the diffusion coefficient within a specific region. Gamma-ray emission coincident with interstellar clouds is typically invoked as evidence of hadronic particle acceleration. Lack of detectable emission may either disfavour hadronic scenarios, or may be due to gamma-ray absorption in surrounding radiation fields, leading to opaque sources (see e.g. Celli et al. (2017)). However, it is not necessarily the case that a gamma-ray flux detectable by current instruments is produced, even if PeV CRs are present and the source is transparent to high energy radiation, as the flux is highly dependent on the specific properties of the interacting systems.
Previous modelling of SNRs and clouds of target material have enabled hadronic emission to be confirmed through the characteristic pion-bump signature in several SNRs, including IC 443, W 44, W 28 and W 51C (Ackermann et al. 2013;Jogler & Funk 2016;Cui et al. 2018). For this study, we adopt a baseline model adapted from Aharonian & Atoyan (1996) and Kelner et al. (2006), to describe the evolution of the particle flux during transport through the ISM and the production of gamma rays upon reaching the cloud respectively. The model is described in detail in section 2.
Using this model, we quantify which combinations of accelerator and cloud properties achieve detectable levels of gamma-ray flux under standard assumptions. We present this information through plots demonstrating the favourable phase space for interstellar clouds, in an analogous manner to the well-known phase space for astrophysical sites of CR acceleration (Hillas 1984). In particular, this can be indicated with curves of constant gamma-ray flux, such as the 50 hour sensitivity of the future Cherenkov Telescope Array (CTA) at a gamma-ray energy of 100 TeV (Cherenkov Telescope Array Consortium et al. 2019).
CTA will outperform the current generation of Imaging Atmospheric Cherenkov Telescope (IACT) arrays, such as H.E.S.S., MAGIC and VERITAS, towards 100 TeV; yet data from these experiments can already be used to place constraints on the flux above 10 TeV Aleksić et al. 2015;Meagher & VER-ITAS Collaboration 2015). Water Cherenkov Detector (WCD) based facilities, such as HAWC, provide a complementary view of the very high energy gamma-ray sky -with improved sensitivity around 100 TeV at the cost of a degraded angular resolution (Abeysekara et al. 2017b). Operating as wide field-of-view survey instruments, WCD based experiments cover large areas of the sky reducing the bias in source detection. The future facilities LHAASO (already producing first results in a partial configuration) and SWGO (planned for the Southern hemisphere) may detect multiple such interstellar clouds illuminated by nearby SNRs (Bai et al. 2019;Albert et al. 2019).
Using current catalogues of known SNRs and interstellar clouds, we explore which combinations of accelerators and clouds are the most promising targets to look for evidence of PeVatron activity with current instruments, such as H.E.S.S. or HAWC, and for future detectability with CTA, LHAASO or SWGO. By applying our model to different plausible combinations of SNRs and nearby interstellar clouds, we produce a ranked shortlist of the most promising candidate targets.

MODEL
In this study, we model the injection of particles (treated as a pure proton flux) from SNRs. Figure 1 illustrates the scenario under consideration schematically. Three key aspects comprise the model: SNR evolution and particle escape; particle transport through the ISM; and particle interactions with interstellar clouds to produce gamma rays. We adopt a dynamical description for the injection of particles into the ISM, which allows us to treat SNRs as impulsive accelerators. In fact, particles of different energies are released at different times, with high energy particles being able to leave the shock before the low energy ones (Celli et al. 2019b). While details of the escape process are still poorly understood, expectations are such that the turbulence which confines particles inside the SNR decays as the shock expands, hence the most energetic particles with larger mean free path are able to leave the system earlier than lower energy ones. As the SNR expands, this also modifies the distance travelled by the particles to reach the target interstellar cloud. Less energetic particles, released later in the SNR evolution, will have a reduced distance to travel to the target cloud.
The propagation of particles from accelerator to target and subsequent modification to the particle spectrum arriving at the interstellar Schematic showing the geometry assumed in the model. Particles are released from the time-dependent SNR radius SNR ( ), propagate through the ISM and penetrate nearby clouds to produce gamma rays. For a nearby cloud located at a distance 1 , the cosmic rays fully traverse the cloud; whereas for a larger cloud located further away at a distance 2 , the cosmic rays are only able to partially traverse the cloud within the available time. A cell based approach is used to treat this case, as described in section 2.
cloud are treated following Aharonian & Atoyan (1996). The gammaray flux produced due to particle interactions in the interstellar clouds are then calculated from Kelner et al. (2006). This enables us to explore the phase space of different accelerator and interstellar cloud properties and the resulting dependence of the gamma-ray flux at a given energy on these properties, assuming the same initial particle spectrum.

Particle spectrum reaching the target
For an impulsive source with a power law acceleration spectrum of slope , we use the following expression to determine the number of protons of energy at a distance from the centre of the SNR at the age SNR of the SNR: where the function ( , , ) describes the probability density function of the particles. The normalisation 0 over the particle spectrum is obtained from: with CR the total energy budget available to cosmic rays from the supernova explosion, which we take to be CR = 10 50 erg for a canonical supernova explosion with kinetic energy SN = 10 51 erg.
Considering those particles which have escaped the SNR and are diffusing into the ISM, the probability density function is (Aharonian & Atoyan 1996): where ( , , ) has units of cm −3 GeV −1 , and the diffusion radius is given by: being the slope with energy of the diffusion coefficient. The variables and correspond to the time spent and distance travelled in the ISM respectively, after CRs escape from the SNR. We note here that equation (3) from Aharonian & Atoyan (1996) is a Green's function solution which already takes the 3D isotropic spherical expansion into account. However, it does not include the energy losses that the particles undergo while the remnant expands. Such adiabatic losses are not relevant for the highest energy particles (at about 1 PeV) that we are interested in, as they leave the system when the remnant still has a well-contained size, as shown in equation (11). This expression (4) can be used for diffusion both through the ISM and within the denser target material. However, within the ISM p-p collisions are irrelevant due to the low ambient density ( ) and the simplification ≈ 2 √︁ ( ) can be used. The proton lifetime against particle interactions is given by the expression: where is the number density of the ambient gas and ≈ 0.45 is the interaction inelasticity in the energy range of interest, which is approximately energy independent at GeV to TeV energies. The cross-section ( ) for particle interactions follows an energy dependence described in Kafexhiu et al. (2014).
where is the proton kinetic energy and ℎ = 0.2797 GeV the threshold kinetic energy for pion production. For the diffusion coefficient we consider a power-law dependence on energy as: with index , suppression factor (which relates to the level of turbulence within the propagation region) and reference diffusion coefficient 0 at 1 GeV, given in table 1. Note that two values for the normalization of the diffusion coefficient are considered, with the higher value corresponding to fast diffusion, i.e. a small diffusion time over a distance ; ( ) = 2 / ( ). Such a value is consistent with Boron over Carbon measurements, namely with the average time spent by CRs in our Galaxy. This should be considered the reference situation, with = 1 within the ISM. A reduced diffusion coefficient will be eventually considered (when explicitly mentioned), in order to investigate the case for enhanced magnetic turbulence around CR sources, as induced e.g. by the presence of CRs themselves, which can modify the diffusion (D'Angelo et al. 2018; Inoue 2019). The Table 1. Parameters of the model and their values, unless otherwise specified. The 'slow case' value of D 0 is taken from ).

Note
Parameter Value Slope of particle spectrum ( ) 2 Slope in energy of ( ) 0.5 (Inverse) Slope in momentum of esc ( ) 2.5 Suppression factor of D(E) in clouds 0.05 Normalization of ( ) (fast case) 0 (1 GeV) 3 × 10 27 cm 2 s −1 Normalization of ( ) (slow case) 0 (1 GeV) 3 × 10 26 cm 2 s −1 ISM density 1cm −3 Sedov time sed 1.6 kyr magnetic field ( ) has a dependence on the density of the cloud where (Crutcher et al. 2010): We note that the estimation of the dependence of the magnetic field on cloud density relies on the Zeeman splitting effect which is significant only for dense clouds. Although a constant value of is assumed for < 300 cm −3 , in practice the magnetic field is likely to vary with , albeit a parameterisation of this relation is highly uncertain at the present time. Measurements of in lower density clouds must use other approaches, such as Faraday tomography (Van Eck et al. 2017;Crutcher 2012). For the ISM, a constant value of = 3 is assumed (Jansson & Farrar 2012). As the original formulation in Aharonian & Atoyan (1996) refers to particle released from a point (such as the centre of the SNR), the normalisation factor 0 (see Eq. (3)) is used to account for release of particles at a time dependent radius SNR ( ). This was derived from the requirement that ∫ ( , , ) d = ∫ ( , , ) 4 2 d = 1, leading to the solution: A list of the benchmark values used for constants in the model (unless otherwise specified) is given in table 1. Results will be presented for the case of slow diffusion as default, with fast diffusion occasionally shown as indicated for comparison. To find the particle spectrum reaching the target interstellar cloud, we first need to know the time at which particles escaped the SNR, esc ( ), such that the time spent in the ISM is = − esc with the time elapsed since the supernova event occurred (i.e. the remnant age). We consider the following equation to describe the momentumdependence of escape time esc ( ): where sed ∼ 1.6 kyr is the Sedov time for core collapse supernovae of 10 solar mass ejecta expanding into a uniform medium of density 1 proton/cm 3 ; and is the maximum momentum of the particles at the Sedov time, set to 1 PeV/c in order to simulate the PeVatron phase of the SNR (Celli et al. 2019b). Note that for type Ia supernovae a shorter Sedov time is expected sed ∼ 250 yr. We adopt the normalisation for core collapse supernovae as default. The parameter shown in Eq. (10) is one of the main factors of uncertainty in the model: in fact, its value depends on the temporal evolution of the magnetic turbulence, which in turn is affected by the physical mechanisms from which it originates. Throughout this work we fixed = 2.5, consistently with gamma-ray observations of middle-aged SNRs (Celli et al. 2019b). The particle momentum is related to the kinetic energy through the standard relativistic energy-momentum relation. From equation (10) we can find the time available for the particles to travel, namely the difference between the system age and particle escape time esc ( ). The distance to the interstellar cloud is calculated between the centres of both the SNR and the cloud. Meanwhile, the SNR expands to a non-negligible size; so we find the radius of the SNR, SNR at the time : assuming that the SNR is undergoing adiabatic expansion within the Sedov-Taylor phase, which typically lasts for ∼ 40 kyr (Truelove & McKee 1999;Reynolds 2008). The radius at which the particles of energy escape from the SNR is therefore given by: such that the distance travelled within the ISM is = − esc with measured to the centre of the SNR. As this study focuses on the PeV particles, which are released at the start of the Sedov-Taylor phase, we do not take further evolutionary stages into consideration. Consequently, the distance travelled by the particles to reach the target also varies with energy and is found by the total distance less the SNR radius at esc as given by equation (12).
In exploring the properties of cloud and accelerator which are ideal for tracing PeVatrons, we will focus on PeV particles and the corresponding 100 TeV gamma-ray flux. The escape time for PeV particles (corresponding to ∼ 100 TeV gamma rays) is 1.6 kyr for core collapse supernovae, corresponding to the Sedov time as in equation (10) (Celli et al. 2019b). 1

Gamma-ray flux produced: Two-step model
The gamma-ray flux produced is now evaluated in a two-step model. The expression for the particle flux, equation (3), is applied to propagation of the particles through the ISM to reach the target material. Treatment of particle propagation through and interaction with the interstellar cloud takes place in a second step.

Step 1: Interstellar Medium
The time available for particle propagation through the ISM is given by the accelerator age SNR minus the escape time esc for particles of a given energy, as determined by equation (10). Similarly, the distance that particles of a given energy must travel in order to reach the cloud is given by the total distance from the stated central position of the SNR and cloud, less the cloud radius and less the radius of the SNR at the escape time of those particles, esc ( ), equation (12).
For this propagation step, an average density = 1 cm −3 is assumed for the ISM. This first step is evaluated using equation (3) Figure 2. Illustration of the cell based approach used to calculate the total flux from a cloud in cases when particles do not fully traverse the cloud. The particles arrive from the left hand side and penetrate the cloud to a constant depth (with constant density assumed throughout the cloud). The flux shown is integrated along the line of sight depth, with spherical geometry assumed.

with
= − − esc through the ISM and is either the available time SNR − esc or the time needed to travel through the ISM ISM = 2 / ( ), whichever is the smaller.

Step 2: Interstellar Cloud
To calculate the subsequent gamma-ray flux produced by particles reaching the target cloud, we use the treatment of Kelner et al. (2006), with the particle spectrum obtained in the previous section; ( , , ) from equation (3). The produced gamma-ray emissivity at fixed space and time position Φ ( , , ), in ph cm −3 s −1 TeV −1 , is given by: where ( , , ) is the particle density arriving at the cloud, as given by equation (3). Clouds are assumed to have a constant density profile, although a spatially varying profile could be considered in future work. For this second step, we wish to calculate the time for which the particles interact with the target material. If the particles have sufficient time to reach the cloud, i.e. ISM ( ) < SNR − esc ( ), then the remaining time spent in the cloud is Δ = SNR − [ esc ( ) + ISM ( )]. The depth, , to which particles penetrate the cloud within the available time, Δ , is an energy dependent quantity given by equation (4). If the cloud can be fully traversed, then the flux observed from the cloud ( , ), in ph cm −2 TeV −1 s −1 , is: where the gamma-ray emissivity was integrated over the volume of the cloud (assumed to be spherical) and is the distance from the cloud to the Earth. However, in cases where the particles do not fully traverse the   cloud, < 2 , we divide the cloud into a grid of cells and evaluate the total flux from the cloud by summing the contributions from all of the cells reached by the particles. Note that due to the energy dependent diffusion process, particles of different energy will reach different locations within the cloud. The volume in equation (14) is replaced by the sum of the cell volumes, , where is the 2D projected area of the ith cell face and is the depth of the ith cell along the line of sight. Figure 2 demonstrates this for a case where particles penetrate to a fixed depth.

Resulting gamma-ray spectrum
Figures 3 to 5 show how the produced gamma-ray spectrum varies with properties of the accelerator-cloud system; the accelerator age, interstellar cloud density and distance between accelerator and cloud. The cloud radius is fixed to 10 pc and the distance of the system to Earth is fixed to 1 kpc throughout this section. For each variable, ten logarithmically spaced values were tested between 1-1000 cm −3 for the density; 1-500 pc for separation distance and 1-500 kyr for accelerator age.
In figure 3, the gamma-ray flux is higher for more dense target clouds as expected (see Eq. (13)). As the accelerator age and separation distance are fixed, the change in density affects mainly the flux normalisation with an otherwise consistent shape to the spectra. Lower energy particles take longer to reach the target material, such that the proton spectrum reaching the cloud shows a low-energy cut-off, which is seen in gamma rays in the form of a smoothed suppression of the flux. For the slow diffusion case, the flux normalisation is higher and the curve smoother with a less significant drop towards higher energies, as the particles typically spend a longer amount of time within the cloud.
The trend in figure 4 also matches that expected, with higher flux for lower separation distances. Again, a drop towards lower energies can be seen. Larger separation distances can be reached by lower energy particles in the case of fast diffusion. Similarly, the energy at which the flux peaks also moves towards higher energies for larger separation distances, as the higher energy particles will reach the cloud earlier and interact for longer than less energetic particles. The largest separation distances reached are 500 pc or 125 pc for fast or slower diffusion respectively; although with a comparatively low predicted flux ∼ 10 −13 − 10 −16 TeV cm −2 s −1 . The shortest distance shown in figure 4 is 31.5 pc due to the requirement that the particles must first traverse the ISM (i.e. that the cloud does not overlap the SNR). At the escape time for 1 PeV particles, esc = 6 pc from equation (11),whilst the cloud radius is fixed to 10 pc, such that the minimum separation distance must be greater than 16 pc. 2 Figure 5 shows how the age of the accelerator affects the predicted gamma-ray flux; the peak flux shifts towards lower energies with increasing age. Note that in the case of slower diffusion, a higher normalisation is achieved as particles spend a longer time interacting within the cloud. This also emphasises how the shape of gamma-ray spectra can be used to infer accelerator ages, due to the energy losses of the underlying particle spectrum as well as the time needed for particles to travel through the ISM to reach the target cloud.
For fixed and , the gamma-ray flux is governed by the amount of target material available for interactions -i.e. by the density or size of the cloud. Once the particles have had sufficient time to reach and fully traverse the cloud, the gamma-ray flux produced is fixed by properties of the cloud. The finite size of the target cloud therefore imposes a limit on the flux.
These figures 3 to 5 therefore demonstrate that properties of both the SNR and target cloud need to be favourable in order to generate a detectable gamma-ray flux.

Accelerator Target Distribution
In principle, particles released from an accelerator may either travel through the ISM before reaching an interstellar cloud; or could be injected directly into a neighbouring cloud. The latter case may occur where the cloud is in close proximity and the SNR expands such that the forward shock encounters the cloud. Indeed, there are prominent examples where this is thought to have occurred, indicated by high levels of turbulence within a cloud and significant offset of leptonic gamma-ray emission away from the cloud direction, such as in the case of the pulsar wind nebula HESS J1825-137 and its progenitor supernova (Voisin et al. 2016). However, we do not consider the case of particle injection into a directly adjacent cloud for two reasons.
Firstly, the result of the impact of a shock on a cloud strongly depends on the system properties, namely shock speed and density contrast between the cloud and the surrounding medium. Diffuse clouds are typically destroyed by the shock, which means that the system would lose the spherical symmetry assumed here. In fact, for this model we consider evolution of the SNR radius with time and corresponding particle release using equations (10) and (11) with spherically symmetric geometry; these approximations are no longer a valid description for the non-symmetric cases. In turn, more dense clumps might be able to survive the shock passage, however in this case the interaction between the shock and the clump might generate additional magnetic turbulence that would modify the transport of particles all around the clump, thus affecting their ability to interact with target gas there (Inoue et al. 2012;Celli et al. 2019a). The extent to which the particle transport is modified by magnetic turbulence depends on the length scale of the fluctuations, and may be limited for high energy particles (Roh et al. 2016).
Secondly, the likelihood of such a situation occurring is low in comparison to the case of an intervening step through the ISM first. The volume of the Milky Way is approximately 10 66 cm 3 , assuming a disk with radius 15 kpc and height 0.3 kpc (Rix & Bovy 2013). To estimate the volume occupied by SNRs in the Sedov-Taylor phase, we first assume an approximately uniform distribution of SNRs throughout our Galaxy, with a uniform distribution of SNR ages. On average, a supernova explosion occurs approximately three times per century. As the Sedov-Taylor phase typically lasts about 40 kyr, this suggests a total of 600 hypothetical SNRs younger than 40 kyr for this calculation, with a radius corresponding to their age as determined by the Sedov-Taylor self-similar expansion given in equation (11). Assuming that each SNR forms a perfect sphere, the total volume occupied by SNRs with an age < 40 kyr is 6.7 × 10 62 cm 3 .
The total volume of interstellar clouds within our Galaxy can be simply obtained from radii reported in the corresponding catalogues, yet assuming a spherical volume. From the Rice et al. (2016) catalogue, as will be used in section 4, the total volume is estimated at 1. volume, whilst clouds occupy 0.25% (or 3.3%, using the more complete catalogue of Miville-Deschênes et al. (2017)). Hence, SNRs occupy a factor ∼ 100 less Galactic volume than the most prominent clouds, implying that direct injection into a neighbouring cloud is unlikely. However, within 40 kyr particles escaping from an SNR could reasonably travel to ∼ 100 pc beyond the SNR radius (corresponding to particle energies of 100 GeV and 10 TeV for the fast and slow diffusion cases respectively). Let us consider instead the reachable volume with a radius of 100 pc, attainable by more energetic particles within a shorter timescale. This increases the volume reachable by SNRs by a factor ∼ 350 and is roughly 3.7% of the total Milky Way volume. This latter estimate is comparable to the volume occupied by clouds.
Considering that the true distribution is not uniform and both clouds and SNRs are more likely to be found in denser regions of the Galaxy, such as along spiral arms, this improves the likelihood of cloud-SNR encounters. Nevertheless, we can conclude that the two step approach with intervening propagation within the ISM is a much more likely scenario than that of direction injection of energetic particles by an SNR into an immediately adjacent cloud.
To further support this conclusion, we note that OH masers are an indicator of SNR shocks interacting directly with interstellar clouds. Of the over 800 maser sites recently catalogued, only ∼20 were associated with SNRs, which even after considering selection effects, may indicate that direct injection is a less common scenario (Beuther et al. 2019).
Lastly, the duration of >100 TeV gamma-ray emission from an adjacent cloud would be quite short-lived, as the highest energy particles will fully traverse the cloud early in the evolution of the SNR.

DEPENDENCE OF GAMMA-RAY FLUX ABOVE 100 TEV ON SYSTEM PROPERTIES
In this section, we explore how different system properties affect the gamma-ray flux produced. This is shown for our two-step model with initial propagation in the ISM to reach the target cloud. Again, we assume a power law acceleration spectrum with slope = 2, as in table 1. The following properties are fixed throughout this section: fast value for 0 , SNR radius at 30 pc, cloud radius at 10 pc and distance from Earth to the SNR at 1 kpc.

Variation in System Properties
To explore the phase space, densities from 1 to 1000 cm −3 were chosen as test values, along with separation distances from 1 to 500 pc and ages from 1 to 500 kyr. All possible combinations of these three variables were evaluated and are shown in figure 6.
The largest flux is seen for low separation distances and high densities as expected. Provided the time is sufficient for particles to fully traverse the cloud, then the resulting flux is subsequently independent of the age; this is seen particularly at low separation distances. Regions of the phase space at large distances and young ages do not have a detectable gamma-ray flux above 100 TeV, because the time elapsed is insufficient for particles of the corresponding energy (∼PeV) to travel the distance to the target material. Additionally, the flux is reduced towards lower target material densities, as the probability for p-p interactions is reduced.

Diffusion Coefficient Variation
Equation (7) describes the diffusion coefficient dependence on the particle energy , the ambient magnetic field and the suppression factor , which is related to the level of magnetic field turbulence ). It should be noted that is assumed constant within the ISM and in this section we consider variation in within the cloud only, with particles fully traversing the cloud. Figure 7 shows how the same gamma-ray flux level is reached by variation in the interstellar cloud properties, keeping the accelerator age and distance to the cloud fixed. In this case, the diffusion suppression coefficient, , is plotted against the cloud density , which relates to the amount of target material. With smaller values of , the diffusion coefficient is more strongly suppressed, potentially corresponding to higher levels of magnetic field turbulence. This shows that the same gamma-ray flux can be achieved for a range of interstellar clouds with different properties (roughly analogous to the phase space for particle accelerators (Hillas 1984)). The lines of constant flux show that the same flux level can be achieved with lower cloud densities if diffusion there is more efficient, namely particles remain trapped for longer times. For fixed , higher densities lead to higher gamma-ray fluxes, as expected from equation (13) and also seen in figure 3. For fixed density, lower values lead to higher gamma-ray fluxes. Suppression of diffusion coefficient might be achieved in the presence of a high ion-to-neutral ratio within the cloud (Nava et al. 2016). Conversely, a high presence of neutral particles would act in the opposite direction, namely by suppressing magnetic turbulence in the cloud because of ion-neutral damping (Kulsrud & Pearce 1969), hence increasing the diffusion coefficient.
As the density of an interstellar cloud can be inferred from catalogues and the separation distance and accelerator age determined for known systems, this relation between and could be used to constrain the diffusion properties within the cloud, given a gamma-ray flux measurement. For this purpose, gamma-ray upper limits could also be used to constrain the diffusion properties of the system.

Predicting gamma rays above 100 TeV from SNR-cloud pairs
To find suitable target clouds for gamma-ray detection, we use the Green's catalogue of SNRs (Green 2019) and the interstellar cloud catalogue of Rice et al. (2016) using 12 CO data from the DAME survey (Dame et al. 2001), which includes a distance estimator for the interstellar clouds. Rice et al. (2016) covers Galactic latitudes ranging from | | ≤ 1 • to | | ≤ 5 • , depending on the Galactic quadrant. Longitudes within ∼ 13 • of the Galactic Centre ( = 0 • ) are omitted by the survey. We scan the catalogues for pairings of interstellar clouds and SNRs that either overlap along the line of sight or are within 100 pc distance of each other, as evaluated from the angular separation at a given distance.
Only a limited number of SNRs have age and distance estimates from either Green (2019), Ferrand & Safi-Harb (2012) or Vukotić et al. (2019); where this information is available, we impose a distance cut of the cloud distance being within 20% of the SNR distance along the line of sight for a plausible pair. In the majority of cases, however, no distance estimate is provided for the SNR and we assume for the purposes of the flux calculation that the SNR and interstellar cloud are located at the same distance along the line of sight (i.e. that of the interstellar cloud). For several SNRs, this results in calculated fluxes for multiple interstellar clouds at different assumed distances. For SNRs without an age estimate, the SNR age is evaluated by inverting equation (11), given the SNR radius in pc (at the assumed distance).
An acceleration particle spectrum of ( , , ) ∝ −2 is assumed at the source location, as in section 2. Here we take the true size of the interstellar cloud, again with spherical symmetry assumed.
A total of 550 potential SNR-cloud pairs are found; however, for many of these pairs the age of the system is either not sufficient to allow particles to reach the cloud, or is long enough that the highest energy particles have fully traversed the cloud. Less than 10% of these pairings (∼ 40) result in a detectable very high energy gammaray flux above 10 TeV produced through particle interactions with the cloud under the slow diffusion scenario. Figure 7 shows that density dominates over diffusion in determining the gamma-ray flux. The gamma-ray flux was predicted using both fast and slow values for 0 and fixed (see table 1). In both cases, few combinations result in a gamma-ray flux above 100 TeV > 10 TeV Results presented in the remainder of this section adopt the slow diffusion case in all figures. As the spectral shape predicted by this model starts to drop at the highest energies at early times, a larger number of clouds are predicted to result in a detectable gamma-ray flux above 10 TeV. We therefore provide results both for the integral flux above 100 TeV, corresponding to the most promising candidates, and for the integral flux above 10 TeV as providing a larger selection of candidates to which current instruments, including H.E.S.S. and HAWC, are sensitive. This will enable us to compare our model results to existing data in Appendix B.
The predicted gamma-ray flux above 10 TeV for SNR-cloud pairs detectable by CTA-S within 50 hours 4 is shown in colour scale in the phase space of system properties in figure 8. Non-detectable systems are also plotted in figure 8 according to their properties as smaller circles in grey.
The median properties of systems detectable > 100 TeV are listed in table 2. In general, the gamma-ray detection prospects are enhanced for systems that are closer to Earth, for younger SNRs and for larger amounts of target material. The median cloud density for detectable systems is higher than that of non-detectable, as well as the detectable clouds being on average more massive. Gamma-ray spectra for a system with the median properties of both detectable and non-detectable systems shown in figure 9. WCDs such as HAWC, LHAASO and SWGO have improved sensitivity over IACTs towards 100 TeV, and may therefore be sensitive to a wider variety of systems.

Total gamma-ray spectrum per cloud
The total predicted gamma-ray flux from a given interstellar cloud is here overestimated, under the assumption that all SNRs contributing to the total flux are within a distance comparable to the cloud along the line of sight. This obviously cannot be true in reality, as SNRs without an intrinsic distance estimate have different assumed distances for various cloud pairings. Therefore, if one cloud is as bright  Table 2. Median properties of SNR-cloud pair systems predicted to have a gamma-ray flux above 10 TeV or 100 TeV that is detectable (det.) or non-detectable (non-det.) by CTA-S, under the slow diffusion scenario. as predicted, then another cloud associated with the same SNR yet at a different distance along the line of sight must be correspondingly dimmer in gamma-ray flux.
In practise, when observing the sky in gamma rays, it will be challenging to disentangle the contributions of different SNRs (or other nearby accelerators) to the total gamma-ray flux emanating from an interstellar cloud. We therefore sum the predicted contributions of all SNRs to obtain a prediction for the maximum total gamma-ray flux from a given interstellar cloud. Figure 10 shows the total gamma-ray spectrum from an arbitrary sub-sample of target clouds (approx. 5%), not all of which would be detectable by CTA. The right panel of figure 10 shows the fractional contribution of individual SNRs to the total flux per cloud; it is clear that in the majority of cases the emission is dominated by a single SNR, with other contributions only minor in comparison.
In table 3, the interstellar clouds are listed according to the predicted gamma-ray flux above 100 TeV. The predicted gamma-ray fluxes are compared to flux upper limits taken from the H.E.S.S. Galactic Plane Survey (HGPS) with a 0.2 • correlation radius. As the HGPS coverage is from = 250 • to 65 • and | | ≤ 3 • , no upper limits were evaluated for clouds located outside of this range. These flux upper limits evaluated using a 0.2 • correlation radius were scaled to the true angular size of the cloud assuming a scaling behaviour as: where psf = 0.2 • in this case. The flux upper limits from the HGPS assume a power law gamma-ray spectrum with an index of −2.3, enabling the integral flux upper limits and energy threshold to be converted into a comparable integral energy flux upper limit above 100 TeV.
Although the predicted flux from this model exceeds the upper limit from the HGPS in several cases, we note that for most of these clouds, multiple SNRs contribute to the total flux prediction for a given cloud. Figure 10 shows that in the majority of cases, the fractional contribution from a single SNR dominates the total flux (> 90%) whilst other SNRs only providing a minor contribution (< 10%). In table 3, the fractional contributions from different SNRs are explicitly given; in many cases, one SNR dominates and the majority of the individual contributions to the total flux from different SNR lie below the upper limit from the HGPS.
As noted above, in reality the corresponding SNRs will be located at different distances along the line of sight and cannot simultaneously contribute to the gamma-ray flux from a single cloud. In fact, flux upper limits corresponding to the locations of illuminated clouds could be used to constrain the line-of-sight distance to SNRs; in cases where a detectable flux is expected should the SNR be at a comparable distance to the cloud.
For several cases, there is no known distance estimate to the SNR, which was assumed to be located at the same line-of-sight distance as the cloud for the purpose of this model. Therefore, if the predicted flux were to exceed the flux upper limit from H.E.S.S., this would imply that the SNR is not located at a distance comparable to the cloud. Those SNRs for which a distance estimate was available are indicated in bold in table 3.
Further uncertainties that could contribute to the predicted gammaray flux are discussed in section 4.4.

Including contributions from 'invisible' SNRs traced by pulsars
Given that on average three SNR are expected every century, the total number of known SNR according to Green's catalogue, 294, is rather low. Out of the 126 SNR with age estimates, 86 are estimated to be younger than 20 kyr and therefore within the Sedov-Taylor expansion phase; compared with the 600 SNR expected from the average rate. This could be due to the fraction of the Galaxy that has been observed; surveys such as THOR and GLEAM are showing promising results  However, there are many more pulsars known, with a total of 2360 detected within our Galaxy (Manchester et al. 2005). 5 Every pulsar is born in a supernova explosion, due to the collapse of the stellar core as the outer layers are expelled outwards into the ISM. As such, the presence of pulsars can be used as an indicator for SNRs that are otherwise 'invisible'; dimly radiating or no longer actively accelerating particles. The ATNF catalogue was used to select the most promising pulsar candidates, with 63 listed as located within our Galaxy (distance < 25 kpc); characteristic age younger than 30 kyr and are not millisecond pulsars (rotation period > 10 ms).
The age of the progenitor SNR is then assumed to correspond to 5 Millisecond pulsars are excluded from this estimate.
the pulsar characteristic age = /2 ; however this assumes magnetic dipole radiation (with braking index = 3) and is known to be a rough approximation to the true pulsar behaviour (Gaensler & Slane 2006). Additionally, several pulsars have a measurable proper motion; during the original SNR explosion, the pulsar may receive a kick velocity propelling it to travel away from the birth location (Gaensler & Slane 2006). In these cases, the birth location of the pulsar was determined by extrapolating the proper motion of the pulsar in the opposite direction for a distance corresponding to the characteristic age -this birth location was assumed to correspond to the central position of the corresponding 'invisible' SNR. Similarly, the hypothetical SNR radius could be calculated based on the age , and additional plausible pairs of hypothetical SNRs and nearby interstellar clouds were identified. The predicted gamma-ray fluxes arising from these clouds due to invisible SNRs as traced by pulsars were then added to the gamma-ray flux predictions calculated above. This then provided a revised estimate for the predicted gamma-ray flux from each cloud. Fractional contributions to the gamma-ray flux above 100 TeV from hypothetical SNRs are given in table C1 of Appendix C.
Many more systems are predicted to result in a detectable flux above 10 TeV, which nevertheless implies the presence of 0.1 PeV protons. Up to ∼ 40 of these may be detectable by current IACT facilities. The full list is provided in table A1 of Appendix A. Maps of the predicted flux above 10 TeV from clouds in the Galactic Plane in figure 11; and above 100 TeV in figure 12 show the total predicted gamma-ray flux, including contributions from both catalogued and hypothetical SNRs.

Systematic Uncertainties due to model assumptions
In order to estimate some of the systematic uncertainties inherent due to assumptions in the model, we tested the variation in predicted flux arising from varying the spectral slope of proton injection.
As described in section 2, we assumed a power law acceleration spectrum with slope = 2. This was found to have a strong influence on the predicted integral flux above 100 TeV. Adopting a hard value of = 1.8 increased the predicted flux by up to a factor ∼ 30, whilst a softer slope of = 2.2 decreased the flux marginally with respect to = 2. Note that observationally most of the gamma-ray detected SNRs show steeper spectra than −2 , especially middle-aged ones. An alternative spectral shape for the initial proton spectrum, such as a broken power law or power law with exponential cut-off, may also be expected to reduce the predicted flux at the highest energies. The resulting variation in cloud detectability is indicated in table D1 in the Appendix.
Similarly, the flux upper limits from the HGPS were derived assuming a spectral index of 2.3, whilst the mean spectral index of sources listed in the HGPS is 2.4 ± 0.3. Therefore, we varied the spectral index for the derived upper limit and found that an index of 2.1 relaxes the upper limit by a factor 3, whilst an index of 2.5 tightens the upper limit by a factor 3 (with respect to the quoted value for an assumed index of 2.3).
Flux predictions provided by this model can hence be considered an order of magnitude estimate. In table 3 , the predicted fluxes that currently lie below the detectable limit may be enhanced by an order of magnitude within the limitations of the model.
As mentioned in section 2, the Sedov time occurs at characteristically different ages for different supernova types. Although a distinction between core collapse and type Ia events can be made for a small number of SNRs based on their immediate environment, for the majority of cases this is uncertain. The core collapse scenario was adopted as default, due to the progenitor massive stars being somewhat younger than the progenitor white dwarfs of type Ia events and therefore more likely to be in environments harbouring rich molecular clouds. Nevertheless, we explored the influence of a type Ia Sedov time on the results, finding that only two of the clouds predicted to be detectable by CTA remain so under the type Ia scenario.
The variation in predicted flux under type IIb and type Ia scenarios implies that flux constraints on clouds may also be used to investigate intrinsic properties of the SNR.

DISCUSSION
As explained in section 4.1, although the integral gamma-ray flux above 100 TeV indicates the most promising candidates for evidence of PeVatron activity, the integral gamma-ray flux above 10 TeV provides a broader range of candidates to which current generation instruments are sensitive. Bright gamma-ray flux above 10 TeV nevertheless indicates the presence of energetic particles within an order of magnitude of the CR knee. A summary of the clouds predicted by our model to be detectable by H.E.S.S. above 10 TeV (under the slow diffusion scenario) is provided in Appendix A.

Potentially bright clouds overlapping HGPS sources
Appendix B shows an enlarged view of the clouds in the Galactic plane predicted to harbour a gamma-ray flux above 10 TeV detectable by H.E.S.S. (as in figure 11). Figures B1 to B4 show significance contours at 3, 5 and 15 from the H.E.S.S. Galactic Plane Survey (HGPS) overlaid, together with the name and nominal location of each known source. Due to the limited coverage of the surveyed clouds by Rice et al. (2016) and Dame et al. (2001), the Galactic latitude is restricted to ≤ 2.5 • . Similarly, due to the limited coverage of the HGPS, the longitude is restricted to between = 64 • to = 262 • running through the Galactic Centre.
In several regions, potentially detectable clouds are seen to overlap with known H.E.S.S. sources, although caution should be exercised over the assumption of spherical symmetry. We discuss here a few selected regions in more detail in turn.

49 • > > 43 •
At ∼ 46 • , there are two clouds overlapping with a gamma-ray significance contour not currently associated with a H.E.S.S. source. This position is coincident with the unidentified TeV source discovered by HAWC, 2HWC J1914+117. The larger cloud at (45.86,0.01) is illuminated by the SNR G043.9+1.6, with which a shell-like Fermi source is associated, whilst the small cloud at this location is associated with the SNR G046.8-0.3. Two pulsars are located nearby, yet without estimates: PSR J1915+1144 and PSR J1915+1149 at distances of 7.2 kpc and 14 kpc respectively. The predicted gamma-ray flux is compatible with the upper limit derived from the HGPS. A further cloud coincident with HESS J1911+090 is illuminated by particles escaping from the SNR G042.8+0.6, which itself is tentatively associated with the nearby source HESS J1908+063. HESS J1911+090 is identified as W49B, corresponding to the SNR G043.3-0.2 already known to be interacting with molecular clouds.

40 • > > 37 •
In this region, the cloud with the highest predicted flux, at (39.33,-0.32), is well compatible with the upper limit from the HGPS, with contributions from the SNRs G040.5-0.5 and G038.7-1.3, with a distances constrained at 3.2 kpc and 4 kpc respectively. It is unlikely that both associated SNR contribute to the total flux; the prediction from either SNR separately is approximately half of the total predicted flux for this cloud.

37 • > > 34 •
This region is shown in the second panel of figure B1, where a group of clouds overlap with HESS J1857+026 and HESS J1858+020. Both of these sources remain unidentified, although it is known that there is substantial molecular material in the vicinity (Paredes et al. 2014). HESS J1857+026 was best described by two components that were merged to a single source in the HGPS, whereas the MAGIC collaboration favour a two source interpretation, associating MAGIC J1857.2+0263 with the pulsar PSR J1856+0245, whilst MAGIC J1857.6+0297 remains unidentified (H.E.S.S. Collaboration et al. 2018;MAGIC Collaboration et al. 2014). A higher resolution view of the source reveals that the gamma-ray peak may coincide with a wind-blown cavity, suggesting a leptonic origin of the radiation; however there is considerable overlap with molecular gas and a hadronic origin is not excluded by MAGIC Collaboration et al. (2014).
The predicted gamma-ray flux from the brighter cloud at (36.100,-0.140) is provided by escaped particles from the SNR G035.6-00.4, for which the distance is constrained. MAGIC J1857.6+0297 is positionally coincident with this cloud, whilst G035.6-00.4 is itself associated with the TeV source HESS J1858+020. Detailed gammaray studies of this region with high angular resolution are required to disentangle the multiple components.
Similarly, for HESS J1858+020, (Paredes et al. 2014) note the rich molecular environment and suggest that the non-thermal emission may arise from interactions between a SNR and nearby clouds. The predictions of our model that clouds should be detectable in gammarays in this region is therefore consistent with the data.
The large and rather bright cloud at (34.99,-0.96) is illuminated predominantly by the SNR G036.6-0.7, for which no TeV detection has yet been made, despite a GeV detection by Fermi-LAT (Ferrand & Safi-Harb 2012). However, as the distance to the SNR is known the flux prediction, which lies below the HGPS upper limit, is uncertain.

34 • > > 31 •
The cloud overlapping with HESS J1852-000 -at (32.73,-0.42) -has a predicted flux with contributions from both the SNR G033.2-0.6 and G033.6+0.1, as well as lesser contributions from several other SNRs. The flux predicted by our model is consistent with the upper limit derived from the HGPS. However, it has been previously suggested that this unidentified TeV source is associated with a molecular cloud interaction, although a PWN associated with pulsars in the vicinity remains a plausible alternative scenario (H.E.S.S. Collaboration et al. 2018).
A few clouds coincide with contours around HESS J1848-018; a complex source without a firm identification yet potentially corresponding to an SNR interacting with adjacent molecular clouds. The location is adjacent to several clouds, such as the cloud located at (28.770, -0.090) and illuminated predominantly by the SNR G028.8+1.5. This prediction is fairly reliable, as the SNR has a known distance of 2.8-4 kpc. Similarly, the adjacent cloud at (27.64,0.07) is predicted to have a gamma-ray flux produced by particles escaping from the SNR G028.6-0.1, of known distance and potentially associated to the source HESS J1843-033.

27 • > > 25 •
Several clouds overlap with the TeV gamma-ray sources HESS J1841-055 (unidentified) and HESS J1837-069, which is firmly identified as a pulsar wind nebula, powered by the energetic pulsar PSR J1838-0655 (H.E.S.S. Collaboration et al. 2018). Nonetheless, the SNR associated with clouds around = 25 • − 26 • in our model is G024.7+0.6, from which MAGIC has reported the detection of a source not listed in the HGPS; MAGIC J1835-069 located at (24.94, 0.37) (MAGIC Collaboration et al. 2019). The predicted fluxes are, however, compatible with derived upper limits from the HGPS. As G024.7+0.6 has a known distance of 3.2 -3.7 kpc, this provide some confidence in the possibility that escaped particles may indeed be interacting with these clouds; yet the resulting flux is likely below that of the dominant gamma-ray source HESS J1837-069, and therefore not currently identifiable. Note that the adjacent cloud at (23.80,1.58) is also associated with SNR G024.7+0.6.

23 • > > 20 •
In this region, there is a large, bright cloud nestled between three known H.E.S.S. sources; the gamma-ray binary HESS J1832-093; the pulsar wind nebula HESS J1833-105 (associated with the pulsar PSR J1833-104 and the SNR G021.5-00.9); and the dark source HESS J1828-099, with no known associations. The cloud in question, 4 at (21.970, -0.290), derives its predicted flux above 10 TeV from a combination of the SNRs G021.6-00.8, about which little is known, and from G021.5-0.9. The predicted integral flux above 10 TeV is consistent with the derived upper limit, although as discussed in section 4.4, uncertainties of around an order of magnitude are inherent in our model assumptions.
From figure B1, it can be seen that there appear to be significant contours at the location of this cloud. Indeed, there is a Gaussian component HGPSC 072 reported here with a TS (test statistic) of 59, well above the detection threshold of TS= 30, and a reported flux above 1 TeV of 8.1 × 10 −13 cm −2 s −1 (H.E.S.S. Collaboration et al. 2018). This component was, however, not reported as a gamma-ray source due to the lack of a significant detection in the independent cross-check analysis (see section 4.9 of H.E.S.S. Collaboration et al. (2018)). We conclude that the lack of a currently known TeV gammaray source at the location of this cloud does not rule out the existence of detectable gamma-ray emission emanating from this cloud due to interactions with energetic particles originating from the SNRs G021.6-00.8 and G021.5-0.9.

20 • > > 14 •
This is a complex region in which the large, bright pulsar wind nebula HESS J1825-137 dominates, exhibiting strong energy dependent morphology with extended emission reaching as far as the gammaray binary LS 5039 (HESS J1826-148) and rendering the source HESS J1826-130 identifiable only at the highest energies (H.E.S.S. Collaboration et al. 2019). The presence of molecular material in this region is, however, well established; the offset of the leptonic gammaray emission (i.e. the parent electron population) from the pulsar PSR J1826-1334 is thought to be due to shock interactions between the PWN and molecular material towards HESS J1826-130. Supporting this scenario is the identification of high turbulence within the nearby cloud at (18.16, -0.34) (Voisin et al. 2016). It is worth noting that emission from this region has recently been confirmed to extend beyond 100 TeV, as reported for the source eHWC J1825-134, which overlaps multiple H.E.S.S. sources due to the limited angular resolution of HAWC (Abeysekara et al. 2020).
Adjacent to HESS J1826-130, the cloud at (18.81,-0.46) is illuminated by the SNR G018.1-0.1 which is already potentially associated to the HESS source itself; the distance to both cloud and SNR is estimated at around 5.6-6.8 kpc.
Many clouds are seen to overlap with HESS J1825-137; however, clouds with predicted fluxes lower than that measured from the PWN are likely dwarfed by its gamma-ray emission. The larger cloud at (16.970, 0.530) slightly offset from the PWN has two SNRs contributing to the total predicted flux; G019.1+0.2 and G018.9-01.1, whereby only the latter has an estimated distance. The contribution from this more certain SNR-cloud association is at the level of ∼ 40% of the total predicted flux (see table A1). For this system, the total predicted flux is compatible with the derived upper limits from the HGPS. Note that the bright cloud at (16.97, 0.53) is additionally illuminated by a hypothetical SNR, corresponding to the pulsar J1826-1256; which is the pulsar potentially powering the HESS J1826-130 PWN (see table  C1).
The brighter cloud at (16.24,-1.02) has a predicted flux that exceeds the HGPS upper limit and is almost evenly split between two SNRs, G016.4-0.5 and G016.0-0.5 (see tables 3 and 4), neither of which has an available distance estimate. As the region has been extensively observed in recent years, we can surmise that it is likely the primary contributing SNR, G016.4-00.5 (flux >10 TeV of 5.81 × 10 −13 TeV cm −2 s −1 , table A1), is not located at ∼ 2 kpc and hence not associated with this cloud.
The smaller cloud at (14.14,-0.59) has a predicted flux exceeding the HGPS upper limit by ∼ 50%, arising solely from the SNR G015.1-1.6, which has a known distance of ∼ 2.1−2.2 kpc well compatible with the distance to the cloud of 2.09 kpc. In this case, one or more of the model assumptions concerning propagation through the intervening medium, age of the SNR or spectral properties of the accelerated particles, must be incorrect.

347 • > > 341 •
Two of the known TeV sources in this region, HESS J1708-410 and HESS J1702-420, remain without any clear counterparts or multiwavelength associations. Several of the smaller clouds coincident with significant emission in the region are predictions for particles accelerated by the SNR G344.7-0.1, which is known to be interacting with clouds.
An extended cloud at (343.64,-0.54) is illuminated primarily by the SNR G343.1-2.3, for a total flux at a level of a quarter of the upper limit obtained from the HGPS; for this SNR again the distance estimate of 2 kpc is based on an associated pulsar.
The small bright cloud at (341.34,0.21) is illuminated by the SNR G341.2+0.9 for a predicted flux slightly above the corresponding HGPS upper limit, although again the distance to this SNR is unknown, suggesting that it is plausibly not located at a distance comparable to the cloud at ∼ 2.4 kpc Several of the cloud predictions at Galactic longitudes slightly lower than HESS J1702-420, derive a significant fraction of their predicted flux from the SNR G342.0-0.2, for which no distance or age estimates are provided. Given the absence of significant gammaray emission, this may indicate that the SNR is not located at the same distance as the clouds (of ∼ 2.5 kpc).

341 • > > 336 •
Westerlund 1, with the designation HESS J1646-458, is a massive stellar cluster with highly complex gamma-ray morphology comprising multiple components (Abramowski et al. 2012). Massive stellar clusters are thought to be a suitable complementary (or even alternative) PeVatron candidate to SNRs, accelerating cosmic rays through stellar winds and multiple shock interactions from several supernovae (Aharonian et al. 2019). There is evidence for at least one young SNR 10 kyr in Westerlund 1, established by the presence of a magnetar associated with the cluster (Muno et al. 2006). Although for this paper we focus on SNRs, a dedicated study would be required to model the complex Westerlund 1 region in detail. Here, we merely note that the predicted fluxes for clouds at Galactic longitudes 338 • are in general compatible with derived upper limits from the HGPS, whilst the reported flux has been quoted with an uncertainty of a factor 2 (H.E.S.S. Collaboration et al. 2018).
There are several clouds situated around HESS J1640-465 and HESS J1641-463, where the former is a composite SNR and the latter is (similarly to HESS J1826-130) a hard spectrum source of unknown origin that reveals its presence towards higher energies (Abramowski et al. 2014). Molecular material coincident with the emission 'bridge' between the two H.E.S.S. sources has been identified (Lau et al. 2017b), however we note that the gamma-ray emission contours do not consistently correspond to all of these clouds. For the cloud at (337.84,-0.40) located at ∼ 2.91 kpc, the SNR G337.2-0.7 has an estimated distance spanning a wide range from 2 − 9.3 kpc such that the true distance may indeed be different from ∼ 3 kpc. The other contributing SNRs do not have firm distance estimates, such that these model results are in any case intrinsically uncertain.
Next, we turn to clouds overlapping at ∼ 337 • with HESS J1634-472 -an unidentified source known to overlap dense CO gas and with plausible SNR associations (H.E.S.S. Collaboration 2010). In this case, it seems that our model may be verified and the emission is associated with escaped particles from nearby SNRs. At least some of the distance estimates are broadly consistent with the large 14 kpc distance -for several clouds in the vicinity, the reported distance is ∼ 10 − 14 kpc.
Lastly, we note the brighter cloud at (336.73,-0.98) which also derives its flux predominantly from the aforementioned SNR, G337.2-0.7, which has an estimated distance in the wide range 2 -9.3 kpc, comptable with this cloud at 3.24 kpc.

334 • > > 330 •
Previous studies have established the presence of diffuse gas overlapping with HESS J1616-508 (a likely PWN) and HESS J1614-518 (an SNR candidate with shell-like morphology) along the line of sight, yet without being a plausible association (Lau et al. 2017a).
The neighbouring bright cloud at (333.46,-0.31) derives its predicted flux almost evenly split between the SNR G332.4-0.4 and a hypothetical SNR -whilst the contributions are separately comparable to the derived upper limit, the summation is in violation of it. As the SNR G332.4-0.4 has an estimated distance of 2.7 -3.3 kpc and is already known to be interacting with molecular clouds, it is the contribution from the hypothetical SNR which is considered less reliable and can be neglected.

330 • > > 325 •
All of the predictions for illuminated clouds in this region are compatible with upper limits obtained from the HGPS. The brightest prediction, for the cloud at (328.91,-0.13) and a distance of 1.51 kpc, originates from an association with the SNR G329.7+0.4, for which neither age nor distance are constrained. 5.1.13 324 • > > 318 • Firstly, we note that the most prominent cloud in this region at (322.51,0.17); as listed in Table A1, remains nevertheless compatible with the upper limit derived from the HGPS. This cloud is assumed to be illuminated by CRs from two SNRs, G321.9-1.1 and G320.6-1.6, about which little is known and no distance estimate is currently available. As the total flux is dominated by the former SNR, tighter upper limits would constrain the distance to the former SNR more than the latter.
Further clouds can be seen to overlap with HESS J1457-593, associated with the SNR G318.2+0.1 (distance estimated between 3.5 kpc and 9.2 kpc); a source that was first announced with the HGPS (H.E.S.S. Collaboration et al. 2018). This SNR is suggested by our model to illuminate the cloud at (319.05,0.61), adjacent to HESS J1503-582, whilst clouds overlapping with HESS J1457-593 are illuminated by G318.9+0.4, for which no distance is available.
We note that a cloud located at (318.07, -0.21) was recently analysed by Aharonian et al. (2020) using data from Fermi-LAT to probe the CR sea by testing for agreement between predicted and measured flux. This cloud was among those exhibiting a higher gamma-ray flux than expected due to the CR sea alone, whilst our model does not predict an additional contribution due to nearby SNRs.

316 • > > 308 •
At ∼ 311 • two clouds overlap with significance contours not associated to any known TeV source. Both are predicted to be illuminated by CRs from the SNR G309.8+0.0 about which little is known. A dedicated analysis to search for emission from these clouds is worthwhile to either confirm our model predictions, or place constraints on properties of the corresponding SNR.
The rather bright cloud at (309.2,-0.96) results from the sum of CR contributions from three SNRs; whilst the individual contributions are well compatible with upper limits, the summation likely overpredicts the total flux associated to this cloud, resulting in a total flux incompatible with the HGPS upper limit in the slow diffusion case. Of the two contributing SNRs, only G308.8-0.1 has a constrained distance and can be considered the more reliable prediction, with the flux prediction due to this SNR alone lying approximately 20% below the HGPS upper limit (the other is G309.8+0.0 mentioned above).

287 • > > 282 •
Westerlund 2 (HESS J1023-575) is another massive stellar cluster; as for Westerlund 1, in such a rich environment (of massive stars and molecular material) we can expect energetic particles illuminating nearby clouds to be a frequent occurrence (Dame 2007). The resulting morphology is complex and may be comprised of multiple CR accelerators; with the improved angular resolution of CTA in conjunction with multi-wavelength information, the origins of gammaray emission in this region may be revealed (H.E.S.S. Collaboration et al. 2011).
The cloud at (284.12,-0.75) does, however, have a distance which is nevertheless compatible with the kinematic distance of molecular material in the region, ∼ 6 kpc as determined in previous studies, although the distance to Westerlund 2 itself is thought to be closer to ∼ 8 kpc (Dame 2007;Furukawa et al. 2014). This suggests that such a scenario may indeed contribute to, but not dominate the flux in the region.
The environments of massive stellar clusters may be worth further dedicated modelling efforts, to establish if these environments are indeed suitable PeVatron candidates (Aharonian et al. 2019). We note that the larger cloud at (286.62,-0.36) which does not overlap with any known TeV source is predicted to be illuminated by the SNR G286.5-1.2, generating a gamma-ray flux at a level approximately half that of the corresponding HGPS upper limit, although neither age nor distance are constrained for this SNR. With a continued lack of detectable emission with increased exposure in this region, this could constrain the distance to this SNR as lying outside the range of 2-2.5 kpc, where the cloud is situated.

Potentially bright clouds observable by HAWC
Although H.E.S.S. is based in Namibia observing the Southern sky and HAWC, based in Mexico, observes the Northern sky, there is a region of the Galactic Plane from 60 • 10 • which is observable by both. A dedicated study will compare the two instruments data in this overlapping region (Jardin-Blicq et al. 2019;Abdalla et al. 2021). For clouds within this common region, we do not repeat discussion of H.E.S.S. sources that have also been observed by HAWC, with the exception of sources for which gamma rays > 100 TeV have been detected (Abeysekara et al. 2020). These include eHWC J1907+063 and eHWC J1825-134; the location of the former is consistent with the source HESS J1908+063, whilst the latter is associated with the complex region discussed in section 5.1.8. From this study, there are no bright clouds > 100 TeV predicted that are coincident with eHWC J1907+063, which is consistent with the current scenario of the observed TeV sources above 100 TeV being associated with energetic pulsars (Abeysekara et al. 2020).
The centroid of the highest energy emission >100 TeV from eHWC J1825-134 is located between the two known sources HESS J1825-137 and HESS J1826-130. A recent detailed study attempts to disentangle contributions from multiple components in the region, with indications for an independent emission region emerging at the highest energies, potentially associated with molecular material (Albert et al. 2021).
There are several new TeV sources discovered by HAWC in the common region -with work ongoing to confirm their detection in the H.E.S.S. data (Jardin-Blicq et al. 2019). We now briefly discuss two of these regions.

54 • > > 52 •
Adjacent to HESS J1930+188 and nearby clouds lies the source HAWC J1928+178, newly discovered by HAWC and confirmed in subsequent HAWC catalogues, although it was not detected in an initial follow-up study by VERITAS and Fermi-LAT (Abeysekara et al. 2018). The dedicated comparison study of HAWC and H.E.S.S. analysis techniques does, however, appear to reveal HAWC J1928+178 in the HGPS dataset.
A dedicated study of this source suggested an association with the energetic pulsar PSR J1928+1746, implying a leptonic, PWN origin for the gamma-ray emission (Lopez-Coto et al. 2017). More recently, alternative scenarios including an illuminated molecular cloud or a gamma-ray binary were shown to remain plausible by Mori et al. (2020). They also state that further X-ray observations are necessary to clarify a PWN association for the emission.

47 • > > 45 •
The TeV source HAWC J1914+118 is located at (46.00, 0.25) and overlaps with the clouds located at (45.92,-0.45) and (45.86,0.01). As yet, there is no firm identification for HAWC J1914+118, with several HII regions, pulsars and X-ray sources either coincident with or nearby to the gamma-ray emission Abeysekara et al. (2017c). Our model predicts illumination of the former cloud by the SNR G046.8-0.3, and of the latter cloud by G043.9+1.6, both of which have constrained distance (Ferrand & Safi-Harb 2012). We suggest that further investigations are necessary to firmly establish the nature of this new TeV source.

Regions outside the HGPS: 260 • > > 60 •
Although there is a large portion of the Galactic Plane that is not covered by the HGPS, there are only a few candidate clouds identified by our model as potential targets for CRs from nearby SNRs. The most prominent of these is the cloud at (110.43, 1.89) illuminated by the potentially young (7.7kyr) SNR G114.3+0.3, at a known nearby distance of 0.7 kpc, compatible with the distance to the cloud of 0.6 kpc. The angular size of the cloud is quite large at 1.30 • , which will make it more challenging to detect with point source searches and will likely require a dedicated analysis (for any pointing gammaray instrument). Although this region is observable by MAGIC and VERITAS, observations of this region may have not yet taken place and/or dedicated searches not yet performed. An additional challenge comes from the sensitivity of HAWC degrading towards this region of the sky, with declination > +60 • Abeysekara et al. (2017c). According to our model, this could be a promising region to study for evidence of PeVatron activity from SNRs, suggesting that dedicated observations and analysis searching for extended sources in this region would be worthwhile.

Limitations of the model
Despite our efforts to construct a realistic model for particle acceleration, transport and interaction, there are nevertheless several additional known limitations of the modelling. Firstly, we have assumed that all SNRs are in the Sedov-Taylor phase, when evaluating the SNR age based on its current radius or equivalently in calculating the spatial evolution of the SNR to establish the radius at the escape time for particles of a given energy. This is typically valid for SNRs with ages between ∼ 300 yr and ∼ 20 kyr, and is therefore valid for the majority of the systems we consider. However, SNRs generally undergo three phases during their evolution; an initial free-expansion phase, prior to adiabatic expansion during the Sedov-Taylor phase; followed by a period of radiative expansion up to approximately ∼ 500 kyr. An extension to the modelling could therefore be to incorporate the radiative phase into the evolution. In the radiative phase, particles that remain attached to the shock (have not yet escaped) become less energetic with time (Chevalier 1999).
However, the transition time between different phases strongly depends on the system properties, including details of the supernova explosion and pre-existing conditions in the circumstellar medium, which are often unknown. Similarly, throughout the modelling of both the SNRs and the interstellar clouds we have assumed the morphology to be spherically symmetric, whereas in practise the morphologies can be significantly more complex.
The largest uncertainty of the model resides in the diffusion properties of both the ISM and the interstellar clouds, which remain so far elusive. We have assumed a suppression of the diffusion coefficient within the clouds of = 0.05 with respect to the ISM. As figure 7 shows, for the same cloud density , different values can lead to variations in gamma-ray flux. Also, the average diffusion coefficient within the ISM is uncertain and unlikely to be uniform.
The normalisation of 0 assumed in table 1 is a typical Galactic value as derived from CR diffusion time in our Galaxy . Recent measurements have shown, however, that the diffusion coefficient may be suppressed within regions of the ISM with turbulent activity due to nearby accelerators (Abeysekara et al. 2017a). Transport of energetic particles within the intervening ISM between an accelerator and a nearby cloud may therefore proceed somewhat slower than our assumed transport here. Comparing the results using two different assumed values of 0 provides an indication of this effect.
To estimate the magnetic field strength within interstellar clouds, we adopted in equation (8) a parameterisation based on the cloud density from Crutcher et al. (2010), using estimates from Zeeman splitting which is significant only for dense clouds. Although a fixed value comparable to the ISM field strength has been used for clouds with < 300 cm −3 , the magnetic field strength is likely to vary further. Alternative methods for estimating the magnetic field must then be used, such as Faraday tomography, which could revise the flux prediction (Van Eck et al. 2017).
To describe the energy dependence of the p-p interaction crosssection, we adopted the most recent parameterisation from Kafexhiu et al. (2014) in equation (6). This differs quite substantially from the parameterisation used in Kelner et al. (2006) at around 1 GeV, however the agreement improves and the two expressions are quite compatible from a few tens of GeV towards higher energies. A minimum energy of 10 GeV was used in the modelling, above which these two parameterisations are in reasonable agreement.
Additionally, Green's catalogue (Green 2019) has limited availability of distance and age estimates for the SNRs; where these are provided (for SNRs indicated in bold in tables 3, 4 and A1), uncertainties are typically large with a range of possible values given. Distance estimates to interstellar clouds may suffer from a distance ambiguity between 'near' and 'far' estimates based on the observed velocites, although Rice et al. (2016) removed this effect as far as possible, using the kinematic distance model of Reid et al. (2009).

Applications of the results
Despite the uncertainties surrounding several input parameters of the model, the consequence is that observations can in turn be used to constrain the system parameters. For example, in many cases there was no distance estimate available for a particular SNR, such that it was assumed to be located at a comparable distance to the nearby cloud. However, if the resulting flux prediction exceeds a measured upper limit, then the distance to the SNR can be constrained (as unrelated to the cloud). Lack of detectable gamma-ray emission at the highest energies may, however, also be due to gamma-ray absorption in surrounding radiation fields. For clouds at distances 10 kpc, this effect should be corrected for (Porter et al. 2018).
Similarly, as shown in section 3.2, the flux produced depends on the highly uncertain diffusion properties of the cloud. For a given density, the measured flux (or flux upper limit) can be used to constrain the diffusion properties of the cloud, provided the SNR -cloud association is secure.
As discussed in section 1, particles accelerated to PeV energies and beyond will escape the source environment, such that it is considerably more likely that evidence for PeVatron activity will be found from target clouds located nearby an accelerator, than at the accelera-tors themselves (Celli et al. 2020). We consider the predictions made here for the clouds listed in tables 3, 4 and A1 to be suitable targets for future observing programs with CTA, searching for evidence of PeVatron activity from SNRs.

Further work and next steps
This study can be extended to other source classes, such as stellar clusters and pulsars, considering not only impulsive but also continuous acceleration scenarios. Indeed, it has been shown that pulsars are capable of accelerating hadronic particles to energies beyond 1 PeV and injecting these into the pulsar wind; albeit at a rate dwarfed by the production of electron -positron pairs (Lemoine et al. 2015;Kotera et al. 2015). Massive stellar clusters have also been established as prime PeVatron candidates, that warrant dedicated modelling efforts Aharonian et al. (2019). Additionally, detailed studies of individual sources and specific regions can be made using the modelling framework established here.

CONCLUSIONS
Interstellar clouds are a suitable tool that can be exploited to search for evidence of PeVatron activity within our Galaxy. The case of target material for hadronic CR interactions being located near to but separated from an astrophysical accelerator is reasonably common. Under the hypothesis that energetic CRs escaping an SNR and illuminating nearby clouds could provide an alternative signature for PeV particle acceleration to the gamma-ray spectrum at the accelerator itself, we constrain the necessary phase space of system properties.
In general, a higher target cloud density and lower separation distances are preferred for a detectable gamma-ray flux at very high energies. The predicted flux is dominated by the amount of target material, however there is a slight preference for older systems, to allow sufficient time for particles to reach the cloud. Faster diffusion, consistent with measurements from the Boron to Carbon ratio in the ISM and corresponding to a lower magnetic field within the cloud, leads to lower flux levels due to fewer interactions within the cloud. Slower diffusion within the ISM, however, may prevent particles from reaching the cloud at all for larger separation distances.
The four brightest clouds, for which a detectable flux above 10 TeV is predicted under a range of model scenarios, are located at (l,b) = (333.46,-0.31), (16.97,0.53), (110.43,1.89) and (336.73,-0.98). We remark that these clouds are spatially extended with respect to the angular resolution that CTA will achieve. However, no significant degradation of the minimum flux detectable by the instrument is expected above 100 TeV for sources whose extension is less than 1 deg, which applies to the majority of the clouds in this study (Ambrogi et al. 2018).
On average, a detectable gamma-ray flux is more likely for more massive clouds; systems with lower separation distance between the SNR and cloud; and for slightly older SNRs. We provide a further key target list for observations with future facilities, and suggest that even gamma-ray upper limits can provide useful information, constraining the distances to SNRs or diffusion properties of the ISM.
Currently operational IACT facilities such as H.E.S.S., MAGIC and VERITAS have already demonstrated the capabilities of IACTs to perform detailed studies of specific regions, untangling the origins of the gamma-ray emission thanks to their good angular and energy resolution. The forthcoming CTA will further improve on these; this work provides a suitable list of target areas for further investigation. Facilities such as HAWC, LHASSO and the proposed SWGO employ WCD technology; these are survey instruments with a large effective area particularly towards the highest gamma-ray energies, ideally complementing the IACT observations (Abeysekara et al. 2017b;Bai et al. 2019;Albert et al. 2019). Continuous observations by WCD facilities will rapidly establish locations of the highest energy gamma-ray emission, suitable for dedicated follow-up studies with IACTs. We encourage detailed studies of the promising regions highlighted by our model and look forward to further developments in the ongoing search for PeVatron candidates. Figure 11. Predicted integral gamma-ray flux above 10 TeV from interstellar clouds illuminated by escaped particles from nearby SNRs, using interstellar clouds identified in (Rice et al. 2016). Clouds detectable by CTA are shown in colour scale, whilst non-detectable are shown in grey. The total includes contributions from SNRs as listed in Green's catalogue (green circles) as well as contributions from hypothetical SNRs as traced by the energetic pulsar population (Green 2019). APPENDIX A: PREDICTED GAMMA-RAY FLUX FROM CLOUDS ABOVE 10 TEV Table A1. As for Table 3, for the predicted flux and upper limits above 10 TeV. All clouds in this table are predicted to have a flux above the H.E.S.S. 50 hour sensitivity threshold in the slow diffusion case.
Cloud ID Cloud coordinates Cloud size Cloud distance > 10 TeV total > 10 TeV # SNRs SNR flux contribution (l, b) deg deg kpc TeV cm −2 s −1 TeV cm −2 s −1 TeV cm −2 s −1    Table 3, for the predicted flux and upper limits above 100 TeV due to cosmic rays escaped from hypothetical SNRs, as traced by pulsars. In this case, the J2000 name of the pulsar is given, and all pulsars have an associated distance estimate. Note that only fractional contributions from hypothetical SNRs are shown; fractional contributions to the total flux from known SNRs are not given here (see tables 3 and 4). All clouds in this table are predicted to have a flux above the CTA 50 hour sensitivity threshold in the slow diffusion case.

APPENDIX D: SUMMARY OF PREDICTED CLOUD DETECTABILITY
This paper has been typeset from a T E X/L A T E X file prepared by the author. Table D1. Summary of cloud detection prospects based on the gamma-ray flux predicted by this model under different assumptions. These are evaluated for the integral gamma-ray flux above 100 TeV (compared to CTA), or 10 TeV (compared to H.E.S.S.) respectively. Unless otherwise specified, the faster value for the diffusion coefficient 0 from table 1, a particle spectrum index = 2 and core collapse supernovae (as opposed to type Ia, see section 2) are assumed.