# The population of Galactic supernova remnants in the TeV range

Rowan Batzofin, Pierre Cristofari, Kathrin Egberts, Constantin Steppa, Dominique Meyer

► To cite this version:

Rowan Batzofin, Pierre Cristofari, Kathrin Egberts, Constantin Steppa, Dominique Meyer. The population of Galactic supernova remnants in the TeV range. Astronomy & Astrophysics - A&A, 2024, 687, pp.A279. 10.1051/0004-6361/202449779. hal-04570582

HAL Id: hal-04570582

https://hal.science/hal-04570582v1

Submitted on 16 Aug 2025

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# The population of Galactic supernova remnants in the TeV range

Rowan Batzofin<fnref n="1" />, Pierre Cristofari<fnref n="2" />, Kathrin Egberts<fnref n="1" />, Constantin Steppa<fnref n="1" />, and Dominique M.-A. Meyer<fnref n="3" />

Received 28 February 2024 / Accepted 26 April 2024

## ABSTRACT

*Context.* Supernova remnants (SNRs) are likely to be significant sources of cosmic rays up to the knee of the local cosmic-ray (CR) spectrum. They produce gamma rays in the very-high-energy (VHE) (*E* > 0.1 TeV) range mainly via two mechanisms: hadronic interactions of accelerated protons with the interstellar medium and leptonic interactions of accelerated electrons with soft photons. Observations with current instruments have lead to the detection of about a dozen SNRs emitting VHE gamma rays and future instruments should significantly increase this number. Yet, the details of particle acceleration at SNRs and of the mechanisms producing VHE gamma-rays at SNRs remain poorly understood.

*Aims.* We aim to study the population of SNRs detected in the TeV range and its properties and confront it to simulated samples in order to address fundamental questions concerning particle acceleration at SNR shocks. Such questions concern the spectrum of accelerated particles, the efficiency of particle acceleration, and the gamma-ray emission being dominated by hadronic or leptonic interactions.

*Methods.* By means of Monte Carlo methods, we simulated the population of SNRs in the gamma-ray domain and confronted our simulations to the catalogue of sources from the High Energy Stereoscopic System (H.E.S.S.) Galactic Plane Survey (HGPS).

*Results.* We systematically explored the parameter space defined in our model, including for example, the slope of accelerated particles *α*, the electron-to-proton ratio *K*<sub>ep</sub>, and the efficiency of particle acceleration *ξ*. In particular, we found possible sets of parameters for which ≥90% of Monte Carlo realisations are found to be in agreement with the HGPS. These parameters are typically found at 4.2 ≥ *α* ≥ 4.1, 10<sup>-5</sup> ≤ *K*<sub>ep</sub> ≤ 10<sup>-4.5</sup>, and 0.03 ≤ *ξ* ≤ 0.1. We are able to strongly argue against some regions of the parameter space describing the population of Galactic SNRs in the TeV range, such as *α* ≤ 4.05, *α* ≥ 4.35, or *K*<sub>ep</sub> ≥ 10<sup>-3</sup>.

*Conclusions.* Our model is so far able to explain the SNR population of the HGPS. Our approach, when confronted with the results of future systematic surveys, such as the Cherenkov Telescope Array Observatory, will help remove degeneracy from the solutions and to better understand particle acceleration at SNR shocks in the Galaxy.

**Key words.** acceleration of particles – astroparticle physics – ISM: supernova remnants

## 1. Introduction

Cosmic rays (CRs) were discovered thanks to the pioneering work of Domenico Pacini and Victor Hess in the 1910s, and they have been extensively studied ever since (Hess 1912). Remarkably, the most fundamental question of CR physics remains unanswered; that is, we do not know where are they produced. The standard paradigm for the origin of Galactic CRs is that supernova remnants (SNRs) accelerate the bulk of CRs via the first-order Fermi mechanism called diffusive shock acceleration (DSA; Ptuskin 2010). DSA can efficiently energise protons (and ions), and electrons. Subsequently, the accelerated protons and electrons can interact with the interstellar medium (ISM) to produce gamma rays in the very-high-energy (VHE) domain, mainly through two mechanisms: (1) the production and decay of neutral pions, in the collision of accelerated protons with hadrons of the ISM (hadronic mechanism), and (2) the inverse Compton scattering of accelerated electrons on soft photons (cosmic microwave background, infrared, optical); leptonic mechanism. Numerous SNR shocks have now been detected in the gamma-ray domain, with instruments such as the Fermi Large Area Telescope (Fermi-LAT; Acero et al. 2016)

in the GeV range and, with Cherenkov instruments such as the Very Energetic Radiation Imaging Telescope Array System (VERITAS; Acciari et al. 2011), Major Atmospheric Gamma Imaging Cherenkov (MAGIC; Albert et al. 2007) and the High Energy Stereoscopic System (H.E.S.S.; H.E.S.S. Collaboration 2018a) in the TeV range. While detailed studies of individual objects have been performed extensively in the past (for example: RXJ1713−3946 H.E.S.S. Collaboration 2018e, Cas A van Veelen et al. 2009, Vela Jr Sushch et al. 2018, SN 1006 Pais & Pfrommer 2020), the study of the SNR population itself can unveil valuable information.

In particular, H.E.S.S., an array of imaging atmospheric Cherenkov telescopes located in the Khomas highlands in Namibia, performed a systematic survey of the Galactic plane. The H.E.S.S. Galactic Plane Survey (HGPS) covers the Galactic plane from longitudes of 250° to 65° and latitudes of −3° to 3°. This systematic survey found 78 sources in the TeV domain; eight were clearly detected SNRs, and 47 sources were classified as unidentified.

In this work, we intend to simulate the population of TeV-emitting SNRs and confront it to the population detected in the HGPS. We followed the framework presented in Cristofari et al.

<footnote n="1">¹ Universität Potsdam, Institut für Physik und Astronomie, Campus Golm, Haus 28, Karl-Liebknecht-Str. 24/25, 14476 Potsdam-Golm, Germany
e-mail: rowan.batzofin@uni-potsdam.de</footnote>

<footnote n="2">² Observatoire de Paris, PSL Research University, 61 Avenue de l'Observatoire, Paris, France</footnote>

<footnote n="3">³ Institute of Space Sciences (ICE, CSIC), Campus UAB, Carrer de Can Magrans s/n, 08193 Barcelona, Spain</footnote>
(2013, C13) and built on it with a refined modelling and a comparison with the HGPS SNR catalogue based on the method introduced in Steppa & Egberts (2020).

For this purpose, we simulated populations of Galactic SNRs using a Monte Carlo method. We simulated the explosion times and positions of the supernovae and calculated the gamma-ray emission from the SNRs at the current time. We constrained the simulations with measurements from the HGPS by identifying the detectable sources in the simulations based on the longitude, latitude, and angular-extension-dependent HGPS sensitivity and comparing it with the detected sample of SNRs in the HGPS. The inclusion of the complex dependencies of the HGPS sensitivity is crucial for drawing conclusions on the SNR population based on this data comparison. Due to the large number of unidentified sources (47) in the HGPS, the SNRs that are identified as such in the HGPS (8) are treated as a lower limit. We used the sum of firmly identified SNRs (8), firmly identified composite objects (8), and unidentified sources associated (by position) with a known SNR (12) as a stringent upper limit (28). For a strict upper limit (63), we include all the sources in the stringent upper limit as well as the rest of the unidentified sources (35).

With respect to C13, the present work notably includes several refined aspects with regard to the modelling; these include (1) the inclusion of a connection between the ejecta masses and explosion energies for the SNRs; (2) a refined description of the magnetic-field amplification and corresponding maximum energy for protons and electrons; (3) the inclusion of diffusive shock reacceleration at SNR shocks; and (4) a discussion of several prescriptions for the spatial distribution of SNRs in the Galaxy.

The paper is structured as follows: In Sect. 2, we describe the population model; this description includes the physics of individual SNRs and the distributions used to place them in the Milky Way. In Sect. 3, we discuss the methods used to compare our simulations to the HGPS, and we present our results in Sect. 4. We provide our conclusions in Sect. 5.

## 2. SNR population model

Following C13, in order to simulate the population of Galactic SNRs emitting in the TeV range, we need (a) a model for particle acceleration at SNRs and the subsequent gamma-ray emission (from hadronic and leptonic interactions of accelerated particles with the surrounding medium), which should be computationally inexpensive; (b) a description of the spatial distribution of SNRs in the Galaxy; and (c) a description of the gas density distribution in which the SNR shocks are expanding. These three ingredients are discussed in the following subsections.

### 2.1. Modelling the gamma-ray emission of SNRs

#### 2.1.1. Particle acceleration at SNR shocks

We relied on the usual assumption that a fraction $\xi_{CR}$ of the SNR shock ram pressure $\rho u_{sh}^2$ is converted into accelerated particles (loosely referred to as CRs in the following). To this end, we assume that the CR particle distribution at a time $t$ and at a position $r$ inside the shock follows a power law in momentum, as expected with DSA, such that $f_{CR}(R_{sh}, p, r) = A(r, t) \left(\frac{p}{m_{p}c}\right)^{-\alpha}$, where $p$ is the particle momentum and $\alpha$ is treated as a free parameter in the range $\alpha \ge 4$. We can find the normalisation $A$ by requiring the cosmic-ray pressure to be equal to a given fraction, $\xi_{CR}$, of the shock ram pressure at the shock location,

$$ \rho u_{sh}^2: \quad \frac{1}{3} \int_{p_{\min}}^{p_{\max}} dp \, 4\pi p^2 f_{CR}(p) \rho v(p) = \xi_{CR} \rho u_{sh}^2, \quad (1) $$

where

$$ A = \frac{3}{4\pi} \frac{\xi_{CR} \rho u_{sh}^2}{m_{p}^4 c^5 I(\alpha)} \quad (2) $$

and

$$ I(\alpha) = \int_{\frac{p_{\min}}{m_{p}c}}^{\frac{p_{\max}}{m_{p}c}} dx, \quad \frac{x^{4-\alpha}}{(1+x^2)^{1/2}}. \quad (3) $$

The spatial distribution of cosmic rays inside the SNR can be computed by solving the transport equation:

$$ \frac{\partial f}{\partial t} + u \nabla f - \nabla D \nabla f - \frac{p}{3} \nabla u \frac{\partial f}{\partial p} = 0, \quad (4) $$

where $D$ is the momentum-dependent diffusion coefficient for cosmic rays. Since particles with momenta $<p_{max}$ are expected to be well confined within the SNR, the diffusion term $(\nabla D \nabla f)$ is much smaller than the advection term $(u \nabla f)$. Eq. (4) can then be solved using the method of characteristics and using the following boundary condition: $f(R_{sh}, p, t) = f_0(p, t)$.

#### 2.1.2. Particle reacceleration at SNR shocks

To determine the contribution of gamma-ray emission from the reaccelerated particles, we followed Cristofari & Blasi (2019). The contribution from the reaccelerated particles is given by

$$ f_{0}^{\text{seed}}(p) = \alpha \int_{p_{0}}^{p} \frac{dp'}{p'} \left(\frac{p'}{p}\right)^{\alpha} f_{\infty}(p'), \quad (5) $$

where $\alpha$ is the same spectral index as before, $f_{\infty}(p)$ is the distribution function at upstream infinity of the seeds to be reaccelerated (in this case Galactic CRs), and $p_0$ is the minimum momentum of reaccelerated particles. The exact value of $p_0$ is not critical for Galactic CRs, as discussed in Blasi (2017); in this case we used $p_0 = 10^{-2}$ GeV. We assume that the spectrum of Galactic protons and electrons is the same as the local interstellar spectrum and use the parameterisations as in Bisschoff et al. (2019), which described the Voyager I (Cummings et al. 2016) and Payload for Antimatter Matter Exploration and Light-nuclei Astrophysics (PAMELA; Adriani et al. 2011) data well.

#### 2.1.3. Acceleration of electrons

Just as with protons, DSA is also efficiently energising electrons. The main difference is that the electrons suffer synchrotron and inverse Compton losses, thus affecting the spectra of accelerated particles. At low energies, where these losses can be neglected, the spectra of protons and electrons have the same shape. We introduced a parameter, $K_{ep}$, to describe the ratio of electrons to protons at these low energies. The maximum energy of electrons at the shock is found by equating the acceleration time at the shock to the synchrotron energy-loss time. Following the approach of Vannoni et al. (2009), the maximum energy of electrons is

$$ E_{\max}^{e} \approx 7.3 \left( \frac{u_{sh}}{1000 \text{ km s}^{-1}} \right) \left( \frac{B_{\text{down}}}{100 \text{ µG}} \right)^{-1/2} \text{ TeV}. \quad (6) $$
The electrons are accelerated on timescales much shorter than the synchrotron loss time. After being accelerated, the electrons are advected downstream of the shock, where they continue to lose energy through synchrotron radiation. The characteristic time for synchrotron radiation is

$$ \tau_{\text{syn}} \approx 1.8 \times 10^3 \left( \frac{E_e}{\text{TeV}} \right)^{-1} \left( \frac{B_{\text{down}}}{100 \, \mu\text{G}} \right)^{-2} \text{yr}, \quad (7) $$

where $E_e$ is the electron energy. The energy-loss time decreases with electron energy, and there could be an energy ($E_{\text{break}}^e$) beyond which the loss time is less than the age of the SNR. Above this energy the electron spectrum is shaped by radiative losses and steepens by one power in momentum (Morlino & Caprioli 2012). The SNR is expanding, so there are also adiabatic losses at a rate of $\tau_{\text{ad}} = \frac{R_{\text{sh}}}{u_{\text{sh}}}$. Hence, $E_{\text{break}}^e$ can be found by solving the following equation (Finke & Dermer 2012):

$$ \tau_{\text{age}}^{-1} = \tau_{\text{syn}}^{-1} + \tau_{\text{ad}}^{-1}. \quad (8) $$

### 2.1.4. Maximum energy of accelerated particles

The question of the maximum energy of accelerated particles at SNR shocks is fundamentally connected to the question of the amplification of magnetic fields at SNR shocks. For fast SNR shocks, the amplification of magnetic fields is expected to be determined by the growth of non-resonant streaming instabilities upstream of the SNR shock (Bell et al. 2013; Schure & Bell 2013; Cristofari 2021). As the SNR shock slows down, resonant streaming instabilities become relevant (Schure & Bell 2012). We account for the magnetic field amplification as in Cristofari et al. (2021a). For the resonant streaming, the upstream magnetic field is estimated as

$$ B_{\text{up}}^2 \approx \frac{1}{M_A} \xi_{\text{CR}} \nu_{\text{sh}}^2, \quad (9) $$

where $M_A = v_{\text{sh}}/v_{\text{A,0}}$ is the Alfvénic Mach number.

In the non-resonant case, the magnetic field stops growing when the energy density in the form of escaping particles equals that in the amplified field:

$$ \frac{\delta B^2}{4\pi} \approx 3 \frac{v_{\text{sh}}}{c} \frac{\xi_{\text{CR}} \rho v_{\text{sh}}^2}{(\alpha - 3)I(\alpha)} \left( \frac{p_{\text{min}}}{m_{\text{p}}c} \right)^{4-\alpha}. \quad (10) $$

Assuming the turbulent field upstream of the shock is roughly isotropic and the perpendicular components are compressed at the shock, the downstream magnetic field is given by $B_{\text{down}} \approx B_{\text{up}} \sqrt{\frac{1+2r^2}{3}}$, where r is the compression factor (set to 4 in this work). The amplified magnetic field for the SNR is given by: $\sqrt{B_{\text{resonant}}^2 + B_{\text{non-resonant}}^2 + B_0^2}$, where $B_0$ is the magnetic field in the ISM (3 µG in this work; Ferrière 2001).

Two different methods were used to estimate the maximum energy of the protons accelerated at the shock. The first method (which we refer to as Bell) follows the work of Bell et al. (2013) and Cristofari et al. (2020), where the maximum energy of the protons is estimated using the following condition: $\int_{0}^{t'} dt' \gamma_{\text{max}}(t') \approx 5$, where $\gamma_{\text{max}}$ is the maximum non-resonant hybrid growth rate, which leads to

$$ \left( \frac{p_{\text{max}}}{m_{\text{p}}c} \right)^{\alpha-3} = \frac{e R_{\text{sh}} B_{\text{up}}^2}{20 \sqrt{\pi \rho} v_{\text{sh}} m_{\text{p}} c^2} \left( \frac{p_{\text{min}}}{m_{\text{p}}c} \right)^{\alpha-4}. \quad (11) $$

The second method (which we refer to as Hillas) follows C13, where the maximum momentum of the accelerated protons, $p_{\text{max}}$, is determined by the following equation:

$$ l_d = \frac{D(p_{\text{max}})}{u_{\text{sh}}} \approx \zeta R_{\text{sh}}, \quad (12) $$

where D represents the momentum-dependent diffusion coefficient for the CRs upstream of the shock, $l_d$ is the diffusion length of the shock, and $\zeta$ is some fraction; we adopted $\zeta = 0.1$ here. The diffusion coefficient depends on the magnetic field's strength and structure. We assume that the CR diffusion coefficient is of the Bohm type and $D = \frac{R_L c}{3}$, where $R_L = \frac{\text{pc}}{qB}$ is the particle Larmor radius, q is the elementary charge, c is the speed of light, and B is the amplified magnetic field.

We used a diffusion coefficient for the cosmic rays at the shock (Zirakashvili & Ptuskin 2012):

$$ D = D_B \left( 1 + \frac{v_d^2}{u_{\text{sh}}^2} \right)^g, \quad (13) $$

where $D_B$ is the Bohm diffusion coefficient, $v_d$ is the velocity that defines the importance of wave damping in limiting the field of amplification, and the parameter $g$ depends on the nature of the dominant damping mechanism. We used $g = 1.5$ and, following Zirakashvili & Aharonian (2010), adopted the following expression for $v_d$:

$$ v_d = \left( \frac{\sigma^2 B_{\text{up}}^2}{8\pi\xi_{\text{B}}\varrho_{\text{up}}} \right)^{1/2}, \quad (14) $$

where $\varrho_{\text{up}}$ represents the gas density upstream of the shock and $\xi_{\text{B}}$ is the fraction of the shock ram pressure, $\varrho_{\text{up}} u_{\text{sh}}^2$, converted into magnetic field. Völk et al. (2005) found that $\xi_{\text{B}} \approx 3.5\%$ by looking at X-ray data from young SNRs, and we used that value in this work. Plugging all this information into Eq. (12) we can solve for the maximum energy of protons.

### 2.1.5. Dynamics of the SNR shocks

In the description of particle acceleration and of the maximum energy of accelerated particles, the dynamics of the SNR's shocks, radius ($R_{\text{sh}}$), and speed ($u_{\text{sh}}$) plays an important role.

We followed the approach of Ptuskin & Zirakashvili (2003, 2005), where it is assumed that there is very efficient acceleration and the cosmic ray pressure inside the supernova governs its expansion. For thermonuclear SNRs, the shock speed and radius in the ejecta-dominated phase is described by Chevalier (1982), Ptuskin & Zirakashvili (2005)

$$ R_{\text{sh}} = 5.3 \left( \frac{\epsilon_{51}^2}{n_0 M_{\text{ej},\odot}} \right)^{1/7} t_{\text{kyr}}^{4/7} \text{pc} \quad (15) $$

and

$$ u_{\text{sh}} = 3.0 \times 10^3 \left( \frac{\epsilon_{51}^2}{n_0 M_{\text{ej},\odot}} \right)^{1/7} t_{\text{kyr}}^{-3/7} \text{km s}^{-1}, \quad (16) $$

where $\epsilon_{51}$ is the explosion energy in $10^{51}$ erg, $n_0$ is the ambient gas density in cm$^{-3}$, $M_{\text{ej},\odot}$ is the mass of the ejecta in solar-mass units and $t_{\text{kyr}}$ is the time after the explosion in kiloyears.
The shock speed and radius in the adiabatic phase can be calculated with the following equations (Truelove & McKee 1999; Ptuskin & Zirakashvili 2005):

$$R_{\text{sh}} = 4.3 \left(\frac{\epsilon_{51}}{n_0}\right)^{1/5} t_{\text{kyr}}^{2/5} \left(1 - \frac{0.06 M_{\text{ej,}\odot}^{5/6}}{\epsilon_{51}^{1/2} n_0^{1/3} t_{\text{kyr}}}\right)^{2/5} \text{pc} \quad (17)$$

and

$$u_{\text{sh}} = 1.7 \times 10^3 \left(\frac{\epsilon_{51}}{n_0}\right)^{1/5} t_{\text{kyr}}^{-3/5} \left(1 - \frac{0.06 M_{\text{ej,}\odot}^{5/6}}{\epsilon_{51}^{1/2} n_0^{1/3} t_{\text{kyr}}}\right)^{-3/5} \text{km s}^{-1}, \quad (18)$$

respectively. Equations (15) and (16) connect smoothly to Eqs. (17) and (18) at the transfer time $t_{\text{transfer}} \approx 260 (M_{\text{ej,}\odot}^{5/6} \epsilon_{51}^{-1/2} n_0^{-1/3})$ yr.

Core-collapse supernova shocks propagate in the wind-blown bubble generated by the wind of the progenitor star (Das et al. 2022). The wind-blown bubble is divided into two regions: a dense red-super-giant wind and a tenuous hot bubble inflated by the wind of the progenitor star. We assume the gas-density distribution of the wind-blown region is spherically symmetric and that the stellar wind has a velocity $u_w = 10^6$ cm s$^{-1}$. The mass-loss rate is of the order of $\dot{M} = 10^{-5} \dot{M}_{-5}$, where $\dot{M}_{-5}$ is the wind mass-loss rate (Meyer et al. 2021). The density of the gas in the wind is $n_w = \frac{\dot{M}}{4\pi m_a u_w r^2}$, where $u_w$ is the stellar wind velocity, $r$ is the radius of the wind and $m_a = \mu m_p$ is the mean interstellar atom mass ($\mu$ is assumed to be 1.4 here). The density of the hot bubble is (Weaver et al. 1977): $n_b = 0.01 (L_{36}^6 n_0^{19} t_{\text{Myr}}^{-22})^{1/35}$ cm$^{-3}$, where $L_{36}$ is the star wind power in units of $10^{36}$ erg s$^{-1}$ ($L_{36} \sim 1$) and $t_{\text{Myr}}$ is the wind lifetime in units of megayears ($t_{\text{Myr}} \sim 1$) (Longair 1994). The radius of the wind is set by equating the wind pressure to the thermal pressure of the hot bubble: $r_w = \sqrt{\frac{\dot{M} u_w}{4\pi k n_b T_b}}$, where $k$ is the Boltzmann constant and $T_b = 1.6 \times 10^6 (n_0^2 L_{36}^8 t_{\text{Myr}}^{-6})^{1/35}$ K (Castor et al. 1975). The radius of the hot bubble is $r_b = 28 (L_{36}/\mu m_0)^{1/5} t_{\text{Myr}}^{3/5}$ pc.

For core-collapse SNRs, the ejecta-dominated phase is described by the following equations (Chevalier 1982; Ptuskin & Zirakashvili 2005):

$$R_{\text{sh}} = 7.7 \left( \frac{\epsilon_{51}^{7/2} \frac{u_w}{10^6}}{\dot{M}_{-5} M_{\text{ej,}\odot}^{5/2}} \right)^{1/8} t_{\text{kyr}}^{7/8} \text{pc} \quad (19)$$

and

$$u_{\text{sh}} = 6.6 \times 10^3 \left( \frac{\epsilon_{51}^{7/2} \frac{u_w}{10^6}}{\dot{M}_{-5} M_{\text{ej,}\odot}^{5/2}} \right)^{1/8} t_{\text{kyr}}^{-1/8} \text{km s}^{-1}. \quad (20)$$

To describe the evolution of core-collapse SNRs in the adiabatic phase, we adopt the thin-shell approximation, implying that the gas swept up by the shock is concentrated in a thin layer behind the shock front. Using this we can describe the shock speed and the age of the SNR as functions of the shock radius using the following equations (Ptuskin & Zirakashvili 2005):

$$u_{\text{sh}}(R_{\text{sh}}) = \frac{\gamma_{\text{ad}} + 1}{2} \left[ \frac{12(\gamma_{\text{ad}} + 1)\epsilon}{(\gamma_{\text{ad}} - 1)M^2(R_{\text{sh}})R_{\text{sh}}^{6(\gamma_{\text{ad}} - 1)/(\gamma_{\text{ad}} + 1)}} \times \int_{0}^{R_{\text{sh}}} dr \, r^6 (\gamma_{\text{ad}}^{\gamma_{\text{ad}} - 1})^{-1} M(r) \right]^{1/2}, \quad (21)$$

and

$$t(R_{\text{sh}}) = \int_{0}^{R_{\text{sh}}} \frac{dr}{u_{\text{sh}}(r)}, \quad (22)$$

where $\gamma_{\text{ad}}$ is the gas adiabatic index and $M$ is the total gas mass (ejecta and swept up) inside the shock:

$$M(R) = M_{\text{ej}} + 4\pi \int_{0}^{R} dr \, r^2 \rho(r), \quad (23)$$

where $M_{\text{ej}}$ is the ejected mass and $\rho$ is the density of ambient gas.

When the mass of the swept-up gas becomes larger than the mass of the ejecta, the SNR enters the adiabatic (Sedov–Taylor; ST) phase (in this work, we assumed that this is when the mass of the swept up gas is twice the mass of the ejecta). Equations (19) and (20) are multiplied by normalisation constants such that at the transfer time they match the shock radius and shock speed calculated using Eqs. (21) and (22).

This description of the SNR shock holds until the end of the ST phase, which is assumed to be when the acceleration of particle stops being efficient. The end of the ST phase can typically be estimated by equating the cooling time to the age of the SNR (Blondin et al. 1998):

$$t_{\text{cool}} \approx 10^3 \left(\frac{n_0}{1 \text{ cm}^{-3}}\right)^{-1} \left(\frac{v_{\text{sh}}(t)}{10^8 \text{ cm s}^{-1}}\right)^3 \text{kyr.} \quad (24)$$

### 2.1.6. Gamma rays

The gamma-ray emission from accelerated (and reaccelerated) protons and electrons can be computed by making use of the NAIMA Python package (Zabalza 2015). The hadronic emission is calculated using the pion decay model (Kafexhiu et al. 2014), and the proton spectrum is described by a power law with exponential cut-off. The normalisation is calculated using Eq. (2); the cut-off energy is $p_{\text{max}}$, which is calculated using either Eq. (11) or Eqs. (12)–(14); and the power-law index is $\alpha - 2$ since $\alpha$ is the power law in momentum. The leptonic emission is calculated using the inverse Compton model (Khangulyan et al. 2014) of the NAIMA package, with only the cosmic microwave background as a radiation field. The electron spectrum is described by a broken power law with exponential cut-off. This is calculated similarly to the protons except the normalisation is multiplied by the electron-to-proton ratio ($K_{\text{ep}}$) and the break energy is calculated by solving Eq. (8); the power-law index after the break is $\alpha - 1$ and the cut-off energy is calculated with Eq. (6). The total simulated flux for an individual SNR is the sum of the hadronic emission and the leptonic emission from both the accelerated and reaccelerated protons and electrons.

## 2.2. Distribution of Galactic SNRs

### 2.2.1. Distribution of the physical parameters

The simulated SNRs were drawn relying on a Monte Carlo approach: the age of the SNRs was taken following a uniform distribution and we assumed a constant supernova explosion rate of three per century. We accounted for SNRs from both thermonuclear (32%) and core-collapse (68%) supernovae (SNe; Smartt et al. 2009). For SNRs from thermonuclear SNe, the total explosion energy is $10^{51}$ erg and the mass of the ejecta ($M_{\text{ej}}$) is $1.4 M_{\odot}$. For SNRs from core-collapse SNe, the masses of the ejecta depend on the zero-age main-sequence mass of the star,
the explosion energies depend on the mass of the ejecta. The progenitor's last mass-loss rate is assumed to be that of a red supergiant star, that is $\dot{M} = 10^{-5} M_{\odot} yr^{-1}$ for 97% of them and $10^{-4} M_{\odot} yr^{-1}$ for the other 3%.

To obtain the distribution of ejecta masses for core-collapse SNRs, we need the distribution for the initial masses of the massive stars. We use the high-mass regime of the initial mass function, that is the distribution function of how massive stars reach the onset of their main sequence in a given star-formation region (Kroupa 2001). In the regime of interest it reads:

$$dN(M_{\star}) \propto M_{\star}^{-\zeta}, \qquad (25)$$

with $\zeta = 2.3$. Hence, the total number of progenitors obeys:

$$N = \int_{8}^{M_{\max}} dN(M_{\star}) \propto \int M_{\star}^{-\zeta} dM_{\star}, \qquad (26)$$

with $M_{\max} = 120 M_{\odot}$. By normalising Eq. (25) to unity, one can invert Eq. (26) to have a mass for each progenitor $M_{\star}(N_i)$. We used this initial mass distribution for high-mass stars to obtain a progenitor star mass for each of our synthetic SNRs.

Furthermore, the relation between the zero-age main-sequence masses of high-mass stars and their ejecta masses is estimated on the basis of tabulated evolution models from the GENEVA library<fnref n="1" /> (Eggenberger et al. 2008; Ekström et al. 2012). This code calculates the stellar structure in one dimension including convection and diffusion in the stellar interior, determined with the Schwarzschild criterion for the mixing length theory and the Zahn recipes for diffusion coefficient to the viscosity caused by horizontal turbulence and for the transport of angular momentum. The physics of the mass-loss rate is described in the work of Georgy (2010). There are similar libraries of stellar evolution models produced by other evolutionary codes, such as the BOOST database (Szécsi et al. 2022), which is produced by the Bonn stellar evolution code (Brott et al. 2011a,b), or the MIST library (Choi et al. 2016; Dotter 2016), which is calculated using the MESA code (Paxton et al. 2015, 2018, 2019). Such multi-purpose models are widely used in many aspects of stellar astrophysics such as population synthesis (Fragos et al. 2023), hydrodynamical simulations (Meyer et al. 2022, 2023; Meyer & Meliani 2022), and particle acceleration simulations (Das et al. 2022).

We used the SYCLIST interface to the GENEVA library to generate stellar-evolution models for non-rotating massive stars at solar metallicity ($Z = 0.014$), with a zero-age main-sequence mass $M_{\star}$ spanning from 10 to 120 $M_{\odot}$. The initial mass range is discretised by a representative number of non-rotating stars covering it with a mass increment of 5 $M_{\odot}$. The ejecta mass $M_{ej}$ of these core-collapse supernova progenitors is estimated by subtracting the mass of a neutron star from the pre-supernova mass of the star as calculated by GENEVA. It reads

$$M_{ej} = M_{\star} - \int_{t_{ZAMS}}^{t_{SN}} \dot{M}(t) dt - M_{NS}, \qquad (27)$$

where $t_{ZAMS}$ and $t_{SN}$ are the zero-age main-sequence and supernova times, respectively; $\dot{M}$ is the mass-loss rate of the star at a given time, $t_{ZAMS} \le t \le t_{SN}$ of its evolution; and $M_{NS} = 1.4 M_{\odot}$ is the mass of the neutron star produced after the death of the massive progenitor. The obtained zero-age mass-sequence versus ejecta-mass law is shown in Fig. 1 (top), and it is

<figure><img src="image_2.png" /><figcaption>Fig. 1. Relation between zero-age main-sequence mass and ejecta mass of non-rotating core-collapse supernova progenitors at solar metallicity in the GENEVA library (Ekström et al. 2012) and the relation between the zero-age main-sequence mass and the explosion energy (Sukhbold et al. 2016). The blue dots show the calculated points, while the orange lines show the interpolation between them.</figcaption></figure>

implemented in our model via a routine interpolating a table containing the data relative to these $M_{\star}-M_{ej}$ relations.

As the explosion energy $E_{ej}$ is not provided by stellar evolutionary calculations but requires additional calculations taking them as initial conditions, we used the results of Sukhbold et al. (2016); they provide an $M_{\star}-E_{ej}$ law in Table 5, which is referred to as model N20K based on neutrino-powered explosions and makes it possible to interpolate the explosion energy from the progenitor mass. These data are displayed in Fig. 1 (bottom).

Using these ejecta-mass and explosion energy distributions described above, as well as the uniform age distribution, the SNRs are placed throughout the Milky Way according to one of the following papers: Faucher-Giguère & Kaspi (2006), Steiman-Cameron et al. (2010), Green (2015), or Reid et al. (2019). The source distributions will be called CAFG, Steiman-Cameron, Green, and Reid from here on and are described in more detail in Sect. 2.2.2. The Steiman-Cameron and Green models are both implemented following the work by Steppa & Egberts (2020). In this work, we used a galactocentric coordinate system with the Sun located at $x = 0, y = 8.5$ kpc, $z = 0$.

### 2.2.2. Spatial distribution of SNRs

**The CAFG model.** The four-spiral-arm model described in Faucher-Giguère & Kaspi (2006) is used to model pulsar positions. We use it by imposing that the birth locations of the supernovae follow the radial distribution described for the pulsars and that the arm be chosen at random. The locations of the supernovae also follow a $z$ distribution, with thermonuclear supernovae following the HII gas density and core-collapse supernovae following the HII gas density. The radial distribution of sources uses the following probability density function:

$$\rho(r) = A \left(\frac{r + R_1}{R_{\odot} + R_1}\right)^a \exp\left[-b \left(\frac{r - R_{\odot}}{R_{\odot} + R_1}\right)\right], \qquad (28)$$

where $A$ is a normalisation constant, $r$ is the galactocentric radius, $R_1 = 0.55$ kpc, $a = 1.64$, and $b = 0.24$.

<footnote n="1">¹ https://www.unige.ch/sciences/astro/evolution/en/database/syclist/</footnote>
<table><caption>**Table 1.** Parameter values for Steiman-Cameron arms.</caption><thead><tr><th rowspan="2">Spiral arm</th><th rowspan="2">β<sub>i</sub></th><th rowspan="2">a<sub>i</sub></th><th rowspan="2">R (kpc)</th><th colspan="2">σ<sub>r</sub> (kpc)</th><th rowspan="2">σ<sub>z,2</sub> (kpc)</th><th rowspan="2">δ (deg)</th><th rowspan="2">A<sub>i</sub></th></tr><tr><th>(r &lt; R)</th><th>(r &gt; R)</th></tr></thead><tbody><tr><td>Sagittarius-Carina</td><td>0.242</td><td>0.246</td><td>2.9</td><td>0.7</td><td>3.1</td><td>0.070</td><td>15</td><td>169</td></tr><tr><td>Scutum-Crux</td><td>0.279</td><td>0.608</td><td>2.9</td><td>0.7</td><td>3.1</td><td>0.070</td><td>15</td><td>266</td></tr><tr><td>Perseus</td><td>0.249</td><td>0.449</td><td>2.9</td><td>0.7</td><td>3.1</td><td>0.070</td><td>15</td><td>339</td></tr><tr><td>Norma-Cygnus</td><td>0.240</td><td>0.378</td><td>2.9</td><td>0.7</td><td>3.1</td><td>0.070</td><td>15</td><td>176</td></tr></tbody></table>

<table><caption>**Table 2.** Parameter values for arms in the Reid distribution of sources.</caption><thead><tr><th rowspan="2">Spiral arm</th><th rowspan="2">β range (deg)</th><th rowspan="2">β<sub>kink</sub> (deg)</th><th rowspan="2">R<sub>kink</sub> (kpc)</th><th colspan="2">ψ (deg)</th></tr><tr><th>(β &lt; β<sub>kink</sub>)</th><th>(β &gt; β<sub>kink</sub>)</th></tr></thead><tbody><tr><td>Sagittarius-Carina</td><td>2 → 97</td><td>24</td><td>6.04</td><td>17.1</td><td>1.09</td></tr><tr><td>Scutum-Centaurus</td><td>0 → 104</td><td>23</td><td>4.91</td><td>14.1</td><td>12.1</td></tr><tr><td>Perseus</td><td>-23 → 115</td><td>40</td><td>8.87</td><td>10.3</td><td>8.7</td></tr><tr><td>Norma</td><td>5 → 54</td><td>18</td><td>4.46</td><td>-1.0</td><td>19.5</td></tr><tr><td>Local</td><td>-8 → 34</td><td>9</td><td>8.26</td><td>11.4</td><td>11.4</td></tr><tr><td>3 kpc</td><td>15 → 18</td><td>15</td><td>3.52</td><td>-4.2</td><td>-4.2</td></tr><tr><td>Outer</td><td>-16 → 71</td><td>18</td><td>12.24</td><td>3.0</td><td>9.4</td></tr></tbody></table>

**The Steiman-Cameron model.** A four-spiral-arm model based on the interstellar medium is described in Steiman-Cameron et al. (2010). The model matches the data obtained by the FIRAS instrument of the Cosmic Background Explorer for the [CII] 158 µm and [NII] 205 µm lines of the interstellar medium. The distribution of sources is described by the following probability density equation:

$$ \rho(r, \phi, z) = \sum_{i=1}^{4} A_i \exp \left( -\frac{1}{\delta^2} \left( \phi - \frac{\ln\left(\frac{r}{a_i}\right)}{\beta_i} \right)^2 \right) \times \exp \left( -\frac{|r-R|}{\sigma_r} \right) \exp \left( -\frac{z^2}{2\sigma_{z,2}^2} \right), \quad (29) $$

where $\sigma_r$ is the scale length of the radius, $R$ is the local maximum radius, $\delta$ is the azimuthal scale angle, $\beta_i$ is the pitch angle of the arm, $a_i$ is the orientation of the arm and $\sigma_{z,2}$ is the scale height. The parameters for the arms can be seen in Table 1.

**The Green model.** An azimuthally symmetric model based on SNR data is described in Green (2015). The spatial distribution we use is described by the following probability density function:

$$ \rho(r, z) = \left(\frac{r + r_{\text{off}}}{R_{\odot} + r_{\text{off}}}\right)^{\alpha} \exp\left(-\beta \frac{r - R_{\odot}}{R_{\odot} + r_{\text{off}}}\right) \exp\left(-\frac{|z|}{z_0}\right), \quad (30) $$

with $R_{\odot}$ being the distance of the Sun to the Galactic centre, $z_0$ the scale height of the Galactic disc, $\alpha$ the shape parameter, and $\beta$ the rate parameter. The parameter $r_{\text{off}}$ accounts for a non-zero density at $r = 0$. Following Steppa & Egberts (2020), we set $z_0 = 0.083$ kpc, $r_{\text{off}} = 0$ kpc, $\alpha = 1.09$, and $\beta = 3.87$.

**The Reid model.** A four-spiral-arm model based on data from massive stars with maser parallaxes is described in Reid et al. (2019). This model also includes some additional features: the Local Arm; a 3 kpc arm around the centre of the Galaxy; and an outer arm, which connects with the Norma Arm. The Local Arm is described as an arm segment and does not link

with any of the other arms. The 3 kpc arm is associated with the Galactic bar and might not be a true spiral arm. The spiral arms are described by Eq. (31) and Table 2:

$$ \ln\left(\frac{R}{R_{\text{kink}}}\right) = -(\beta - \beta_{\text{kink}}) \tan\psi, \quad (31) $$

where $R$ is the galactocentric radius at azimuth $\beta$, $R_{\text{kink}}$ is the radius of the kink in the arm at azimuth $\beta_{\text{kink}}$, and $\psi$ is the pitch angle.

### 2.3. The gas distribution

For the matter distribution we need to know the hydrogen density of the region we place our sources into. In this work, we used the matter distribution described by Shibata et al. (2010), which includes an empirical model for the matter distribution. They describe the gas density of interstellar hydrogen (HI and H$_2$) as a function of galactocentric distance and height above or below the Galactic plane. This empirical model closely matches the one used in the GALPROP code (Strong & Moskalenko 1998).

## 3. Confronting our simulated populations with available H.E.S.S. data

Following the modelling described in Sect. 2, we can simulate the Galactic population of SNRs. The typical parameters of interest are:

- $K_{\text{ep}}$: the electron-to-proton ratio.

- $\alpha$: the spectral index of the accelerated particles.

- $T_{\text{ST}}$: the duration of the ST phase.

- $\eta$: the CR acceleration efficiency.

- The spatial distribution of sources.

We typically compute 100 realisations for a given set of parameters to ensure the stability of our results. We note that the positions of simulated SNRs change for each of the 100 realisations for a given set of properties, but these 100 different position sets are constant when changing the population properties other than
<figure><img src="image_3.png" /><figcaption>Fig. 2. Example of one realisation of a simulation of the population of Galactic SNRs. The grey shaded region is the HGPS detectability range for point-like sources with a luminosity of $5 \times 10^{33}$ photons s$^{-1}$ (equivalent to a differential flux at 1 TeV of $\sim 4 \times 10^{-11}$ TeV$^{-1}$ cm$^{-2}$ s$^{-1}$ at a distance of 1 kpc). The blue plus symbols mark the simulated sources that are detectable, and the black dots mark the simulated sources that are not detectable. 3.7% of the emitting SNRs are detectable in this realisation. The location of the Sun is marked by the red cross.</figcaption></figure>

the source distribution model. This allows us to closely examine how each of the properties affects the population.

For each source, using the gamma-ray luminosity, size, and position we are able to determine whether it would have been detected in the HGPS. To account for the selection bias in the H.E.S.S. catalogue we followed the work in Steppa & Egberts (2020); in particular, we accounted for the fact that the HGPS sensitivity varies with position and source extent. For a simulated source to be detectable it must exceed the threshold of $5\sigma$ above the background. This condition reads

$$F_{\min}(\theta_{\text{source}}) = \begin{cases} F_{\min,0} \sqrt{\frac{\theta_{\text{source}}^2 + \theta_{\text{PSF}}^2}{\theta_{\text{PSF}}^2}}, & \theta_{\text{source}} \le 1^{\circ}, \\ \infty, & \theta_{\text{source}} > 1^{\circ} \end{cases} \quad (32)$$

where $F_{\min,0}$ is the point-source sensitivity, $\theta_{\text{source}}$ is the source extent (radius) and $\theta_{\text{PSF}}$ is the size of the H.E.S.S. point-spread function. The limited field of view of H.E.S.S. combined with the applied technique of deriving background measurements means that sources $\ge 1^{\circ}$ are not detectable. It is worth noting that the source in the HGPS with the largest extent, Vela Junior, was the subject of a dedicated analysis rather than the analysis method used in the HGPS. Vela Junior has an extent of $1^{\circ}$; the largest extent of any source analysed using the method in the HGPS is Vela X, with an extent of $0.58^{\circ}$. Another well-known SNR, RXJ1713.7-3946 was also the subject of a dedicated analysis and has an extent of $0.5^{\circ}$. In Fig. 2, the HGPS detectability range for point-like sources with a luminosity of $5 \times 10^{33}$ photons s$^{-1}$ is displayed, demonstrating the large variations in sensitivity as a function of Galactic longitude. This variation in sensitivity shows why it is not sufficient to have a cut-off in the integrated flux to determine if a source is detectable, and thus why we need to include the sensitivity of H.E.S.S. when testing each simulated source.

Because many of the HGPS sources are unidentified and at least some of those should be SNRs, we create an upper and lower limit for the HGPS sources we compare to. The lower limit

is the firmly identified SNRs in the HGPS (8), and the upper limit includes all the unidentified sources (47) and the sources identified as composite SNRs (8); this gives us an upper limit of 63 sources to compare to. Moreover, only 12 of the unidentified sources have been found with a position compatible with a known SNR: this allows us to consider a more stringent upper limit of 28 sources. Of the eight firmly detected SNRs, which have been modelled individually (H.E.S.S. Collaboration 2011, 2015, 2018b,c,d,e; Nguyen-Dang et al. 2023; Cui et al. 2018), four of them are clearly found with $E_{\max} > 10$ TeV, two of them with $E_{\max} \sim 10$ TeV and, two of them with $E_{\max} < 10$ TeV. We use this as an additional constraint on the populations.

## 4. Results

We compare our synthetic populations to the HGPS sample. For illustrative purposes, we define a reference set of parameters as follows: the Steiman-Cameron source distribution, the Bell method to estimate $p_{\max}$, $\alpha = 4.3$, $K_{\text{ep}} = 10^{-3}$, $\eta = 7\%$, and $T_{\text{ST}} = 20$ kyr. It is important to note that this particular set of properties is not the only set that can describe the data, but it is a useful starting point for making comparisons.

The effect of varying $\alpha$ from 4.0 to 4.4 is shown in Fig. 3. In this case, only one of the sets of parameters (with a spectral index of 4.3) is found in agreement with the data. With a spectral index of 4.3, almost all (96%) the populations fit between the eight firmly detected SNRs (blue line) and the stringent upper limit (28), which only includes sources associated with SNRs (green line). All the populations for this parameter set fit between the eight firmly detected SNRs (blue line) and the stringent upper limit, which includes all the unidentified sources (orange line). For $\alpha > 4.3$, there are too few detectable SNRs, while $\alpha < 4.3$ overproduces the H.E.S.S. results. The log N-log S distribution of fluxes for the simulated SNRs has a very similar shape to that of the H.E.S.S. SNRs. Harder spectra ($\alpha$ close to 4) increase the number of detectable SNRs, which is expected since harder spectra correspond to an increase in the flux in the TeV range, thus increasing the horizon of detectability. The total number of SNRs remains the same for all the values of the spectral index (right panel Fig. 3).

Changing the efficiency, $\eta$, from 1% to 10% directly impacts the number of detectable SNRs, as shown in Fig. 4. With our benchmark case, we find that $\eta < 5\%$ and $\eta > 9\%$ seem to under- or over-produce the H.E.S.S. data, respectively.

Exploring values of $K_{\text{ep}}$ in the $10^{-2}$ to $10^{-5}$ range, we illustrate in Fig. 5 that a value of $10^{-3.5}$ leads to results compatible within the firmly detected H.E.S.S. SNRs and our stringent upper limit. With the initial benchmark set of parameters we considered, a value of $K_{\text{ep}} = 10^{-2.5}$ leads to a simulated population of SNRs close to the upper limit of all the unknown H.E.S.S. sources; thus, this is a possible but extreme and quite unlikely case.

Figure 6 shows the mean hadronic ratio (ratio of the number of SNRs with gamma-ray luminosity dominated by the hadronic component to the total number of SNRs) against the integrated flux (>1 TeV) of the detectable SNRs using the same cumulative binning as in Figs. 3–5 and its strong dependence on the parameters $K_{\text{ep}}$, $\eta$, and $\alpha$. The dependence on $K_{\text{ep}}$ is easily understood since this parameter directly affects the normalisation of the electron spectrum, and thus an increase in $K_{\text{ep}}$ leads to a decrease of the hadronic ratio (increase in the number of leptonic-dominated SNRs), which can be seen in Fig. 6 (top left). At higher integrated fluxes, we have fewer detectable SNRs. For $K_{\text{ep}} = 10^{-4}$ and $K_{\text{ep}} = 10^{-5}$ we have, on average, <1 of the brightest SNRs. As seen in Fig. 5, this brightest detectable
<figure><img src="image_4.png" /><figcaption>**Fig. 3.** log N–log S distribution of detectable SNRs (left) and all emitting SNRs (right) based on 100 Monte Carlo realisations following a Steiman-Cameron source distribution, $E_{max}$ estimation considering Bell instabilities, $K_{ep} = 10^{-3.0}$, and $\eta = 7\%$. The spectral index of the accelerated particles $\alpha$ varies from 4 to 4.4. Shaded areas indicate $\pm$ one standard deviation. Thick solid lines correspond to the HGPS limits (blue: identified SNRs; green: identified and unidentified sources associated with SNRs; orange: all unidentified sources). The hatched region shows the ideal region for the simulated SNRs: between the identified SNRs and the sources associated with SNRs.</figcaption></figure>

SNR is always dominated by hadronic emission, and this is why there is no statistical error at the edge of Fig. 6. An increase in $\eta$, a decrease in $\alpha$, or both also leads to a decrease of the hadronic ratio (increase of the number of leptonic dominated SNRs). This can be understood as the result of several combined effects: by decreasing $\alpha$, increasing $\eta$, or both, in a given horizon, we allow for a part of the SNR population (older SNRs, or those located in less dense parts) to become detectable, and thus we unveil a part of the SNR population that is typically more leptonic dominated (see Table 3), on top of the already detected hadronic SNRs. This leads to a decrease in the hadronic ratio.

We additionally investigate the influence of the duration of the active acceleration phase, treating the duration of the ST phase $T_{ST}$ as a free parameter that we are artificially changing. Changes within a factor of a few typically do not affect our results substantially, as shown in Fig. 7. This is somewhat expected since mostly young $\le 5–10$ kyr old SNRs account for the brightest TeV SNRs.

Diffusive shock reacceleration, as discussed in Sect. 2.1.2, is found to be especially relevant at old (≥10 kyr) SNR shocks (Bell 1978; Blasi 2004, 2017) and for small $\alpha$; this does not substantially impact our results in the TeV range.

The different spatial distributions of SNRs are compared in Fig. 8. The importance of the source distribution is highlighted in the fact that the differences in the number of the detectable sub-sample of sources between the spatial distributions (left panel) is larger than when considering all emitting sources (right panel). At the same time, this shows that a thorough discussion on the spatial model requires taking into account the refined description of the sensitivity of the instrument, including its spatial dependence, as done in this work. Interestingly, the Steiman-Cameron and the Green distributions produce very similar detectable simulated SNRs despite having very different source distributions.

All spatial distributions considered can account for the HGPS data when adjusting the other parameters of the SNR population.

**Systematic exploration of the parameter space.** In addition to varying each parameter at a time for the benchmark case discussed above, we performed a systematic exploration of the parameter space with $4.0 \le \alpha \le 4.4$, $10^{-5} \le K_{ep} \le 10^{-2}$ and $0.01 \le \eta \le 0.1$, considering $E_{max}$ from the Bell and Hillas estimates. For these simulated properties, the maximum length of the ST phase was set to 20 kyr, and the Steiman-Cameron source distribution were used.

Exploring the parameter space, we constrained the results of our simulations by the HGPS catalogue, especially investigating which parameters can lead to realisations compatible with the detected Galactic population.

Moreover, we check that the maximum energy of accelerated particles, $E_{max}$, at the simulated SNRs is compatible with the current $E_{max}$ estimated from the measured gamma-ray differential spectrum at each Galactic SNR. Schematically, out of the eight firmly detected SNRs, four of them are clearly found with $E_{max} > 10$ TeV (HESS J0852–463, HESS J1442–624, HESS J1713.7–3946 and HESS J1731–347), two of them with $E_{max} \sim 10$ TeV (HESS J1534–571 and HESS J1718–347). We should mention that the cases of HESS J1801–233 (W28; Cui et al. 2018) and HESS J1911+090 (W49B; H.E.S.S. Collaboration 2018d) are rather tricky, as these two sources are thought to be SNRs interacting with molecular clouds, so it is not easy to ascribe the gamma-ray emission to the SNR shock, to clouds, or to a combination of the two. In addition, the question of whether there is an actual exponential suppression in the gamma-ray flux is yet not clearly answered. It remains that the gamma-ray luminosity is steeply decreasing in the ~1–10 TeV range, and we considered that $E_{max} < 10$ TeV in this work.
<figure><img src="image_5.png" /><figcaption>Fig. 4. As in Fig. 3, fixing $\alpha = 4.3$ and varying $\eta$ in the 1%–10% range.</figcaption></figure>

<figure><img src="image_6.png" /><figcaption>Fig. 5. As in Fig. 3, fixing $\alpha = 4.3$ and varying the electron-proton ratio in the 10$^{-2}$–10$^{-5}$ range.</figcaption></figure>

For all realisations of our simulations, we checked the num-
ber of realisations for which (1) results are found between
the firmly detected H.E.S.S. SNRs (8) and all the detected
H.E.S.S. sources associated with SNRs (28); and (2) the $E_{\max}$
of the simulated SNRs that are detectable in our sample is
compatible with the number of SNRs with $E_{\max} > 10$ TeV in
the HGPS (≥4). This is illustrated in Fig. 9, where, assum-
ing an $E_{\max}$ following Bell, one realisation of a simulated
population is shown in the parameter space of the SNR age
and $E_{\max}$. The population parameters are chosen such that
they lead to the best agreement with HGPS data ($\alpha = 4.2$,

$K_{ep} = 10^{-5}$, $\eta = 9\%$, with ~97% of the realisations compat-
ible with data). Interestingly, in this situation, the number of
SNR PeVatrons (accelerating particles up to the PeV range) is
close to zero, which would be compatible with the idea that
PeV CRs are not accelerated at typical SNRs (Cristofari et al.
2020). The prescription of Hillas naturally leads to possible
higher values of $E_{\max}$, but it does not substantially increase
the chances of detecting a SNR PeVatron; this can be seen in
Fig. 10.

The results are illustrated in Figs. 11 and 12. In Fig. 11, the
results are displayed assuming $E_{\max}$ follows the Bell description.
<figure><img src="image_7.png" /></figure>

<figure><img src="image_8.png" /></figure>

<figure><img src="image_9.png" /><figcaption>**Fig. 6.** Mean ratio of hadronic-dominated detectable SNRs versus integrated flux, as in Fig. 3. Top left panel illustrates the importance of $K_{ep}$ (fixing $\alpha = 4.3$ and $\eta = 7\%$); top right panel illustrates the effect of varying $\eta$ (fixing $K_{ep} = 10^{-3}$ and $\alpha = 4.3$); and bottom panel illustrates the effect of varying $\alpha$ (fixing $\eta = 7\%$ and $K_{ep} = 10^{-3}$). The rest of the properties are as in Fig. 3.</figcaption></figure>

In Fig. 12, we show the results with $E_{max}$ estimated using the Hillas criterion.

It follows from Fig. 11 that one can rule out the following parameters: $\alpha \ge 4.35$ and $K_{ep} > 10^{-3}$. One can only find a set of parameters for which >80% of the realisations are within the associated SNR limits if $4.1 \le \alpha < 4.3$, $K_{ep} < 10^{-3.5}$, and $\eta > 0.02$. However, Fig. 12 illustrates that when the $E_{max}$ is estimated with the Hillas criterion, if $\alpha < 4.15$, or if $K_{ep} > 10^{-3}$ no region of the parameter space can account for the HGPS data. The lower limit placed on $\alpha$ is a direct result of including the emission from reaccelerated particles in our SNR model. This emission is larger at smaller $\alpha$ and causes the populations with those parameters to have more detectable SNRs than the HGPS. In both plots, the correlations between parameters can be seen: an increase in $\eta$ is compensated by either a decrease in $K_{ep}$ or larger values of $\alpha$. This intertwining makes it arduous to identify a preferred region of the parameter space.

The systematic exploration of the parameter space allows us to identify the set of parameters that lead to a maximum number of realisations in agreement with the HGPS data (see for example Figs. 11 and 12), but also to examine the properties of these synthetic populations. Assuming an $E_{max}$ following Bell, the typical characteristics of the simulated population for which a maximum number of Monte Carlo realisations are in agreement with the HGPS are shown in Table 3.

A few remarks are in order. The set of parameters leading to the highest fraction of realisations in agreement with the HGPS is $\alpha = 4.2$, $K_{ep} = 10^{-5}$, $\eta = 0.09$. For all populations leading to 90% of compatibility with the HGPS, we found that $4.2 \ge \alpha \ge 4.1$ and $10^{-4.5} \ge K_{ep} \ge 10^{-5}$ for a CR efficiency $0.1 \ge \eta \ge 0.03$. Interestingly, the electron-to-proton ratio $K_{ep}$ is found to be on the lower side of typically inferred values, and at the same time, typically about 40%-60% of the gamma-ray emission of the simulated SNRs detectable in the survey is dominated by leptonic interactions (see hadronic ratio). This is somehow counterintuitive, as the idea that the VHE gamma-ray emission of SNRs is leptonic usually implies that the acceleration of leptons is efficient. In our case, low values of $K_{ep}$ are compensated by sufficiently hard spectra and efficiencies, to allow for leptonic-dominated gamma rays to be detectable. The leptonic-dominated SNRs are typically older than the hadronic-dominated SNRs (the mean age of leptonic ones is typically twice the mean age of hadronic ones), which is consistent with the consensual picture (Pohl 1996; Reynolds 2008; Acero et al. 2015), for instance because of the reduced efficiency of magnetic-field amplification at slower shocks, thus reducing losses suffered by electrons. As illustrated in Fig. 6, the efficiency of particle acceleration is indeed found to strongly impact the fraction of leptonic sources. Indeed, the spectra of both accelerated protons and electrons are directly proportional to the acceleration efficiency, but an increase of the efficiency allows for the dimmer leptonic-dominated sources to reach the sensitivity of the HGPS, and thus unveils the leptonic-dominated SNRs. With $\alpha = 4.3$, and $K_{ep} = 10^{-3}$, varying the CR efficiency from 0.01 to 0.1 changes the fraction of hadronic sources from ≈1 to ≈0.2. It should be noted that this outcome was achieved assuming that the maximum energy at SNRs is determined by the development of non-resonant streaming instabilities, and ensuring that the count of simulated SNRs detected above approximately 10 TeV aligns with the HGPS data. Considering a maximum energy determined by the Hillas criterion, or relaxing the constraint
<figure><img src="image_10.png" /><figcaption>**Fig. 7.** As in Fig. 3, fixing $\alpha = 4.3$ and varying $T_{\text{ST}}$ from 5 kyr to 100 kyr.</figcaption></figure>

<figure><img src="image_11.png" /><figcaption>**Fig. 8.** As in Fig. 3, fixing $\alpha = 4.3$ and varying the source distribution models: Stemian-Cameron, CAFG, Green, and Reid.</figcaption></figure>

on the number of SNRs detected above about 10 TeV, opens up a wider range of possibilities where parameter combinations with an $\alpha$ value around 4.4 or a $K_{\text{ep}}$ value around $10^{-2}$ can be identified. Both population (Ellison et al. 2007; Cristofari et al. 2013) and theoretical studies (Corso et al. 2023) have not yet led to stringent constraints for the values of $K_{\text{ep}}$ or the ratio of hadronic-dominated sources. Mandelartz & Becker Tjus (2015) interpreted 21 out of 24 SNRs in the GeV range as hadronic dominated, but the work of Zeng et al. (2019) led to more complex discussions.

**5. Conclusion**

Relying on a Monte Carlo approach, we simulated populations of the Galactic SNRs in the TeV range (built on a physically motivated model of the gamma-ray emission from accelerated protons and electrons at SNR shocks) and confronted the results of the simulated populations to available data from the systematic survey of the Galactic plane performed by H.E.S.S.

We performed a systematic exploration of the parameter space defined by $K_{\text{ep}}$, which is the electron-to-proton ratio;
<figure><img src="image_12.png" /><figcaption>**Fig. 9.** Maximum energy of accelerated protons versus age of each simulated SNR for a given Monte Carlo realisation. Parameters adopted are $\alpha = 4.2$, $K_{ep} = 10^{-5}$, $\eta = 9\%$, and with $E_{max}$ following the Bell prescription. The horizontal line at 10 TeV illustrates the criterion required to account for the HGPS data.</figcaption></figure>

<figure><img src="image_13.png" /><figcaption>**Fig. 10.** As in Fig. 9, but with $E_{max}$ following the Hillas prescription.</figcaption></figure>

$\alpha$, the power-law spectral index of the accelerated particles at the SNR shocks; $T_{ST}$, the duration of the ST phase; and $\eta$, the CR acceleration efficiency. We did this considering different prescriptions for the description of the maximum energy of accelerated particles $E_{max}$ and different prescriptions for the spatial distribution of sources. Our conclusions are listed below.

1. We find that the Galactic population of TeV SNRs can be explained by multiple regions of the parameter space. Given the current population in the HGPS, only some parts of the parameter space seem to be excluded. For instance, $\alpha \lesssim 4.05$ seems to be disfavoured in all cases. $K_{ep} \gtrsim 10^{-2.5}$, is excluded and $K_{ep} \sim 10^{-3}$ requires $\alpha \gtrsim 4.35$ and $\eta \lesssim 0.02$.

2. Despite the limitations of our approach, the systematic exploration of the parameter space allows us to look for the region of the parameter space producing the most realisations

within the limits of the HGPS sample. In addition to the requirement concerning the total number of sources, the additional constraint on the maximum energy (ensuring that $E_{max}$ at SNRs is compatible with values estimated at individual SNRs) allows us to further reduce the allowed parameter space. One possible solution rendering ~97% populations in the limits of the HGPS is $\alpha = 4.2$, $K_{ep} = 10^{-5}$, $\eta = 9\%$, assuming $E_{max}$ from the Bell prescription. In this case, the mean ratio of hadronic-to-leptonic-dominated SNRs is 0.62, which is compatible with the current idea that the TeV emission can be leptonic or hadronic for the different known SNRs. We note that the relatively high efficiency $\eta = 9\%$ allows for a non-negligible fraction of leptonic SNRs, which are older or located in low-density, regions to be detectable.

3. The sets of parameters that have $\ge 90\%$ realisations compatible with the HGPS are obtained with $4.1 \le \alpha \le 4.2$, $10^{-5} \le K_{ep} \le 10^{-4.5}$ and $0.03 \le \eta \le 0.1$. Interestingly, this leads to a fraction of $\approx 40-60\%$ of leptonic-dominated SNRs in the TeV range, along with a reduced fraction of electrons injected at the SNR shocks. This means that our study naturally highlights situations in which we can account for a significant (~50%) fraction of the sources being dominated by the leptonic emission, with a low value of $K_{ep}$ (~$10^{-5}$). For most SNRs, the dominant processes at play (leptonic vs. hadronic) are still debated, and thus the overall situation for the Galactic population is also unknown (see, for example discussion in Zeng et al. 2019). In fact, even for what is probably the best-studied SNR of the TeV range RX J1713-3946, numerous works have still not managed to conclude on the dominant mechanism at play (Ellison et al. 2010; Zirakashvili & Aharonian 2010; Morlino & Caprioli 2012; Gabici 2017; Ohira & Yamazaki 2017; Cristofari et al. 2021b; Fujita et al. 2022). In other words, the simulated populations are thus directly compatible with a situation where half of the TeV-detected SNRs are leptonic dominated, and at the same time the fraction of injected electrons at SNR shocks is substantially lower than the ~1/100 (~1/50) electron-to-proton ratio measured in the local CR spectrum approximately in the giga-electron-volt range. This would thus be an additional argument in favour of SNRs not being dominant sources of CR electrons (Gabici et al. 2019; Lipari 2019), and at the same time, that a substantial fraction of SNRs are detected in the TeV range with a dominant leptonic component. We note that this result was obtained under the assumption that the maximum energy at SNRs is set by the growth of non-resonant streaming instabilities, and ensuring that the number of simulated SNRs detected above $\ge 10$ TeV is compatible with the HGPS data. Considering a maximum energy set by the Hillas criterion or relaxing the condition on the number of SNRs detected above $\ge 10$ TeV leads to more diverse possibilities where sets of parameters with $\alpha \sim 4.4$ or $K_{ep} \sim 10^{-2}$ can be found.

4. Some parameters do not substantially affect our results. For $T_{ST}$ this is coherent with the idea that most TeV-emitting SNRs are young ($\lesssim$ a few thousand years). The various prescriptions for the spatial distribution of sources do not substantially affect our results.

Any additional constraint on one of the parameters explored in our study would naturally narrow down the parameter space compatible with the HGPS. Moreover, an increase in the sample of the firm SNR detections (8, our lower limit) would help our approach. This increased number of detections is clearly
<figure><img src="image_14.png" /><figcaption>**Fig. 11.** Two-dimensional histogram plots showing the percentage of populations that have at least as many detectable simulated SNRs as firmly detected H.E.S.S. SNRs, do not exceed the number of H.E.S.S. sources associated with SNRs, and they have at least as many detectable SNRs with $E_{max} \ge 10$ TeV as HGPS SNRs. The source distribution is Steiman-Cameron, and $E_{max}$ follows the Bell description.</figcaption></figure>

<figure><img src="image_15.png" /><figcaption>**Fig. 12.** As in Fig. 11, but with $E_{max}$ estimated with the Hillas criterion.</figcaption></figure>

expected with next-generation instruments such as the CTAO (Cristofari et al. 2017; Cherenkov Telescope Array Consortium 2019; CTA Consortium 2023) or LHAASO (Liu et al. 2016). Cristofari et al. (2017), CTA Consortium (2023) discussed that CTAO could detect up to several hundreds of SNRs in the systematic exploration of the sky.

Additionally, joint observations of SNRs with next-generation radio observatories (Jonas & MeerKAT Team 2016; Tingay et al. 2013; Dewdney et al. 2009) will provide data of unprecedented value concerning the content of non-thermal electrons accelerated at SNR shocks and help disentangle the hadronic and leptonic content, thus allowing us to better understand the details of DSA at SNRs.
<table><caption>**Table 3.** Typical characteristics of detectable SNRs in the simulated populations.</caption><thead><tr><th></th><th>Population parameters</th><th>% of realisations compatible with HGPS</th><th>Hadronic ratio</th><th>Mean no. detectable SNRs</th><th>Mean Had. age (kyr)</th><th>Mean Lep. age (kyr)</th><th>Mean Had. dist. (kpc)</th><th>Mean Lep. dist. (kpc)</th></tr></thead><tbody><tr><td>α = 4.2</td><td>K<sub>ep</sub> = 10<sup>-5.0</sup></td><td>η = 0.09</td><td>97.0</td><td>0.62</td><td>16.84</td><td>2.15</td><td>4.86</td><td>5.65</td><td>4.88</td></tr><tr><td>α = 4.2</td><td>K<sub>ep</sub> = 10<sup>-4.5</sup></td><td>η = 0.07</td><td>96.0</td><td>0.43</td><td>16.14</td><td>1.94</td><td>4.36</td><td>5.64</td><td>4.9</td></tr><tr><td>α = 4.15</td><td>K<sub>ep</sub> = 10<sup>-5.0</sup></td><td>η = 0.05</td><td>95.0</td><td>0.51</td><td>16.41</td><td>2.06</td><td>5.21</td><td>5.62</td><td>4.79</td></tr><tr><td>α = 4.2</td><td>K<sub>ep</sub> = 10<sup>-5.0</sup></td><td>η = 0.08</td><td>93.0</td><td>0.66</td><td>13.6</td><td>2.0</td><td>4.88</td><td>5.63</td><td>5.06</td></tr><tr><td>α = 4.1</td><td>K<sub>ep</sub> = 10<sup>-5.0</sup></td><td>η = 0.03</td><td>92.0</td><td>0.37</td><td>19.56</td><td>2.05</td><td>5.7</td><td>5.61</td><td>4.63</td></tr><tr><td>α = 4.2</td><td>K<sub>ep</sub> = 10<sup>-5.0</sup></td><td>η = 0.1</td><td>92.0</td><td>0.6</td><td>20.64</td><td>2.32</td><td>4.92</td><td>5.66</td><td>4.76</td></tr><tr><td>α = 4.25</td><td>K<sub>ep</sub> = 10<sup>-4.0</sup></td><td>η = 0.09</td><td>87.0</td><td>0.27</td><td>19.34</td><td>1.11</td><td>4.34</td><td>5.71</td><td>5.06</td></tr><tr><td>α = 4.15</td><td>K<sub>ep</sub> = 10<sup>-5.0</sup></td><td>η = 0.06</td><td>86.0</td><td>0.48</td><td>22.83</td><td>2.31</td><td>5.2</td><td>5.65</td><td>4.63</td></tr><tr><td>α = 4.15</td><td>K<sub>ep</sub> = 10<sup>-4.5</sup></td><td>η = 0.04</td><td>85.0</td><td>0.37</td><td>15.26</td><td>1.87</td><td>4.65</td><td>5.66</td><td>4.83</td></tr><tr><td>α = 4.2</td><td>K<sub>ep</sub> = 10<sup>-4.5</sup></td><td>η = 0.08</td><td>85.0</td><td>0.38</td><td>22.67</td><td>2.03</td><td>4.6</td><td>5.64</td><td>4.88</td></tr><tr><td>α = 4.25</td><td>K<sub>ep</sub> = 10<sup>-4.0</sup></td><td>η = 0.08</td><td>83.0</td><td>0.32</td><td>13.97</td><td>0.96</td><td>4.13</td><td>5.73</td><td>5.03</td></tr><tr><td>α = 4.2</td><td>K<sub>ep</sub> = 10<sup>-4.0</sup></td><td>η = 0.05</td><td>76.0</td><td>0.29</td><td>15.11</td><td>1.04</td><td>4.19</td><td>5.74</td><td>5.03</td></tr><tr><td>α = 4.2</td><td>K<sub>ep</sub> = 10<sup>-5.0</sup></td><td>η = 0.07</td><td>74.0</td><td>0.68</td><td>10.41</td><td>1.82</td><td>4.89</td><td>5.63</td><td>5.22</td></tr><tr><td>α = 4.25</td><td>K<sub>ep</sub> = 10<sup>-4.5</sup></td><td>η = 0.1</td><td>74.0</td><td>0.58</td><td>10.87</td><td>1.72</td><td>4.14</td><td>5.66</td><td>5.08</td></tr><tr><td>α = 4.2</td><td>K<sub>ep</sub> = 10<sup>-4.0</sup></td><td>η = 0.06</td><td>73.0</td><td>0.21</td><td>24.56</td><td>1.19</td><td>4.54</td><td>5.71</td><td>5.01</td></tr><tr><td>α = 4.1</td><td>K<sub>ep</sub> = 10<sup>-4.0</sup></td><td>η = 0.02</td><td>70.0</td><td>0.2</td><td>19.06</td><td>1.41</td><td>4.74</td><td>5.78</td><td>4.88</td></tr><tr><td>α = 4.2</td><td>K<sub>ep</sub> = 10<sup>-4.5</sup></td><td>η = 0.06</td><td>69.0</td><td>0.49</td><td>11.26</td><td>1.65</td><td>4.3</td><td>5.69</td><td>4.99</td></tr><tr><td>α = 4.3</td><td>K<sub>ep</sub> = 10<sup>-3.5</sup></td><td>η = 0.1</td><td>69.0</td><td>0.19</td><td>18.19</td><td>0.78</td><td>4.05</td><td>5.82</td><td>5.14</td></tr><tr><td>α = 4.15</td><td>K<sub>ep</sub> = 10<sup>-4.5</sup></td><td>η = 0.05</td><td>68.0</td><td>0.31</td><td>25.19</td><td>2.17</td><td>4.76</td><td>5.66</td><td>4.81</td></tr><tr><td>α = 4.25</td><td>K<sub>ep</sub> = 10<sup>-4.0</sup></td><td>η = 0.1</td><td>66.0</td><td>0.24</td><td>25.46</td><td>1.2</td><td>4.56</td><td>5.68</td><td>5.02</td></tr><tr><td>α = 4.15</td><td>K<sub>ep</sub> = 10<sup>-5.0</sup></td><td>η = 0.04</td><td>61.0</td><td>0.57</td><td>10.67</td><td>1.7</td><td>5.24</td><td>5.67</td><td>5.1</td></tr><tr><td>α = 4.15</td><td>K<sub>ep</sub> = 10<sup>-4.0</sup></td><td>η = 0.03</td><td>60.0</td><td>0.26</td><td>14.63</td><td>1.13</td><td>4.28</td><td>5.78</td><td>4.98</td></tr><tr><td>α = 4.25</td><td>K<sub>ep</sub> = 10<sup>-3.5</sup></td><td>η = 0.07</td><td>50.0</td><td>0.14</td><td>26.48</td><td>0.78</td><td>4.35</td><td>5.82</td><td>5.12</td></tr></tbody></table>

**Notes.** These are the simulated populations found to be the most compatible with the HGPS results, as described in Fig. 11. Column 3: the hadronic ratio is the ratio of the number of SNRs with gamma-ray luminosity dominated by the hadronic component to the total number of SNRs. In Columns 5–8, Had. and Lep. refer to simulated SNRs dominated by hadronic emission and simulated SNRs dominated by leptonic emission respectively.

**Acknowledgements.** This work has been funded by the Deutsche Forschungsgemeinschaft (DFG, German Research Foundation) with the grant 500120112 and supported by the grant PID2021-1245810B-I00 funded by MCIN/AEI/10.13039/501100011033 and 2021SGR00426 of the Generalitat de Catalunya. This work was also supported by the Spanish program Unidad de Excelencia María de Maestu CEX2020-001058-M. P.C. acknowledges support from the PSL Starting Grant GALAPAGOS. The authors gratefully acknowledge the computing time made available to them on the high-performance computer “Lise” at the NHR Center NHR@ZIB. This center is jointly supported by the Federal Ministry of Education and Research and the state governments participating in the NHR (https://www.nhr-verein.de/unsere-partner).

## References

Acciari, V. A., Aliu, E., Arlen, T., et al. 2011, ApJ, 730, L20

Acero, F., Lemoine-Goumand, M., Renaud, M., et al. 2015, A&A, 580, A74

Acero, F., Ackermann, M., Ajello, M., et al. 2016, ApJS, 224, 8

Adriani, O., Barbarino, G. C., Bazilevskaya, G. A., et al. 2011, Science, 332, 69

Albert, J., Aliu, E., Anderhub, H., et al. 2007, ApJ, 664, L87

Bell, A. R. 1978, MNRAS, 182, 443

Bell, A. R., Schure, K. M., Reville, B., & Giacinti, G. 2013, MNRAS, 431, 415

Bisschoff, D., Potgieter, M. S., & Aslam, O. P. M. 2019, ApJ, 878, 59

Blasi, P. 2004, Astropart. Phys., 21, 45

Blasi, P. 2017, MNRAS, 471, 1662

Blondin, J. M., Wright, E. B., Borkowski, K. J., & Reynolds, S. P. 1998, ApJ, 500, 342

Brott, I., de Mink, S. E., Cantiello, M., et al. 2011a, A&A, 530, A115

Brott, I., Evans, C. J., Hunter, I., et al. 2011b, A&A, 530, A116

Castor, J., McCray, R., & Weaver, R. 1975, ApJ, 200, L107

Cherenkov Telescope Array Consortium (Acharya, B. S., et al.) 2019, Science with the Cherenkov Telescope Array (World Scientific Publishing Co.)

Chevalier, R. A. 1982, ApJ, 258, 790

Choi, J., Dotter, A., Conroy, C., et al. 2016, ApJ, 823, 102

Corso, N. J., Diesing, R., & Caprioli, D. 2023, ApJ, 954, 1

Cristofari, P. 2021, Universe, 7, 324

Cristofari, P., & Blasi, P. 2019, MNRAS, 489, 108

Cristofari, P., Gabici, S., Casanova, S., Terrier, R., & Parizot, E. 2013, MNRAS, 434, 2748

Cristofari, P., Gabici, S., Humensky, T. B., et al. 2017, MNRAS, 471, 201

Cristofari, P., Blasi, P., & Amato, E. 2020, Astropart. Phys., 123, 102492

Cristofari, P., Blasi, P., & Caprioli, D. 2021a, A&A, 650, A62

Cristofari, P., Niro, V., & Gabici, S. 2021b, MNRAS, 508, 2204

CTA Consortium 2023, ArXiv e-prints [arXiv:2310.02828]

Cui, Y., Yeung, P. K. H., Tam, P. H. T., & Pühlhofer, G. 2018, ApJ, 860, 69

Cummings, A. C., Stone, E. C., Heikkila, B. C., et al. 2016, ApJ, 831, 18

Das, S., Brose, R., Meyer, D. M. A., et al. 2022, A&A, 661, A128

Dewdney, P. E., Hall, P. J., Schilizzi, R. T., & Lazio, T. J. L. W. 2009, IEEE Proc., 97, 1482

Dotter, A. 2016, ApJS, 222, 8

Eggenberger, P., Meynet, G., Maeder, A., et al. 2008, Ap&SS, 316, 43

Ekström, S., Georgy, C., Eggenberger, P., et al. 2012, A&A, 537, A146

Ellison, D. C., Patnaude, D. J., Slane, P., Blasi, P., & Gabici, S. 2007, ApJ, 661, 879

Ellison, D. C., Patnaude, D. J., Slane, P., & Raymond, J. 2010, ApJ, 712, 287

Faucher-Giguère, C.-A., & Kaspi, V. M. 2006, ApJ, 643, 332

Ferrière, K. M. 2001, Rev. Mod. Phys., 73, 1031

Finke, J. D., & Dermer, C. D. 2012, ApJ, 751, 65

Fragos, T., Andrews, J. J., Bavera, S. S., et al. 2023, ApJS, 264, 45

Fujita, Y., Yamazaki, R., & Ohira, Y. 2022, ApJ, 933, 126

Gabici, S. 2017, AIP Conf. Ser., 1792, 020002

Gabici, S., Evoli, C., Gaggero, D., et al. 2019, Int. J. Mod. Phys. D, 28, 1930022

Georgy, C. 2010, PhD Thesis, University of Geneva, Switzerland

Green, D. A. 2015, MNRAS, 454, 1517

Hess, V. F. 1912, Phys. Z., 13, 1084

H.E.S.S. Collaboration (Abramowski, A., et al.) 2011, A&A, 531, A81

H.E.S.S. Collaboration (Abramowski, A., et al.) 2015, A&A, 574, A100

H.E.S.S. Collaboration (Abdalla, H., et al.) 2018a, A&A, 612, A3

H.E.S.S. Collaboration (Abdalla, H., et al.) 2018b, A&A, 612, A7

H.E.S.S. Collaboration (Abramowski, A., et al.) 2018c, A&A, 612, A4

H.E.S.S. Collaboration (Abdalla, H., et al.) 2018d, A&A, 612, A5

H.E.S.S. Collaboration (Abdalla, H., et al.) 2018e, A&A, 612, A6

Jonas, J., & MeerKAT Team 2016, MeerKAT Science: On the Pathway to the SKA, 1

Kafexhiu, E., Aharonian, F., Taylor, A. M., & Vila, G. S. 2014, Phys. Rev. D, 90, 123014

Khangulyan, D., Aharonian, F. A., & Kelner, S. R. 2014, ApJ, 783, 100

Kroupa, P. 2001, MNRAS, 322, 231
Lipari, P. 2019, Phys. Rev. D, 99, 043005

Liu, Y., Cao, Z., Chen, S., et al. 2016, ApJ, 826, 63

Longair, M. S. 1994, High Energy Astrophysics. Volume 2: Stars, the Galaxy and the Interstellar Medium (Cambridge: Cambridge University Press)

Mandelartz, M., & Becker Tjus, J. 2015, Astropart. Phys., 65, 80

Meyer, D. M. A., & Meliani, Z. 2022, MNRAS, 515, L29

Meyer, D. M. A., Mignone, A., Petrov, M., et al. 2021, MNRAS, 506, 5170

Meyer, D. M. A., Velázquez, P. F., Petruk, O., et al. 2022, MNRAS, 515, 594

Meyer, D. M. A., Pohl, M., Petrov, M., & Egberts, K. 2023, MNRAS, 521, 5354

Morlino, G., & Caprioli, D. 2012, A&A, 538, A81

Nguyen-Dang, N. T., Pühlhofer, G., Sasaki, M., et al. 2023, A&A, 679, A48

Ohira, Y., & Yamazaki, R. 2017, J. High Energy Astrophys., 13, 17

Pais, M., & Pfrommer, C. 2020, MNRAS, 498, 5557

Paxton, B., Marchant, P., Schwab, J., et al. 2015, ApJS, 220, 15

Paxton, B., Schwab, J., Bauer, E. B., et al. 2018, ApJS, 234, 34

Paxton, B., Smolec, R., Schwab, J., et al. 2019, ApJS, 243, 10

Pohl, M. 1996, A&A, 307, L57

Ptuskin, V. 2010, ApJ, 718, 31

Ptuskin, V. S., & Zirakashvili, V. N. 2003, A&A, 403, 1

Ptuskin, V. S., & Zirakashvili, V. N. 2005, A&A, 429, 755

Reid, M. J., Menten, K. M., Brunthaler, A., et al. 2019, ApJ, 885, 131

Reynolds, S. P. 2008, ARA&A, 46, 89

Schure, K. M., & Bell, A. R. 2013, MNRAS, 435, 1174

Schure, K. M., Bell, A. R., & O’C Drury, L., & Bykov, A. M., 2012, Space Sci. Rev., 173, 491

Shibata, T., Ishikawa, T., & Sekiguchi, S. 2010, ApJ, 727, 38

Smartt, S. J., Eldridge, J. J., Crockett, R. M., & Maund, J. R. 2009, MNRAS, 395, 1409

Steiman-Cameron, T. Y., Wolfire, M., & Hollenbach, D. 2010, ApJ, 722, 1460

Steppa, C., & Egberts, K. 2020, A&A, 643, A137

Strong, A. W., & Moskalenko, I. V. 1998, ApJ, 509, 212

Sukhbold, T., Ertl, T., Woosley, S. E., Brown, J. M., & Janka, H. T. 2016, ApJ, 821, 38

Sushch, I., Brose, R., & Pohl, M. 2018, A&A, 618, A155

Szécsi, D., Agrawal, P., Wünsch, R., & Langer, N. 2022, A&A, 658, A125

Tingay, S. J., Goeke, R., Bowman, J. D., et al. 2013, PASA, 30, e007

Truelove, J. K., & McKee, C. F. 1999, ApJS, 120, 299

van Veelen, B., Langer, N., Vink, J., García-Segura, G., & van Marle, A. J. 2009, A&A, 503, 495

Vannoni, G., Gabici, S., & Aharonian, F. A. 2009, A&A, 497, 17

Völk, H. J., Berezhko, E. G., & Ksenofontov, L. T. 2005, A&A, 433, 229

Weaver, R., McCray, R., Castor, J., Shapiro, P., & Moore, R. 1977, ApJ, 218, 377

Zabalza, V. 2015, Proc. Int. Cosm. Ray Conf., 2015, 922

Zeng, H., Xin, Y., & Liu, S. 2019, ApJ, 874, 50

Zirakashvili, V. N., & Aharonian, F. A. 2010, ApJ, 708, 965

Zirakashvili, V. N., & Ptuskin, V. S. 2012, Astropart. Phys., 39, 12