# Prospects for a survey of the Galactic plane with the Cherenkov Telescope Array

S Abe, J Abhir, A Abhishek, F Acero, A Acharyya, R Adam, A Aguasca-Cabot, I Agudo, A Aguirre-Santaella, J Alfaro, et al.

► To cite this version:

S Abe, J Abhir, A Abhishek, F Acero, A Acharyya, et al.. Prospects for a survey of the Galactic plane with the Cherenkov Telescope Array. JCAP, 2024, 10, pp.081. 10.1088/1475-7516/2024/10/081. hal-04246716

HAL Id: hal-04246716

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

Submitted on 27 Jan 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
PAPER • OPEN ACCESS

# Prospects for a survey of the galactic plane with the Cherenkov Telescope Array

To cite this article: S. Abe et al JCAP10(2024)081

View the [article online](#) for updates and enhancements.

You may also like

- Multi-messenger Observations of a Binary Neutron Star Merger

B. P. Abbott, R. Abbott, T. D. Abbott et al.

- First Measurement of the Hubble Constant from a Dark Standard Siren using the Dark Energy Survey Galaxies and the LIGO/Virgo Binary-Black-hole Merger GW170814

M. Soares-Santos, A. Palmese, W. Hartley et al.

- Fast b-tagging at the high-level trigger of the ATLAS experiment in LHC Run 3

G. Aad, B. Abbott, K. Abeling et al.
# Prospects for a survey of the galactic plane with the Cherenkov Telescope Array

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

cherenkov
telescope
array

## The CTA Consortium<fnref n="1" />

E-mail: fabio.acero@cea.fr, masha.chernyakova@dcu.ie,  
barbara.olmi@inaf.it, quentin.remy@mpi-hd.mpg.de,  
luigi.tibaldo@irap.omp.eu

**ABSTRACT:** Approximately one hundred sources of very-high-energy (VHE) gamma rays are known in the Milky Way, detected with a combination of targeted observations and surveys. A survey of the entire Galactic Plane in the energy range from a few tens of GeV to a few hundred TeV has been proposed as a Key Science Project for the upcoming Cherenkov Telescope Array Observatory (CTAO). This article presents the status of the studies towards the Galactic Plane Survey (GPS). We build and make publicly available a sky model that combines data from recent observations of known gamma-ray emitters with state-of-the-art physically-driven models of synthetic populations of the three main classes of established Galactic VHE sources (pulsar wind nebulae, young and interacting supernova remnants, and compact binary systems), as well as of interstellar emission from cosmic-ray interactions in the Milky Way. We also perform an optimisation of the observation strategy (pointing pattern and scheduling) based on recent estimations of the instrument performance. We use the improved sky model and observation strategy to simulate GPS data corresponding to a total observation time of 1620 hours spread over ten years. Data are then analysed using the methods and software tools under development for real data. Under our model assumptions and for the realisation considered, we show that the GPS has the potential to increase the number of known Galactic VHE emitters by almost a factor of five. This corresponds to the detection of more than two hundred pulsar wind nebulae and a few tens of supernova remnants at average integral fluxes one order of magnitude lower than in the existing sample above 1 TeV, therefore opening the possibility to perform unprecedented population studies. The GPS also has the potential to provide new VHE detections of binary systems and pulsars, to confirm the existence of a hypothetical population of gamma-ray pulsars with an additional TeV emission component, and to detect bright sources capable of

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

<footnote n="1">¹Full author list at p. 61.</footnote>
accelerating particles to PeV energies (PeVatrons). Furthermore, the GPS will constitute a pathfinder for deeper follow-up observations of these source classes. Finally, we show that we can extract from GPS data an estimate of the contribution to diffuse emission from unresolved sources, and that there are good prospects of detecting interstellar emission and statistically distinguishing different scenarios. Thus, a survey of the entire Galactic plane carried out from both hemispheres with CTAO will ensure a transformational advance in our knowledge of Galactic VHE source populations and interstellar emission.

KEYWORDS: gamma ray experiments, gamma ray theory

arXiv ePRINT: 2310.02828
Contents

<table><thead><tr><th></th><th></th></tr></thead><tbody><tr><td>1 Introduction</td><td>2</td></tr><tr><td>2 Sky model</td><td>4</td></tr><tr><td>2.1 Known sources</td><td>5</td></tr><tr><td>2.2 Synthetic source populations</td><td>6</td></tr><tr><td>2.3 Interstellar emission</td><td>11</td></tr><tr><td>3 Observation plans</td><td>11</td></tr><tr><td>3.1 Instrument description, simulation strategy and software</td><td>12</td></tr><tr><td>3.2 Pointing pattern</td><td>13</td></tr><tr><td>3.3 Observation scheduling</td><td>16</td></tr><tr><td>3.4 Observation simulation</td><td>17</td></tr><tr><td>4 Source catalogues</td><td>17</td></tr><tr><td>4.1 Analysis outline</td><td>17</td></tr><tr><td>4.2 Diagnostics and results</td><td>18</td></tr><tr><td>5 Population studies: PWNe and SNRs</td><td>21</td></tr><tr><td>6 Dedicated analyses of other source classes</td><td>23</td></tr><tr><td>6.1 Gamma-ray binaries</td><td>23</td></tr><tr><td>6.2 Pulsars</td><td>26</td></tr><tr><td>6.3 PeVatrons</td><td>28</td></tr><tr><td>7 Diffuse emission</td><td>30</td></tr><tr><td>7.1 Unresolved sources</td><td>30</td></tr><tr><td>7.2 Interstellar emission</td><td>33</td></tr><tr><td>8 Summary and perspectives</td><td>35</td></tr><tr><td>A Details on the construction of the sky model</td><td>40</td></tr><tr><td>A.1 Pulsar selection and models</td><td>40</td></tr><tr><td>A.2 Other dedicated models</td><td>41</td></tr><tr><td>A.3 Spectral cutoffs</td><td>41</td></tr><tr><td>A.4 Removal of bright synthetic sources</td><td>44</td></tr><tr><td>B Details on observation scheduling</td><td>46</td></tr><tr><td>C Details on catalogue A production</td><td>47</td></tr><tr><td>C.1 Object detection</td><td>47</td></tr><tr><td>C.2 Model fitting</td><td>47</td></tr><tr><td>D Details on catalogue B production</td><td>48</td></tr></tbody></table>
<table><tr><td>D.1 Datasets</td><td>48</td></tr><tr><td>D.2 Object detection</td><td>48</td></tr><tr><td>D.3 Classification, spatial model selection, and candidate filtering</td><td>49</td></tr><tr><td>D.4 Model fitting</td><td>50</td></tr><tr><td>D.5 Models refinement</td><td>51</td></tr><tr><td><strong>E Maps of unresolved point sources</strong></td><td><strong>51</strong></td></tr><tr><td><strong>The CTA Consortium</strong></td><td><strong>61</strong></td></tr></table>

# 1 Introduction

Current gamma-ray observations in the very-high-energy (VHE) domain (above a few tens of GeV) have revealed a wide range of particle acceleration and interaction phenomena in Galactic astrophysical objects including pulsar wind nebulae (PWNe, e.g., [1]), supernova remnants (SNRs, e.g., [2]), and gamma-ray binaries [e.g. 3], as well as many sources whose nature is still unknown. Important results include the prevalence of extended sources powered by pulsars in the VHE sky [e.g. 4], the measurement of soft gamma-ray spectra for most SNRs which challenges the standard paradigm for the origin of Galactic cosmic rays (CRs) up to the knee [e.g. 5], and the existence of particle confinement or slow diffusion around many of these sources [e.g. 6]. Despite these results, the mechanisms of particle acceleration and their efficiencies, the maximum achievable energies, the nature of the accelerated particles and the radiation processes, the role played by the source's environment and the evolution with the source's age, and the contributions to the population of Galactic CRs remain uncertain in most cases.

Surveys, which provide large samples of sources with a small selection bias, have enabled major advances in astronomy in general, and in the field of VHE astronomy in particular. VHE surveys were carried out with imaging atmospheric Cherenkov telescopes (IACTs) which were the first instruments with sufficient sensitivity to detect the brightest TeV sources [7], and with ground-based air shower arrays that only recently reached sensitivity to detect larger numbers of sources. Air shower arrays have a worse angular resolution and worse sensitivity for short observations, but they reach higher energies and have a larger field of view so that they naturally survey the entire observable sky and are better suited to detect very extended sources and Galactic diffuse emission.

Historically, the first survey of a quarter of the Milky Way (with Galactic longitude range $-2^{\circ} < l < 85^{\circ}$) at very high energies was performed by the HEGRA system of imaging atmospheric Cherenkov telescopes (IACTs), with no significant detection of gamma-ray sources [8]. The largest IACT survey to date was performed by H.E.S.S. [4], and covered the longitude range $-110^{\circ} < l < 65^{\circ}$ and latitudes $|b| \le 3^{\circ}$, with a total of 78 sources detected. Additionally, VERITAS performed a survey of the nearby Cygnus star-forming region, leading to the detection of a new source and the detailed study of three other VHE sources [9].

Among air-shower arrays, Milagro observed the Galactic plane at longitudes $30^{\circ} < l < 216^{\circ}$ and detected TeV emission towards 14 previously known Galactic gamma-ray sources [10]. ARGO-YBJ surveyed the Milky Way at longitudes $-10^{\circ} < l < 70^{\circ}$ detecting four gamma-ray excesses with statistical significance greater than five standard deviations
(> 5σ) toward the Galactic plane [11]. The third catalog of sources detected by the HAWC
observatory [12], which covers Galactic longitudes 0° < l < 110° and 150° < l < 240°,
includes 65 sources located mostly towards the Galactic plane. The LHAASO collaboration
have recently released their first catalog of gamma-ray sources [13] with 90 sources detected
with high significance (> 5σ) above 1 TeV, including 32 new TeV sources, as well as 43 sources
with lower significance (> 4σ) above 100 TeV. LHAASO has detected photons from Galactic
sources up to a maximum of 1.4 PeV [14].

Furthermore, due to its wide field of view and decade-long exposure, the Large Area Telescope (LAT) aboard the *Fermi Gamma-ray Space Telescope* has reached into the VHE domain. The LAT has surveyed the entire sky, detecting 189 high-significance sources at energies above 30 GeV at $|b| < 10°$ [15]. Of those, 61 are VHE emitters known from previous surveys.

Diffuse TeV gamma-ray emission was also observed towards the Galactic plane [16–20]. While at GeV energies Galactic diffuse emission is interpreted as the result of CR interactions with interstellar gas and radiation fields, at TeV it is likely a superposition of interstellar emission and collective emission from populations of sources that cannot be detected individually given the sensitivities of the current instruments [e.g. 6, and references therein].

The Cherenkov Telescope Array (CTA) is the next-generation IACT system [21, 22]. CTA will consist of a southern array located at the Paranal Observatory (Chile) and a northern array located at Roque de los Muchachos Observatory (Spain), allowing observations of almost the entire sky. Its energy range will extend from 20 GeV to beyond 200 TeV, with a sensitivity improvement of an order of magnitude with respect to current IACT systems. With a large field of view, reaching > 8° diameter at the highest energies,<fnref n="1" /> and an angular resolution of a few arcminutes, CTA is well-suited to perform large surveys.

A survey of the entire Galactic plane was proposed as a Key Science Project for the
CTA Observatory (CTAO) [23, 24] with the goals of:

1. providing a census of Galactic VHE emitters, such as SNRs and PWNe, through the detection of hundreds of new sources and, therefore, enabling population studies;

2. identifying a list of promising candidates for follow-up observations, such as new gamma-ray binaries, PeVatron candidates (objects capable of accelerating particles to PeV, i.e., up to the knee of the CR spectrum);

3. determining the properties of TeV diffuse emission from the Galactic plane;

4. searching for new and unexpected high-energy phenomena in the Milky Way;

5. providing a legacy dataset to the astronomical and astroparticle physics communities.

The results from the Galactic Plane Survey (GPS) during the first years of operations will also
be instrumental in seeding other observational programmes, including several Key Science

<footnote n="1">¹For an IACT array the effective field of view is determined by the camera’s geometrical field of view, and also by a decreasing detection efficiency at increasing angular distances from the centre of the cameras. The webpage https://www.cta-observatory.org/science/ctao-performance/#1472563544190-020879e1-468f shows how the CTA sensitivity to a point-like source varies as a function of angular distance from the centre of the field of view and of photon energy.</footnote>
Projects and the general observer programme. Therefore, periodic releases of data and scientific results to the public during the execution of the GPS are proposed [24]. We note that other dedicated Key Science Projects are proposed to study transient and flaring sources, including Galactic sources and extragalactic sources observed towards the Galactic plane [24], for which GPS observations can also be exploited. Galactic transients are discussed in a companion publication [25].

In this article we provide an updated and improved assessment of prospects for the GPS compared to [24]. We describe the proposed implementation and expected scientific returns, as well as the development of analysis pipelines. In section 2 we describe the sky model used in our simulations, based on recent gamma-ray source catalogues and on observations of individual sources, as well as state-of-the-art models of the main populations of VHE gamma-ray sources (SNRs, PWNe, binaries) and interstellar gamma-ray emission. In section 3 we describe the implementation of the survey in terms of pointing pattern and observation scheduling and the simulations of the observations used for this paper. In section 4 we describe how we created catalogues of sources in the entire Galactic plane based on the results of the analyses of the simulated data. In section 5 we discuss the properties of the two dominant source classes detected in the survey: PWNe and SNRs. In section 6 we present analyses of the simulated data aimed at studying gamma-ray binaries, pulsars, and PeVatrons. In section 7 we discuss Galactic diffuse emission, including the emission from unresolved sources and from interstellar processes. Finally, in section 8, we present a summary of the results and discuss perspectives for the future.

## 2 Sky model

Gamma-ray emission from the Milky Way in the VHE domain is dominated by individual emitters, which include SNRs, pulsars, PWNe and pulsar halos, gamma-ray binaries, star-forming regions (SFRs), and sources of still unknown nature [e.g., 4]. A large-scale diffuse component is also detected and interpreted as a mix of interstellar emission from the interactions of CRs with interstellar gas and radiation fields and of collective emission from populations of sources that cannot be detected individually given the sensitivities of the current instruments [e.g., 17].

The sky model used in this paper consists of three main components that are described in this section:

1. a set of known sources modelled on the basis of observations from current and past instruments;

2. synthetic population models tuned to existing observations for the three main classes of Galactic VHE emitters (SNRs, PWNe, and gamma-ray binaries) to model the sources that will be detected by CTA as well as those that will remain undetected and will then constitute an unresolved diffuse emission component;

3. models for interstellar emission.
## 2.1 Known sources

We model the population of known sources with a compilation of existing catalogues from IACTs, the Fermi-LAT space telescope, and air-shower arrays. For the known VHE sources detected by IACTs we used the compilation provided by gamma-cat<fnref n="2" /> that we enhanced based on recent results not yet included in the official gamma-cat release. The completeness of this source list was checked against TeVCat<fnref n="3" /> [26].

Galactic variable sources detected at very high energies (i.e. binaries and pulsars) are also included in our sky model. For pulsars, we include dedicated spectral and temporal models for 38 objects preselected among those detected by Fermi LAT [27, 28]. We arbitrarily implement two spectral classes, either extending from GeV to TeV energies as simple steep power laws extrapolated from Fermi-LAT measurements (Crab-like pulsars), or with an additional TeV component modelled as a power-law with an exponential cutoff whose brightness is randomly generated assuming a TeV to GeV flux ratio ranging from $5 \times 10^{-5}$ to $10^{-2}$, while complying at the same time with upper-limits available in the literature (Vela-like pulsars). More details and the full list of pulsars simulated are provided in appendix A.1. Gamma-ray binaries are known to be variable on orbital and sometimes superorbital timescales [see, e.g., 3, for a review]. Some of them also display spectral variability as a function of the orbital phase. Therefore, for six known gamma-ray binaries and one candidate gamma-ray binary within the sky area covered by the GPS, we introduce dedicated temporal and spectral models as detailed in table 7.

We then add sources detected by Fermi LAT [3FHL catalogue, 28] and by HAWC [2HWC catalogue,<fnref n="4" /> 29]. Sources are added only if they are not associated with already included TeV sources. Furthermore, we exclude 3FHL and 2HWC sources for which a source among those discussed above lies within their $3\sigma$ position uncertainty. Finally, to limit the number of very soft sources of limited interest for CTA, we include only 3FHL sources detected at energies > 50 GeV with a significance > $3\sigma$, and associated with at least one photon of measured energy > 100 GeV.

In gamma-cat, 2HWC, and 3FHL, the morphology of extended gamma-ray sources is modelled using simple geometrical templates (disc or Gaussian shapes). In order to test source detection and characterisation algorithms on more realistic source models, we include 14 sources modelled by dedicated templates, that either replace the corresponding catalogue model or are added to the sky model for a few recent detections. Templates are also used to model some diffuse features not captured by the large-scale interstellar emission models described below. The list and a description of the templates are provided in table 8.

Furthermore, the gamma-ray catalogues describe the source spectra using simple analytical functions. Often the limited statistics above a few TeV prevent current instruments from detecting spectral curvature, therefore the catalogues report the spectra as simple power laws. However, based on the maximum energy achievable by the accelerators, radiative cooling of the particles, and the physics of the gamma radiation mechanisms, a power law

<footnote n="2">²https://gamma-cat.readthedocs.io.</footnote>

<footnote n="3">³http://tevcat2.uchicago.edu/.</footnote>

<footnote n="4">⁴The 3HWC catalogue, based on a 30% longer integration time, and the first catalogue of LHAASO sources became available only when this work was already in an advanced state.</footnote>
is not expected to be a good approximation of the gamma-ray spectrum up to several hundreds of TeV for many sources (see e.g. [30]). In order to assess the capabilities of CTA to measure spectral curvature in a more realistic scenario, we add a physically-motivated exponential cutoff to the power-law spectra with a hard spectral index (< 2.4) as described in appendix A.3. The cutoff is applied to the spectrum of known SNRs, PWNe, active galactic nuclei<fnref n="5" /> (AGNs), and unidentified sources, which are treated as if they belong to the dominant source class according to the energy band of the observations (PWNe for gamma-cat and 2HWC, AGNs for 3FHL). We do not modify the spectrum of the two PeVatron candidates Westerlund 1 [31] and HESS J1641−463 [32], which are therefore included with the currently measured power-law spectrum without any curvature.

## 2.2 Synthetic source populations

CTA is expected to extend our horizon of detectability and to observe a number of presently undetected faint sources. Therefore, our sky model needs to include a population of synthetic sources representing fainter, and yet undetected objects. The following subsections describe the methodology we use to create synthetic populations of the three main classes of VHE sources, namely SNRs, PWNe, and binaries.

### 2.2.1 Supernova Remnants

SNRs are simulated using the approach described in [33, 34]: a supernova (SN) event is randomly placed in the Galaxy in time and position, assuming a constant rate of 3 explosions per 100 years [35] over the last $10^5$ years.

Four types of SNe are considered: thermonuclear (TN), high-ejecta-mass core-collapse (CC-HEM), low-ejecta-mass core-collapse (CC-LEM) and high-explosion-energy core-collapse (CC-HEE). The following relative rates for the various SN types are assumed [36]: 0.32, 0.44, 0.22, and 0.02 respectively. Each type is characterised by a typical range of physical parameters: the mass of the ejecta (from $1.4 M_{\odot}$ to $20 M_{\odot}$), the mass loss-rate and the velocity of the wind (as in [36]). The explosion energy is fixed to $10^{51}$ erg. For a better description of core-collapse systems associated with PWNe, we consider the mass of the ejecta to have a Gaussian distribution with peak at $13 M_{\odot}$ and $\sigma = 3 M_{\odot}$, truncated at a minimum value of $5 M_{\odot}$ and a maximum of $20 M_{\odot}$ [37].

The Galaxy is modelled with four logarithmic spiral arms as in [38]; core-collapse SNe follow the distribution of pulsars [39] characterised by the Galactocentric radial distribution described in [38], while thermonuclear SNe are placed considering a flat distribution for the inner region [40], plus a radial distribution in the outer zone [41, 42]. The height above the Galactic plane follows the distribution of the molecular gas.

The evolution of SNRs requires us to follow the time variation of the position of the forward shock ($R_{sh}$), marking the shell boundary in the ambient medium, modelled following the three-dimensional distribution of the atomic and molecular hydrogen presented in [43, 44], with density varying between $10^{-5}$ and $10 \text{ cm}^{-3}$.

<footnote n="5">⁵AGNs are background sources for Galactic studies. Those detected by current instruments within the sky area covered by the GPS are included in the sky model.</footnote>
Gamma-ray emission is produced by particles accelerated at the SNR shock, both via neutral-pion decay (in case of nuclei) and inverse-Compton (IC) scattering (in case of electrons). Their spectra are modelled as a power-law in the particle momentum p, following [33, 34]: $f(R_{\text{sh}}, p, t) = A(t)p^{-\alpha}$, where $\alpha$ is considered as a free parameter, ranging over $4.1 \le \alpha \le 4.4$ (e.g. [45–50]). The normalisation $A(t)$ is determined by requiring that the CR pressure at the shock is a fraction of the shock ram pressure (fixed to 0.1, [47]). The maximum energy is determined by the assumption that particles escape the acceleration region when their diffusion length equals 1/10 of the shock radius.

Electrons are accelerated at the same rate as protons, and the two types of particles have identical spectral shapes at low enough energies, where radiation losses affecting electrons can be neglected. The ratio between the electron and proton spectra is assumed to be $10^{-4}$ [33, 34, 51, 52]. The electron spectrum is shaped by radiation losses above a threshold energy (where the characteristic energy-loss time becomes shorter than the SNR age), and it steepens by 1.0 in spectral index with respect to the injection spectrum [52].

The spatial distribution of the accelerated particles inside the SNR is determined by solving a transport equation, assuming a 20 µG magnetic field inside the SNR shell, a momentum-dependent Bohm-like diffusion coefficient, advection and radiative losses (details can be found in [33]).

The gamma-ray emission from each SNR of the population is finally obtained by calculating: (i) the hadronic component from proton-nuclei interactions (see [53] and [54] for nuclei heavier than hydrogen); (ii) the (dominant) leptonic component from IC scattering of the accelerated electrons on CMB photons [55].

### 2.2.2 Interacting Supernova Remnants (iSNRs)

In recent years, SNRs interacting with molecular clouds (iSNRs) have emerged as a new class of gamma-ray emitters with observational characteristics which distinguish them from other members of the SNR class. In these systems, the gamma-ray emission originates from inelastic interactions of CRs accelerated by the SNR shock with the matter in the cloud (through neutral-pion decay). These systems are characterised by emission that is spatially coincident with the molecular clouds [e.g. 56, and references therein].

In order to produce a synthetic population of such systems, three ingredients are needed: a synthetic population of molecular clouds, the properties of CRs inside these clouds, and the probability that a molecular cloud has an SNR close enough to interact. A synthetic population of Galactic molecular clouds is created using the spiral arm distribution from [38]. A mass is assigned randomly to each cloud, assuming a power-law mass function of spectral index $-1.6$ for clouds with masses > $10^3 M_{\odot}$ [57]. A broken power-law model is assumed for the CRs inside the cloud. The spectral parameters (i.e. the slopes and the energy breaks) as well as the CR densities inside the clouds are extracted randomly from the distributions obtained for the observed SNR-molecular cloud systems and assigned to each molecular cloud [58]. The gamma-ray spectrum is then calculated assuming hadronic emission via the proton-proton interaction. The probability of interaction with an SNR is assumed to be 1.5% in order to reproduce the observed flux distribution for known systems. A possible leptonic component associated with these systems is not considered.
### 2.2.3 Pulsar Wind Nebulae

Synthetic PWNe are produced by associating a pulsar to every core-collapse SN from the population described in section 2.2.1. The properties of each pulsar are drawn following the prescription by [59], based on gamma-ray emitting pulsars. Each source is evolved in time assuming a pure dipole braking for the pulsar spin-down (i.e. a braking index $n = 3$). We note that no important variations are expected in the population if the index is varied over the range $2.3 \le n \le 3$ [60]. A three-dimensional kick velocity is associated with each pulsar, following a double-sided exponential distribution with mean value $\sim 380 \text{ km s}^{-1}$ [38].

Previous attempts to provide a population of Galactic PWNe, and to estimate its gamma-ray contribution, were based on the assumption of an “average” PWN, with properties tuned to match observations (see e.g. [1]). Here instead, we use the same approach as in [61], where the (approximated) dynamical and radiative evolution of each source in the population is reproduced using a one-zone model. This provides a physically-informed description of the intermediate phase of the PWN evolution, when the nebula starts to interact with the shocked ejecta of the SNR (usually known as the *reverberation* phase). This phase may cause important variations in the PWN spectral energy distribution at all energies, due to the system compression, and so must be properly treated in order not to introduce biases in the population study [60, 62].

Leptons responsible for the emission are injected in each PWN following a broken power-law in energy: $Q(E, t) = Q_0(t)(E/E_b)^{\alpha_i}$, with the break energy randomly drawn from a lognormal distribution centred at $E_b = 0.28 \text{ TeV}$ and with a $0.12 \text{ TeV}$ spread. The injection indices range over $1.0 \le \alpha_1 \le 1.7$ for $E \le E_b$, and $2.0 \le \alpha_2 \le 2.7$ otherwise. The normalisation $Q_0(t)$ is determined by requiring that the power injected in particles equals a fraction $(1-\eta)$ of the pulsar luminosity $L(t)$, with the magnetic fraction $\eta$ ranging over $0.02-0.2$. The radiation properties of each PWN are then computed by evolving the injection function in time, considering both adiabatic and radiation losses (synchrotron and IC) plus the possible escape through diffusion. The magnetic field is modelled as in [63], considering the variation of the magnetic energy due to both the adiabatic expansion/contraction of the nebula, and the energy injection from the pulsar. Gamma-ray emission through IC is finally evaluated using the Klein-Nishina cross-section for the photon-electron interaction [55], considering different local photon fields: photons from the CMB, synchrotron photons, and near and far IR photons from the stellar background, described by a normal distribution centred on the value estimated from the GALPROP model at the specific position of each source in the Galaxy [64].

A considerable fraction of the pulsars escape their SNR on a time scale smaller than the final age of the simulation. The escaping pulsars leave a relic nebula at their original position, which remains an active source of gamma-ray emission through IC radiation as long as the pulsar-injected particles survive against losses and adiabatic expansion. Escaping pulsars form then a bow shock nebula that inject high-energy particles in the ambient medium [65] and are possibly associated with the formation of TeV halos. The number of expected halos in the Galaxy is still debated — ranging from a few to hundreds [66–68] — and a conclusive physical model is still lacking. For this reason, here we have only included two *ad hoc* models for the halos surrounding Monogem and Geminga [69], based on their observational properties (table 8).
### 2.2.4 Gamma-Ray Binaries

The synthetic population of binaries consists of a list of simulated systems giving the location in the Galaxy, the orbital period, and the orbital lightcurve above 1 TeV. The population is modelled following [70], to which we refer for details and references. The binary is assumed to be a 1.4 M☉ neutron star in orbit with a 30 M☉ massive star, i.e. a generic pulsar-driven gamma-ray binary. The location in the Galaxy is randomly drawn from the radial distribution of OB stars and following a model of the spiral arms [38, 71]. The gamma-ray emission assumes that high-energy electrons injected close to the neutron star inverse Compton scatter UV photons from the massive star to VHE energies. The model also takes into account the absorption of the VHE emission due to pair production with the UV photons from the star following [72], as well as physical eclipses of the VHE emission region (assumed pointlike) by the massive star. The orbital period is drawn from a flat distribution in $\log P_{\text{orb}}$ from 1 to $10^4$ days while the eccentricity is drawn from a thermal eccentricity distribution [73] corrected for circularisation at short orbital periods. To normalise the lightcurves, the injected power in VHE-emitting electrons is assumed to be ≈ 1% of the total pulsar spin-down power. This was estimated from the VHE observations of the rotation-powered pulsars PSR J2032+4127 and PSR B1259−63, which both orbit a massive star. The radiative efficiency is taken into account by scaling the power to the ratio of escape timescale to radiative timescale. The total available power is randomly drawn from the observed distribution of spin-down power of young pulsars.

The resulting lightcurves display a variety of complex shapes as the relative position of the star, the particles and the observer changes as a function of orbital phase. The shapes range from sinusoidal behaviour when the orbit is circular to flare-like behaviour for highly-eccentric orbits (see [70] for examples). Hence, even though this model was built to reproduce a binary population consisting of young rotation-powered pulsars in orbit around a massive star (gamma-ray binaries), the resulting lightcurves also provide reasonable test cases for the outburst-like behaviour expected of other classes of binary sources, such as microquasars.

### 2.2.5 Comparison of the synthetic source populations with the known gamma-ray sky and assembly of the source model

We have described in the previous subsections the construction of a sky model including known sources and populations of synthetic sources. In order to reduce computing time and resources we eliminate the low-flux tails of the populations and only retain for the following synthetic sources above an integral photon flux threshold: $6 \times 10^{-16} \text{ cm}^{-2} \text{ s}^{-1}$ in the 100 GeV-1 TeV energy interval for steady sources, $2 \times 10^{-17} \text{ cm}^{-2} \text{ s}^{-1} > 1 \text{ TeV}$ at the lightcurve peak for binaries. Figure 1 compares source counts as a function of integral flux $> 1 \text{ TeV}$ for the different populations. To generate the synthetic populations and, in particular, to tune the properties of PWNe, we make the simplifying hypothesis that all currently unidentified sources are PWNe. Therefore, as expected, our source model is largely dominated by PWNe.

The first LHAASO catalogue [13] became available when this work was already advanced and was not included in our sky model. However, our synthetic PWN model naturally accounts for a comparable number of sources extending beyond 100 TeV. Indeed, the PWN population contains 72 sources distributed over the entire sky with flux $> 100 \text{ TeV}$ similar to the 43 sources detected by LHASSO in the Northern hemisphere only.
<figure><img src="image_4.png" /><figcaption>**Figure 1.** Cumulative number of sources for Galactic latitude < 6° as a function of integral source flux above 1 TeV, showing known emitters and the synthetic source populations from our sky model for different source classes. For known emitters the vertical bars show the estimated Poisson fluctuations. The shaded blue vertical band shows the target sensitivity of the GPS defined in [24].</figcaption></figure>

In order to build the sky model used as input for realistic simulations, we need to remove sources from the synthetic populations that have equivalents in the known source part of our model to avoid double-counting some sources. Indeed, as the synthetic source models account for the entire populations, the brightest sources have most likely already been detected with existing telescopes. For each source already detected (SNR, iSNR, PWN, composite SNR-PWN system, binary), therefore, we exclude the most similar synthetic source belonging to the same class. The similarity is established based on source position, angular extension, and integral flux > 1 TeV. The exact procedure used to identify the most similar source is detailed in appendix A.4. After removing those sources, we obtain the global population model of Galactic sources that will be used for the rest of this paper. This model is also made publicly available to the community.<fnref n="6" /> We stress that, for sources below the detection threshold for current instruments, the model is based on a single realisation of the synthetic populations and not on the mean expectation.

<footnote n="6">⁶The model is available at https://doi.org/10.5281/zenodo.8402519.</footnote>
## 2.3 Interstellar emission

Uncertainties in the predictions for interstellar emission across the CTA energy range are remarkably large [e.g. 74, 75]. Therefore we use a small set of models to assess the uncertainties in our results. The reference interstellar emission model (IEM) used to produce the simulated CTA dataset (hereafter “base model”, or IEM-base) is based on the DRAGON code that evaluates the propagation of charged CRs in the Galaxy [76, 77]. The key assumption behind the base model is that CR transport is homogeneous and isotropic, while the Galaxy is axisymmetric. The interstellar emission, originating from neutral-pion decay, IC scattering on the diffuse low-energy photon background, and bremsstrahlung, is computed with the HERMES code [78]. The hadronic component is based on a gas model composed of a set of column density maps in (l, b) Galactic coordinates for atomic and molecular gas, associated with Galactocentric rings.<fnref n="7" /> For IEM-base we model the CO-to-H₂ conversion factor, X<sub>CO</sub>, using the functional form from [80] with parameters tuned to be consistent with recent observations of the Milky Way [e.g. 81]:

$$X_{\text{CO}} = 0.6 \times 10^{0.4 \left( \frac{R}{5\,\text{kpc}} - 1 \right)} \times 10^{20} \, \text{cm}^{-2} (\text{K km s}^{-1})^{-1}, \quad (2.1)$$

where R is the Galactocentric radius.

The base model provides a minimal expectation for interstellar emission in the TeV energy range. A first alternative model (IEM-varmin) features an inhomogeneous CR transport including a linear gradient in the index of the diffusion coefficient. This model is discussed in a recent paper [74], and predicts fluxes in the multi-TeV domain larger than the base model by a factor of a few to ten. The normalisation of this model is adjusted mainly by tuning the masses of target molecular gas via X<sub>CO</sub> over large regions of the sky. The resulting X<sub>CO</sub> factors are 2–5 times higher than the values usually measured in the Milky Way. Therefore, we introduce a second alternative model (IEM-varmin rescaled) with the same CR transport setup, but with a rescaling of the X<sub>CO</sub> values and gamma-ray emissivity per gas nucleon at 8 GeV that match Fermi-LAT measurements [82]. This model predicts an intermediate flux between the base and varmin models at TeV energies.

## 3 Observation plans

The GPS will consist of observations of the entire Galactic plane using both the southern and northern arrays [24]. The target sensitivity for isolated point-like sources is, at integral photon fluxes above 1 TeV, ≈ 5 × 10<sup>-14</sup> cm<sup>-2</sup> s<sup>-1</sup>, a factor deeper than existing surveys, and coherently covering the full Milky Way. The survey will be graded so that regions with a higher expected number of sources, especially the inner part of the Galaxy, will receive significantly more observation time from a total allocated budget of 1620 hours over ten years, partitioned between the two arrays and different regions, as described in table 6.3 of [24] and summarised in table 1. These observing times will be used as a working hypothesis for the rest of this paper.

<footnote n="7">⁷The details of the gas model are provided in [79] and https://fermi.gsfc.nasa.gov/ssc/data/analysis/soft ware/aux/4fgl/Galactic_Diffuse_Emission_Model_for_the_4FGL_Catalog_Analysis.pdf.</footnote>
<table><caption>**Table 1.** Observing times in different regions of the Galactic plane. The total times are split in a short-time programme (STP) and a long-time programme (LTP), see main text for more details. Abridged from table 6.3 of [24].</caption><thead><tr><th>Region</th><th>STP (h)</th><th>LTP (h)</th><th>Total (h)</th></tr></thead><tbody><tr><td><strong>SOUTH</strong></td><td></td><td></td><td></td></tr><tr><td>300°-60°, Inner region</td><td>300</td><td>480</td><td>780</td></tr><tr><td>240°-300°, Vela, Carina</td><td></td><td>180</td><td>180</td></tr><tr><td>210°-240°</td><td></td><td>60</td><td>60</td></tr><tr><td><strong>NORTH</strong></td><td></td><td></td><td></td></tr><tr><td>60°-150°, Cygnus, Perseus</td><td>180</td><td>270</td><td>450</td></tr><tr><td>150°-210°, anticentre</td><td></td><td>150</td><td>150</td></tr></tbody></table>

## 3.1 Instrument description, simulation strategy and software

The results presented in the paper are based on the full-scope array described in [83], also known as the CTA baseline or Omega configuration. This includes 4 Large-Sized Telescopes (LSTs), 25 Medium-Sized Telescopes (MSTs), and 70 Small-Sized Telescopes (SSTs) in the Southern hemisphere, plus 4 LSTs and 15 MSTs in the Northern hemisphere. We refer the reader to [83] for a detailed description of the telescopes and their characteristics.

The data discussed and analysed in the rest of the paper are produced through simulations of the sky model described in section 2. Simulated data in the context of this work consist of science-ready event lists, that is, lists of candidate photons with their estimated arrival directions and energies. Such lists will be produced from an analysis of the Cherenkov images recorded by the telescopes. They are one of the main products expected from CTAO (dubbed Data Level 3, or DL3) and are planned to be openly distributed to the community for the scientific exploitation of the data [see, e.g. 84].

We simulate the event lists based on instrument response functions (IRFs) extracted from detailed Monte Carlo simulations of the interactions of gamma rays (and background CRs) with the atmosphere and the subsequent collection and recording of their Cherenkov light signal by the telescopes, followed by image analysis, shower reconstruction, and event selection. The IRFs include a description of the gamma-ray effective area, point spread function (PSF), and energy dispersion, as well as a description of the residual CR background as a function of position with respect to the camera centre and measured energy. We used the IRFs prod3v2 [85] optimised for observation durations of 50 hours.<fnref n="8" /> For the rest of the paper we will always assume that the IRFs are perfectly known.

During the preparation of this work, an initial array configuration with a reduced number of telescopes (the Alpha configuration) was approved for construction [22, 84]. While the

<footnote n="8">⁸Specifically, the South_z20_50h, South_z40_50h and South_z60_50h IRFs were used for the Southern array, while the North_z20_50h and North_z40_50h IRFs were used for the Northern array. The IRF optimisation takes into account a trade-off between background rejection and gamma-ray efficiency that depends on the duration of the observations. In practise, for the observations described in this article, the IRFs optimised for 50 hours are always the most appropriate among those currently available.</footnote>
telescope reduction will impact the prospects of the survey outlined in this paper to some degree, this is expected to be mitigated by ongoing improvements in the event-level analysis tools. In particular, we expect to maintain the off-axis sensitivity which is critical for survey observations. A comparison between current performance curves for the Alpha configuration and those corresponding to the IRFs used in our paper is presented in [86].

For easier reading we also list here the simulation and analysis software tools used in the rest of the paper with the corresponding versions and references.

* `ctools`, version 1.6.3 or higher [87];

* `Gammapy` version 0.17 or higher [88].

Simulations were performed with `ctools`. Catalogues were produced with both `ctools` and `Gammapy` independently. For other parts of the analysis one of the two software packages was conveniently chosen.

## 3.2 Pointing pattern

The GPS pointing pattern has been discussed in [23] and in [24]. Two pointing patterns were considered, either with pointing directions aligned along the Galactic equator (single row), or alternating above and below the equator following an equilateral triangle pattern (double row). Here we revisit the GPS pointing strategy by considering more pattern candidates (shown in figure 2), the IRFs prod3b-v2, and a lower energy threshold of 70 GeV (above which CTA observations of a typical duration of 50 h become competitive with 10 years of Fermi-LAT observations for steady point-like sources). Besides the aforementioned single-row and double-row patterns we consider two more possibilities:

* a triple-row pattern with pointing directions alternating on and above/below the Galactic equator, aimed at obtaining more exposure at higher latitudes;

* a non-equilateral double-row pattern, with pointing directions alternating above and below the equator but with different spacing in latitude and longitude (see figure 2 for an illustration); this is motivated by the need to split the total observing times into short individual observations (with a typical duration of 30 minutes) in order to sample the lightcurves of periodic sources (see below) and, at the same time, to allow to use a smaller spacing in longitude (which yields a more uniform coverage) while still preserving good coverage at higher latitudes.

For consistency with [23, 24] we refer to the distance between adjacent pointing directions of a single-row or equilateral pointing pattern as the pattern step $s$ (see figure 2). As we will discuss later, the latitude spacing $h$ is one of the most critical parameters. For easier comparison with the other patterns we define the step of the non-equilateral double row pattern as $s = \sqrt{4/3}h$, i.e., the step of an equilateral double row pattern with the same latitude spacing $h$. To characterise the non-equilateral double row pattern we also introduce a parameter $w$ that corresponds to the longitude distance between adjacent pointings in alternate rows.
<figure><img src="image_5.png" /><figcaption>**Figure 2.** Schematic view of the pointing patterns considered. Dashed lines represent the Galactic equator. Circles represent the CTA field of view (FoV) for an individual pointing. Red arrows show the pattern step (*s*), i.e., the distance between pointing directions of adjacent pointings for single row and equilateral patterns. Green arrows show the latitude spacing (*h*) for patterns with multiple rows. The orange arrow shows twice the longitude spacing (*w*) for the non-equilateral double row pattern. The figure is not to scale. Notably the actual CTA field of view, in relation to the chosen step size, is larger than shown in this schematic view.</figcaption></figure>

We calculate exposure and effective (PSF) maps for all the patterns under consideration and for energy thresholds of 70 GeV, 125 GeV, and 1 TeV, using the IRFs described above<fnref n="9" /> and for observation zenith angles of 20° and 40°. The calculations described here have been performed using the *ctexpcube* and *ctpsfcube* tools from *ctools*. For this step and the rest of the pattern optimisation procedure the time per pointing is always adjusted to result in a fixed time spent per unit longitude.

Steps $\lesssim 3^\circ$ yield exposure fluctuations as a function of longitude $< 7\%$ for latitudes $|b| < 2^\circ$ and, for multiple-row patterns, a spill-over of exposure $< 10\%$ at high latitudes $|b| > 5^\circ$ (where the number density of Galactic sources and the intensity of interstellar Galactic emission are very small). The instrument PSF varies as a function of photon energy and inclination with respect to the pointing direction. The effective PSF, calculated for each position in the sky by taking into account the pointing pattern, is improved for smaller steps, but the effect in terms of variations of the 68% containment angle is $< 5\%$ for the patterns considered here with $s \lesssim 5^\circ$.

Then we investigate the sensitivity to an isolated point-like source for all the configurations described above. We calculated the sensitivity using *cssens* (from *ctools*) assuming 6.5 hours

<footnote n="9">⁹Results are unchanged in terms of optimal pattern if we use IRFs optimised for 5 hours of observations, while the latter yield slightly worse sensitivities overall.</footnote>
<figure><img src="image_6.png" /><figcaption>**Figure 3.** Integral flux sensitivity to an isolated point-like source at energies $E > 125$ GeV averaged over two latitude bands $|b| < 2^\circ$ and $|b| < 5^\circ$ for observations with the southern array at 40° zenith angle. We assume an observing time of 6.5 hours per degree in longitude. The dashed line highlights the step chosen for the non-equilateral double row pattern used in the rest of the paper.</figcaption></figure>

of time spent per degree in longitude (foreseen for the inner Galaxy) and $w = 0.1^\circ$ as a representative value of longitude spacing for the non-equilateral double-row pattern.<fnref n="10" />

In figure 3 we show the sensitivity at energies > 125 GeV for observations with the southern array at 40° zenith angle for all the pointing patterns considered in two latitude bands: $|b| < 2^\circ$ and $|b| < 5^\circ$. The non-equilateral double-row pattern yields the best sensitivity in the Galactic plane, with a broad minimum for steps $s$ between 2° and 3° and the least dependence upon the exact choice of step. At higher latitudes, the single-row pattern yields much worse sensitivity than multiple-row schemes. Multiple-row patterns yield a sensitivity that improves as a function of increasing step size and tends to a plateau for $s > 2^\circ$. The non-equilateral double-row pattern is only marginally worse than the triple-row scheme for the same step. Results are similar for the other energy ranges considered, for other zenith angles, for the northern array, and for sources with moderate angular extensions ($< 0.2^\circ$).

Based on these considerations, for the rest of the paper we adopt the non-equilateral double-row pattern with a step $s = 2.25^\circ$ ($h = 1.95^\circ$). We split the total observing times per longitude range given in table 1 in individual pointings of 30 minutes each. The longitude spacing $w$ in each longitude range is set to the ratio of the longitude span over the total number of pointings in that range within each year of the programme (see the next section on scheduling and appendix B). Figure 4 shows the observing time and resulting sensitivity for a point-like source as a function of position in the Galactic plane and in three broad energy ranges. Including observations at an offset $< 5^\circ$ for a given direction, the total observing time at any given position on the Galactic equator varies from 15 h to 60 h as a function of longitude, with an average offset with respect to the field of view centre of 2.5°. We note that the differences between the two arrays in the number of telescopes and their

<footnote n="10">¹⁰This value corresponds to the approximate spacing to be used for the inner Galaxy region based on the observing time allocated for 30-minute observations.</footnote>
<figure><img src="image_7.png" /></figure>

**Figure 4.** Top panel: observing time as a function of position in the sky, including all observations at an offset < 5° for a given direction. Other panels: sensitivity to a point-like source with power-law spectrum of index 2.34 in three energy bands. The sensitivity is defined as the minimum energy flux integrated over the energy range of interest to achieve a 5σ detection. The sensitivity estimates assume that the longitude range 210°-60° is observed using the southern array, the rest of the plane using the northern array.

optimisation for different energy ranges are clearly visible when comparing the observing time to the achieved sensitivity as a function of longitude. Regions observed from the southern hemisphere (210°-60°) benefit from better performance, especially at the highest energies.

## 3.3 Observation scheduling

The scheduling of the observations affects the detectability of variability in source fluxes. We implement a realistic scheduling strategy that follows the overall plan outlined in [24]. Specifically, the simulated GPS observations were spread over ten years, with a short-term programme (STP) for which 480 h of observing time were allocated over the first two years, and a long-term programme (LTP) for which 1140 h of observing time were allocated over the following eight years (table 1). The simulated observations were scheduled from 1 January 2021 in dark conditions only, requiring that the Sun is at least 15° below the horizon and that the Moon is also below the horizon. Observations were scheduled so that a given pointing is observed as close as possible to its minimum zenith angle. Other factors affecting the observability such as weather and hardware failures were not taken into account for simplicity.
In addition, observations were distributed over the year so that a given location of the sky is revisited at different time intervals, enabling the detection of periodic flux variations that cover periodicities between a few days and up to a few years. In this way, the periods of all known, and most expected, gamma-ray binaries should be accessible using the GPS data. Practically this was achieved by tentatively distributing the observations as a function of time, and then altering the time distribution algorithm until we reached an exposure as a function of time difference between pointings at all longitudes as uniform as possible and without any gaps in the interval between 5 days and about 1 year. A detailed description of the algorithm is provided in appendix B.

## 3.4 Observation simulation

For each pointing, ctobssim (from ctools) was used to generate mock event lists, based on the sky model described in section 2 and the IRFs (taking into account the instrument energy dispersion). To cope with the limited number of zenith angles for which IRFs were available, pointings with a zenith angle smaller than 30° were then assigned the 20° IRF, pointings with a zenith angle between 30° and 50° were assigned the 40° IRF, and other pointings were assigned the 60° IRF. Only a few pointings were finally scheduled which required the 60° IRF. Most of the results presented below, unless stated otherwise, are based on a single realisation of the simulated event lists.

# 4 Source catalogues

The source catalogues are planned to be one of the major products delivered by CTAO to the community. In the following we present the developments done to prepare a catalogue production workflow and discuss its potential outcome in term of detectable sources from the CTA GPS survey.

## 4.1 Analysis outline

In figure 5 we show the excess counts above the true CR background for the full GPS survey at latitudes $|b| < 6°$ and for energies between 0.07 and 200 TeV. Using the simulated data and the true IRFs we build catalogues of sources in the entire Galactic plane in this energy range.

The first step of the catalogue production is to build, in a short amount of computational time, a list of candidate objects from the structures found in the excess or significance maps. The candidate objects are then fitted with different models in order to determine the best-fit model and its optimal parameters. The fitted candidates are filtered such that, for each object, the test statistic $TS = 2 \Delta \ln(L) > 25$ with $\Delta \ln(L)$ being the difference in log-Likelihood between the best-fit model including the source and the model for the null hypothesis (no source). An additional step would be needed to determine the optimal threshold in order to minimise both false positive and false negatives rates. Furthermore, the threshold should be set in terms of significance, since the statistical meaning of the TS value depends on the number of degrees of freedom associated with the source model. However, these aspects are not critical for the results we present in the following as we can compare them to the true sky model.
<figure><img src="image_8.png" /><figcaption>**Figure 5.** Excess counts above the CR background in the 0.07–200 TeV energy range from the entire CTA GPS. The spatial bin width is 0.06° and a Gaussian smoothing with σ = 0.03° is applied.</figcaption></figure>

The catalogue production was performed independently using two different pipelines. In the following we will refer to the output catalogues as A and B, respectively. The main differences between the output catalogues result from the analysis strategies used for the initial object detection (different algorithms), and the model fitting refinement (order and type of models tested). Catalogue A is based on the work presented in [89] and summarised in appendix C. Catalogue B partly relies on the techniques and ideas discussed in [90], but the complete work-flow is detailed in appendix D. In the next section we will comment on how the differences between the two strategies affect the results.

## 4.2 Diagnostics and results

In order to match the detected objects in the catalogue with the sources in the true sky model we test for spatial coincidence using two criteria. For each object we first select true sources within a given angular separation to the object centre defined as $d_{\text{inter-center}} < d_{\min} + f_R \times R_{\text{object}}$, where $R_{\text{object}}$ is the 68% containment radius of the fitted model convolved by the PSF. We set $d_{\min} = 0.1°$ which corresponds to about twice the 68% containment radius of the PSF at 1 TeV and $f_R = 30\%$ which is close to the dispersion in relative error on the radius guess from the first step of the catalogue production. After this first selection we report the association maximising the surface overlap fraction, defined as the ratio of the intersection to the union of the candidate object and true source sur-
<figure><img src="image_9.png" /><figcaption>**Figure 6.** Integral flux of the sources in the true sky model versus integral flux of the associated objects from catalogue B in two energy bands 0.07–1 TeV (left) and 1–200 TeV (right). The dashed black line corresponds to a one-to-one match. The “Known” label corresponds to a compilation of sources detected by the current generation telescopes at GeV or TeV energies (see section 2.1).</figcaption></figure>

<table><caption>**Table 2.** Number of detectable sources (TS > 25 for the true sky model, labelled True detectable) and detected objects (TS > 25 in the catalogue) in the 0.07–200 TeV energy range. Unmatched corresponds to detected objects without direct match with any simulated sources, ΔF/F is the relative error on the total flux of the detectable sources, $f_{\text{match}}$ is the fraction of the detected objects matching a true source, and $f_{\text{reco}}$ is the fraction of the true detectable sources matching a detected object. The first columns correspond to the synthetic populations and the “Known” column to sources already detected by the current generation telescopes (at GeV or TeV energies). Further details on the content of the true sky model are given in section 2 and appendix A.</caption><thead><tr><th>Name</th><th>PWN</th><th>SNR</th><th>iSNR</th><th>Binaries</th><th>Known</th><th>Unmatched</th><th>Total</th><th>ΔF/F</th><th>f<sub>match</sub></th><th>f<sub>reco</sub></th></tr></thead><tbody><tr><td>True detectable</td><td>294</td><td>37</td><td>24</td><td>10</td><td>134</td><td>—</td><td>499</td><td>—</td><td>—</td><td>—</td></tr><tr><td>Catalogue A</td><td>241</td><td>16</td><td>20</td><td>10</td><td>111</td><td>169</td><td>567</td><td>-12.5%</td><td>0.70</td><td>0.80</td></tr><tr><td>Catalogue B</td><td>257</td><td>31</td><td>14</td><td>10</td><td>122</td><td>36</td><td>470</td><td>3.8%</td><td>0.92</td><td>0.87</td></tr></tbody></table>

faces: $SF_{\text{overlap}} = (S_{\text{object} \cap \text{source}}) / (S_{\text{object} \cup \text{source}})$, where the surfaces are delimited by the iso-contours at 68% containment in flux of the model convolved by the PSF. We choose to report only the association that maximises the surface overlap fraction in order to limit the possible associations for extended objects. We also enforce that each detected object can be associated with only one true source and vice-versa. Moreover, only the associations with $SF_{\text{overlap}} > 0.25$ are reported in order to limit spurious associations. We checked on a sample of mock catalogues that this threshold maximizes the balanced accuracy defined as the average of true positive and true negative rates for associations.

In order to test the overall quality of the catalogue produced we introduce the matching fraction, $f_{\text{match}} = N_{\text{match}}/N_{\text{objects detected}}$, defined as the fraction of the detected objects matching a true source; and the reconstruction fraction, $f_{\text{reco}} = N_{\text{match}}/N_{\text{sources detectable}}$, defined as the fraction of the true detectable sources that match a detected object. A source
<figure><img src="image_10.png" /><figcaption>**Figure 7.** Left panel: detection ratio of sources above a given integrated flux (for E = 0.07–200 TeV). The grey curve corresponds to the expected detections with TS > 25 for the true sky model, the coloured curves correspond to the detections reported in the catalogues at the same threshold. Note that the detection ratio of a catalogue can exceed 100% because of the confusion bias or modelling bias (see details in the text). Right panel: histogram of the nearest neighbour distance between detected objects.</figcaption></figure>

is defined as detectable if it has a TS > 25 for the true sky model (the same threshold as for the detected objects in the catalogue). $f_{\text{match}}$ and $f_{\text{reco}}$ are indicators of purity and completeness, respectively.

Table 2 gives the number of sources simulated that are detectable for the different synthetic source populations and for the known sources. Based on the matching criterion previously defined, we also report the potential detections from the catalogues associated to the same populations. These results show that we may detect up to 500 sources in the 0.07–200 TeV energy range from the CTA GPS which is more than 6 times the number of objects in the HESS-GPS [4] or the 3HWC [91] catalogues.

Overall, the relative error on the integral flux of the detectable sources is about 4% (for catalogue B); the good agreement in flux of those objects spatially associated with true sources can be seen in figure 6. The larger deviation for the fainter sources can be explained by source confusion, in particular in the composite systems where the PWN and/or the SNR are not significant enough to be individually detected. Outliers are expected for the known sources as some of these are simulated with complex templates or multiple models, rather than with a single parametric model, which makes the one-to-one associations more ambiguous. In a few cases, sources simulated with multiple components of similar spatial extent but with different spectra can be associated with the wrong fitted component, as the association criterion is purely spatial. In the most common case, one complex source is fitted with multiple simpler models. This can be seen as a modelling bias (further discussed below). Finally, some of these deviations may be peculiar to the individual realisation of the simulated dataset considered.

We also define the detection ratio as the number of objects detected with TS > 25 above a given integrated flux, divided by the total number of simulated sources above the same flux (figure 7, left). The detection ratio of a catalogue can exceed 100% because of confusion or modelling biases. In the case of the confusion bias, the emission from the unresolved sources biases the flux of the sources near the threshold upwards and so enhances their detection probability. In the case of the modelling bias, we see that sources simulated
with more complex models (shell, elliptical Gaussian, or template) than those considered in the catalogue construction (Gaussian, disc, and point-like) can be fragmented into multiple smaller objects of lower flux. Another indicator of this effect is the increase of unmatched object detections with a small distance to their nearest neighbour, as shown in the right panel of figure 7 (with reference to catalogue A).

The fragmentation of complex sources into multiple sub-structures explains most of the discrepancy observed for catalogue A in figure 7. It also explains the large number of objects detected but the low matching fraction, as shown in table 2. At this stage, the filtering and merging of the sub-structures is more a classification than a statistical problem. Thus, the solution introduced in the production of catalogue B to solve this issue is to use pattern recognition techniques: (i) to determine *a priori* the best-suited morphological model between a shell or a generalized Gaussian (see appendix D.4); (ii) to identify *a posteriori* the groups of objects that could be merged together or replaced by a different model (see appendix D.5). The use of a larger variety of spatial models and the identification of multiple objects as a single entity allows a better description of the complex sources in the true sky model. This leads to a better estimate of the number of source detections and their individual fluxes which explains the closer agreement of the detection ratio in catalogue B seen in figure 7. The results of catalogue B in terms of detection ratio and matching fraction show that we can produce a catalogue close to the limit of detections expected from the true sky model.

For the following sections one of the two catalogues, A or B, was conveniently chosen. Differences between the quality of the catalogues are such that they do not affect the main conclusions discussed.

# 5 Population studies: PWNe and SNRs

In this section we will discuss the properties of the two dominant source classes detected in the simulated survey: PWNe and SNRs. We caution the reader that the absolute detection numbers depend on the individual realisations of the population models and the CTA configuration considered. However, this will not affect the main conclusions.

PWNe are the dominant source class at TeV energies in the Galactic plane. About 250 synthetic PWN detections in addition to the known ones are included in the catalogues described above and a brief overview of their properties is discussed below. We will focus on the population study of the synthetic sources, as opposed to real sources observed with current instruments, in order to have a uniform sample for which all physical parameters are known (Ê, age, Galactic coordinates, distance, etc.).

Figure 8 (left panel, upper row) shows the spatial distribution of the entire synthetic PWN population generated, along with the objects detected in catalogue B and associated with a synthetic PWN. The CTA GPS sensitivity makes it possible to detect a large number of objects even at the far edge of our Galaxy. As expected, a selection bias is observed in the Galactic distribution of detected sources as a function of heliocentric distance. Only the most luminous objects are detected at the farthest edges of the Galaxy.

About one third of the detected synthetic PWNe are found to be spatially extended at a statistical significance greater than 3σ which is an important feature to help identify the source category (i.e. likely PWN or SNR origin). The completeness of this survey is illustrated
<figure><img src="image_11.png" /></figure>

**Figure 8.** Properties of the synthetic PWN population and corresponding detections. Left panel, upper row: galactic distribution of the synthetic (shown with blue circles) and detected objects (circles with colour coding in terms of luminosity). The Sun is represented with a yellow star at a distance of 8.5 kpc from the Galactic centre. Note that the lack of detections in the solar neighbourhood is due to nearby synthetic sources having similar properties to known sources being removed from the population (see section 2). Right panel, upper row: distribution of the synthetic and detected PWNe in the luminosity-$\dot{E}$ space. The dotted lines indicate that the survey could detect 50% of PWN with $\dot{E} > 10^{36}$ erg s$^{-1}$ and 70% of sources with luminosity $> 10^{33}$ s$^{-1}$. Bottom row: integral flux histogram of the CTA GPS detections compared to the entire synthetic PWN population and the HESS GPS PWN catalogue [1].

in figure 8 (right panel, upper row) in terms of the intrinsic properties of the sources. We can see that this survey is able to detect half of all the PWNe in the Galaxy currently powered by an energetic pulsar with spin-down luminosity $\dot{E} > 10^{36}$ erg s$^{-1}$. Concerning detection limits, 95% of the sources we have detected have an $\dot{E} > 2.2 \times 10^{34}$ erg s$^{-1}$ and a luminosity $> 2.7 \times 10^{32}$ erg s$^{-1}$. Note that the $\dot{E}$ value used here is the present $\dot{E}$ value (not the value at birth). This, combined with environmental factors, explains the scatter in the Luminosity-$\dot{E}$ correlation as PWNe at TeV energies act as calorimeters and reflect the integral of power injected into the system since birth.
The comparison in figure 8 (bottom row) of our current population of known PWNe [HESS PWN catalogue: 1] with the members of our synthetic population detected in the simulation and the entire synthetic population emphasises the transformational jump in population size that CTA will bring to the field of PWN population studies (about 7 times the current sample, or 2.5 times if we consider that most of the unidentified sources are PWNe as was done in the construction of the population model). In terms of sensitivity, 95% of the HESS PWNe have an integral flux above 1 TeV > $2.2 \times 10^{-13} \text{ cm}^{-2} \text{ s}^{-1}$ compared to > $2.4 \times 10^{-14} \text{ cm}^{-2} \text{ s}^{-1}$ for the CTA sample, an order of magnitude improvement. For a flux > $2.4 \times 10^{-14} \text{ cm}^{-2} \text{ s}^{-1}$, the completeness factor is ~55%, meaning that CTA is expected to detect more than half of the PWN population in our Galaxy above that flux. We note a small discrepancy at the high-flux end in figure 8 (bottom row) where the population of the brightest simulated PWNe is not entirely detected. This is mainly due to a bias in the catalogue analysis. Some extended PWNe can end up being split into multiple sources with smaller extension by the catalog pipeline, leading to a lower reconstructed flux. In addition, in composite SNRs cross talk between the SNR and the PWN can lead to an inaccurate reconstruction of the PWN flux.

The second-most numerous class detected in this survey is SNRs. Focusing only on the synthetic shells and interacting SNRs, 45 (31 and 14 respectively) sources were detected in catalogue B. This suggests that the CTA GPS may be able to increase, by a factor larger than two, the number of SNRs observed at TeV energies. About half of the synthetic SNRs detected are significantly spatially extended, which is a valuable feature for the identification of sources with multiwavelength counterparts. Among the sources simulated as shells (young SNRs), 19 are effectively detected as shells. The distribution of SNRs in integral true simulated flux and distance in figure 9 shows that new sources can be detected out to the other side of our Galaxy (up to 20 kpc in this realisation). The flux sensitivity is 5–10 times better than the current TeV SNR sample, with sources being detected down to an integral flux of a few $10^{-14} \text{ cm}^{-2} \text{ s}^{-1}$. The newly detected shell-type SNRs have an age range from 0–10 kyrs with a detection efficiency of 15–30% in each age bin and an average of 3 SNRs per 1 kyr age bin. This wider distribution of ages than current samples will pave the way to a more detailed population study.

This is a major step forward which allows to explore the population of Galactic SNRs by discovering and measuring the spatial extension of new SNRs. This is possible for CTA even with the relatively limited observing times (< 50 h at an average offset of 2.5°) provided by the GPS. Conversely, the current population of SNRs was established through dedicated deep (> 100 h) observations for the faintest objects.

# 6 Dedicated analyses of other source classes

## 6.1 Gamma-ray binaries

### 6.1.1 Known sources

Six known gamma-ray binaries and one candidate gamma-ray binary (see table 7) were included in the simulation. For all of these we derived a phase-folded lightcurve. In the phase-resolved analysis we included all sources detected in catalogue A within 3.5° from the source of interest. The flux and spectral parameters were left free to vary only for the
<figure><img src="image_12.png" /><figcaption>**Figure 9.** Comparison of the known SNRs and the population of synthetic shell and interacting SNRs in the integral true simulated flux-distance parameter space. New objects can be detected up to a distance of 20 kpc and down to an integral flux of a few $10^{-14} \text{ cm}^{-2} \text{ s}^{-1}$. Note that for a given flux, the detectability of a SNR depends on its extension.</figcaption></figure>

<figure><img src="image_13.png" /><figcaption>**Figure 10.** Phase-folded lightcurve above 100 GeV of LS 5039 on a 30-minute time scale (one point per individual observation) as reconstructed from simulated data. The flux normalisation $dN/dE$ is given at 1 TeV. The dash-dotted black line shows the simulated profile and the black dashed line shows the mean differential flux. See A.2 for a description of the input to the simulation.</figcaption></figure>

source of interest and for nearby sources (within 1°), while they were frozen to the catalogue
A values for more distant ones.
<figure><img src="image_14.png" /><figcaption>**Figure 11.** Variability of the differential flux at 1 TeV of source bin095 on timescales of a single observation (30 minutes, left) and a month (right). Red dots show > 5σ detections and green stars show 95% c.l. upper limits when a significant detection was not achieved. The blue line shows the average flux.</figcaption></figure>

All gamma-ray binaries included in the simulations are bright and are clearly detected by CTA in a 30-minute observation during the orbital phases with differential fluxes at 100 GeV higher than ~2.5×10⁻¹³ cm⁻² s⁻¹ TeV⁻¹ (the exact flux threshold depends on the spectral index of the source). The simulated orbital and spectral variability is clearly detected using the GPS observations (see e.g. the phase-folded lightcurve of LS 5039 in figure 10). In phase bins corresponding to 5–10 hours of observations, spectral parameters of the binaries were reconstructed to better than 10% level.

### 6.1.2 Blind search for variable sources

Catalogue A contains 73 sources, that either were classified by the automatic pipeline as point sources, or associated during the subsequent cross-correlation with point-like sources in the true sky model. For each of these we created a lightcurve using events with energies between 0.1 and 100 TeV from the observations pointed within 3° of the source. For this purpose we fixed the spectral shape of each source to the catalogue results. The flux was left free for the source of interest and nearby sources (located within 1°), and was fixed to the catalogue value for more distant sources. A simple $\chi^2$ test was used to analyse the lightcurves and, if the probability to have a constant flux was less than 5% (without accounting for trials), then the source was classified as being variable. The lightcurves were generated for temporal intervals corresponding to one individual observation (30 minutes), one week, and, for sources simulated with periods longer than a week, one month and one year (see figure 11 for an example).

Table 3 shows the results for sixteen sources simulated as variable on a timescale longer than 30 minutes. Twelve of them were classified as variable in the analysis, according to the criterion described above. On the other hand, four of them were not detected as variable, either because their mean flux is too low or because they have a flat orbital flux profile.

Table 3 also includes two sources that showed a probability of constant flux below 5% in spite of being simulated as constant. In the case of 3FHL J1855.3+0751 we verified that the false variability detection is related to the imperfect modelling of a bright extended shell SNR
<table><caption>**Table 3.** List of variable sources. Sources simulated as variable are given in the upper section, while sources simulated as constant but detected as variable are given in the bottom section. Sources detected as variable and the probabilities *P* for the corresponding time binning are marked in bold. *F*<sub>mean</sub> is the mean differential flux at 1 TeV in units of 10<sup>-12</sup> cm<sup>-2</sup> s<sup>-1</sup>TeV<sup>-1</sup>. *σ*<sub>obs</sub> = (F<sub>max</sub> - F<sub>mean</sub>)/√(ΔF<sup>2</sup><sub>max</sub> + ΔF<sup>2</sup><sub>mean</sub>) is the significance of the maximum flux deviation on the 30 minutes time scale.</caption><thead><tr><th>Source Name</th><th>σ<sub>obs</sub></th><th>F<sub>mean</sub></th><th>P<sub>Obs</sub></th><th>P<sub>Day</sub></th><th>P<sub>Week</sub></th><th>P<sub>Month</sub></th><th>P<sub>Year</sub></th><th>Period<br>d</th></tr></thead><tbody><tr><td>LS I +61° 303</td><td>15.6</td><td>0.97±0.03</td><td><b>2×10<sup>-295</sup></b></td><td><b>8×10<sup>-302</sup></b></td><td><b>5×10<sup>-273</sup></b></td><td>—</td><td>—</td><td>26.5</td></tr><tr><td>PSR B1259-63</td><td>11.1</td><td>0.33±0.01</td><td><b>3×10<sup>-251</sup></b></td><td><b>4×10<sup>-121</sup></b></td><td><b>9×10<sup>-128</sup></b></td><td><b>10<sup>-288</sup></b></td><td><b>10<sup>-165</sup></b></td><td>1241</td></tr><tr><td>bin040</td><td>5.3</td><td>0.06±0.01</td><td>6.5×10<sup>-1</sup></td><td>6.3×10<sup>-1</sup></td><td>1.6×10<sup>-1</sup></td><td><b>5×10<sup>-3</sup></b></td><td><b>5.7×10<sup>-6</sup></b></td><td>3358</td></tr><tr><td>1FGL J1018.6-5856</td><td>4.3</td><td>0.41±0.02</td><td><b>9×10<sup>-22</sup></b></td><td><b>5×10<sup>-25</sup></b></td><td><b>7×10<sup>-29</sup></b></td><td>—</td><td>—</td><td>16.6</td></tr><tr><td>LS5039</td><td>4.0</td><td>4.15±0.03</td><td><b>10<sup>-37</sup></b></td><td><b>10<sup>-37</sup></b></td><td>—</td><td>—</td><td>—</td><td>3.9</td></tr><tr><td>bin095</td><td>3.9</td><td>0.23±0.01</td><td>9.6×10<sup>-2</sup></td><td><b>4.9×10<sup>-2</sup></b></td><td><b>1.8×10<sup>-3</sup></b></td><td><b>3.13×10<sup>-5</sup></b></td><td>—</td><td>200</td></tr><tr><td>HESS J1832-093</td><td>3.5</td><td>2.28±0.02</td><td><b>6.6×10<sup>-7</sup></b></td><td><b>1.8×10<sup>-6</sup></b></td><td><b>1.5×10<sup>-5</sup></b></td><td><b>1.2×10<sup>-4</sup></b></td><td>—</td><td>365</td></tr><tr><td>PSR J2032+4127</td><td>3.4</td><td>3.95±0.05</td><td><b>4.1×10<sup>-26</sup></b></td><td><b>1.2×10<sup>-26</sup></b></td><td><b>5.1×10<sup>-27</sup></b></td><td><b>8.23×10<sup>-29</sup></b></td><td><b>5.36×10<sup>-32</sup></b></td><td><b>1.8×10<sup>4</sup></b></td></tr><tr><td>HESS J0632+057</td><td>3.0</td><td>0.30±0.06</td><td><b>10<sup>-8</sup></b></td><td><b>10<sup>-8</sup></b></td><td><b>7×10<sup>-9</sup></b></td><td><b>6×10<sup>-10</sup></b></td><td>—</td><td>315</td></tr><tr><td>bin074</td><td>2.8</td><td>0.05±0.01</td><td>9.9×10<sup>-1</sup></td><td>9.9×10<sup>-1</sup></td><td>9.3×10<sup>-1</sup></td><td><b>6×10<sup>-3</sup></b></td><td>4.9×10<sup>-1</sup></td><td>840</td></tr><tr><td>bin159</td><td>2.6</td><td>0.11±0.01</td><td>1.6×10<sup>-1</sup></td><td><b>1.6×10<sup>-2</sup></b></td><td><b>4.7×10<sup>-2</sup></b></td><td>—</td><td>—</td><td>5.2</td></tr><tr><td>bin123</td><td>2.3</td><td>0.09±0.01</td><td>3.9×10<sup>-1</sup></td><td>5×10<sup>-1</sup></td><td>3.7×10<sup>-1</sup></td><td><b>1.8×10<sup>-2</sup></b></td><td><b>2×10<sup>-4</sup></b></td><td>522</td></tr><tr><td>bin162</td><td>1.9</td><td>0.10±0.01</td><td>9.9×10<sup>-1</sup></td><td>9.9×10<sup>-1</sup></td><td>9.3×10<sup>-1</sup></td><td>6×10<sup>-2</sup></td><td>4.9×10<sup>-1</sup></td><td>1387</td></tr><tr><td>bin154</td><td>1.9</td><td>0.05±0.01</td><td>9.9×10<sup>-1</sup></td><td>9.7×10<sup>-1</sup></td><td>8.2×10<sup>-1</sup></td><td>6.1×10<sup>-1</sup></td><td>—</td><td>35.8</td></tr><tr><td>bin093</td><td>1.4</td><td>0.06±0.01</td><td>9.9×10<sup>-1</sup></td><td>9.9×10<sup>-1</sup></td><td>9.9×10<sup>-1</sup></td><td>—</td><td>—</td><td>7</td></tr><tr><td>bin146</td><td>1.0</td><td>0.07±0.01</td><td>9.9×10<sup>-1</sup></td><td>9.9×10<sup>-1</sup></td><td>—</td><td>—</td><td>—</td><td>1.5</td></tr><tr><td>pwn2059</td><td>2.5</td><td>2.30±0.17</td><td><b>2.3×10<sup>-5</sup></b></td><td><b>4.8×10<sup>-3</sup></b></td><td><b>2.8×10<sup>-5</sup></b></td><td>—</td><td>—</td><td>N/A</td></tr><tr><td>3FHL J1855.3+0751</td><td>1.1</td><td>1.55±0.37</td><td>9.9×10<sup>-1</sup></td><td><b>10<sup>-5</sup></b></td><td><b>2.7×10<sup>-4</sup></b></td><td><b>1.1×10<sup>-3</sup></b></td><td>—</td><td>N/A</td></tr></tbody></table>

nearby that affects differently the flux determination depending on epoch/pointing direction. By selecting events only within 3° of the source for the analysis we obtain probabilities to have a constant flux > 5%. In the case of the synthetic source pwn2059, the variability detection is determined by a single 30-minute observation for the individual realisation of the simulated dataset considered. A detailed investigation based on 10000 realisations of the observation simulations showed that the low-flux point that drives the variability detection stems from a combination of downward statistical fluctuations in the number of events from the source and the background at the source position. The said fluctuations consist of variations < 10 events, therefore the simple χ<sup>2</sup> test based on Gaussian statistics employed here may not be fully appropriate in the context of blind searches for variable sources with CTA and further developments are needed.

Nevertheless, this study shows that the GPS can be used to find variable sources on different time scales, so that follow-up observations of individual sources to study in detail their properties can be performed.

## 6.2 Pulsars

We analysed data for all known Fermi-LAT pulsars included in the simulations. To this end we selected observations in which the pulsar has a maximum offset of 5° from the centre of the field of view. No relativistic effects or delay in the generation of event time stamps
were considered in the simulations. Therefore, no barycentric corrections are included in
the analysis chain.

As a first step, we searched for pulsed emission at the known positions of the pulsars.
Events were extracted from a region of 5° around the pulsar and weighted by the probability
to originate from the pulsar position. This probability was estimated by using the PSF model
of the instrument as a function of angular distance to the source and energy of each event.
The significance when testing for a periodic signal based on the known ephemerides was
estimated using the weighted H-test [27, 92, 93].

Among the 25 Crab-like pulsars, pulsations from PSR J1833-1034 and PSR J1838-0537 were detected at significance levels (H-test values) of 5.1σ (37) and 7.2σ (70.3), respectively. We note, however, that the quoted significance levels are subject to caution due to a possible contribution of the PWN in the tens of GeV range, owing to the use of phase-averaged spectra to model the pulsed emission in the true sky model (see appendix A.1). Some of the brightest Fermi-LAT pulsars are not detected due to the limitation in latitude coverage of the CTA GPS (see figure 4). Indeed, the latter provides scarce or no observing time/sensitivity at the positions of the Crab and Geminga pulsars, PSR J0007+7303, PSR J1057-5226, and PSR J1514-4946 (for instance, the closest pointing in the simulated GPS lies at 4.8° from the Crab pulsar, located at a Galactic latitude of -5.78°). These objects are however expected to be well-studied through dedicated observations independent of the GPS.

The Vela pulsar, as well as several Vela-like pulsars covered with at least few hours
of observing time by the GPS were significantly detected. Figure 12 shows an example of
the reconstructed phase profile for the Vela pulsar for both GeV and TeV components, as
compared to the input Fermi-LAT phase profile.

In a second step, we performed a maximum likelihood fit to the data without using the
timing information in order to quantify the CTA sensitivity to phase-averaged emission of
pulsars. Events with energies from 30 GeV to 100 TeV were selected from a square region of
4° size around each pulsar. The model spectrum was assumed to be either a simple power-law
function (PWL) or a power law with exponential cutoff (PLEC). The spectral parameters of
nearby sources known from other observations and/or present in Catalogue A were left free.

For all Vela-like pulsars with significant periodic emission we also detected the phase-
averaged TeV component (TS$\ge$25), whereas the GeV component was bright enough to be
detected only for the Vela pulsar. Although the observing time for the latter pulsar is small
(Galactic latitude of -2.79°), we were able to reconstruct with fair precision its spectrum.

None of the Crab-like emission tails are significantly detected in the phase-averaged
analysis. For the two Crab-like pulsars with a pulsation detection we obtained a TS of 0 in the
case of PSR J1838-0537, and a TS of 1.7 in the case PSR J1833-1034. The low TS values
are due to source-confusion effects for these two sources surrounded by multiple overlapping
brighter components, the modelling of which is imperfect in the catalogues. This illustrates
the challenges in searching for faint sources in complex regions of the Galactic plane. The
results on the pulsation search and phase-averaged emission are summarised in table 4.

The detection of pulsations from two Crab-like pulsars illustrates the potential of the
GPS data to be used as a path-finder to probe and constrain the possible extension of some
pulsars' GeV component into the VHE range, and thereby to help optimising the follow-up
<figure><img src="image_15.png" /><figcaption>**Figure 12.** Phasogram of the Vela pulsar from simulated GPS data for two energy ranges: 1–20 TeV (top panel, 0.1° around the pulsar position) and 30–200 GeV (middle panel, 0.25° around the pulsar position). The *Fermi*-LAT profile above 10 GeV used as template in the simulations is shown for comparison in the bottom panel. The grey shaded areas on the two upper panels show the estimated level of the background evaluated in the off-phase interval [0.7–1.0].</figcaption></figure>

observations with CTA. On the other hand, the clear detection of the TeV component in the Vela [94] and hypothetical Vela-like pulsars (with diverse assumptions on their emission spectra) shows that if such a population existed, the GPS would be able to discover it. Thanks to the high sensitivity of CTA GPS in the multi-TeV range, this could even be the case for pulsars which remain still undetected in the GeV range.

## 6.3 PeVatrons

We assess the potential of the GPS to find PeVatrons (or candidates) following the approach presented in [95]. For this analysis we have considered catalogue B from which we extracted a list of sources with $TS > 25$ in the 1–200 TeV energy range and with a number of predicted events from the fitted spectral model at energies higher than 50 TeV greater than 1. In this way we obtained a sample of 231 sources showing spectra extending up to ~100 TeV energies. These sources have been subsequently re-analysed in the full energy range (0.07 to 200 TeV). The analysis setup is the same as for the catalogue generation, with morphological parameters of the source of interest and all parameters of nearby sources and the background models fixed to the catalogue values. The only free parameters of the fit are therefore the spectral parameters of the source of interest.
<table><caption>**Table 4.** Pulsars with evidence of emission from simulated GPS data. The significance levels predicted are subject to caution given possible contamination from the PWN emission in the true sky model.</caption><thead><tr><th rowspan="2">Name</th><th colspan="2">Periodicity search</th><th colspan="2">Phase-averaged analysis</th></tr><tr><th>H-test</th><th>significance</th><th>TS</th><th>spectral model</th></tr></thead><tbody><tr><td colspan="5"><em>Crab-like pulsars</em></td></tr><tr><td>PSR J1833-1034</td><td>37</td><td>5σ</td><td>1.7</td><td>PWL</td></tr><tr><td>PSR J1838-0537</td><td>70.3</td><td>7.2σ</td><td>0</td><td>PWL</td></tr><tr><td colspan="5"><em>Vela pulsar</em></td></tr><tr><td>PSR J0835-4510</td><td>984.9</td><td>≫8σ</td><td>301</td><td>PWL+PLEC</td></tr><tr><td colspan="5"><em>Hypothetical Vela-like pulsars</em></td></tr><tr><td>PSR J1413-6205</td><td>855.2</td><td>≫8σ</td><td>863</td><td>PLEC</td></tr><tr><td>PSR J1709-4429</td><td>135.6</td><td>≫8σ</td><td>110</td><td>PLEC</td></tr><tr><td>PSR J1732-3131</td><td>258.7</td><td>≫8σ</td><td>60</td><td>PLEC</td></tr><tr><td>PSR J1813-1246</td><td>74.8</td><td>7.4σ</td><td>249</td><td>PLEC</td></tr><tr><td>PSR J1952+3252</td><td>578.4</td><td>≫8σ</td><td>109</td><td>PLEC</td></tr><tr><td>PSR J2021+3651</td><td>435.8</td><td>≫8σ</td><td>139</td><td>PLEC</td></tr><tr><td>PSR J2021+4026</td><td>55.6</td><td>6.3σ</td><td>406</td><td>PLEC</td></tr></tbody></table>

To identify a source as a PeVatron we used the PeVatron Test Statistic [PTS, 95] which is a likelihood ratio test with the null hypothesis corresponding to a spectral model with proton energy cutoff equal to 1 PeV. Following the methodology in [95], a source is identified as a PeVatron if one can exclude the null hypothesis with a statistical significance of 5σ which corresponds to a threshold value of the PeVatron Test Statistic (PTS) of 35.6, taking into account the trials factor in the analysis of 231 sources [95, 96]. If a source is not identified as a PeVatron (PTS < 35.6), it is considered as a PeVatron candidate when it presents a 95% CL lower limit on the proton energy cutoff ($LL_p$) higher than 1 PeV.

We found 32 PeVatron candidates and we ranked them according to $LL_p$ (see table 5). Most of these candidates are known sources that have been simulated with a power-law model (see appendix A.3), with the exception of HESS J1813-178 which has been simulated with an exponential cutoff power-law model with a gamma-ray energy cutoff at 127 TeV. We identified three PeVatron sources ($PTS > 35.6$, reported in bold in table 5) which are among the brightest sources in the list of the PeVatron candidates.

We note that the two known sources Westerlund 1 and HESS J1641-463 were simulated with a hard power-law spectrum but were not identified as PeVatron candidates in our analysis. Both sources are located in complex regions, therefore the results are affected by source confusion. Furthermore, Westerlund 1 was modelled in the catalogue analysis as multiple objects that are not associated to the complex template simulated according to the strict association criteria that we adopted. We conclude that a more sophisticated analysis method is needed to identify PeVatron candidates with complex morphologies and affected by source confusion. This is beyond the scope of this paper and left for future work.
We note that our test to define a PeVatron candidate is based on the derived 95% CL lower limit on the underlying proton cutoff energy. However, the test does not differentiate between hadronic and leptonic emission mechanisms. We notice, indeed, the presence of seven synthetic PWNe appearing as PeVatron candidates. The presence of leptonic sources as PeVatron candidates is not surprising since leptonic accelerators can have spectra extending to very high energies, as suggested also by recent LHAASO measurements [14, 97]. Therefore, it is clear that the PeVatron identification methodology discussed in [95] and used in our work needs to be expanded in order to properly distinguish between leptonic and hadronic PeVatrons. However, this is beyond the scope of this paper and left for future work.

Our results show that bright PeVatrons can be detected in the GPS dataset itself, while the survey of the entire Galactic disk provides an ideal pathfinder to pinpoint candidates to follow-up with deeper observations using CTA.

# 7 Diffuse emission

The study of Galactic diffuse gamma-ray emission with IACTs is challenging as, over large scales, the signal is largely dominated by the CR background. Moreover, diffuse emission is comprised of a contribution from unresolved source populations and from interstellar emission from CR interactions in the Galaxy. In this section we will show that: 1) we can extract from the source catalogue an estimate of the contribution from unresolved sources based on minimal modelling assumptions, and 2) CTA will make it possible to detect TeV interstellar emission and to statistically distinguish between different scenarios.

## 7.1 Unresolved sources

Unresolved sources, i.e. sources that individually fall below the detection threshold of the measurement, make a cumulative contribution to the measurable diffuse emission signal. For the determination of the amount of this contribution, a model of the entire gamma-ray population must first be developed based on the catalogue of detected sources. In a second step, source populations can be divided into detectable and unresolved sources by applying the detection threshold of the measurement. For the development of a source model, two different approaches can be followed. In the first approach, a single source class or several source classes are modelled based on the underlying physics and the model is verified by comparison with the detected sources (as was done in the construction of our sky model in section 2). In the second, more data-driven approach, fitting based on minimal modelling assumptions of the detected sources results in a model for a generic gamma-ray source population. Here, we adopted the latter approach to describe the gamma-ray source population based on the catalogue B described in section 4. This allows us to assess our capability to characterise emission from unresolved sources independently from our prior knowledge of the true sky model (described in section 2). Specifically, we adopt the simple source population model from [98].

The detection threshold in flux for extended sources as function of sky position, $F_{\text{min}}(l, b, \sigma_{\text{source}})$, is derived by scaling the detection threshold in flux for point sources (figure 4), $F_{\text{min},0}$, such as:

$$F_{\text{min}}(l, b, \sigma_{\text{source}}) = F_{\text{min},0}(l, b) \sqrt{\frac{\sigma_{\text{source}}^2 + \sigma_{\text{min}}^2}{\sigma_{\text{min}}^2}} \quad (7.1)$$
<table><caption>**Table 5.** List of the best PeVatron candidates ($LL_p > 1$ PeV). The differential flux values are given in cm$^{-2}$ s$^{-1}$ TeV$^{-1}$. The spectral index is obtained from the gamma-ray PLEC fit and is generally in very good agreement with the simulated one. $LL_p$ is the lower limit on the cut-off energy of the proton spectrum. The PTS is the test statistic used to identify PeVatron candidates as explained in the text. Bold text highlights the three identified PeVatron sources (with PTS> 35.6).</caption><thead><tr><th>Source Name</th><th>Differential flux at 1 TeV</th><th>Spectral Index</th><th>LL<sub>p</sub> [TeV]</th><th>PTS</th></tr></thead><tbody><tr><td><strong>HESS J1841-055</strong></td><td>1.364 ± 0.008 × 10<sup>-11</sup></td><td>2.426 ± 0.004</td><td>22931</td><td>119</td></tr><tr><td>HESS J1800-240C</td><td>6.20 ± 0.16 × 10<sup>-13</sup></td><td>2.608 ± 0.018</td><td>9639</td><td>18</td></tr><tr><td><strong>2HWC J1837-065</strong></td><td>2.40 ± 0.01 × 10<sup>-11</sup></td><td>2.905 ± 0.002</td><td>7397</td><td>42</td></tr><tr><td>HESS J1708-410</td><td>5.83 ± 0.05 × 10<sup>-12</sup></td><td>2.605 ± 0.005</td><td>7079</td><td>31</td></tr><tr><td>2HWC J1902+048*</td><td>4.17 ± 0.04 × 10<sup>-12</sup></td><td>3.290 ± 0.004</td><td>5592</td><td>15</td></tr><tr><td>HESS J1834-087</td><td>2.72 ± 0.03 × 10<sup>-12</sup></td><td>2.62 ± 0.008</td><td>5423</td><td>23</td></tr><tr><td><strong>HESS J1614-518</strong></td><td>8.47 ± 0.04 × 10<sup>-12</sup></td><td>2.425 ± 0.004</td><td>3660</td><td>44</td></tr><tr><td>pwn1772</td><td>2.59 ± 0.16 × 10<sup>-13</sup></td><td>2.50 ± 0.05</td><td>3123</td><td>2</td></tr><tr><td>SNR G323.7-1.0</td><td>2.81 ± 0.04 × 10<sup>-12</sup></td><td>2.491 ± 0.011</td><td>3101</td><td>12</td></tr><tr><td>pwn438</td><td>1.39 ± 0.11 × 10<sup>-13</sup></td><td>1.93 ± 0.07</td><td>3006</td><td>1</td></tr><tr><td>HESS J1018-589 A</td><td>4.19 ± 0.19 × 10<sup>-13</sup></td><td>2.13 ± 0.03</td><td>2438</td><td>8</td></tr><tr><td>pwn521</td><td>5.92 ± 0.15 × 10<sup>-13</sup></td><td>2.733 ± 0.019</td><td>2337</td><td>6</td></tr><tr><td>pwn813</td><td>1.78 ± 0.02 × 10<sup>-12</sup></td><td>2.307 ± 0.009</td><td>2217</td><td>14</td></tr><tr><td>Westerlund 2</td><td>2.90 ± 0.06 × 10<sup>-12</sup></td><td>2.631 ± 0.015</td><td>2064</td><td>5</td></tr><tr><td>2HWC J1914+117*</td><td>2.13 ± 0.03 × 10<sup>-12</sup></td><td>2.853 ± 0.008</td><td>2007</td><td>7</td></tr><tr><td>pwn2733</td><td>1.97 ± 0.15 × 10<sup>-13</sup></td><td>2.26 ± 0.05</td><td>1994</td><td>5</td></tr><tr><td>HESS J1503-582</td><td>2.37 ± 0.06 × 10<sup>-12</sup></td><td>2.399 ± 0.019</td><td>1974</td><td>6</td></tr><tr><td>pwn2252</td><td>1.09 ± 0.14 × 10<sup>-13</sup></td><td>1.85 ± 0.07</td><td>1959</td><td>1</td></tr><tr><td>HESS J1119-614</td><td>1.60 ± 0.04 × 10<sup>-12</sup></td><td>2.637 ± 0.017</td><td>1930</td><td>5</td></tr><tr><td>W 51C</td><td>1.02 ± 0.02 × 10<sup>-12</sup></td><td>2.606 ± 0.018</td><td>1913</td><td>5</td></tr><tr><td>HESS J1844-030</td><td>4.49 ± 0.15 × 10<sup>-13</sup></td><td>2.57 ± 0.02</td><td>1900</td><td>5</td></tr><tr><td>2HWC J1819-150*</td><td>1.596 ± 0.008 × 10<sup>-11</sup></td><td>2.89 ± 0.07</td><td>1780</td><td>12</td></tr><tr><td>pwn2913</td><td>8.3 ± 1.1 × 10<sup>-13</sup></td><td>1.953 ± 0.005</td><td>1616</td><td>4</td></tr><tr><td>isnr99</td><td>9.92 ± 0.06 × 10<sup>-12</sup></td><td>2.412 ± 0.005</td><td>1494</td><td>13</td></tr><tr><td>2HWC J1852+013*</td><td>5.21 ± 0.05 × 10<sup>-12</sup></td><td>2.906 ± 0.005</td><td>1364</td><td>5</td></tr><tr><td>HESS J1832-093</td><td>5.92 ± 0.17 × 10<sup>-13</sup></td><td>2.63 ± 0.02</td><td>1313</td><td>4</td></tr><tr><td>pwn2934</td><td>1.37 ± 0.13 × 10<sup>-13</sup></td><td>2.11 ± 0.05</td><td>1232</td><td>3</td></tr><tr><td>HESS J1804-216</td><td>6.09 ± 0.07 × 10<sup>-12</sup></td><td>2.721 ± 0.007</td><td>1230</td><td>4</td></tr><tr><td>CTB 37B</td><td>6.42 ± 0.19 × 10<sup>-13</sup></td><td>2.675 ± 0.019</td><td>1218</td><td>3</td></tr><tr><td>2HWC J1921+131</td><td>1.766 ± 0.015 × 10<sup>-12</sup></td><td>2.743 ± 0.009</td><td>1186</td><td>4</td></tr><tr><td>HESS J1813-178</td><td>2.81 ± 0.02 × 10<sup>-12</sup></td><td>2.133 ± 0.008</td><td>1049</td><td>4</td></tr><tr><td>HESS J1846-029</td><td>6.89 ± 0.15 × 10<sup>-13</sup></td><td>2.356 ± 0.017</td><td>1018</td><td>3</td></tr></tbody></table>
<figure><img src="image_16.png" /><figcaption>**Figure 13.** Number of sources as a function of integral source flux in the 3–100 TeV energy range. We show the distributions for the true sky model, for detected sources (catalogue B), for the population reconstructed blindly from catalogue B, and the corresponding subsample of sources that should have been detected. The values shown for the reconstructed population consist of the mean from 1000 realisations with error bars spanning the 25% to 75% quantiles of the distribution.</figcaption></figure>

for source sizes $\sigma_{\text{source}} \le 2^\circ$ in radius. Above an angular extension of $2^\circ$ the source is assumed to be undetectable given the applied background subtraction technique. The quantity $\sigma_{\min}$ represents the minimum resolvable source size due to the instrument resolution and characteristics of the observations. For this exercise it is fixed to a value of $0.052^\circ$, which is the smallest angular extension of any source found in catalogue B.

The free parameters of the source population were then determined by a likelihood fitting procedure (based on Poisson statistics) to the sources detected in catalogue B. For all sources in the true sky model it is assumed that the source distance is known and therefore the flux and angular extent can be converted into the source luminosity and radius, respectively.

Simulation of gamma-ray observations of the source population so characterised makes it possible to apply the detection threshold so as to divide the population into detectable sources and unresolved sources for the two energy ranges below and above 1 TeV considered in catalogue B (section 4). To estimate the flux of unresolved sources and its dispersion, we simulate 1000 realisations of the source population. The final maps of unresolved sources in the two energy ranges are obtained by taking a bin-wise average of the multiple realisation. An example map is shown in appendix E. The fluxes and distributions of the detectable sources agree well with the sources from the catalogue. We show, as an example, the longitude distribution of the integral flux > 1 TeV in figure 13.
In spite of the relatively simple assumptions, which are different from those used in the generation of the true sky model, one can see that the reconstructed population agrees well with the catalogue results, and that we obtain from the data a reasonable description of sub-threshold sources in agreement with the true sky model. We cannot reproduce all details of the complex true sky model, in particular in the range of very low fluxes ($< 10^{-15} \text{ cm}^{-2} \text{ s}^{-1}$) and small source numbers, but this has a negligible impact on the properties of the collective emission from unresolved sources. Conversely, we see that, for the realisation considered, catalogue B does not perfectly capture the flux distribution of the brightest sources in the true sky model due to a combination of a few bright sources lying at the edge of the survey region, notably the Crab Nebula, and because of fragmentation effects for sources with complex morphology discussed earlier. Nevertheless, the methodology compensates, at least partially, for these effects, and the reconstructed population matches the flux distribution of the true sky model at high fluxes.

The estimated contribution of unresolved sources is shown in figure 14, together with the contribution of the true-source model, the three alternative interstellar models, and the CR background. We can see that, for the models considered, the integral flux of the unresolved sources is smaller by an order of magnitude than interstellar emission. Overall the level of the predicted flux matches that of the true sources below the catalogue threshold (TS<25). The differences between the true sky model below the detection threshold and the unresolved source template extracted from the catalogue are expected as the mean template derived from multiple realisations cannot match exactly a single observed sky. These deviations are larger than statistical fluctuations on the number of detected photons. They will need to be taken into account when interpreting real observations, along with systematic uncertainties not dealt with in the present work.

## 7.2 Interstellar emission

In this section we summarise our results on interstellar emission. For the production of the catalogue the sources' parameters were fitted to the data together with several adjustable parameters for the true CR background and IEM models. Namely, the Galactic plane was split in overlapping regions, each spanning $10^\circ$ in longitude and $12^\circ$ in latitude, with their centres shifted in steps of $\pm 5^\circ$. In each of these regions the CR background and IEM were each adjusted via a free normalisation and a power-law correction to the spectrum with free spectral index. The comparison of the different emission components in the fitted and true models for catalogue B is shown in figure 15. We note that the relative deviation from the true flux is larger for the interstellar emission than for the sources, because most sources are detectable on smaller scales and are therefore easier to spatially distinguish from the CR background.

In order to illustrate how the measurement of the interstellar emission is affected by the confusion with the CR background, we show in figure 16 the Pearson correlation coefficient between the predicted counts of the CR background and the IEM (for two different models and energy bands), the TS of the IEM, and the fitted IEM normalisations as a function of the Galactic longitude.<fnref n="11" /> In the outer Galaxy ($|l| > 90^\circ$), where the interstellar emission is fainter and extends over a larger scale height, the TS is lower and the correlation coefficient

<footnote n="11">¹¹We note that all models are fit to data simulated based on the IEM-base model.</footnote>
<figure><img src="image_17.png" /><figcaption>**Figure 14.** Integral flux profile in the 3–100 TeV energy range as a function of Galactic longitude for the CR background, for different variations of the interstellar emission model, for the true source models above and below the detection threshold of the catalogue, and for the unresolved sources template reconstructed from catalogue B (see section 7.1). Fluxes are integrated over latitudes of ±6° and over a 6° sliding window in longitude.</figcaption></figure>

is larger (as there are less counts and they are more sparsely distributed). In this case the normalisation of the IEM is poorly constrained. If we include energies below 1 TeV we find an increasing bias of the IEM normalisation with longitude associated with a steeper decrease of the value of the TS with longitude and a global increase of the correlation coefficient. This can be explained by the degradation of the instrument performance at lower energies. Indeed the higher CR background lowers the TS, while the broader PSF and energy dispersion increase the confusion. In the $|l| < 60°$ longitude range and $> 1$ TeV energy range, the interstellar emission can be significantly detected and the normalisation of the model can be correctly determined from the data for the IEM-base model used in the simulation.

The normalisation of the IEM-varmin model, fitted to the data simulated based on IEM-base model, is compatible with the flux ratio of the two models in the region where interstellar emission is significantly detected. Moreover, the difference in flux between IEM-base and the two IEM-varmin variants is larger than the statistical error on the fitted model and its deviation from the true model (figure 15 and 16), which suggests that CTA should be able to distinguish between such scenarios.

However, the systematic uncertainties on the CR background modelling may further complicate the exercise. Even if the error on the CR background does not exceed a few percent, it can lead to several tens of percent error on the IEM. We performed a simple test
<figure><img src="image_18.png" /><figcaption>**Figure 15.** Flux distribution in Galactic longitude from different source and background components. Fluxes are integrated over latitudes of ±6° and over a 6° sliding window in longitude for the 3–100 TeV energy range. The fitted models are displayed as solid lines and the simulated models as dotted lines.</figcaption></figure>

by flipping the CR background template in latitude, which leads to a mean relative error on the background of 2%, with a standard deviation of 40%. In this case the longitude range where the fitted normalisation of the interstellar model remains compatible with the true values is reduced, but the results for the innermost region ($|l| < 60°$) remain similar.

Several questions on the effect of systematic uncertainties on interstellar emission studies are still open. Notably, the CR background normalisation usually has to be re-fitted for each observation outside of the region of interest, which is not well-defined for interstellar emission as it spans over the whole field-of-view. Alternatively, the CR background of each observation could be jointly fitted together with the interstellar emission and source models but this would be computationally intensive. Moreover one could investigate the effect of potential systematic errors in the IRFs and CR background modelling. We defer further investigations on these questions to future studies. In this domain we expect important synergies between CTA and neutrino observatories [99–102]. Synergies with HAWC, LHAASO and Fermi-LAT will also be crucial as these instruments have a much better handle on the CR background modelling.

# 8 Summary and perspectives

A survey of the entire Galactic plane has been proposed as a Key Science Project for CTAO. In this article we presented a snapshot of the current status of the project preparation and predictions of the expected results. With respect to previous assessments of the CTA GPS we employed an improved sky model (made publicly available along with this paper)
<figure><img src="image_19.png" /><figcaption>**Figure 16.** Longitude profiles showing from top to bottom: 1) Pearson correlation coefficient (PCC) of the predicted counts from interstellar emission and CR background; 2) $\sqrt{TS}$ of the interstellar emission models; 3) normalisation of the interstellar emission models fitted to data simulated using the IEM-base model. All these quantities are computed in 8° by 12° windows in Galactic longitude and latitude, respectively.</figcaption></figure>
that combines data from recent observations of known gamma-ray emitters with physically-
driven models of synthetic populations of the three main classes of established Galactic VHE
sources (PWNe, young and interacting SNRs, and compact binary systems), as well as of
interstellar emission from CR interactions in the Milky Way. The article also illustrates
the ongoing optimisation of the observation strategy (pointing pattern and scheduling), and
the development and testing of the methods and software tools that will be later used to
build source catalogues from data, to characterise specific source classes (binaries, pulsars,
PeVatrons), and to study diffuse emission.

The approach used has several limitations. First of all, for most of the work we considered
only one single realisation of a given synthetic source population model and of the simulated
data. The effects of the variance of the synthetic populations and the gamma-ray event lists
are not assessed in the paper. Furthermore, uncertainties in the population modelling may
lead to variations in the number of sources detected and the properties of the sources. We
also assumed that the IRFs were perfectly known and used the true background model in
the analysis with simple scaling factors, while with real data the systematic uncertainties
on both are certainly going to complicate the analysis. Note also that in this exercise, the
association of a detected object with its counterpart is done by cross-matching the catalogue
results with the true sky model. For real data, assigning a detected object to a source class
will require the use of multiwavelength data and will be noticeably more difficult. A fraction
of the detected sources may therefore remain unidentified.

Moreover, as discussed in section 3, we presented results for the CTA baseline/Omega
configuration, while so far only the reduced Alpha configuration was approved for construction.
For this reason, and in order to take into account other changes in the instrument configuration
and performance yet to come, as well as advances in the field made possible by other
instruments, the details of the programme proposal are likely to evolve in the coming years
with respect the current snapshot presented in this paper.

Nevertheless, the results presented provide a plausible ballpark estimate of what can be
expected from the CTA GPS for the anticipated observation time of 1620 hours.

* We show that, under our model assumptions and for the realisation considered, the GPS has the potential to increase by a factor five the number of known Galactic VHE emitters.

* In particular, we expect to be able to detect over two hundred PWNe and several tens of SNRs. PWNe should be detected across the entire Milky Way, at average integral fluxes > 1 TeV one order of magnitude lower than in the existing sample, and with good coverage for the parent pulsar intrinsic spindown luminosity > 10³⁶ erg s⁻¹ (50% completeness for the model considered). Similar increases in the number of objects and a reduction in the typical integral flux observable are obtained for SNRs. For about half of the newly-detected SNRs a significant angular extension could be measured, which is valuable for source identification and for constraining physical models of the emission.

* The GPS also has the potential to provide several new VHE detections of gamma-ray binary systems and pulsars, to confirm the existence of a hypothetical population of
gamma-ray pulsars with an additional TeV emission component, and to detect bright PeVatrons with differential fluxes at 1 TeV $\gtrsim 10^{-11}$ cm$^{-2}$ s$^{-1}$ TeV$^{-1}$. Furthermore, the GPS data will provide a pathfinder for deeper follow-up observations of these classes of sources.

* The GPS data will make it possible to detect interstellar emission for all models considered in the $|l| < 60^\circ$ longitude range and $> 1$ TeV energy range, and to statistically distinguish different scenarios. We can also extract from the GPS source catalogue an estimate of the contribution to diffuse emission from unresolved sources based on minimal modelling assumptions.

In addition to improvements in the source detection and characterisation methods, a major development still needed for the scientific exploration of the GPS is to gather multiwavelength/multimessenger data and to develop methods for the statistical association of sources. Indeed, the large number of sources detected will make it impractical to associate objects detected by CTA to emitters known from other wavelengths via visual inspection of each source. Advanced methods already exist in other energy bands [e.g., 103–105], but they must be adapted to the specific attributes of the VHE domain, notably spatial extension and confusion, and exploit as much as possible our understanding of the physical mechanisms underlying the multiwavelength emission.

The GPS will be accompanied by deep observations of individual Galactic sources, either from the population already known or discovered by the GPS itself. Such deep observations will provide fine morphological, spectral, and temporal characterisation of the sources, possibly complemented by a richer set of information on the objects extracted from multiwavelength/multimessenger data and/or detailed modelling [e.g., 106] and will be highly complementary to the population studies made possible by the GPS. Other proposed CTAO Key Science Projects will also address Galactic sources and interstellar emission through observations of the Galactic centre region, transient Galactic sources, the Large Magellanic Cloud, and star-forming systems [24]. The combination of these programs will enable a transformational advance in our understanding of gamma-ray source populations and of the physics of particle acceleration and transport in galaxies.

## Acknowledgments

We gratefully acknowledge financial support from the following agencies and organisations: State Committee of Science of Armenia, Armenia; The Australian Research Council, Astronomy Australia Ltd, The University of Adelaide, Australian National University, Monash University, The University of New South Wales, The University of Sydney, Western Sydney University, Australia; Federal Ministry of Education, Science and Research, and Innsbruck University, Austria; Conselho Nacional de Desenvolvimento Científico e Tecnológico (CNPq), Fundação de Amparo à Pesquisa do Estado do Rio de Janeiro (FAPERJ), Fundação de Amparo à Pesquisa do Estado de São Paulo (FAPESP), Fundação de Apoio à Ciência, Tecnologia e Inovação do Paraná — Fundação Araucária, Ministry of Science, Technology, Innovations and Communications (MCTIC), Brasil; Ministry of Education and Science, National RI Roadmap Project DO1-153/28.08.2018, Bulgaria; The Natural Sciences and Engineering Research Council of Canada and the Canadian Space Agency, Canada; CONICYT-Chile
grants CATA AFB 170002, ANID PIA/APOYO AFB 180002, ACT 1406, FONDECYT-Chile grants, 1161463, 1170171, 1190886, 1171421, 1170345, 1201582, Gemini-ANID 32180007, Chile, W.M. gratefully acknowledges support by the ANID BASAL projects ACE210002 and FB210003, and FONDECYT 11190853; Croatian Science Foundation, Rudjer Boskovic Institute, University of Osijek, University of Rijeka, University of Split, Faculty of Electrical Engineering, Mechanical Engineering and Naval Architecture, University of Zagreb, Faculty of Electrical Engineering and Computing, Croatia; Ministry of Education, Youth and Sports, MEYS LM2015046, LM2018105, LTT17006, EU/MEYS CZ.02.1.01/0.0/0.0/16_013/0001403, CZ.02.1.01/0.0/0.0/18_046/0016007 and CZ.02.1.01/0.0/0.0/16_019/0000754, Czech Republic; Academy of Finland (grant nr.317636 and 320045), Finland; Ministry of Higher Education and Research, CNRS-INSU and CNRS-IN2P3, CEA-Irfu, ANR, Regional Council Ile de France, Labex ENIGMASS, OCEVU, OSUG2020 and P2IO, France; The German Ministry for Education and Research (BMBF), the Max Planck Society, the German Research Foundation (DFG, with Collaborative Research Centres 876 & 1491), and the Helmholtz Association, Germany; Department of Atomic Energy, Department of Science and Technology, India; Istituto Nazionale di Astrofisica (INAF), Istituto Nazionale di Fisica Nucleare (INFN), MIUR, Istituto Nazionale di Astrofisica (INAF-OABRERA) Grant Fondazione Cariplo/Regione Lombardia ID 2014-1980/RST ERC, Italy; ICRR, University of Tokyo, JSPS, MEXT, Japan; Netherlands Research School for Astronomy (NOVA), Netherlands Organization for Scientific Research (NWO), Netherlands; University of Oslo, Norway; Ministry of Science and Higher Education, DIR/WK/2017/12, the National Centre for Research and Development and the National Science Centre, UMO-2016/22/M/ST9/00583, Poland; Slovenian Research Agency, grants P1-0031, P1-0385, I0-0033, J1-9146, J1-1700, N1-0111, and the Young Researcher program, Slovenia; South African Department of Science and Technology and National Research Foundation through the South African Gamma-Ray Astronomy Programme, South Africa; The Spanish groups acknowledge the Spanish Ministry of Science and Innovation and the Spanish Research State Agency (AEI) through the government budget lines PGE2021/28.06.000X.411.01, PGE2022/28.06.000X.411.01 and PGE2022/28.06.000X.711.04, and grants PID2022-139117NB-C44, PID2019-104114RB-C31, PID2019-107847RB-C44, PID2019-104114RB-C32, PID2019-105510GB-C31, PID2019-104114RB-C33, PID2019-107847RB-C41, PID2019-107847RB-C43, PID2019-107847RB-C42, PID2019-107988GB-C22, PID2021-124581OB-I00, PID2021-125331NB-I00; the "Centro de Excelencia Severo Ochoa" program through grants no. CEX2019-000920-S, CEX2020-001007-S, CEX2021-001131-S; the "Unidad de Excelencia María de Maeztu" program through grants no. CEX2019-000918-M, CEX2020-001058-M; the "Ramón y Cajal" program through grants RYC2021-032552-I, RYC2021-032991-I, RYC2020-028639-I and RYC-2017-22665; the "Juan de la Cierva-Incorporación" program through grants no. IJC2018-037195-I, IJC2019-040315-I. They also acknowledge the "Atracción de Talento" program of Comunidad de Madrid through grant no. 2019-T2/TIC-12900; the project "Tecnologiás avanzadas para la exploración del universo y sus componentes" (PR47/21 TAU), funded by Comunidad de Madrid, by the Recovery, Transformation and Resilience Plan from the Spanish State, and by NextGenerationEU from the European Union through the Recovery and Resilience Facility; the La Caixa Banking Foundation, grant no. LCF/BQ/PI21/11830030; the "Programa Operativo" FEDER
2014-2020, Consejería de Economía y Conocimiento de la Junta de Andalucía (ref. 1257737), PAIDI 2020 (ref. P18-FR-1580) and Universidad de Jaén; “Programa Operativo de Crecimiento Inteligente” FEDER 2014-2020 (ref. ESFRI-2017-IAC-12), Ministerio de Ciencia e Innovación, 15% co-financed by Consejería de Economía, Industria, Comercio y Conocimiento del Gobierno de Canarias; the “CERCA” program and the grant 2021SGR00426, both funded by the Generalitat de Catalunya; and the European Union’s Horizon 2020 GA:824064 and NextGenerationEU (PRTR-C17.11); Swedish Research Council, Royal Physiographic Society of Lund, Royal Swedish Academy of Sciences, The Swedish National Infrastructure for Computing (SNIC) at Lunarc (Lund), Sweden; State Secretariat for Education, Research and Innovation (SERI) and Swiss National Science Foundation (SNSF), Switzerland; Durham University, Leverhulme Trust, Liverpool University, University of Leicester, University of Oxford, Royal Society, Science and Technology Facilities Council, U.K.; U.S. National Science Foundation, U.S. Department of Energy, Argonne National Laboratory, Barnard College, University of California, University of Chicago, Columbia University, Georgia Institute of Technology, Institute for Nuclear and Particle Astrophysics (INPAC-MRPI program), Iowa State University, the Smithsonian Institution, V.V.D. is funded by NSF grant AST-1911061, Washington University McDonnell Center for the Space Sciences, The University of Wisconsin and the Wisconsin Alumni Research Foundation, U.S.A.

The research leading to these results has received funding from the European Union's Seventh Framework Programme (FP7/2007-2013) under grant agreements No 262053 and No 317446. This project is receiving funding from the European Union's Horizon 2020 research and innovation programs under agreement No 676134.

This research made use of ctools, a community-developed gamma-ray astronomy science analysis software. ctools is based on GammaLib, a community-developed toolbox for the scientific analysis of astronomical gamma-ray data. This research made use of gammapy,<fnref n="12" /> a community-developed core Python package for TeV gamma-ray astronomy.

# A Details on the construction of the sky model

## A.1 Pulsar selection and models

We preselected from among more than 250 gamma-ray pulsars detected with the Fermi-LAT [27, 28] those whose integral energy flux exceeds 1% of that of the brightest pulsar, Vela (PSR J0835-4510), in the energy range from 1–10 GeV. This amounted to 38 pulsars, including the four brightest GeV gamma-ray pulsars detected with IACTs: Vela, Geminga (PSR J0633+1746), Crab (PSR J0534+2200) and PSR J1709–4429.

Two classes of gamma-ray emitting pulsars were constructed using the information on the VHE spectra of the Crab and Vela pulsars known to date. The first class (“Crab-like” sources, 25 objects; see table 6) comprises the sources with a spectrum consistent with a power law (PWL) in the high-energy (HE) range from one to a few hundred GeV (called hereafter the “GeV component” for the sake of simplicity). Their VHE spectrum above 100 GeV was modelled simply extrapolating the “GeV component” without spectral cutoff, as observed in the Crab pulsar (see e.g. [107]). This approach is similar to the one adopted

<footnote n="12">¹²https://www.gammapy.org.</footnote>
in [108]. The spectral parameters of Crab-like pulsars were taken from 3FHL whereas, for the Crab pulsar, we adopted the power-law spectrum of the peak emission (P1+P2) calculated for the joint Fermi-LAT and MAGIC spectral data as reported in [107].

The second class (“Vela-like” sources, 13 objects; see table 6) comprises the sources with a Fermi-LAT HE spectrum consistent with a power-law with an exponential cutoff (PLEC) model. Their VHE spectra were modelled by summing the “GeV component” with an additional multi-TeV pulsed component (called “TeV component” hereafter), as observed in the Vela pulsar [94]. The parameters of the GeV component were taken from [27] (2PC), except for the Vela and Geminga pulsars for which we used the phase-averaged spectra from [109] and [110], respectively. Although no other pulsar with multi-TeV emission is known to date besides Vela, a TeV component was included in the spectral models for the 12 remaining objects in the same class. For these additional components the spectral index was fixed according to the measurements below 1 GeV with the Fermi-LAT, and the high-energy cutoff was randomly set following a normal distribution centred at 7 TeV with a σ of 3 TeV. The flux of the TeV component was randomly generated assuming a TeV to GeV flux ratio η<sub>TeV/GeV</sub> ranging from 5 × 10<sup>-5</sup> to 10<sup>-2</sup>, while complying at the same time with upper-limits available in the literature [111, 112].

The pulse shapes of both classes of pulsars were generated using templates at energies above 10 GeV made available by the Fermi-LAT team [113]. No evolution with energy of the phasograms was considered. We modelled the frequency as a function of time as a third-order polynomial. The coordinates and timing properties (reference epoch, frequency and its derivatives) for most of the pulsars were taken from [114].<sup>13</sup> For 4 pulsars that are missing in [114] (PSRs J1119−6127, J1648−4611, J1838−0537 and J2215+5135) we used the position and timing parameters reported in the Fermi 2PC [27].

## A.2 Other dedicated models

For gamma-ray binaries, periods and relevant references are listed in table 7. When the sky model was built the orbital period of HESS J1832−093 was unknown and in our sky model we assumed it to be equal to 1 year similar to HESS J0632+057. Now a periodicity of 86 days has been reported for its X-ray counterpart [115]. The binary nature of one more gamma-ray source, 4FGL J1405.1−6119, was discovered when our work was already in an advanced state and thus it is not included in the sky model.

For some well studied extended sources, a spatial template providing a more complex morphology is used in the simulation. Templates are also used to model some diffuse features not captured by the large-scale interstellar emission models. The templates are based either on multi-wavelength observations or models as described in table 8. Some of the templates present an energy-dependent morphology.

## A.3 Spectral cutoffs

For a number of sources the spectrum is not well constrained above 30 TeV with current generation telescopes. We impose an exponential cutoff for known sources described in

<footnote n="13"><sup>13</sup>See also https://confluence.slac.stanford.edu/display/GLAMCOG/LAT+Gamma-ray+Pulsar+Timing+Models.</footnote>
<table><caption>**Table 6.** Pulsars included in the CTA GPS simulations. The spectral properties of the GeV and TeV band are chosen to be either a power-law (PWL), power-law with exponential cutoff (PLEC) or a PWL with super exponential cutoff (super PLEC).</caption><thead><tr><th>Pulsar Name</th><th>GeV model</th><th>TeV model</th></tr></thead><tbody><tr><td colspan="3"><em>Crab-like pulsars</em></td></tr><tr><td>PSR J0534+2200 (Crab)</td><td>PWL</td><td>not present</td></tr><tr><td>PSR J0633+0632</td><td>PWL (3FHL)</td><td>not present</td></tr><tr><td>PSR J1016-5857</td><td>PWL (3FHL)</td><td>not present</td></tr><tr><td>PSR J1028-5819</td><td>PWL (3FHL)</td><td>not present</td></tr><tr><td>PSR J1048-5832</td><td>PWL (3FHL)</td><td>not present</td></tr><tr><td>PSR J1119-6127</td><td>PWL (3FHL)</td><td>not present</td></tr><tr><td>PSR J1459-6053</td><td>PWL (3FHL)</td><td>not present</td></tr><tr><td>PSR J1509-5850</td><td>PWL (3FHL)</td><td>not present</td></tr><tr><td>PSR J1514-4946</td><td>PWL (3FHL)</td><td>not present</td></tr><tr><td>PSR J1620-4927</td><td>PWL (3FHL)</td><td>not present</td></tr><tr><td>PSR J1648-4611</td><td>PWL (3FHL)</td><td>not present</td></tr><tr><td>PSR J1747-2958</td><td>PWL (3FHL)</td><td>not present</td></tr><tr><td>PSR J1803-2149</td><td>PWL (3FHL)</td><td>not present</td></tr><tr><td>PSR J1809-2332</td><td>PWL (3FHL)</td><td>not present</td></tr><tr><td>PSR J1826-1256</td><td>PWL (3FHL)</td><td>not present</td></tr><tr><td>PSR J1833-1034</td><td>PWL (3FHL)</td><td>not present</td></tr><tr><td>PSR J1838-0537</td><td>PWL (3FHL)</td><td>not present</td></tr><tr><td>PSR J1907+0602</td><td>PWL (3FHL)</td><td>not present</td></tr><tr><td>PSR J1954+2836</td><td>PWL (3FHL)</td><td>not present</td></tr><tr><td>PSR J1958+2846</td><td>PWL (3FHL)</td><td>not present</td></tr><tr><td>PSR J2032+4127</td><td>PWL (3FHL)</td><td>not present</td></tr><tr><td>PSR J2111+4606</td><td>PWL (3FHL)</td><td>not present</td></tr><tr><td>PSR J2215+5135</td><td>PWL (3FHL)</td><td>not present</td></tr><tr><td>PSR J2229+6114</td><td>PWL (3FHL)</td><td>not present</td></tr><tr><td>PSR J2238+5903</td><td>PWL (3FHL)</td><td>not present</td></tr><tr><td colspan="3"><em>Vela-like pulsars</em></td></tr><tr><td>PSR J0007+7303</td><td>PLEC (2PC)</td><td>PLEC</td></tr><tr><td>PSR J0633+1746 (Geminga)</td><td>super PLEC</td><td>PLEC</td></tr><tr><td>PSR J0835-4510 (Vela)</td><td>super PLEC</td><td>PLEC</td></tr><tr><td>PSR J1057-5226</td><td>PLEC (2PC)</td><td>PLEC</td></tr><tr><td>PSR J1413-6205</td><td>PLEC (2PC)</td><td>PLEC</td></tr><tr><td>PSR J1418-6058</td><td>PLEC (2PC)</td><td>PLEC</td></tr><tr><td>PSR J1709-4429</td><td>super PLEC (2PC)</td><td>PLEC</td></tr><tr><td>PSR J1732-3131</td><td>PLEC (2PC)</td><td>PLEC</td></tr><tr><td>PSR J1813-1246</td><td>PLEC (2PC)</td><td>PLEC</td></tr><tr><td>PSR J1846+0919</td><td>PLEC (2PC)</td><td>PLEC</td></tr><tr><td>PSR J1952+3252</td><td>super PLEC (2PC)</td><td>PLEC</td></tr><tr><td>PSR J2021+3651</td><td>super PLEC (2PC)</td><td>PLEC</td></tr><tr><td>PSR J2021+4026</td><td>super PLEC (2PC)</td><td>PLEC</td></tr></tbody></table>
<table><caption>**Table 7.** Known and candidate gamma-ray binaries within the GPS footprint.</caption><thead><tr><td>Name</td><td>Period</td><td>References</td></tr></thead><tbody><tr><td>PSR B1259-63</td><td>3.4 yr</td><td>[116]</td></tr><tr><td>PSR J2032+4127</td><td>50 yr</td><td>[117]</td></tr><tr><td>LSI +61°303</td><td>26.5 d</td><td>[118]</td></tr><tr><td>LS 5039</td><td>3.9 d</td><td>[119]</td></tr><tr><td>1FGL J1018.6-5856</td><td>16.6 d</td><td>[120]</td></tr><tr><td>HESS J0632+057</td><td>315 d</td><td>[121]</td></tr><tr><td>HESS J1832-093</td><td>1 yr</td><td>[122]</td></tr></tbody></table>

<table><caption>**Table 8.** Sources modelled by dedicated spatial templates and references.</caption><thead><tr><td>Name</td><td>Source class</td><td>Template type</td><td>References</td></tr></thead><tbody><tr><td>IC 443</td><td>iSNR</td><td>gamma-ray map</td><td>[123]</td></tr><tr><td>SNR G78.2+2.1</td><td>SNR</td><td>gamma-ray map</td><td>[9, 124]</td></tr><tr><td>Vela X</td><td>PWN</td><td>radio and X-ray maps</td><td>[125-128]</td></tr><tr><td>W 28</td><td>iSNR</td><td>gamma-ray map</td><td>[129]</td></tr><tr><td>Pup A</td><td>SNR</td><td>X-ray map</td><td>[130, 131]</td></tr><tr><td>Vela Junior</td><td>SNR</td><td>X-ray map</td><td>[130]</td></tr><tr><td>HESS J1800-240</td><td>iSNR</td><td>gamma-ray map</td><td>[129]</td></tr><tr><td>RX J1713.7-3946</td><td>SNR</td><td>X-ray map</td><td>[132]</td></tr><tr><td>HESS J1825</td><td>PWN</td><td>interpolation of gamma-ray profiles</td><td>[133]</td></tr><tr><td>SS 433 lobes</td><td>microquasar lobes</td><td>two Gaussians</td><td>[134]</td></tr><tr><td>Geminga halo</td><td>pulsar halo</td><td>electron diffusion</td><td>[135]</td></tr><tr><td>PSR B0656+14 halo</td><td>pulsar halo</td><td>electron diffusion</td><td>[135]</td></tr><tr><td>Cygnus cocoon</td><td>SFR</td><td>1/r nuclei profile + gas map</td><td>[31, 136]</td></tr><tr><td>Westerlund 1</td><td>SFR</td><td>1/r nuclei profile + gas map</td><td>[31, 137]</td></tr><tr><td>Galactic centre ridge</td><td>diffuse feature</td><td>gamma-ray map</td><td>[138]</td></tr><tr><td>Fermi bubbles</td><td>diffuse feature</td><td>gamma-ray map</td><td>[139]</td></tr></tbody></table>

catalogues (gamma-cat, 2HWC, 3FHL) using a power-law spectrum with hard spectral index (< 2.4), with the exception of two PeVatron candidates, Westerlund 1 [31] and HESS J1641-463 [32]. This cutoff is added to prevent generating too many artificial PeVatrons due to a power-law extrapolation. The cutoff energy is computed based on physically motivated principles described in the following text for sources associated with SNRs, PWNe, and AGNs. For each catalogue considered, unidentified sources are treated as if they belonged to the dominant source class (PWNe for gamma-cat and 2HWC, AGNs for 3FHL).
SNR properties are extracted from SNRcat [140]. Thermal composite SNRs and SNRs interacting with molecular clouds are considered to be produced by type II supernovae and their gamma-ray emission to be dominated by hadronic processes. Other SNRs are randomly assigned to be type I with 20% probability or type II with 80% probability, and their emission is assumed to be dominated by leptonic processes. When an age estimate is not available the age is inferred from the size of the SNR based on the evolutionary model described in [141]. If the age estimated is unrealistically large (> 10 kyr), we assume an age drawn randomly between 0.5 kyr and 1 kyr. The maximum energy of accelerated particles is then estimated based on equation 9 from [141]. For leptonic emitters the age-limited maximum energy for electrons is corrected to account for synchrotron energy losses. To do so the strength of the amplified magnetic field downstream is computed by assuming that the non-resonant instability reaches saturation upstream of the shock. The gamma-ray cutoff energy is then estimated to be 10% of the maximum proton energy or electron energy (in the deep Klein-Nishina regime).

For PWNe we assume that the associated pulsar is the one at the smallest angular distance with spin-down power > $10^{34}$ erg/s found in the ATNF catalogue [142]. Based on the pulsar properties, the maximum electron energy is computed as the minimum between the maximum potential drop available in the pulsar magnetosphere and the synchrotron-loss-limited maximum energy based on a model for the termination shock magnetic field evolution [143]. As for leptonic SNRs, the gamma-ray cutoff energy is then estimated as 10% of the maximum electron energy.

For AGNs the gamma-ray cutoff energy is set as a function of redshift using an analytical approximation of the maximum energy shown in figure 17 of [28], consistent with expectations for absorption by the extragalactic background light. For AGNs of unknown redshift or unassociated 3FHL sources, a redshift value is randomly drawn following the distribution of 3FHL AGNs with known redshift.

The cutoff energies are shown in figure 17. They mostly lie below 1 TeV for 3FHL sources and below 100 TeV for gamma-cat and 2HWC sources.

## A.4 Removal of bright synthetic sources

To remove bright synthetic sources we consider the following source properties:

* Galactic longitude, *l* (of the centroid for extended sources);

* Galactic latitude, *b* (of the centroid for extended sources);

* angular extension radius, *r* (maximum radius or 99% containment radius for models with non-zero intensities up to arbitrarily large distances from the centre such as Gaussians);

* integral photon flux > 1 TeV, *f*.

Each detected object is compared to synthetic sources belonging to the same class: SNRs, iSNRs, PWNe, composites, binaries. We note that, for this purpose, we identify a subsample of young SNRs and PWNe that are deemed to be observed as composite SNR/PWN systems. We do not use the information that SNR and PWN originate from the same progenitor but,
<figure><img src="image_20.png" /><figcaption>**Figure 17.** Cutoff energies estimated for sources described in catalogues (gamma-cat, 2HWC, 3FHL) using a power-law spectrum with hard spectral index (< 2.4). The histogram is colour-coded according to the source class: pulsar wind nebula (PWN), supernova remnants (SNR), active galactic nucleus (AGN) or unidentified (UNID). See text for details on how the cutoff energies are calculated (appendix A.3).</figcaption></figure>

as for real objects, we consider as composite systems the pairs of SNR and PWN that show (partial) overlap on the sky. The longitude/latitude of the composite system is assumed to be the flux-weighted average of the longitudes/latitudes of the individual objects. The extension radius is the radius of the smallest circle that encloses both objects, and the flux is the sum of the fluxes.

For each detected object we compute a figure of merit to compare it to synthetic objects of the same class. For the $i^{th}$ synthetic source the figure of merit is defined as:

$$F_i = \sqrt{\left(\frac{l_i - l}{u}\right)^2 + \left(\frac{b_i - b}{v}\right)^2 + \left(\frac{1}{w} \frac{r_i - r}{r_i + r}\right)^2 + \frac{1}{z^2} \log_{10}^2\left(\frac{f_i}{f}\right)} \quad (A.1)$$

with $u = 180^{\circ}$, $v = 10^{\circ}$, $w = 2$, $z = 0.6$ if $f_i < f$ or $z = 0.3$ if $f_i > f$. The parameter values were determined empirically so that there was reasonable consistency between the spatial, angular extension, and integral flux distribution for the synthetic source populations and the merged population. We note that:

* longitude and latitude are considered independently to preserve the population distributions for both quantities;

* values of $r_i$ or $r$ lower than 0.05° are set to 0.05° for the calculation of $F_i$ to account for the limited resolution of existing observations, based on the minimum resolved size in [4];

* the asymmetric value of $z$ around $f_i = f$ ensures that brighter synthetic sources, which could have been detected more easily, are more likely to get removed.
<table><caption>**Table 9.** Fractions of the total observing times as a function of year separated for STP and LTP.</caption><thead><tr><th>Programme</th><th colspan="2">STP</th><th colspan="9">LTP</th></tr><tr><th>Year</th><th>1</th><th>2</th><th>3</th><th>4</th><th>5</th><th>6</th><th>7</th><th>8</th><th>9</th><th>10</th></tr></thead><tbody><tr><th>Fraction</th><td>0.6</td><td>0.4</td><td>0.15</td><td>0.2</td><td>0.2</td><td>0.15</td><td>0.1</td><td>0.1</td><td>0.05</td><td>0.05</td></tr></tbody></table>

<table><caption>**Table 10.** Length of the gap in days between nights over which GPS observations are performed.</caption><thead><tr><th>Year</th><th>1</th><th>2</th><th>3</th><th>4</th><th>5</th><th>6</th><th>7</th><th>8</th><th>9</th><th>10</th></tr></thead><tbody><tr><th>South</th><td>1.5</td><td>3</td><td>3</td><td>1.5</td><td>2.5</td><td>5</td><td>7</td><td>10</td><td>12</td><td>14</td></tr><tr><th>North</th><td>4</td><td>5</td><td>5</td><td>4</td><td>4</td><td>5</td><td>6</td><td>8</td><td>12</td><td>13</td></tr></tbody></table>

The synthetic source with the smallest $F_i$ is considered as the most similar to the detected source and therefore removed.

## B Details on observation scheduling

We use the total observing times in several longitude ranges as provided in table 6.3 of [24]. In order to distribute the STP and the LTP observations over 10 years, we somewhat arbitrarily chose the fractions illustrated in table 9.

The pointing positions were then determined independently for the Northern and Southern array for each year. For each array, the number of pointings to be scheduled was determined by dividing the total duration multiplied by the fraction of time used in a given year by the duration of one pointing (30 min). Given the resulting number of pointings, the longitude step size was then determined by dividing each longitude range by the number of pointings.

For each array, the first pointing was scheduled at a given start time. The start time was computed by adding 362.25 days for each year to the start time of the first year, which was taken to be January 1 2021, midnight. An additional offset of 29.5 days was added for the first year, and offsets of $(10 - \text{year}) \times 59 - 365.25$ days were added from the third year on, where we count the number of years starting from 0. In this way, the scheduling within a given year started with one- or two-month date shifts, assuring that similar pointing positions were not always observed during the same period in the year.

All pointings within a given year and for a given array were then scheduled independently. Before scheduling the first pointing, and after scheduling each pointing, the start time was incremented by adding 32 min, assuming 30 min exposure time and 2 min slew time. If, for the resulting time, the Sun had a zenith angle of less than 105° (i.e. the Sun started to rise), several days were added to advance in time. In this way, subsequent observation nights were separated by a few days, so that the GPS observations covered the entire year. The number of days that were added depended on the year of scheduling and the array site, and are specified in table 10.

After this addition, and as long as the Sun had a zenith angle of less than 105° and the Moon had a zenith angle of less than 90°, time was advanced in steps of 4 minutes. In this way, pointings were only scheduled when the Sun was 15° below horizon, and when the Moon
was close to or below the horizon. The suitability of observations under Moonlight for the Galactic plane survey was not investigated. Other factors that limit the time available and affect the observability such as weather and hardware failures were not taken into account. They will need to be addressed in future work.

Once the time of the next pointing was determined, the pointing with the smallest zenith angle was determined in the list of pointings to be scheduled for that year and site. Toggling between pointings with negative and positive Galactic latitudes was enforced. Furthermore it was enforced that the next pointing differed by at least 2° and at maximum 10° in Galactic longitude, assuring some displacement while limiting the slew time between pointings. Pointings that had a zenith angle larger than 5° from the best achievable zenith angle within a given year were excluded. In this way, pointings were always scheduled near their smallest possible zenith angle. In the case that none of the unscheduled pointings satisfied all criteria, the constraints in Galactic longitude were removed and the pointing satisfying the remaining criteria was selected. If no pointing was found after relaxing the constraints, time was advanced according to the procedure described above, and a pointing was searched for the new time. This procedure led to the scheduling of all pointings in a given year.

This algorithm made it possible to reach an exposure as a function of time difference between pointings within longitudes ±5° from every position in the Galactic plane quite uniform and without any gaps in the interval between 5 days and about 1 year. This covers the range of observed periods of known gamma-ray binaries and a large fraction of the binary population that was simulated for the GPS.

## C Details on catalogue A production

The catalogue production strategy from [89] features an iterative process with two steps: detecting significant excesses (objects) and fitting a model for these excesses to the data.

### C.1 Object detection

The significance map used for the object detection was computed using a Poisson likelihood comparing the observed counts to predicted counts. The observed and predicted counts were binned into maps with pixel widths of 0.05° and smoothed by a Gaussian with σ = 0.05°. The predicted counts map initially includes only the CR and interstellar background models, while subsequent iterations also include objects that have been detected and fitted (see next section) in previous steps.

Relevant structures in the significance map were identified using a technique based on the SExtractor algorithm [144]. First, the map is filtered to keep only pixels above a threshold of 3σ, then groups of more than 3 pixels above threshold are isolated. Each group is then decomposed into multiple substructures through a deblending procedure that allows to better disentangle overlapping sources. The positions of the deblended objects are used as initial guess for the models fitted in the next step.

### C.2 Model fitting

The Galactic plane is divided into 180 regions of 3° by 12° in longitude and latitude (with 1° of overlap in longitude). Each analysis region is binned using a spatial pixel width of 0.02°
and 20 logarithmically spaced energy bins from 0.07 TeV to 200 TeV. Each candidate object is fitted independently via a binned maximum-likelihood analysis in three dimensions. The spatial and spectral parameters of the candidate object, along with those of the background models, are fitted simultaneously. Spatial and spectral parameters of nearby sources that contribute to the analysis region are fixed at their previously fitted values. Only sources that are at least marginally significant (TS ≥ 10) are kept in subsequent iterations. Initially, candidate objects from the detection step are represented as point sources with power-law spectra. The spatial models of the sources are sequentially refined, first by comparing a point-source model with a disk model. If the disk model is preferred (TS<sub>disk</sub> − TS<sub>point</sub> ≥ 10) then a Gaussian model is also tested, and chosen if its TS is larger than that of the disk model. Note that the source angular extension is restricted to a maximum value of 3° to minimise the odds of incorrectly detecting background emission as a source.

Once the best-fit spatial model has been determined for all sources, each source is then
tested for spectral curvature. The default power-law model (PL) is compared with both
log-parabola and exponentially cut-off spectral models. The curved spectral model with the
largest TS value is chosen if TS<sub>curved</sub> − T<sub>PL</sub> ≥ 10. Finally, all sources are refitted to find
the best-fit global values. This is done to ensure that each source model accounts for the
best fit spatial and spectral parameters of all other sources.

D Details on catalogue B production

D.1 Datasets

For each observation we select events within a 3° and 5° offset from the pointing position,
below and above 1 TeV, respectively. We impose a smaller offset cut at low energy as the
instrument effective field-of-view is smaller. All observations are stacked together and binned
into three-dimensional maps. The energy axis contains 10 bins per decade in logarithmic
scale from 0.07 TeV to 200 TeV. The spatial binning is 0.06° and 0.03°, below and above
1 TeV, respectively. This energy-dependent spatial binning is motivated by the coarser
angular-resolution at low energy, and allows to improve fitting performance while saving
time and memory.

D.2 Object detection

The goal is to build, in a short amount of computational time, a list of objects without
prior case-specific morphological assumptions. This list provides source candidates, with
robust guesses on their position and morphological parameters, to be tested subsequently
by a conventional template-fitting analysis.

The first step is to compute the significance of the excesses above a given background model
(instrumental and interstellar). The significance maps are produced for different correlation
radii, Rcorr = 0.06°, 0.1°, 0.2° and for the two energy ranges considered E = 0.07 – 1 TeV
and E = 1 – 200 TeV. This separation is useful as different populations of sources may have
different optimal energy ranges for detection. The significance maps are filtered by hysteresis
thresholding [145] using the implementation provided by scikit-image [146]. We use two
thresholds: first, pixels above the higher threshold are selected, and then pixels between the
two thresholds are preserved only if they are continuously connected to a pixel above the high threshold. The low and high thresholds are set to $2\sigma$ and $4\sigma$, respectively.

Subsequent object detection combines three methods:

* Peak detection: we identify local maxima above $5\sigma$ in the hysteresis-filtered significance maps.

* Circle detection: the contours of each group of pixels isolated by the hysteresis filtering are fitted as a circle. If less than 80 % of the pixels are included in the circle, the object is discarded.

* Edge detection and Hough circle detection: for details see [90] and references therein.

Peak detection is performed on significance maps with $R_{corr} = 0.06^\circ, 0.1^\circ, 0.2^\circ$ for each of the two energy ranges. The circle and Hough circle detection are performed only for $R_{corr} = 0.1^\circ$ in the 3–100 TeV energy range in order to take advantage of the better angular-resolution. Removal of likely duplicate objects is performed when the results of the different detection methods are combined. Groups of objects with an inter-center angular separation less than $0.1^\circ$ and a difference in radius less than $0.25^\circ$ are replaced by a single object by averaging their position and radii.

### D.3 Classification, spatial model selection, and candidate filtering

Extended source morphologies are usually more complex than the parametric models we use to describe them, so a single extended source can be detected as a cluster of objects. To address this, we start by classifying the objects depending on their degree of overlap. For every object of radius $R$ in the list we calculate the inter-center angular separation to any other object of radius $R_{other}$. The object is then classified as:

* “isolated”, if $d > R + R_{other}$ for all other objects;

* “outer sub-structure”, if there exists another object such that $d \le R + R_{other}$ and $R < R_{other}$;

* “inner sub-structure”, if there exists another object such that $d \le R$ and $R < R_{other}$;

* “parent”, if there are any other objects such that $d \le R + R_{other}$ and for all $R_{other} \le R$;

In short: isolated objects are non-overlapping, while parents are large objects partially overlapping with smaller ones. These smaller objects are considered as sub-structures of their parent.

For each object, a baseline spatial model is determined using the Pearson correlation coefficient (PCC) of a 5-point radial profile in flux. To this purpose, we integrated the flux map (with $R_{corr} = 0.1^\circ$) in 5 rings of equal area between the centre of the object and its radius. Then we computed the PCC between the radius and flux values. The default spatial model is a generalised Gaussian (see next section). Alternatively, a shell is considered first if the candidate object has $PCC > 1/3$, or if it has $|PCC| < 1/3$ and overlaps with objects classified as inner sub-structures.
For each object we compute a test statistic (TS) defined as the squared significance of the residual excess integrated within a correlation radius equivalent to the object radius. The candidate objects are filtered by requiring TS>10. We then perform outlier detection using the isolation forest algorithm [147] implemented in scikit-learn [148]. The object properties considered when applying the outlier detection are: the mean distance of the 5 nearest neighbours, the distance of a sub-structure to its parent-object, and the PCC. The distribution of objects associated with known sources from existing catalogues<fnref n="14" /> in this parameter space informs us on the expected density of objects with a given morphology. This information is used to train the outlier detection classifier and to set its selection score. Sub-structures below a threshold in selection score are discarded in order to reduce spurious detections in complex sources (those that still remain after the removal of likely duplicate objects). The remaining objects are then ranked according to their selection score, which is used to determine their fitting order in the following (see next section).

Furthermore, after the TS filtering and the outlier detection the selected candidates
are divided into two lists:

* primary: objects associated with known sources, isolated objects, parent objects, and sub-structures more significant than their parent with a difference in TS larger than 25;

* secondary: unassociated sub-structures less significant than their parent.

D.4 Model fitting

For model fitting, the Galactic plane is divided into 10° wide regions separated by 5° (half-overlapping). A 3° border is added to each sub-region to account for the spill-over into the analysis region of photons from the sources outside the region due to the instrument PSF. Regions containing less than five sources are merged with their neighbour in order to limit the number of regions fitted. The 56 sub-regions obtained are then fitted independently. The objects with centres outside of the fit region, but within the 3° border, are merged into a unique background component. So, for each energy range, we have three background components: CR background, IEM, and sources centred outside the fitting region. Note that for the catalogue production, the true model for CR and IEM backgrounds are used and only a normalisation and spectral index correction are fitted.

By default, candidate sources are fitted with a log-parabola<sup>15</sup> as spectral model and a generalised-Gaussian<sup>16</sup> as spatial model. The generalised-Gaussian model has a shape parameter, η, fitted for values in [0.1, 1]. This model is equivalent to a disk when η tends to zero, a Gaussian for η = 0.5 and a Laplace profile for η = 1. The minimum size fitted is 0.06° (slightly larger than the mean PSF radius at 1 TeV). Alternatively, a shell is fitted as the spatial model based on the morphological estimate from the initial detection step (see previous section).

For each region, the candidates in the primary list are fitted first while those in the
secondary list are added only if there is still a significant residual excess after the fitting of

<footnote n="14"><sup>14</sup>The association procedure is described in 4.2. The only difference is that for candidate objects the surface used in the association criterion is computed as a disk using the radius estimated from the detection step.</footnote>

<footnote n="15"><sup>15</sup>https://docs.gammapy.org/1.0.1/user-guide/model-gallery/spectral/plot_logparabola.html.</footnote>

<footnote n="16"><sup>16</sup>https://docs.gammapy.org/1.0.1/user-guide/model-gallery/spatial/plot_gen_gauss.html.</footnote>
the primary candidates ($TS_{resi} > 25$). Once the initial candidate lists are exhausted, more sources are added iteratively (up to 5 per region) at the position of the largest peak above $5\sigma$ in the residual significance map. Finally, for each object, we compute the test statistic for the null hypothesis (no source) and keep only those with $TS_{null} > 25$ in at least one of the energy ranges. Once the fitting of all regions is complete, we assemble the final global model. Models are taken from one region only if they are located within $2.5^\circ$ from the centre of the region (regions are overlapping).

## D.5 Models refinement

At this point, for each source, we compare generalised-Gaussian, shell, and point-like morphological models. The optimal parameters of a model are given by the likelihood maximisation, and the selection among alternative models is performed by minimisation of the Akaike information criterion [AIC, 149] in order to take into account the difference in the number of parameters between the models (which are not necessarily nested).

In order to simplify the global model we also search for groups of sources that could be replaced by an alternative model. In particular, we scan the source distribution to identify linear and circular patterns that could be regrouped into a single elliptical Gaussian or shell, respectively. To do so, we first extract several groups of pixels by applying two thresholds to the flux map of the fitted model: the first (low) threshold (95% percentile of the flux distribution) isolates the different groups and the second (higher) threshold (10% of the group maximum flux) separates the bright peaks from the more diffuse regions. Then, within each of these regions we search for clusters of sources whose positions best match a linear or circular pattern using the RANSAC algorithm [150], implemented in scikit-image [146].

The clusters of sources forming a linear pattern were replaced by an elliptical generalised Gaussian if the difference in AIC between multiple models and one elliptical model was greater than zero. This test resulted in 13 new ellipticals. Similarly, the clusters of sources forming a circular pattern were tested against a shell, but we also considered a shell plus a Gaussian (that could best model a composite system made of a PWN and SNRs from the same progenitor). This test added one composite object to the catalogue. Overall the clustering procedure removed about 30 objects. We also looked for sources surrounded by negative or positive residuals above $1\sigma$ with an ellipticity larger than 0.5, for which we tested an elliptical generalised Gaussian. This test resulted in one additional elliptical source.

Finally, we reduced the catalogue to a list of complexes with similar spectral properties by merging the sub-structures with their parent object if they did not have a different association to a known source or a different spectral shape (if the differences in index and curvature parameters to their mean values within the group were lower than their standard deviations within the group).

## E Maps of unresolved point sources

Figure 18 provides an example of a two-dimensional map of unresolved sources in the energy band $> 1$ TeV reconstructed from the simulated CTA data as described in section 7.1.
<figure><img src="image_21.png" /><figcaption>**Figure 18.** Map of unresolved sources in the energy band > 1 TeV reconstructed from the simulated CTA data as described in section 7.1. We show a zoom of the map in the longitude range $|l| < 90^\circ$ where the emission is more prominent.</figcaption></figure>

## References

[1] H.E.S.S. collaboration, *The population of TeV pulsar wind nebulae in the H.E.S.S. galactic plane survey*, Astron. Astrophys. **612** (2018) A2 [arXiv:1702.08280] [INSPIRE].

[2] H.E.S.S. collaboration, *Population study of galactic supernova remnants at very high γ-ray energies with H.E.S.S.*, Astron. Astrophys. **612** (2018) A3 [arXiv:1802.05172] [INSPIRE].

[3] M. Chernyakova et al., *Overview of non-transient γ-ray binaries and prospects for the Cherenkov Telescope Array*, Astron. Astrophys. **631** (2019) A177 [arXiv:1909.11018] [INSPIRE].

[4] H.E.S.S. collaboration, *The H.E.S.S. galactic plane survey*, Astron. Astrophys. **612** (2018) A1 [arXiv:1804.02432] [INSPIRE].

[5] S. Funk, *Ground- and space-based gamma-ray astronomy*, Ann. Rev. Nucl. Part. Sci. **65** (2015) 245 [arXiv:1508.05190] [INSPIRE].

[6] L. Tibaldo, D. Gaggero and P. Martin, *Gamma rays as probes of cosmic-ray propagation and interactions in galaxies*, Universe **7** (2021) 141 [arXiv:2103.16423] [INSPIRE].

[7] T.C. Weekes et al., *Observation of TeV gamma rays from the Crab nebula using the atmospheric Cerenkov imaging technique*, Astrophys. J. **342** (1989) 379 [INSPIRE].

[8] HEGRA collaboration, *A search for TeV gamma-ray emission from SNRs, pulsars and unidentified GeV sources in the galactic plane in the longitude range between -2 deg and 85 deg*, Astron. Astrophys. **395** (2002) 803 [astro-ph/0209360] [INSPIRE].

[9] VERITAS collaboration, *A very high energy γ-ray survey towards the Cygnus region of the galaxy*, Astrophys. J. **861** (2018) 134 [arXiv:1805.05989] [INSPIRE].

[10] A.A. Abdo et al., *Milagro observations of TeV emission from galactic sources in the Fermi bright source list*, Astrophys. J. Lett. **700** (2009) L127 [Erratum ibid. **703** (2009) L185] [arXiv:0904.1018] [INSPIRE].

[11] ARGO-YBJ collaboration, *TeV gamma-ray survey of the northern sky using the ARGO-YBJ detector*, Astrophys. J. **779** (2013) 27 [arXiv:1311.3376] [INSPIRE].

[12] HAWC collaboration, *3HWC: the third HAWC catalog of very-high-energy gamma-ray sources*, Astrophys. J. **905** (2020) 76 [arXiv:2007.08582] [INSPIRE].

[13] LHAASO collaboration, *The first LHAASO catalog of gamma-ray sources*, Astrophys. J. Suppl. **271** (2024) 25 [arXiv:2305.17030] [INSPIRE].

[14] LHAASO collaboration, *Ultrahigh-energy photons up to 1.4 petaelectronvolts from 12 γ-ray galactic sources*, Nature **594** (2021) 33 [INSPIRE].

[15] FERMI-LAT collaboration, *Incremental Fermi Large Area Telescope fourth source catalog*, Astrophys. J. Supp. **260** (2022) 53 [arXiv:2201.11184] [INSPIRE].
[16] A.A. Abdo et al., *A measurement of the spatial distribution of diffuse TeV gamma ray emission from the galactic plane with Milagro*, *Astrophys. J.* **688** (2008) 1078 [arXiv:0805.0417] [INSPIRE].

[17] H.E.S.S. collaboration, *Diffuse galactic gamma-ray emission with H.E.S.S.*, *Phys. Rev. D* **90** (2014) 122007 [arXiv:1411.7568] [INSPIRE].

[18] HAWC collaboration, *Gamma ray diffuse emission from the galactic plane with HAWC data*, *PoS ICRC2019* (2020) 750 [INSPIRE].

[19] TIBET ASGAMMA collaboration, *First detection of sub-PeV diffuse gamma rays from the galactic disk: evidence for ubiquitous galactic cosmic rays beyond PeV energies*, *Phys. Rev. Lett.* **126** (2021) 141101 [arXiv:2104.05181] [INSPIRE].

[20] LHAASO collaboration, *Measurement of the diffuse gamma-ray emission from galactic plane with LHAASO-KM2A*, *PoS ICRC2021* (2021) 859 [INSPIRE].

[21] CTA CONSORTIUM collaboration, *Design concepts for the Cherenkov Telescope Array CTA: an advanced facility for ground-based high-energy gamma-ray astronomy*, *Exper. Astron.* **32** (2011) 193 [arXiv:1008.3703] [INSPIRE].

[22] W. Hofmann and R. Zanin, *The Cherenkov Telescope Array*, arXiv:2305.12888 [INSPIRE].

[23] CTA CONSORTIUM collaboration, *Surveys with the Cherenkov Telescope Array*, *Astropart. Phys.* **43** (2013) 317 [arXiv:1208.5686] [INSPIRE].

[24] CTA CONSORTIUM collaboration, *Science with the Cherenkov Telescope Array*, World Scientific, Singapore (2018) [DOI:10.1142/10986] [INSPIRE].

[25] CHERENKOV TELESCOPE ARRAY CONSORTIUM collaboration, *Galactic transient sources with the Cherenkov Telescope Array*, arXiv:2405.04469 [INSPIRE].

[26] S.P. Wakely and D. Horan, *TeVCat: an online catalog for very high energy gamma-ray astronomy*, in 30th *International Cosmic Ray Conference 3* (2007), p. 1341 [INSPIRE].

[27] FERMI-LAT collaboration, *The second Fermi Large Area Telescope catalog of gamma-ray pulsars*, *Astrophys. J. Suppl.* **208** (2013) 17 [arXiv:1305.4385] [INSPIRE].

[28] FERMI-LAT collaboration, *3FHL: the third catalog of hard Fermi-LAT sources*, *Astrophys. J. Suppl.* **232** (2017) 18 [arXiv:1702.00664] [INSPIRE].

[29] A.U. Abeysekara et al., *The 2HWC HAWC observatory gamma ray catalog*, *Astrophys. J.* **843** (2017) 40 [arXiv:1702.02992] [INSPIRE].

[30] S. Casanova, *On the search for the galactic PeVatrons by means of gamma-ray astronomy*, *Universe* **8** (2022) 505 [INSPIRE].

[31] F. Aharonian, R. Yang and E. de Oña Wilhelmi, *Massive stars as major factories of galactic cosmic rays*, *Nature Astron.* **3** (2019) 561 [arXiv:1804.02331] [INSPIRE].

[32] H.E.S.S. collaboration, *H.E.S.S. J1641-463*, a very hard spectrum TeV gamma-ray source in the galactic plane, *PoS ICRC2015* (2016) 834 [arXiv:1509.08310] [INSPIRE].

[33] P. Cristofari et al., *Acceleration of cosmic rays and gamma-ray emission from supernova remnants in the galaxy*, *Mon. Not. Roy. Astron. Soc.* **434** (2013) 2748 [arXiv:1302.2150] [INSPIRE].

[34] P. Cristofari et al., *Supernova remnants in the very-high-energy gamma-ray domain: the role of the Cherenkov telescope array*, *Mon. Not. Roy. Astron. Soc.* **471** (2017) 201 [arXiv:1709.01102] [INSPIRE].
[35] W. Li et al., *Nearby supernova rates from the Lick observatory supernova search. III. The rate-size relation, and the rates as a function of galaxy Hubble type and colour*, *Mon. Not. Roy. Astron. Soc.* **412** (2011) 1473 [arXiv:1006.4613] [INSPIRE].

[36] V. Ptuskin, V. Zirakashvili and E.-S. Seo, *Spectrum of galactic cosmic rays accelerated in supernova remnants*, *Astrophys. J.* **718** (2010) 31 [arXiv:1006.0034] [INSPIRE].

[37] S.J. Smartt, *Progenitors of core-collapse supernovae*, *Ann. Rev. Astron. Astrophys.* **47** (2009) 63 [arXiv:0908.0700] [INSPIRE].

[38] C.-A. Faucher-Giguère and V.M. Kaspi, *Birth and evolution of isolated radio pulsars*, *Astrophys. J.* **643** (2006) 332 [astro-ph/0512585] [INSPIRE].

[39] D.R. Lorimer, *The galactic population and birth rate of radio pulsars*, *IAU Symp.* **218** (2004) 105 [astro-ph/0308501] [INSPIRE].

[40] F. Mannucci, M. Della Valle and N. Panagia, *Two populations of progenitors for type IA supernovae*, *Mon. Not. Roy. Astron. Soc.* **370** (2006) 773 [astro-ph/0510315] [INSPIRE].

[41] E. Scannapieco and L. Bildsten, *The type IA supernova rate*, *Astrophys. J. Lett.* **629** (2005) L85 [astro-ph/0507456] [INSPIRE].

[42] N. Prantzos et al., *The 511 keV emission from positron annihilation in the galaxy*, *Rev. Mod. Phys.* **83** (2011) 1001 [arXiv:1009.4620] [INSPIRE].

[43] H. Nakanishi and Y. Sofue, *Three-dimensional distribution of the ISM in the Milky Way galaxy: I. The HI disk*, *Publ. Astron. Soc. Jap.* **55** (2003) 191 [astro-ph/0304338] [INSPIRE].

[44] H. Nakanishi and Y. Sofue, *Three-dimensional distribution of the ISM in the Milky Way galaxy: II. The molecular gas disk*, *Publ. Astron. Soc. Jap.* **58** (2006) 847 [astro-ph/0610769] [INSPIRE].

[45] V.S. Ptuskin and V.N. Zirakashvili, *On the spectrum of high-energy cosmic rays produced by supernova remnants in the presence of strong cosmic-ray streaming instability and wave dissipation*, *Astron. Astrophys.* **429** (2005) 755 [astro-ph/0408025] [INSPIRE].

[46] A.W. Strong, I.V. Moskalenko and V.S. Ptuskin, *Cosmic-ray propagation and interactions in the galaxy*, *Ann. Rev. Nucl. Part. Sci.* **57** (2007) 285 [astro-ph/0701517] [INSPIRE].

[47] V.N. Zirakashvili and V.S. Ptuskin, *Diffusive shock acceleration with magnetic amplification by non-resonant streaming instability in SNRs*, *Astrophys. J.* **678** (2008) 939 [arXiv:0801.4488] [INSPIRE].

[48] Y. Ohira, K. Murase and R. Yamazaki, *Escape-limited model of cosmic-ray acceleration revisited*, *Astron. Astrophys.* **513** (2010) A17 [arXiv:0910.3449] [INSPIRE].

[49] D. Caprioli, E. Amato and P. Blasi, *The contribution of supernova remnants to the galactic cosmic ray spectrum*, *Astropart. Phys.* **33** (2010) 160 [arXiv:0912.2964] [INSPIRE].

[50] S. Gabici, *Cosmic ray escape from supernova remnants*, *Mem. Soc. Ast. It.* **82** (2011) 760 [arXiv:1108.4844] [INSPIRE].

[51] D.C. Ellison, D.J. Patnaude, P. Slane and J. Raymond, *Efficient cosmic ray acceleration, hydrodynamics, and self-consistent thermal X-ray emission applied to SNR RX J1713.7-3946*, *Astrophys. J.* **712** (2010) 287 [arXiv:1001.1932] [INSPIRE].

[52] G. Morlino and D. Caprioli, *Strong evidences of hadron acceleration in Tycho's supernova remnant*, *Astron. Astrophys.* **538** (2012) A81 [arXiv:1105.6342] [INSPIRE].
[53] S.R. Kelner, F.A. Aharonian and V.V. Bugayov, *Energy spectra of gamma-rays, electrons and neutrinos produced at proton-proton interactions in the very high energy regime*, Phys. Rev. D **74** (2006) 034018 [Erratum ibid. **79** (2009) 039901] [astro-ph/0606058] [INSPIRE].

[54] M. Mori, *Nuclear enhancement factor in calculation of galactic diffuse gamma-rays: a new estimate with DPMJET-3*, Astropart. Phys. **31** (2009) 341 [arXiv:0903.3260] [INSPIRE].

[55] G.R. Blumenthal and R.J. Gould, *Bremsstrahlung, synchrotron radiation, and Compton scattering of high-energy electrons traversing dilute gases*, Rev. Mod. Phys. **42** (1970) 237 [INSPIRE].

[56] S. Funk, *High-energy gamma rays from supernova remnants*, in *Handbook of supernovae*, A.W. Alsabti and P. Murdin eds., Springer, Cham, Switzerland (2016) p. 1 [DOI:10.1007/978-3-319-20794-0_12-1].

[57] T.S. Rice et al., *A uniform catalog of molecular clouds in the Milky Way*, Astrophys. J. **822** (2016) 52 [arXiv:1602.02791].

[58] C. Silvia, *Interacting supernova remnants: a population model for the Cherenkov Telescope Array*, arXiv:2306.14359 [INSPIRE].

[59] K.P. Watters and R.W. Romani, *The galactic population of young gamma-ray pulsars*, Astrophys. J. **727** (2011) 123 [arXiv:1009.5305] [INSPIRE].

[60] R. Bandiera et al., *Reverberation of pulsar wind nebulae — II. Anatomy of the ‘thin-shell’ evolution*, Mon. Not. Roy. Astron. Soc. **520** (2023) 2451 [arXiv:2301.04056] [INSPIRE].

[61] M. Fiori et al., *Modelling the γ-ray pulsar wind nebulae population in our galaxy*, Mon. Not. Roy. Astron. Soc. **511** (2022) 1439 [arXiv:2201.02221] [INSPIRE].

[62] R. Bandiera et al., *Reverberation of pulsar wind nebulae (I): impact of the medium properties and other parameters upon the extent of the compression*, Mon. Not. Roy. Astron. Soc. **499** (2020) 2051 [arXiv:2009.10992] [INSPIRE].

[63] J.D. Gelfand, P.O. Slane and W. Zhang, *A dynamical model for the evolution of a pulsar wind nebula inside a non-radiative supernova remnant*, Astrophys. J. **703** (2009) 2051 [arXiv:0904.4053] [INSPIRE].

[64] T.A. Porter, G. Jóhannesson and I.V. Moskalenko, *High-energy gamma rays from the Milky Way: three-dimensional spatial models for the cosmic-ray and radiation field densities in the interstellar medium*, Astrophys. J. **846** (2017) 67 [arXiv:1708.00816] [INSPIRE].

[65] B. Olmi and N. Bucciantini, *On the origin of jet-like features in bow shock pulsar wind nebulae*, Mon. Not. Roy. Astron. Soc. **490** (2019) 3608 [arXiv:1910.01926] [INSPIRE].

[66] T. Sudoh, T. Linden and J.F. Beacom, *TeV halos are everywhere: prospects for new discoveries*, Phys. Rev. D **100** (2019) 043016 [arXiv:1902.08203] [INSPIRE].

[67] G. Giacinti et al., *Halo fraction in TeV-bright pulsar wind nebulae*, Astron. Astrophys. **636** (2020) A113 [arXiv:1907.12121] [INSPIRE].

[68] P. Martin, L. Tibaldo, A. Marcowith and S. Abdollahi, *Population synthesis of pulsar wind nebulae and pulsar halos in the Milky Way — predicted contributions to the very-high-energy sky*, Astron. Astrophys. **666** (2022) A7 [arXiv:2207.11178] [INSPIRE].

[69] HAWC collaboration, *Extended gamma-ray sources around pulsars constrain the origin of the positron flux at earth*, Science **358** (2017) 911 [arXiv:1711.06223] [INSPIRE].

[70] G. Dubus, N. Guillard, P.-O. Petrucci and P. Martin, *Sizing up the population of gamma-ray binaries*, Astron. Astrophys. **608** (2017) A59 [arXiv:1707.05744] [INSPIRE].
[71] I. Yusifov and I. Kucuk, *Revisiting the radial distribution of pulsars in the galaxy, Astron. Astrophys.* **422** (2004) 545 [astro-ph/0405559] [INSPIRE].

[72] G. Dubus, *Gamma-ray absorption in massive X-ray binaries, Astron. Astrophys.* **451** (2006) 9 [astro-ph/0509633] [INSPIRE].

[73] V.A. Ambartsumian, *On the statistics of double stars, Astron. Zh.* **14** (1937) 207.

[74] P.D.T. Luque et al., *Galactic diffuse gamma rays meet the PeV frontier, Astron. Astrophys.* **672** (2023) A58 [arXiv:2203.15759] [INSPIRE].

[75] P.D. Marinos, G.P. Rowell, T.A. Porter and G. Jóhannesson, *The steady-state multi-TeV diffuse γ-ray emission predicted with galprop and prospects for the Cherenkov Telescope Array, Mon. Not. Roy. Astron. Soc.* **518** (2022) 5036 [arXiv:2211.01619] [INSPIRE].

[76] C. Evoli et al., *Cosmic-ray propagation with DRAGON2: I. Numerical solver and astrophysical ingredients, JCAP* **02** (2017) 015 [arXiv:1607.07886] [INSPIRE].

[77] C. Evoli et al., *Cosmic-ray propagation with DRAGON2: II. Nuclear interactions with the interstellar gas, JCAP* **07** (2018) 006 [arXiv:1711.09616] [INSPIRE].

[78] A. Dundovic, C. Evoli, D. Gaggero and D. Grasso, *Simulating the galactic multi-messenger emissions with HERMES, Astron. Astrophys.* **653** (2021) A18 [arXiv:2105.13165] [INSPIRE].

[79] Q. Remy, I.A. Grenier, J.M. Casandjian and FERMI-LAT collaboration, *A 3D view of our galaxy: gas, dust and cosmic rays*, in 8th *International Fermi Symposium*, Baltimore, MD, U.S.A., 14–19 October 2018.

[80] N. Arimoto, Y. Sofue and T. Tsujimoto, *CO-to-H2 conversion factor in galaxies*, *Publ. Astron. Soc. Japan* **48** (1996) 275.

[81] I.A. Grenier, J.H. Black and A.W. Strong, *The nine lives of cosmic rays in galaxies*, *Ann. Rev. Astron. Astrophys.* **53** (2015) 199 [INSPIRE].

[82] FERMI-LAT collaboration, *Development of the model of galactic interstellar emission for standard point-source analysis of Fermi Large Area Telescope data, Astrophys. J. Suppl.* **223** (2016) 26 [arXiv:1602.07246] [INSPIRE].

[83] A. Acharyya et al., *Monte Carlo studies for the optimisation of the Cherenkov Telescope Array layout, Astropart. Phys.* **111** (2019) 35 [arXiv:1904.01426] [INSPIRE].

[84] CTA OBSERVATORY et al. collaborations, *CTA — the world’s largest ground-based gamma-ray observatory, PoS ICRC2021* (2022) 005 [INSPIRE].

[85] CTAO and CTAC collaboration, *CTAO instrument response functions — version prod3b-v2*, Zenodo (2016).

[86] G. Maier, O. Gueta and R. Zanin, *CTAO instrument response functions: comparison of prod5 and prod3b releases*, Zenodo (2023).

[87] J. Knödlseder et al., *GammaLib and ctools: a software framework for the analysis of astronomical gamma-ray data, Astron. Astrophys.* **593** (2016) A1 [arXiv:1606.00393] [INSPIRE].

[88] CTA CONSORTIUM collaboration, *Gammapy — a prototype for the CTA science tools*, PoS ICRC2017 (2018) 766 [arXiv:1709.01751] [INSPIRE].

[89] J. Cardenzana, *A catalog pipeline for sources in the CTA galactic plane survey*, in the proceedings of the SF2A-2018, (2018).
[90] Q. Remy, Y.A. Gallant and M. Renaud, Prospects on high-energy source searches based on pattern recognition: object detection in the H.E.S.S. galactic plane survey and catalogue cross-matches, Astropart. Phys. 122 (2020) 102462 [arXiv:2005.05176] [INSPIRE].

[91] HAWC collaboration, 3HWC: the third HAWC catalog of very-high-energy gamma-ray sources, Astrophys. J. 905 (2020) 76 [arXiv:2007.08582] [INSPIRE].

[92] O.C. de Jager, B.C. Raubenheimer and J.W.H. Swanepoel, A powerful test for weak periodic signals with unknown light curve shape in sparse data, Astron. Astrophys. 221 (1989) 180.

[93] M. Kerr, Improving sensitivity to weak pulsations with photon probability weighting, Astrophys. J. 732 (2011) 38 [arXiv:1103.2128] [INSPIRE].

[94] F.T.H. Djannati-Ataï, H.E.S.S. discovery of 20 TeV gamma-ray pulsations from the Vela pulsar, in 31st Texas symposium on relativistic astrophysics, (2022).

[95] CHERENKOV TELESCOPE ARRAY CONSORTIUM collaboration, Sensitivity of the Cherenkov Telescope Array to spectral signatures of hadronic PeVatrons with application to galactic supernova remnants, Astropart. Phys. 150 (2023) 102850 [arXiv:2303.15007] [INSPIRE].

[96] L. Wasserman, All of statistics: a concise course in statistical inference, Springer, New York, NY, U.S.A. (2004) [doi:10.1007/978-0-387-21736-9].

[97] E. Wilhelmi de Oña, R. López-Coto, E. Amato and F. Aharonian, On the potential of bright, young pulsars to power ultrahigh gamma-ray sources, Astrophys. J. Lett. 930 (2022) L2 [arXiv:2204.09440] [INSPIRE].

[98] C. Steppa and K. Egberts, Modelling the galactic very-high-energy γ-ray source population, Astron. Astrophys. 643 (2020) A137 [arXiv:2010.03305] [INSPIRE].

[99] ICECUBE collaboration, Observation of high-energy neutrinos from the galactic plane, Science 380 (2023) abc9818 [arXiv:2307.04427] [INSPIRE].

[100] G. Schwefer, P. Mertsch and C. Wiebusch, Diffuse emission of galactic high-energy neutrinos from a global fit of cosmic rays, Astrophys. J. 949 (2023) 16 [arXiv:2211.15607] [INSPIRE].

[101] ANTARES collaboration, Hint for a TeV neutrino emission from the galactic ridge with ANTARES, Phys. Lett. B 841 (2023) 137951 [arXiv:2212.11876] [INSPIRE].

[102] D. Gaggero et al., The gamma-ray and neutrino sky: a consistent picture of Fermi-LAT, Milagro, and IceCube results, Astrophys. J. Lett. 815 (2015) L25 [arXiv:1504.00227] [INSPIRE].

[103] FERMI-LAT collaboration, Fermi Large Area Telescope first source catalog, Astrophys. J. Suppl. 188 (2010) 405 [arXiv:1002.2280] [INSPIRE].

[104] M. Salvato et al., Finding counterparts for all-sky X-ray surveys with Nway: a Bayesian algorithm for cross-matching multiple catalogues, Mon. Not. Roy. Astron. Soc. 473 (2018) 4937 [arXiv:1705.10711] [INSPIRE].

[105] B. Barkus et al., The application of ridgelines in extended radio source cross-identification, Mon. Not. Roy. Astron. Soc. 509 (2021) 1 [arXiv:2110.05254].

[106] CTA CONSORTIUM collaboration, Prospects for Cherenkov Telescope Array observations of the young supernova remnant RX J1713.7-3946, Astrophys. J. 840 (2017) 74 [arXiv:1704.04136] [INSPIRE].

[107] MAGIC collaboration, Teraelectronvolt pulsed emission from the Crab pulsar detected by MAGIC, Astron. Astrophys. 585 (2016) A133 [arXiv:1510.07048] [INSPIRE].
[108] A. Burtovoi, T.Y. Saito, L. Zampieri and T. Hassan, Prospects for the detection of high-energy (E > 25 GeV) Fermi pulsars with the Cherenkov Telescope Array, Mon. Not. Roy. Astron. Soc. 471 (2017) 431 [arXiv:1706.07228] [INSPIRE].

[109] H.E.S.S. collaboration, First ground-based measurement of sub-20 GeV to 100 GeV γ-rays from the Vela pulsar with H.E.S.S. II, Astron. Astrophys. 620 (2018) A66 [arXiv:1807.01302] [INSPIRE].

[110] M.L. Ahnen et al., Search for VHE gamma-ray emission from Geminga pulsar and nebula with the MAGIC telescopes, Astron. Astrophys. 591 (2016) A138 [arXiv:1603.00730] [INSPIRE].

[111] H.E.S.S. collaboration, Search for pulsed VHE gamma-ray emission from young pulsars with H.E.S.S., Astron. Astrophys. 466 (2007) 543 [astro-ph/0702336] [INSPIRE].

[112] A. Archer et al., A search for pulsed very high-energy gamma-rays from 13 young pulsars in archival VERITAS data, Astrophys. J. 876 (2019) 95 [arXiv:1904.09329] [INSPIRE].

[113] FERMI-LAT collaboration, The third Fermi Large Area Telescope catalog of gamma-ray pulsars, Astrophys. J. 958 (2023) 191 [arXiv:2307.11132] [INSPIRE].

[114] M. Kerr et al., Timing gamma-ray pulsars with the Fermi Large Area Telescope: timing noise and astrometry, Astrophys. J. 814 (2015) 128 [arXiv:1510.05099] [INSPIRE].

[115] G. Martí-Devesa and O. Reimer, X-ray and γ-ray orbital variability from the γ-ray binary HESS J1832-093, Astron. Astrophys. 637 (2020) A23 [arXiv:2001.02701] [INSPIRE].

[116] H.E.S.S. collaboration, H.E.S.S. and Fermi-LAT observations of PSR B1259-63/LS 2883 during its 2014 and 2017 periastron passages, Astron. Astrophys. 633 (2020) A102 [arXiv:1912.05868] [INSPIRE].

[117] VERITAS and MAGIC collaborations, Periastron observations of TeV gamma-ray emission from a binary system with a 50-year period, Astrophys. J. Lett. 867 (2018) L19 [arXiv:1810.05271] [INSPIRE].

[118] P.C. Gregory, Bayesian analysis of radio observations of the Be X-ray binary LS I +61°303, The Astrophysical Journal 575 (2002) 427.

[119] H.E.S.S. collaboration, 3.9 day orbital modulation in the TeV gamma-ray flux and spectrum from the X-ray binary LS 5039, Astron. Astrophys. 460 (2006) 743 [astro-ph/0607192] [INSPIRE].

[120] H.E.S.S. collaboration, Discovery of variable VHE γ-ray emission from the binary system 1FGL J1018.6-5856, Astron. Astrophys. 577 (2015) A131 [arXiv:1503.02711] [INSPIRE].

[121] VERITAS and H.E.S.S. collaborations, Long-term TeV and X-ray observations of the gamma-ray binary H.E.S.S. J0632+057, Astrophys. J. 780 (2014) 168 [arXiv:1311.6083] [INSPIRE].

[122] H.E.S.S. collaboration, Discovery of the VHE gamma-ray source H.E.S.S. J1832-093 in the vicinity of SNR G22.7-0.2, Mon. Not. Roy. Astron. Soc. 446 (2015) 1163 [arXiv:1411.0572] [INSPIRE].

[123] VERITAS collaboration, The TeV morphology of the interacting supernova remnant IC 443, PoS ICRC2015 (2016) 875 [arXiv:1512.01911] [INSPIRE].

[124] FERMI-LAT collaboration, Search for extended sources in the galactic plane using six years of Fermi-Large Area Telescope pass 8 data above 10 GeV, Astrophys. J. 843 (2017) 139 [arXiv:1702.00476] [INSPIRE].
[125] H.E.S.S. collaboration, Probing the extent of the non-thermal emission from the Vela X region at TeV energies with H.E.S.S., Astron. Astrophys. 548 (2012) A38 [arXiv:1210.1359] [INSPIRE].

[126] D.C.J. Bock, M.I. Large and E.M. Sadler, Sumss: a wide-field radio imaging survey of the southern sky. I. Science goals, survey design and instrumentation, Astron. J. 117 (1999) 1578 [astro-ph/9812083] [INSPIRE].

[127] T. Murphy et al., The Molonglo Galactic Plane Survey (MGPS-2): compact source catalogue, Mon. Not. Roy. Astron. Soc. 382 (2007) 382 [arXiv:0708.3092] [INSPIRE].

[128] W. Voges et al., The ROSAT all-sky survey bright source catalogue, Astron. Astrophys. 349 (1999) 389 [astro-ph/9909315] [INSPIRE].

[129] H.E.S.S. collaboration, Discovery of very high energy gamma-ray emission coincident with molecular clouds in the W28 (G6.4-0.1) field, Astron. Astrophys. 481 (2008) 401 [arXiv:0801.3555] [INSPIRE].

[130] J. Trümper, ROSAT — a new look at the X-ray sky, Science 260 (1993) 1769.

[131] H.E.S.S. collaboration, H.E.S.S. reveals a lack of TeV emission from the supernova remnant Puppis A, Astron. Astrophys. 575 (2015) A81 [arXiv:1412.6997] [INSPIRE].

[132] F. Acero et al., A joint spectro-imaging analysis of the XMM-Newton and H.E.S.S. observations of the supernova remnant RX J1713.7-3946, Astron. Astrophys. 505 (2009) 157 [arXiv:0906.1073] [INSPIRE].

[133] H.E.S.S. collaboration, Particle transport within the pulsar wind nebula H.E.S.S. J1825-137, Astron. Astrophys. 621 (2019) A116 [arXiv:1810.12676] [INSPIRE].

[134] HAWC collaboration, Very high energy particle acceleration powered by the jets of the microquasar SS 433, Nature 562 (2018) 82 [Erratum ibid. 564 (2018) E38] [arXiv:1810.01892] [INSPIRE].

[135] R. López-Coto et al., Effect of the diffusion parameters on the observed γ-ray spectrum of sources and their contribution to the local all-electron spectrum: the EDGE code, Astropart. Phys. 102 (2018) 1 [arXiv:1709.07653] [INSPIRE].

[136] M. Ackermann et al., A cocoon of freshly accelerated cosmic rays detected by Fermi in the Cygnus superbubble, Science 334 (2011) 1103 [INSPIRE].

[137] S. Ohm, J.A. Hinton and R. White, Gamma-ray emission from the Westerlund 1 region, Mon. Not. Roy. Astron. Soc. 434 (2013) 2289 [arXiv:1306.5642] [INSPIRE].

[138] H.E.S.S. collaboration, Discovery of very-high-energy gamma-rays from the galactic centre ridge, Nature 439 (2006) 695 [astro-ph/0603021] [INSPIRE].

[139] L. Herold and D. Malyshev, Hard and bright gamma-ray emission at the base of the Fermi bubbles, Astron. Astrophys. 625 (2019) A110 [arXiv:1904.01454] [INSPIRE].

[140] G. Ferrand and S. Safi-Harb, A census of high-energy observations of galactic supernova remnants, Adv. Space Res. 49 (2012) 1313 [arXiv:1202.0245] [INSPIRE].

[141] M. Cardillo, E. Amato and P. Blasi, On the cosmic ray spectrum from type II supernovae expanding in their red giant presupernova wind, Astropart. Phys. 69 (2015) 1 [arXiv:1503.03001] [INSPIRE].

[142] R.N. Manchester, G.B. Hobbs, A. Teoh and M. Hobbs, VizieR online data catalog: ATNF pulsar catalogue (Manchester+, 2005), VizieR online data catalog B/psr, (2016).
[143] N. Bucciantini et al., The effects of spin-down on the structure and evolution of pulsar wind nebulae, Astron. Astrophys. 422 (2004) 609 [astro-ph/0404328] [INSPIRE].

[144] E. Bertin and S. Arnouts, SExtractor: software for source extraction, Astron. Astrophys. Suppl. Ser. 117 (1996) 393 [INSPIRE].

[145] J. Canny, A computational approach to edge detection, IEEE Trans. Pattern Anal. Machine Intell. 8 (1986) 679.

[146] S. van der Walt et al., Scikit-image: image processing in python, PeerJ 2 (2014) e453.

[147] S. Hariri, M.C. Kind and R.J. Brunner, Extended isolation forest, IEEE Trans. Knowledge Data Eng. 33 (2021) 1479 [arXiv:1811.02141].

[148] F. Pedregosa et al., Scikit-learn: machine learning in python, J. Machine Learning Res. 12 (2011) 2825 [arXiv:1201.0490] [INSPIRE].

[149] H. Akaike, A new look at the statistical model identification, IEEE Trans. Automatic Control 19 (1974) 716 [INSPIRE].

[150] M.A. Fischler and R.C. Bolles, Random sample consensus: a paradigm for model fitting with applications to image analysis and automated cartography, Commun. ACM 24 (1981) 381.
The CTA Consortium

S. Abe¹, J. Abhir², A. Abhishek³, F. Acero⁴,⁵, A. Acharyya⁶, R. Adam⁷,⁸, A. Aguasca-Cabot⁹,
I. Agudo¹⁰, A. Aguirre-Santaella¹¹,¹², J. Alfaro¹³, N. Alvarez-Crespo¹⁴, R. Alves Batista¹²,
J.-P. Amans¹⁵, E. Amato¹⁶, G. Ambrosi¹⁷, F. Ambrosino¹⁸, E. O. Angüner¹⁹, C. Aramo²⁰,
C. Arcaro²¹, L. Arrabito²², K. Asano¹, Y. Ascasíbar¹², J. Aschersleben²³, L. Augusto Stuani²⁴,
M. Backes²⁵,²⁶, C. Balazs²⁷, M. Balbo²⁸, J. Ballet⁴, A. Baquero Larriva¹⁴,²⁹, V. Barbosa Martins³⁰,
U. Barres de Almeida³¹,³², J. A. Barrio¹⁴, I. Batković³³, R. Batzofin³⁴, J. Baxter¹,
J. Becerra González³⁵, G. Beck³⁶, L. Beiske³⁷, R. Belmont⁴, W. Benbow³⁸, E. Bernardini³³,
J. Bernete³⁹, K. Bernlöhr⁴⁰, A. Berti⁴¹, B. Bertucci¹⁷, V. Beshley⁴², P. Bhattacharjee⁴³,
S. Bhattacharyya⁴⁴, B. Bi⁴⁵, N. Biederbeck³⁷, A. Biland², E. Bissaldi⁴⁶,⁴⁷, J. Biteau⁴⁸,⁴⁹,
O. Blanch⁵⁰, J. Blazek⁵¹, F. Bocchino⁵², C. Boisson¹⁵, J. Bolmont⁵³, L. Bonneau Arbeletche²⁴,
G. Bonnoli⁵⁴,⁵⁵, A. Bonollo⁵⁶,⁵⁷, P. Bordas⁹, Z. Bosnjak⁵⁸, E. Bottacini³³, C. Braiding⁵⁹,
E. Bronzini⁶⁰, R. Brose⁶¹, A. M. Brown⁶², F. Brun⁶³, G. Brunelli⁶⁰, N. Bucciantini¹⁶,⁶⁴,⁶⁵,
A. Bulgarelli⁶⁰, I. Burelli⁶⁶, L. Burmistrov⁶⁷, M. Burton⁶⁸,⁶⁹, A. Burtovoi¹⁶, T. Bylund⁴,
P. G. Calisse⁷⁰, A. Campoy-Ordaz⁷¹, B. K. Cantlay⁷²,⁷³, A. Caproni⁷⁴, R. Capuzzo-Dolcetta¹⁸,⁷⁵,
P. Caraveo⁵⁷, S. Caroff⁴³, A. Carosi¹⁸, R. Carosi⁵⁵, E. Carquin⁷⁶, M.-S. Carrasco⁷⁷, E. Cascone⁷⁸,
F. Cassol⁷⁷, N. Castrejon⁷⁹, A. J. Castro-Tirado¹⁰, D. Cerasole⁸⁰, M. Cerruti⁸¹, P. M. Chadwick⁶²,
P. Chambery⁸², S. Chaty⁸¹, A. W. Chen³⁶, M. Chernyakova⁸³, A. Chiavassa⁸⁴,⁸⁵, L. Chytka⁵¹,
A. Cifuentes³⁹, C. H. Coimbra Araujo⁸⁶, V. Conforti⁶⁰, F. Conte⁴⁰, J. L. Contreras¹⁴, J. Cortina³⁹,
A. Costa⁸⁷, H. Costantini⁷⁷, G. Cotter⁸⁸, S. Crestani⁵⁷, P. Cristofari¹⁵, O. Cuevas⁸⁹,
Z. Curtis-Ginsberg⁹⁰, A. D'Ai⁹¹, G. D'Amico⁹², F. D'Ammando⁹³, M. Dadina⁶⁰, M. Dalchenko⁶⁷,
L. David³⁰, F. Dazzi⁹⁴, M. de Bony de Lavergne⁶³, V. De Caprio⁷⁸, F. De Frondat Laadim¹⁵,
E. M. de Gouveia Dal Pino³², B. De Lotto⁶⁶, M. De Lucia²⁰, D. de Martino⁷⁸, R. de Menezes⁸⁴,⁸⁵,
M. de Naurois⁸, E. de Ona Wilhelm¹⁰, V. de Souza²⁴, L. del Peral⁷⁹, A. G. Delgado Giler²⁴,²³,
C. Delgado³⁹, M. Dell'aiera⁴³, M. Della Valle⁷⁸,²⁰, D. della Volpe⁶⁷, D. Depaoli⁴⁰,
T. Di Girolamo⁹⁵,²⁰, A. Di Piano⁶⁰,⁹⁶, F. Di Pierro⁸⁴, R. Di Tria⁸⁰, L. Di Venere⁴⁷, C. Díaz³⁹,
S. Diebold⁴⁵, A. Dinesh¹⁴, A. Djannati-Atai⁸¹, J. Djuvsland⁹², A. Domínguez¹⁴, R. M. Dominik³⁷,
A. Donini¹⁸, J. Dörner⁹⁷, M. Doro³³, R. D. C. dos Anjos⁸⁶, J.-L. Dournaux¹⁵, C. Duangchan⁹⁸,⁷³,
C. Dubos⁴⁸, G. Dubus⁹⁹, S. Duffy⁸³, D. Dumora⁸², V. V. Dwarkadas¹⁰⁰, J. Ebr⁵¹, C. Eckner⁴³,¹⁰¹,
K. Egberts³⁴, S. Einecke⁵⁹, D. Elsässer³⁷, G. Emery⁷⁷, M. Errando¹⁰², C. Escanuela⁴⁰,
P. Escarate¹⁰³,⁷⁶, M. Escobar Godoy¹⁰⁴, J. Escudero¹⁰, P. Esposito⁵⁶,⁵⁷, C. Evoli¹⁰⁵,
D. Falceta-Goncalves¹⁰⁶, A. Fattorini³⁷, S. Fegan⁸, K. Feijen⁸¹, Q. Feng³⁸, G. Ferrand¹⁰⁷,¹⁰⁸,
F. Ferrarotto¹⁰⁹, E. Fiandrini¹⁷, A. Fiasson⁴³, M. Filipovic¹¹⁰, V. Fioretti⁶⁰, M. Fiori¹¹¹, H. Flores¹⁵,
L. Foffano¹¹², L. Font Guiteras⁷¹, G. Fontaine⁸, S. Fröse³⁷, Y. Fukazawa¹¹³, Y. Fukui¹¹⁴, S. Funk⁹⁸,
A. Furniss¹⁰⁴, D. Gaggero⁵⁵, G. Galanti⁵⁷, G. Galaz¹³, Y. A. Gallant²², S. Gallozzi¹⁸, V. Gammaldi¹²,
M. Garczarczyk³⁰, C. Gasbarra¹¹⁵, D. Gasparrini¹¹⁵, M. Gaug⁷¹, A. Ghalumyan¹¹⁶, M. Giarrusso¹¹⁷,
J. Giesbrecht³¹, N. Giglietto⁴⁶,⁴⁷, F. Giordano⁸⁰, R. Giuffrida¹¹⁸,⁵², A. Giuliani⁵⁷, J.-F. Glicenstein⁶³,
J. Glombitza⁹⁸, N. Godinovic¹¹⁹, P. Goldoni¹²⁰, J. M. González¹²¹, J. Goulart Coelho¹²²,
J. Granot¹²³,¹²⁴, D. Grasso⁵⁵, R. Grau⁵⁰, L. Gréaux⁴⁸, D. Green⁴¹, J. G. Green⁴¹, T. Greenshaw¹²⁵,
I. Grenier¹²⁶, G. Grolleron⁵³, M.-H. Grondin⁸², O. Gueta³⁰,¹²⁷, S. Gunji¹²⁸, J. Hackfeld⁹⁷,³⁷,
D. Hadasch¹, W. Hanlon³⁸, S. Hara¹²⁹, V. M. Harvey⁵⁹, T. Hassan³⁹, K. Hayashi¹,¹³⁰,
L. Heckmann⁴¹, M. Heller⁶⁷, G. Hermann⁴⁰, S. Hernández Cadena¹³¹, O. Hervet¹⁰⁴, J. Hinton⁴⁰,
N. Hiroshima¹, B. Hnatyk¹³², R. Hnatyk¹³², W. Hofmann⁴⁰, J. Holder¹³³, M. Holler¹³⁴, D. Horan⁸,
P. Horvath¹³⁵, T. Hovatta¹³⁶, M. Hrabovsky¹³⁵, M. Iarlori¹³⁷, T. Inada¹, F. Incardona⁸⁷, S. Inoue¹⁰⁸, F. Iocco⁹⁵,²⁰, M. Iori¹⁰⁹, M. Jamrozy¹³⁸, P. Janecek⁵¹, F. Jankowsky¹³⁹, C. Jarnot¹⁴⁰, P. Jean¹⁴⁰, I. Jiménez Martínez³⁹, W. Jin¹⁴¹, C. Juramy-Gilles⁵³, J. Jurysek⁵¹, M. Kagaya¹, O. Kalekin⁹⁸, D. Kantzas¹⁰¹, V. Karas¹⁴², H. Katagiri¹⁴³, J. Kataoka¹⁴⁴, S. Kaufmann⁶², D. Kazanas¹⁴⁵, D. Kerszberg⁵⁰, B. Khélifi⁸¹, D. B. Kieda¹⁴⁶, R. Kissmann¹³⁴, T. Kleiner³⁰, G. Kluge¹⁴⁷, W. Kluźniak¹⁴⁸, J. Knödlseder¹⁴⁰, Y. Kobayashi¹, K. Kohri¹⁴⁹, N. Komin³⁶, P. Kornecki¹⁵, K. Kosack⁴, D. Kostunin³⁰, G. Kowal¹⁰⁶, H. Kubo¹, J. Kushida¹⁵⁰, A. La Barbera⁹¹, N. La Palombara⁵⁷, M. Láinez¹⁴, A. Lamastra¹⁸, J. Lapington¹⁵¹, P. Laporte¹⁵, S. Lazarević¹¹⁰, J. Lazendic-Galloway²⁷, M. Lemoine-Goumard⁸², J.-P. Lenain⁵³, F. Leone¹⁵², G. Leto⁸⁷, F. Leuschner⁴⁵, E. Lindfors¹³⁶, M. Linhoff³⁷, I. Liodakis¹³⁶, S. Lombardi¹⁸, F. Longo¹⁵³, R. López-Coto¹⁰, M. López-Moya¹⁴, A. López-Oramas³⁵, S. Loporchio⁴⁶,⁴⁷, J. Lozano Bahilo⁷⁹, F. Lucarelli¹⁸, P. L. Luque-Escamilla¹⁵⁴, E. Lyard²⁸, O. Macias¹⁵⁵, J. Mackey⁶¹, G. Maier³⁰, D. Malyshev⁴⁵, D. Mandat⁵¹, G. Manicò¹¹⁷,¹⁵², A. Marcowith²², P. Marinos⁵⁹, M. Mariotti³³, S. Markoff¹⁵⁵, P. Marquez⁵⁰, G. Marsella¹¹⁸,¹¹⁷, J. Martí¹⁵⁴, P. Martin¹⁴⁰, G. A. Martínez³⁹, M. Martínez⁵⁰, O. Martinez¹⁵⁶,¹⁵⁷, C. Marty¹⁴⁰, A. Mas-Aguilar¹⁴, M. Mastropietro¹⁸, G. Maurin⁴³, D. Mazin¹,⁴¹, S. McKeague⁸³, A. J. T. S. Mello⁸⁶,¹⁵⁸, S. Menchiari¹⁶, S. Mereghetti⁵⁷, E. Mestre¹⁵⁹, J.-L. Meunier⁵³, D. M.-A. Meyer³⁴, D. Miceli²¹, M. Miceli¹¹⁸,⁵², M. Michailidis⁴⁵, J. Michałowski¹⁶⁰, T. Miener⁶⁷, J. M. Miranda¹⁵⁶,¹⁶¹, A. Mitchell⁹⁸, T. Mizuno¹⁶², R. Moderski¹⁴⁸, L. Mohrmann⁴⁰, M. Molero³⁵, C. Molfese⁹⁴, E. Molina³⁵, T. Montaruli⁶⁷, A. Moralejo⁵⁰, D. Morcuende¹⁰, K. Morik³⁷, G. Morlino¹⁶, A. Morselli¹¹⁵, E. Moulin⁶³, V. Moya Zamanillo¹⁴, R. Mukherjee¹⁶³, K. Munari⁸⁷, T. Murach³⁰, A. Muraczewski¹⁴⁸, H. Muraishi¹⁶⁴, S. Nagataki¹⁰⁸, T. Nakamori¹²⁸, R. Nemmen³²,¹⁶⁵, L. Nickel³⁷, J. Niemiec¹⁶⁰, D. Nieto¹⁴, M. Nievas Rosillo³⁵, M. Nikołajuk¹⁶⁶, L. Nikolić³, K. Noda¹⁶⁷,¹, D. Nosek¹⁶⁸, B. Novosyadlyj¹⁶⁹, V. Novotny¹⁶⁸, S. Nozaki⁴¹, M. Ohishi¹, Y. Ohtani¹, A. Okumura¹⁷⁰,¹⁷¹, J.-F. Olive¹⁴⁰, B. Olmi⁵²,¹⁶, R. A. Ong¹⁴¹, M. Orienti⁹³, R. Orito¹⁷², M. Orlandini⁶⁰, E. Orlando¹⁵³, S. Orlando⁵², M. Ostrowski¹³⁸, I. Oya⁷⁰, I. Pagano⁸⁷, A. Pagliaro⁹¹, M. Palatiello⁶⁶, G. Panebianco⁶⁰, D. Paneque⁴¹, F. R. Pantaleo⁴⁷,⁴⁶, R. Paoletti³, J. M. Paredes⁹, N. Parmiggiani⁶⁰, S. R. Patel⁴⁸, B. Patricelli¹⁸,¹⁷³, D. Pavlović¹⁷⁴, M. Pech⁵¹, M. Pecimotika¹⁷⁴,¹⁷⁵, M. Peresano⁸⁵,⁸⁴, J. Pérez-Romero⁴⁴, M. A. Pérez-Torres¹⁰, G. Peron⁸¹, M. Persic¹¹¹,¹⁷⁶, P.-O. Petrucci⁹⁹, O. Petruk⁴², G. Piano¹¹², E. Pierre⁵³, E. Pietropaolo¹⁷⁷, M. Pihet²¹, F. Pintore⁹¹, C. Pittori¹⁸, C. Plard⁴³, F. Podobnik³, M. Pohl³⁴,³⁰, E. Pons⁴³, G. Ponti⁵⁴, E. Prandini³³, G. Principe¹⁵³, C. Priyadarshi⁵⁰, N. Produit²⁸, D. Prokhorov¹⁵⁵, E. Pueschel³⁰, G. Pühlhofer⁴⁵, M. L. Pumo¹⁵²,¹¹⁷, M. Punch⁸¹, F. Queiroz¹⁷⁸,¹⁷⁹, A. Quirrenbach¹³⁹, R. Rando³³, T. Ravel¹⁴⁰, S. Razzaque¹⁸⁰,¹²⁴, M. Regeard⁸¹, P. Reichherzer⁸⁸,⁹⁷, A. Reimer¹³⁴, O. Reimer¹³⁴, Q. Remy⁴⁰, M. Renaud²², T. Reposeur⁸², W. Rhode³⁷, D. Ribeiro¹⁸¹, M. Ribó⁹, T. Richtler¹⁸², J. Rico⁵⁰, F. Rieger⁴⁰, M. Rigoselli⁵⁷, V. Rizi¹⁷⁷, E. Roache³⁸, G. Rodriguez Fernandez¹¹⁵, J. J. Rodríguez-Vázquez³⁹, P. Romano⁵⁴, G. Romeo⁸⁷, J. Rosado¹⁴, A. Rosales de Leon⁵³, G. Rowell⁵⁹, B. Rudak¹⁴⁸, A. J. Ruiter¹⁸³, C. B. Rulten⁶², F. Russo⁶⁰, I. Sadeh³⁰, L. Saha³⁸, T. Saito¹, H. Salzmann⁴⁵, M. Sánchez-Conde¹², P. Sangiorgi⁹¹, H. Sano¹, M. Santander⁶, A. Santangelo⁴⁵, R. Santos-Lima³², V. Sapienza⁵²,¹¹⁸, T. Šarić¹¹⁹, S. Sarkar⁸⁸, F. G. Saturni¹⁸, A. Scherer¹³, F. Schiavone⁸⁰, P. Schipani⁷⁸, B. Schleicher¹⁸⁴,², P. Schovanek⁵¹, J. L. Schubert³⁷, F. Schussler⁶³, U. Schwanke¹⁸⁵, G. Schwefer⁴⁰, M. Seglar Arroyo⁵⁰, I. Seitenzahl¹⁸³, O. Sergijenko¹³²,¹⁸⁶,¹⁸⁷, M. Servillat¹⁵, V. Sguera⁶⁰, P. Sharma⁴⁸, H. Siejkowski¹⁸⁸, C. Siqueira²⁴, P. Sizun¹⁸⁹, V. Sliusar²⁸, A. Slowikowska¹⁹⁰, H. Sol¹⁵, S. T. Spencer⁹⁸,⁸⁸, D. Spiga⁵⁴,
A. Stamerra<sup>18,127</sup>, S. Stanič<sup>44</sup>, R. Starling<sup>151</sup>, Ł. Stawarz<sup>138</sup>, S. Steinmassl<sup>40</sup>, C. Steppa<sup>34</sup>,
T. Stolarczyk<sup>4</sup>, Y. Suda<sup>113</sup>, T. Suomijärvi<sup>48</sup>, H. Tajima<sup>170,171</sup>, R. Takeishi<sup>1</sup>, S. J. Tanaka<sup>191</sup>,
F. Tavecchio<sup>54</sup>, T. Tavernier<sup>51</sup>, Y. Terada<sup>192</sup>, R. Terrier<sup>81</sup>, M. Teshima<sup>41</sup>, W. W. Tian<sup>1</sup>,
L. Tibaldo<sup>140</sup>, O. Tibolla<sup>62</sup>, F. Torradeflot<sup>193,39</sup>, D. F. Torres<sup>159</sup>, N. Tothill<sup>110</sup>, F. Toussenel<sup>53</sup>,
V. Touzard<sup>140</sup>, P. Travnicek<sup>51</sup>, G. Tripodo<sup>118,117</sup>, A. Trois<sup>194</sup>, A. Tsiahina<sup>140</sup>, A. Tutone<sup>91</sup>,
G. Umana<sup>87</sup>, L. Vaclavek<sup>135</sup>, M. Vacula<sup>135,51</sup>, P. Vallania<sup>84,195</sup>, C. van Eldik<sup>98</sup>, V. Vassiliev<sup>141</sup>,
M. L. Vazquez Acosta<sup>35</sup>, M. Vecchi<sup>23</sup>, S. Ventura<sup>3</sup>, S. Vercellone<sup>54</sup>, G. Verna<sup>3</sup>, A. Viana<sup>24</sup>,
N. Viaux<sup>196</sup>, A. Vigliano<sup>66</sup>, J. Vignatti<sup>196</sup>, C. F. Vigorito<sup>84,85</sup>, J. Villanueva<sup>89</sup>, J. Vink<sup>155</sup>,
V. Vitale<sup>115</sup>, V. Vodeb<sup>44</sup>, V. Voisin<sup>53</sup>, S. Vorobiov<sup>44</sup>, G. Voutsinas<sup>67</sup>, I. Vovk<sup>1</sup>, T. Vuillaume<sup>43</sup>,
V. Waegebaert<sup>140</sup>, S. J. Wagner<sup>139</sup>, R. Walter<sup>28</sup>, M. Wechakama<sup>72,73</sup>, R. White<sup>40</sup>,
A. Wierzcholska<sup>160</sup>, D. A. Williams<sup>104</sup>, F. Wohlleben<sup>40</sup>, R. Yamazaki<sup>191</sup>, L. Yang<sup>180,197</sup>,
T. Yoshida<sup>143</sup>, T. Yoshikoshi<sup>1</sup>, M. Zacharias<sup>139,26</sup>, G. Zaharijas<sup>44</sup>, L. Zampieri<sup>111</sup>, R. Zanin<sup>70</sup>,
D. Zavrtanik<sup>44</sup>, M. Zavrtanik<sup>44</sup>, A. A. Zdziarski<sup>148</sup>, A. Zech<sup>15</sup>, V. I. Zhdanov<sup>132</sup>, K. Zietara<sup>138</sup>,
M. Živec<sup>44</sup>, J. Zuriaga-Puig<sup>12</sup>, P. De la Torre Luque<sup>198</sup>, L. Guillemot<sup>199,200</sup>, D. A. Smith<sup>82,201</sup>

<footnote n="1">¹ Institute for Cosmic Ray Research, University of Tokyo, 5-1-5, Kashiwa-no-ha, Kashiwa, Chiba 277-8582, Japan</footnote>

<footnote n="2">² ETH Zürich, Institute for Particle Physics and Astrophysics, Otto-Stern-Weg 5, 8093 Zürich, Switzerland</footnote>

<footnote n="3">³ INFN and Università degli Studi di Siena, Dipartimento di Scienze Fisiche, della Terra e dell'Ambiente (DSFTA), Sezione di Fisica, Via Roma 56, 53100 Siena, Italy</footnote>

<footnote n="4">⁴ Université Paris-Saclay, Université Paris Cité, CEA, CNRS, AIM, F-91191 Gif-sur-Yvette Cedex, France</footnote>

<footnote n="5">⁵ FSLAC IRL 2009, CNRS/IAC, La Laguna, Tenerife, Spain</footnote>

<footnote n="6">⁶ University of Alabama, Tuscaloosa, Department of Physics and Astronomy, Gallalee Hall, Box 870324 Tuscaloosa, AL 35487-0324, U.S.A.</footnote>

<footnote n="7">⁷ Université Côte d'Azur, Observatoire de la Côte d'Azur, CNRS, Laboratoire Lagrange, France</footnote>

<footnote n="8">⁸ Laboratoire Leprince-Ringuet, CNRS/IN2P3, École polytechnique, Institut Polytechnique de Paris, 91120 Palaiseau, France</footnote>

<footnote n="9">⁹ Departament de Física Quàntica i Astrofísica, Institut de Ciències del Cosmos, Universitat de Barcelona, IEEC-UB, Martí i Franquès, 1, 08028, Barcelona, Spain</footnote>

<footnote n="10">¹⁰ Instituto de Astrofísica de Andalucía-CSIC, Glorieta de la Astronomía s/n, 18008, Granada, Spain</footnote>

<footnote n="11">¹¹ Institute for Computational Cosmology and Department of Physics, Durham University, South Road, Durham DH1 3LE, United Kingdom</footnote>

<footnote n="12">¹² Instituto de Física Teórica UAM/CSIC and Departamento de Física Teórica, Universidad Autónoma de Madrid, c/ Nicolás Cabrera 13-15, Campus de Cantoblanco UAM, 28049 Madrid, Spain</footnote>

<footnote n="13">¹³ Pontificia Universidad Católica de Chile, Av. Libertador Bernardo O'Higgins 340, Santiago, Chile</footnote>

<footnote n="14">¹⁴ IPARCOS-UCM, Instituto de Física de Partículas y del Cosmos, and EMFTEL Department, Universidad Complutense de Madrid, E-28040 Madrid, Spain</footnote>

<footnote n="15">¹⁵ LUTH, GEPI and LERMA, Observatoire de Paris, Université PSL, Université Paris Cité, CNRS, 5 place Jules Janssen, 92190, Meudon, France</footnote>

<footnote n="16">¹⁶ INAF — Osservatorio Astrofisico di Arcetri, Largo E. Fermi, 5 — 50125 Firenze, Italy</footnote>

<footnote n="17">¹⁷ INFN Sezione di Perugia and Università degli Studi di Perugia, Via A. Pascoli, 06123 Perugia, Italy</footnote>

<footnote n="18">¹⁸ INAF — Osservatorio Astronomico di Roma, Via di Frascati 33, 00040, Monteporzio Catone, Italy</footnote>

<footnote n="19">¹⁹ TÜBİTAK Research Institute for Fundamental Sciences, 41470 Gebze, Kocaeli, Turkey</footnote>

<footnote n="20">²⁰ INFN Sezione di Napoli, Via Cintia, ed. G, 80126 Napoli, Italy</footnote>

<footnote n="21">²¹ INFN Sezione di Padova, Via Marzolo 8, 35131 Padova, Italy</footnote>

<footnote n="22">²² Laboratoire Univers et Particules de Montpellier, Université de Montpellier, CNRS/IN2P3, CC 72, Place Eugène Bataillon, F-34095 Montpellier Cedex 5, France</footnote>

<footnote n="23">²³ Kapteyn Astronomical Institute, University of Groningen, Landleven 12, 9747 AD, Groningen, The Netherlands</footnote>

<footnote n="24">²⁴ Instituto de Física de São Carlos, Universidade de São Paulo, Av. Trabalhador São-carlense, 400 — CEP 13566-590, São Carlos, SP, Brazil</footnote>