© C. Steppa and K. Egberts 2020

1. Introduction

The past two decades have witnessed the birth and explosive development of teraelectronvolt astronomy. A major breakthrough for the development of the field and especially the Galactic very-high-energy (VHE; 100 GeV ≤ E ≤ 100 TeV) γ-ray sky has been the High Energy Stereoscopic System (H.E.S.S.) Galactic plane survey (HGPS). For 12 years H.E.S.S. has scanned the central part of the Milky Way (extending from Galactic longitudes of l = 250° to 65° and covering latitudes of |b| ≤ 3°) and acquired a data set of nearly 2700 h of good-quality observations (H.E.S.S. Collaboration 2018a). The HGPS has revealed a plethora of γ-ray sources (H.E.S.S. Collaboration 2018a) and a faint component of a so-called diffuse (large-scale unresolved) emission (H.E.S.S. Collaboration 2014). With the brightest and closest sources, the sample of detected γ-ray sources represents only the tip of the iceberg of the overall population of VHE γ-ray emitters. A larger percentage of sources are expected to remain unresolved with the given H.E.S.S. exposure and sensitivity due to being too faint and/or too far away to be significantly detected, thus forming a contribution to the measured large-scale diffuse emission.

Previous studies of the VHE-detected γ-ray source classes that are based on the HGPS characterise the sample of pulsar wind nebulae (PWNe; H.E.S.S. Collaboration 2018b) and supernova remnants (SNRs; H.E.S.S. Collaboration 2018c), but with limited insight into their respective population. A characterisation of the overall population of sources can be achieved by population synthesis with the simulation of synthetic source samples and comparison with observations in the range of detectability of the data set. This procedure is customarily followed for the study of object classes such as pulsars (Gonthier et al. 2018). The study of a generic source population is a slightly different approach. Rather than aiming to derive properties of a specific class of objects, this approach characterises the overall population of sources at a certain wavelength. This procedure allows the prediction of the number of sources detectable with future instruments (e.g. the Cherenkov Telescope Array CTA; Funk et al. 2008) and to estimate the amount of unresolved sources that contribute to the diffuse emission measurements.

In this work, we follow this latter strategy to describe the VHE γ-ray emitters generically, and we derive luminosity as well as radius functions for this generic VHE γ-ray source population. A similar approach has already been applied to data from the Energetic Gamma Ray Experiment Telescope on board the Compton Gamma Ray Observatory and data from the Large Area Telescope on board the Fermi Gamma-ray Space Telescope (Fermi-LAT) for the estimation of unresolved sources in the high-energy (HE; 100 MeV ≤ E ≤ 100 GeV) Galactic diffuse γ-ray emission (Strong 2007; Acero et al. 2015). As already identified by Casanova & Dingus (2008) for the case of the diffuse emission measured by MILAGRO (Abdo et al. 2008), compared to HE the contribution of unresolved sources is expected to rise, and likely dominates, at VHE. This is also reflected in recent measurements of Galactic diffuse emission above 1 TeV by the High-Altitude Water Cherenkov Observatory (HAWC; Nayerhoda et al. 2019) which, like the H.E.S.S. measurements at 1 TeV (H.E.S.S. Collaboration 2014), overshoot predictions considerably. Only an assessment of the entire Galactic source population can disentangle the two components of unresolved γ-ray-source emission and diffuse emission from propagating cosmic rays and allow for the study of cosmic-ray propagation properties in the H.E.S.S. and HAWC data sets.

2. Construction of the model

The VHE source population model presented in this work consists of two distinct components: the spatial distribution of sources and distribution of source properties, that is their luminosities and radii. To determine the spatial distribution, we tested various models based on the assumed source classes and the Galactic structure. We followed a data-driven approach to derive the distribution of source properties. Alternatively, this distribution can be derived from detailed source modelling. However, a population model based on individual source models involves a fair amount of assumptions, for instance about source classes contained in the population, the age of these sources, and their environmental conditions. In contrast, for the data-driven approach we only assume that the source sample is representative for the population in its range of detectability and that sources are distributed according to certain spatial models. In the following, the derivation of each component of the model is described in detail. After assessing the spatial distribution, the derivation of the second component is described, which is based on a combination of the spatial model with observed quantities of the HGPS source sample, namely integrated flux above 1 TeV (henceforth referred to as flux), angular extent, and location in the sky.

2.1. Spatial distribution

The number of detected Galactic VHE γ-ray sources is yet too small to determine the spatial distribution of the entire population. However, it is possible to construct models of the spatial distribution based on few reasonable assumptions. Since most of the known sources are associated with SNRs or PWNe, the corresponding distributions of SNRs and pulsars can be used as templates. Their source densities ρ are well described by an azimuth-symmetric function that only depends on the Galactocentric distance r and the height over the Galactic disc z as follows:

$$\begin{array}{c}\hfill \rho (r,z)={\left(\frac{(r+r_{\mathrm{off}})}{(R_{\odot }+r_{\mathrm{off}})}\right)}^{\alpha }exp\left(-\beta \frac{(r-R_{\odot })}{(R_{\odot }+r_{\mathrm{off}})}\right)exp\left(-\frac{\left|z\right|}{z_0}\right),\end{array}$$(1)

where R_⊙ is the distance of the sun to the Galactic centre, z₀ the scale height of the Galactic disc, shape parameter α, and rate parameter β. The parameter r_off accounts for a non-zero density at r = 0. Based on the assumption that SNRs are the dominant class of γ-ray sources we probe a model (mSNR) by applying Eq. (1) with parameters as given in Green (2015) and Xu et al. (2005). Likewise, we probe a model that is based on the assumption that PWNe are the dominant class (mPWN) using parameters as given in Yusifov & Küçük (2004) and Lorimer et al. (2006). Both parameter sets are listed in Table 1.

While mSNR and mPWN are two examples for azimuth-symmetric source distributions, there is good reason to assume that the spatial distribution of γ-ray sources might deviate from this symmetry. The progenitors of VHE sources, for instance massive stars, typically form in dense regions of gas and dust. Therefore, the distribution of VHE sources might be affected by the spiral structure of the Galaxy, for instance, observed in the distribution of interstellar matter (ISM; Steiman-Cameron et al. 2010). Following the study of Kissmann et al. (2015) on the impact of spiral-arm source distributions on the Galactic cosmic-ray flux we probed three different models of a non-symmetric source distribution. To represent a four-arm distribution, we adopted the model by Steiman-Cameron et al. (2010), denoted as mSp4. Compared to Eq. (1), in this case the source density explicitly depends on the azimuth ϕ and is described as follows:

$$\begin{array}{c}\hfill \begin{array}{cc}\hfill \rho (r,\varphi ,z)=& {\displaystyle \sum _{i=1}^4}{A}_{i}exp(-\frac{1}{{\delta }^2}{(\varphi -\frac{ln\left(\frac{r}{a_i}\right)}{{\beta }_i})}^2)\hfill \\ \hfill & \times exp(-\frac{|r-R|}{{\sigma }_{\mathrm{r}}})exp(-\frac{z^2}{2{\sigma }_{z,2}^2})·\hfill \end{array}\end{array}$$(2)

The radial dependence is defined by the scale length σ_r and a local maximum at R. Likewise, the azimuthal dependence is defined by the scale length δ and two constants, β_i determining the pitch angle of the spiral and a_i giving its orientation. Finally, the z-dependence is governed by the scale height σ_z, 2. For this model we adopted the best-fit values corresponding to the ISM measurement traced by CII cooling lines (see Table 2). We adopted another four-arm model with different, less pronounced arm profiles from Cordes & Lazio (2002), which is based on the free electron density as traced by pulsar dispersion measurements. We refer to this model as mFE. To calculate the source density for this model we made use of the code provided by the authors¹. Finally, we also probed a two-arm model with an additional central bar, whose existence in our Galaxy has been indicated by Spitzer data (Benjamin et al. 2005) and more recently confirmed by Gaia data (Anders et al. 2019). This model, which we refer to as mSp2B, was adopted from Werner et al. (2015). It is based on the model of Steiman-Cameron et al. (2010) but only includes the Scutum-Crux and Perseus arms. The additional component for the Galactic bar is given by

$$\begin{array}{c}\hfill {\rho }_{\mathrm{bar}}(r,\varphi ,z)=\{\begin{array}{cc}{A}_{\mathrm{bar}}exp(-\frac{z^2+r^2{(sin(\varphi )-cos(\varphi )sin(\theta ))}^2}{{\sigma }_{z,1}^2}),\hfill & \mathrm{if}r<l_b\hfill \\ 0,\hfill & \mathrm{otherwise}\hfill \end{array}\end{array}$$(3)

with its radial extent l_b = 3.5 pc, rotation angle relative to the solar-galactic centre line θ = 30.0 deg and scale height σ_z, 1 = 0.31 pc. The normalisation is chosen such that the bar is equally contributing to the source density as the spiral arms (A_bar = 364).

For our study we simulated the source distribution by using a set S_sim of ∼2 × 10⁶ uniformly, randomly distributed points in a box of size 30 kpc × 30 kpc × 10 kpc yielding a mean distance among points of ∼35 pc. Each point is then weighted according to the source density given by the tested model².

2.2. Luminosity and radius distribution

Next to the spatial distribution, each model comprises a distribution function for source properties, namely luminosity (L) and radius (R). Here we assume that the variables L and R are independent and each one follows a power law such that the joint probability density function (PDF) P(L, R) can be written as

$$\begin{array}{c}\hfill P(L,R)=N{\left(\frac{L}{L_0}\right)}^{{\alpha }_{\mathrm{L}}}{\left(\frac{R}{R_0}\right)}^{{\alpha }_{\mathrm{R}}},\end{array}$$(4)

with scaling factors L₀, R₀ and a normalisation factor N that depends on the boundaries set for L and R. The number of detected sources N_det is related to Eq. (4) via

$$\begin{array}{c}\hfill N_{\det }=N_{\mathrm{FoV}}{\displaystyle \int }\mathrm{d}L{\displaystyle \int }\mathrm{d}RC(L,R)P(L,R),\end{array}$$(5)

where the observation bias inherent to the sample of detected sources is accounted for by the correction function C(L, R) and the total number of sources in the probed field of view (FOV) is N_FoV. In order to reconstruct the parameters of the global distribution function from the biased sample of detected sources, we applied a likelihood maximisation as follows. Dividing the probed L × R phase space in equally sized bins of 0.1 × 0.1 on logarithmic scale, we derived the true number of detected sources N_true, i that lie within a bin i of this phase space from the HGPS catalogue. The expected number of sources N_pred, i is approximated via

$$\begin{array}{c}\hfill N_{\mathrm{pred},i}=C(L_i,R_i)N_{\mathrm{FoV}}{{\displaystyle \int }}_{L_{\min ,i}}^{L_{\max ,i}}\mathrm{d}L{{\displaystyle \int }}_{R_{\min ,i}}^{R_{\max ,i}}\mathrm{d}RP(L,R),\end{array}$$(6)

where the correction function is only evaluated at the centre of the respective bin. With this type of counting exercise the distribution of the true number of detected sources per bin is expected to follow a Poissonian P_λ(N_true) with λ = N_pred. Thus, the log likelihood for the maximisation is

$$\begin{array}{c}\hfill logL={\displaystyle \sum _i}log\left(P_{N_{\mathrm{pred},i}}\left(N_{\mathrm{true},i}\right)\right),\end{array}$$(7)

where we sum over the grid of bins. We applied this method to derive the two power-law indices of the joint PDF P(L, R)³. The distinct feature of this procedure is the recognition and inclusion of the observation bias corresponding to the analysed sample of sources and in particular the consideration of its dependency on the radius of sources, which is calculated in Sect. 2.2.2.

2.2.1. Source selection

Although the HGPS catalogue results from the most systematic search for VHE γ-ray sources to date, it suffers from several deviations of its generation from a fully automated pipeline. These deviations are a result of the large extent of many VHE γ-ray sources and their associated complex morphologies and complicate a treatment in a population-synthesis approach. The HGPS combines sources that are detected by a fixed pipeline based on maps of the detection significance for a correlation radius of R_c = 0.1° or R_c = 0.2°, plus sources labelled as external, which are detected and characterised by custom-tailored analyses (H.E.S.S. Collaboration 2018a). Furthermore, the Gaussian components obtained by the automated detection pipeline are then manually merged into γ-ray sources, resulting in a description of their complex morphologies as a combination of various Gaussians. These HGPS procedures render a rigorous treatment of the data in a population synthesis almost impossible. With the approach followed in the study presented in this work, we limited ourselves to extended sources with known distance in order to derive their L − R distribution. Point-like sources were excluded from the sample because they lack the extension information that is necessary for the determination of P(L, R). The criterion of a distance estimation being available guarantees that flux and angular size measurements can be transformed into L and R values.

For simplification we treated extended sources as being observed as symmetric two-dimensional Gaussians in the projected plane on the sky. The angular extent σ_source refers to the 68% containment radius of the measured flux. Sources with complex morphologies, for instance shell-like structures, were treated the same way, in which the shell radius is taken as σ_source. According to the extended nature of these sources the sensitivity map with correlation radius of R_c = 0.2° was used to describe the correction function that accounts for the observation bias. Out of the 78 sources in the HGPS, 64 pass the sensitivity threshold with R_c = 0.2°. Selecting for extended sources reduces this sample to 50 sources. From those 50 sources only 16 (32%) are firmly identified sources with available distance estimates. From this small sample we derived the parameters of the luminosity-radius PDF. The distributions of observable quantities, that is flux, extent, and composition, of the sample of these 16 selected sources do not deviate significantly from the distributions derived for the complete sample of sources except for the missing class of binaries⁴. Thus, it is assumed that the small sample is representative and distance estimates are independent of the luminosities and radii of sources. In addition, the boundaries on L and R for the phase space that we probe are derived from this data set as well, yielding 10^32.1 ph s⁻¹ ≤ L ≤ 10^34.8 ph s⁻¹ and 10^0.4 pc ≤ R ≤ 10^1.7 pc. Thus, the presented models cover a dynamical range of almost three orders of magnitude in the luminosity (cf. with a dynamical range of three orders of magnitude chosen by Strong 2007 and five orders of magnitude chosen by Acero et al. 2015) and one order of magnitude in size. The scale factors in Eq. (4) are set to R₀ = 1 pc and L₀ = 10³⁴ ph s⁻¹. As can be derived from the units, we only considered the number of emitted photons above 1 TeV per second as proxy for the luminosity. The luminosity function in units of erg s⁻¹ can be derived by scaling with a characteristic mean photon energy. To calculate this mean photon energy an additional assumption about the shape of the spectral energy distribution (SED) of sources is required. Assuming the SED follows a power law, a mean spectral index of −2.4 is found for the HGPS source sample, which yields a mean photon energy of 3.52 erg ph⁻¹ in the energy range 1 TeV−10 TeV.

2.2.2. Correction function determination

Derivation of the VHE γ-ray source population properties from observational data needs to account for the strong selection bias in the H.E.S.S. catalogue, which can be expected to distinctly shape the sample of detected γ-ray sources and is based on the HGPS sensitivity. The HGPS sensitivity varies strongly as a function of Galactic longitude and latitude due to the observation pattern, which is a combination of dedicated survey observations and additional follow-up observations and in-depth measurements of detected sources. The inhomogeneity of the HGPS sensitivity as a function of Galactic longitude and latitude is demonstrated in Fig. 1 through the detection horizon of point-like sources with a luminosity of 10³³ ph s⁻¹.

Fig. 1. Distribution of VHE sources in the Milky Way following the four-arm spiral model by Steiman-Cameron et al. (2010) as described in Sect. 2 (grey contours). Additionally, the HGPS detection horizon for point-like sources with a luminosity of 1033 ph s−1 is plotted on top (blue contours).

Besides the direction dependency, the sensitivity is also a function of the angular extent of a source since the number of background events increases with σ_source. To exceed the detection threshold of 5σ above background, the flux of an extended source needs to be greater than

$$\begin{array}{c}\hfill F_{\min }({\sigma }_{\mathrm{source}})=\{\begin{array}{cc}{F}_{\min ,0}\sqrt{\frac{{\sigma }_{\mathrm{source}}^2+{\sigma }_{\mathrm{PSF}}^2}{{\sigma }_{\mathrm{PSF}}^2}},\hfill & \le 1^{°}\hfill \\ \infty ,\hfill & >1^{°}\hfill \end{array}\end{array}$$(8)

(H.E.S.S. Collaboration 2018a), where F_min, 0 is the point-source sensitivity, σ_source the source extent, and σ_PSF the size of the H.E.S.S. point-spread function (PSF). In addition, the limited FOV of the instrument in combination with the applied background subtraction technique of deriving background measurements from within the FOV renders sources ≳1° not detectable. For our purpose we adopted Eq. (8) to account for the fact that we only selected for extended sources. The HGPS does not provide a criterion for the minimal detectable extent of a source. Therefore, we set the threshold to the cited average value of the PSF (σ_PSF = 0.08°), which yields a sensitivity for extended sources that is written as

$$\begin{array}{c}\hfill F_{\min }^{\mathrm{extended}}({\sigma }_{\mathrm{source}})=\{\begin{array}{cc}\infty ,\hfill & {\sigma }_{\mathrm{source}}\le {\sigma }_{\mathrm{PSF}}\hfill \\ F_{\min ,0}\sqrt{\frac{{\sigma }_{\mathrm{source}}^2+{\sigma }_{\mathrm{PSF}}^2}{{\sigma }_{\mathrm{PSF}}^2}},\hfill & {\sigma }_{\mathrm{PSF}}<{\sigma }_{\mathrm{source}}\le 1^{°}\hfill \\ \infty ,\hfill & 1^{°}<{\sigma }_{\mathrm{source}}.\hfill \end{array}\end{array}$$(9)

We define the correction function as the fraction of detectable sources, namely sources within the adopted sensitivity range of the HGPS, in the total amount of sources in the FoV of the HGPS. Since the sensitivity decreases for faint or extended sources, the correction function depends on source properties, that is luminosity and radius C(L, R).

Under the assumption that average source properties are identical throughout the Milky Way, we expect sources of any given properties to follow the same spatial distribution function. To calculate the correction function, we first derive the subset of sources which lie within the FOV of the HGPS from the set of simulated sources S_FoV ⊂ S_sim. Assigning each source the same luminosity L and source radius R we can then calculate the corresponding fluxes, angular extents, and locations in the sky as they would be observed at Earth. Based on these observables and the sensitivity limit (Eq. (9)) of the HGPS, the subset of detectable sources for the given luminosity and radius can be derived S_det(L, R)⊂S_FoV. With this, a two-dimensional correction function C(L_i, R_i) is derived for the same grid of source luminosities and radii mentioned above via

$$\begin{array}{c}\hfill {\displaystyle C(L_i,R_i)=0.32·\sum _{\mathit{x}\in S_{\det }(L_i,R_i)}\rho \left(\mathit{x}\right)/\sum _{\mathit{x}\in S_{\mathrm{FoV}}}\rho \left(\mathit{x}\right),}\end{array}$$(10)

where ρ is the source density corresponding to the assumed model. The additional factor 0.32 accounts for the fact that 68% of the detected sources are disregarded owing to missing distance estimates. The correction function for the spatial distribution model mSp4 is shown in Fig. 2 together with the distribution of HGPS sources that fulfil our selection criterion. Based on this estimation we show, for instance, that the distribution of the luminosity for L <  10³³ ph s⁻¹ cannot be well constrained by the HGPS data. The estimation of the luminosity distribution in this regime is affected by statistical fluctuations of the data and special care has to be taken to explore the range of validity of the model.

Fig. 2. Estimated correction function for mSp4 showing the fraction of sources in the FOV of the HGPS that for a given luminosity and radius are detectable according to Eq. (9). The low values for small sources result from the rejection of point-like sources in the analysis. Additionally those sources of the HGPS that fulfil the selection criterion are shown.

2.2.3. Monte Carlo verification

We studied the capability of reconstructing properties of the parent population of the method presented in this work by the means of Monte Carlo simulations. For this purpose, we simulated source populations with a set of luminosity and radius functions for all spatial models discussed previously. For each simulated population, N spatial coordinates x were randomly drawn following the distribution defined by Eqs. (1) and (2), or Eqs. (2) and (3), according to the spatial model, with model specific parameters, together with N random samples of L and R following Eq. (4) with the given parameters. For each combination of a luminosity and radius function and a source distribution, the subsample of detected sources was determined according to Eq. (9). As for the 50 extended sources detected with the HGPS the statistics of the simulated source population was adjusted to yield on average the same number of detected sources. To account for the fact that only a fraction of 32% of the HGPS sources comes with a distance estimation, all but 16 of the detected simulated sources were randomly discarded. Each data set was then reconstructed using the machinery discussed before and the indices α_L and α_R of luminosity and radius functions were calculated. For each choice of spatial model and luminosity and radius function 600 populations were simulated and reconstructed. We performed these tests for all combinations of the indices α_L = −3, −2, −1 and α_R = −2, −1, 0. For model mSp4 the mean of the reconstructed α values and their standard deviations are listed exemplarily in Table 3. For any combination, the mean of the reconstructed α agrees with the true value within Δα <  0.1. These results are remarkably consistent between spatial models. In particular, for all models and combinations of α_L, true and α_R, true, the reconstructed α_L and α_R are always compatible with the true values of the simulations. An example is given in Fig. 3 showing for all spatial models the distribution of reconstructed α_L and α_R given true values of α_L, true = −2 and α_R, true = −1. The reconstructed values are always centred at their respective true values, with standard deviations around 0.3 and 0.6. We repeated this exercise for varying values to investigate the influence of the boundaries set on the luminosity and radius. No effect on the reconstructed α_L and α_R could be observed.

Fig. 3. Distribution of reconstructed αL (top) and αR (bottom) for 600 toy models using different spatial models and true values of αL, true = −2 and αR, true = −1. True values are denoted by the vertical line.

In Fig. 4, we show one-dimensional luminosity and radius distributions, which are derived from the 600 samples drawn from model mSp4 for α_L, true = −2 and α_R, true = −1. The blue line indicates the bin-wise mean of the luminosity and radius distribution for the whole population averaged over the 600 simulated samples, respectively. The blue shaded regions depicts the standard deviation accordingly. Likewise, the orange area shows the distribution for the selected samples, each comprised of 16 extended sources within the sensitivity range of the HGPS and with known L and R. Besides the large spread of the latter distribution it is obvious that the distribution of detected sources on average does not resemble the global distribution. Therefore, it is inevitable for the reconstruction to properly account for the observation bias. The result of our reconstruction method is indicated in green. The solid line shows a power law with the mean reconstructed index. The power law is normed to the total number of sources of the simulated populations and perfectly matches the input distribution. The green shaded area represents the bin-wise quartile deviation.

Fig. 4. One-dimensional luminosity distribution (top) and radius distribution (bottom) for simulated populations. Original distributions are given in blue; distribution of detected sources in orange and the bias-corrected reconstruction of the original distribution are shown in green. Details are given in the text.

2.2.4. Result

The derived power-law indices of P(L, R) for the parent source population of the HGPS sample under the assumption of different underlying spatial source distributions are listed in Table 4. Again the results are fairly consistent among the spatial distributions with average values of ⟨α_L⟩= − 1.77 and ⟨α_R⟩= − 1.26. Owing to the small sample size used in this study the errors on these reconstructed values are expected to be comparable to those listed in Fig. 3 (Δα_L ∼ 0.3, Δα_R ∼ 0.6). Besides this statistical uncertainty, an additional cause of error is the choice of the boundaries of L and R. In general, deviations of the results due to this choice are found to be less than the stated errors, while the upper bound of the luminosity affects the reconstructed values most. However, this upper bound is well constrained since, according to the HGPS sensitivity, the chance of detecting a source with L ≥ 10^34.8 ph s⁻¹ (disregarding the availability of a distance estimate) is close to one. In comparison, previous studies on the VHE source population not taking into account source sizes or the inhomogeneity of the sensitivity, have estimated the luminosity function to follow a power law with harder index (e.g. −1 >  α >  −1.5; Casanova & Dingus 2008).

3. Comparison with observable quantities

In order to probe the validity of the derived models, for each model we compared the distribution of observable quantities from simulated source populations, namely Galactic longitude and latitude, fluxes, and angular extents, with those from observations. For each model a set M of 3000 synthetic source populations S ∈ M were simulated. Because the distribution functions for source positions and source properties were fixed, we were able to estimate the total number of sources in the population based solely on the number of observed sources. Thus, it was not necessary to limit the analysed sample to extended sources, but we could increase the accuracy of our prediction by determining the expected number of detectable sources according to Eq. (8), including point-like sources, that match the 64 sources in the HGPS that pass this criterion. The numbers of sources per population for the individual models are listed in Table 6 and further discussed in Sect. 4.1. The distributions that are investigated in this section are solely derived from those detectable sources within the sensitivity range of the HGPS S_det ⊂ S according to Eq. (8).

3.1. Flux and angular extent

The observed flux and angular extent of a source both depend on the distance of the source to the observer. In addition to this observational correlation, an intrinsic correlation between the luminosity and radius of sources can shape the observed distribution of fluxes and angular extents. Instrumental selection effects, such as the dependency of the sensitivity on the angular extent, can shape this distribution as well. To account for correlations, we compared the observed two-dimensional flux-extent distribution by means of its PDFs against the model prediction. We derived the PDFs from kernel density estimations for which the optimal bandwidth for a Gaussian kernel was derived individually in the range [10⁻³, 10] via the GridSearchCV method from the Python package scikit-learn. The distributions are very similar among the models, thus only the result for mSp4 is shown in comparison with the HGPS distribution in Fig. 5. Contour lines indicate the 1σ, 2σ, and 3σ containment fraction of the derived PDFs for the observed distribution and the predicted distribution. For all distributions we observe an increase of the angular extent with flux, although this correlation appears less pronounced for the model distributions. Additionally, model distributions are noticeably wider than the observed distribution. This discrepancy is further reflected in the fraction of extended sources in the sample of detected sources. While for the HGPS we yield a fraction of 78% (50 out of 64) of extended sources, for the models we observe on average a fraction between 23% and 38%, the rest being point-like. This discrepancy might be an effect that is intrinsic to the HGPS. More point sources might be present in the data set but are “lost” in extended sources owing to source confusion (Ambrogi et al. 2016) or detected and later merged with an overlapping source. We did not account for these effects in our model. Alternatively, the model assumptions might not reflect reality and the number of point-like sources is overestimated: this might be an effect of the simplified definition of source extent, namely describing distinct and complex morphologies altogether by a single parameter. In addition, it is likely that the assumed independence of radius and luminosity and the power law for the source radius do not capture the true nature of this source property. Given the connection to the observation bias and especially the impact on the detectability of nearby sources, this relation is worth investigating in a follow-up study. Nevertheless, the flux-extent distribution of our models and observations are considered to be close enough to make reasonable predictions in the following.

Fig. 5. Kernel density estimation for the flux-extent distribution of the HGPS sample (grey) and mSp4 (orange). The contour lines indicate the 1σ, 2σ, and 3σ containment fraction. Flux is given in “crab units” (c.u.).

3.2. Spatial distribution of sources

The source distributions in Galactic longitude for the probed models are depicted in Fig. 6. The shaded region shows the standard deviation around the mean value of the different samples S_det while the solid line represents the mean. For comparison, the observed distribution of the HGPS is given by the grey bars. Furthermore, the source density of the respective model within the FOV of the HGPS is shown by the dotted line, which is scaled for better visibility (∫_FoVρ(x) dx = 2N_HGPS). On the top right panel of Fig. 6 the dynamic range of the sensitivity over Galactic longitude is shown, which is expressed by $F_{\min }^{-3/2}$, where F_min is the point-source sensitivity of the HGPS map at b = 0° and the corresponding longitude bin. Here, for a given luminosity $F_{\min }^{-3/2}$ is proportional to the sampled volume. For better visibility this distribution is scaled in the same way as the source densities.

The models are generally in good agreement with observations, although for mSp2B the longitude distribution appears to be somewhat too narrow as it falls off too steeply in the outskirts of the Galactic plane. The result for mSp2B suggests that if the ISM density distribution, which is used as proxy for the distribution of regions with high star formation rates, indeed follows the assumed shape, the VHE source population must comprise at least one source class that evolves outside those regions. In the central regions all models commonly tend to overpredict the actual source distribution. According to the models, this is the region of highest source density. Therefore, we can arguably expect that the detection of sources is also affected by source confusion and an increased background resulting from the existence of bright diffuse emission in the Galactic ridge region (H.E.S.S. Collaboration 2006); both effects lead to a deficit of detected sources in that region. A future iteration of the model, taking both effects into account, will be required to unambiguously test whether this discrepancy in the central region can be attributed to an inaccurate spatial source distribution. Regarding the distribution of the source density we observe that, for model mSp4, peaks in this distribution align well with peaks in the observed source distribution. This might be suggestive that the Galactic population of VHE γ-ray sources indeed follows a similar spiral structure as derived from ISM measurements. However, with the inhomogeneity of the sensitivity, which yields similar distributions for detectable sources independently of the probed source distribution model, it is not feasible to make strong claims. Quantitatively we investigate the compatibility between model predictions and observations by means of the Kolmogorov–Smirnov test statistic d_n as follows:

$$\begin{array}{c}\hfill d_n=\underset{x}{sup}|F_n(x)-F_0(x)|,\end{array}$$(11)

where F(x) is the cumulative distribution over the variable x (i.e. Galactic longitude). In this equation, F₀(x) is derived from the mean distribution of a given model as shown by the solid line in Fig. 6. For each simulated population d_n is calculated with the cumulative distribution F_n(x) of detectable sources, yielding the probability distribution P(d_n). From the cumulative distribution of observed sources we derive d_HGPS and calculate the p-value P(d_n ≥ d_HGPS). The values listed in Table 5 confirm that only the longitude distribution of model mSp2B is incompatible with observations at a level of significance of 5%.

Fig. 6. Comparison of the longitudinal source distribution between HGPS sources (grey bars) and predictions of different models (coloured lines). The dotted lines represent the modelled source densities within the HGPS region, arbitrarily scaled for the sake of visibility. In addition, the dynamic range of the HGPS sensitivity is shown on the top right panel. Details are given in the text.

Source distributions in Galactic latitude are presented in Fig. 7 in the same way as for the longitudinal distributions. While for the latitude distribution all models are compatible with observations according to the statistical test, we can see some obvious deviation in this plot. The number of observed sources is falling rapidly outwards from the Galactic disc. The most prominent feature of the observed distribution is an asymmetry towards the southern sky that is not covered by any of the assumed models. Given that here we average over a broad source distribution over Galactic longitude, source confusion can be assumed to play a less pronounced role and the simulations correctly reflect the observation bias. Thus, this asymmetry appears to be a real feature of the spatial source distribution (cf. e.g. Skowron et al. 2019), which is not accounted for in the models. Besides this, data show a stronger peak of the latitude distribution towards the Galactic equator compared to simulations. This is most notably the case for the model mFE, whose flatter distribution appears to be in tension with observation.

Fig. 7. Comparison of the latitudinal source distribution as in Fig. 6.

4. Global properties of the Galactic source population

While in the previous section it is shown that most models can describe observations reasonably well within the sensitivity range of the HGPS, in the following section these models are used to predict global properties of the Galactic VHE source population, namely the total number of sources in the Milky Way, their contribution to the observed γ-ray flux, and their cumulative luminosity.

4.1. Total number of sources

As described in Sect. 3, we can derive the average number of detectable sources ⟨N_det⟩ according to Eq. (8) for any given total number of sources in a population N_tot: ⟨N_det(N_tot)⟩. With the probability of detecting 64 HGPS sources given by the Poissonian P_⟨Ndet⟩(64), the number of Galactic sources was derived from the maximum of the distribution f(N_tot) = P_{⟨Ndet(Ntot)⟩}(64) and cited errors from the corresponding 68% containment area around this maximum. These numbers vary considerably among the probed models, ranging from 831 sources (mSp4) up to 7038 sources (mFE) (see Table 6).

Although sources are treated generically as VHE γ-ray emitters in this model, that is no source type is explicitly assumed, a source count as high as seen for model mFE is challenging for the paradigm that SNRs and PWNe are the dominant source classes of VHE γ-ray emission. With a Galactic supernova rate of one per 40 yr (Tammann et al. 1994) a source count of 7000 implies a maximum age of emitters of ∼3 × 10⁵ yr. Interestingly, the models mSNR, mSp4, and mSp2B yield very similar results regarding the total number of sources. The similarity between model mSNR and mSp4 is also apparent in the cumulative source distribution over flux (log(N)−log(F)) as shown in Fig. 8, while mSp2B shows distinct differences. In Fig. 8 the distribution of observed sources in the HGPS is given by grey points with Poissonian errors. The mean distribution for the whole Galactic population⁵ according to the different models is given by coloured lines. It appears that model mSp2B on average does not comply with the observed distribution, especially for F >  0.03 c.u. (c.u.: integral flux of the Crab Nebula above 1 TeV) yielding too few sources in this regime. In contrast, the other four models allow for sources of high flux (e.g. > 0.1 c.u.) being undetected by the HGPS. The latter point is more clearly shown in Fig. 9. In this figure the predicted completeness range of the HGPS is shown; more precisely, the figure shows the median fraction of detected sources within the FOV for sources exceeding a given flux level. This number decreases with decreasing flux levels as fainter sources are less likely to be detected. For most models a kink in the distribution is observed between 0.1 c.u. and 1 c.u. That is caused by the limited sensitivity of imaging atmospheric Cherenkov telescopes (IACTs) to extended sources. For flat spatial source distributions the likelihood for close-by and, therefore, bright and extended sources to be found in the Galactic population increases. For model mSp4 we find that 2 ± 1 and for model mPWN 3 ± 2 sources exceeding the threshold of 1° are expected to be found in the Galactic population. To probe this regime either different data analysis techniques or different observation techniques, for example with water Cherenkov telescopes that can make use of a large FOV, can be exploited. Indeed, two extended sources, Geminga (Abdo et al. 2007; Abeysekara et al. 2017) and 2HWC J0700+143 (Abeysekara et al. 2017), have been detected this way, which is in good agreement with both predictions.

Fig. 8. Cumulative source distribution over flux (log(N)−log(F)). The HGPS source distribution is given by grey points with horizontal error bars depicting the bin width and vertical error bars showing the Poissonian error. Coloured lines represent the distribution of the whole Galactic source population averaged over the different realisations of the respective model. As the whole population also includes sources outside the FoV or too extended to be detectable, most models also overshoot the HGPS data in the range of completeness for point-like sources. Only mSp2B is in clear conflict with HGPS data.

Fig. 9. Fraction of detectable sources in the FOV of the HGPS cumulatively over source flux. The shaded regions represent the interquartile range; the lines depict the median.

Taking this one step further, we used these models to predict the number of sources that are detectable with the next generation of IACTs, CTA. Aiming for a point-source sensitivity of 2 mCrab in the longitude range |l| < 60° and latitude range |b| < 2° (Cherenkov Telescope Array Consortium 2019), the predicted numbers of detectable sources lie in the range 295 (mSp4)−457 (mFE). Since most of these sources are expected to appear point-like to CTA this number does not suffer considerably from a degradation of the sensitivity with source extent. Especially, this estimation is not affected by the inaccuracy of the description of source radii inherent to our models. The derived number is valid for the boundaries chosen for the luminosity. The implications for probing a larger dynamical range as it might be possible with CTA are discussed in Sect. 5. With the HGPS providing 53 sources in the region to be observed by the CTA Galactic plane scan, the CTA sample would increase the current source sample substantially by a factor between 5 and 9 according to these models.

4.2. Flux of unresolved sources

Although the HGPS is expected to comprise only 1%−9% of all sources in the Milky Way, these sources can already account for a significant fraction of the measurable flux (total flux given in Table 6). In Fig. 10 we give an illustrative example of one realisation of a synthetic VHE source population for model mSp4 in a face-on view of the Galaxy. For all of the 831 sources the luminosity and radius are encoded in the colour and radius of the circles representing them (radius not to scale), while sources that can be detected with the HGPS sensitivity are additionally denoted by an orange circle. Corresponding sky maps of the fluxes for this realisation of the population are shown in Fig. 11. The middle panel of Fig. 11 shows the sample of detectable sources. In comparison with the HGPS sample shown in the top panel a similarity of the two samples is recognisable, although the HGPS clearly shows a larger fraction of extended sources. The discrepancy in the ratio of extended to point-like sources between observation and model prediction was already discussed in Sect. 3.1. The sky map of the same synthetic population when observed with infinite sensitivity is shown in the bottom panel of Fig. 11. There, a band of faint sources along the Galactic plane, which can contribute to the unresolved, large-scale VHE γ-ray flux, is clearly seen. Focussing only on the region scanned by the HGPS, the flux of all sources detected by the HGPS exceeds the prediction of model mSp2B by ∼24% as already indicated by the log(N)−log(F) distribution (see Fig. 8). The other four models predict that unresolved sources make up about 13%−32% of the total flux stemming from the source population within this region. The H.E.S.S. measurement of large-scale diffuse emission in the HGPS from regions that do not contain any significant γ-ray emission quotes a similar number of ∼28% of the total measured VHE emission in large-scale diffuse emission (H.E.S.S. Collaboration 2014). However, these numbers are not directly comparable since the sky regions they are derived from are not identical. While the model allows us to remove detectable sources easily, complex exclusion regions were applied for the H.E.S.S. measurement to exclude contribution from sources, most of which accumulate at small Galactic latitude values. Still, these numbers are suggestive that unresolved sources might very well be the dominant component of the measured diffuse emission.

Fig. 10. Example realisation of a synthetic population of VHE γ-ray sources for model mSp4. Source luminosities are given by the colour scale; source radii are proportional to circle radii. Detectable sources within the HGPS sensitivity are indicated with orange circles.

Fig. 11. Sky maps of the VHE source population in the HGPS region. Top: from the HGPS catalogue. This map has been constructed by assuming a spherical source shape for better comparability with the simulations. Middle: detectable sources from the synthetic population shown in Fig. 10, based on the mSp4 model. Bottom: all sources from this simulated realisation.

Additionally, the two prominent extended sources that are seen at ∼ − 50° longitude on the sky map at the bottom of Fig. 11, but not in the sky map in the centre, demonstrate the effect of the maximum extent detectable by H.E.S.S., which has been discussed with respect to the catalogue completeness (see Sect. 4.1).

4.3. Luminosity of the Galactic source population

Using the mean photon energy of 3.52 erg ph⁻¹ from Sect. 2.2, the total luminosity of the Galactic VHE source population is estimated to lie in the range (5.07 × 10³⁶ − 2.22 × 10³⁷) erg s⁻¹. Assuming that γ-ray sources are the dominant contribution to the overall VHE γ-ray luminosity of the Milky Way and that the diffuse emission originating from propagating cosmic rays adds only a small contribution, these values can be compared with the total luminosity at megaelectronvolt and gigaelectronvolt energies (∼3 × 10³⁸ erg s⁻¹ and ∼8 × 10³⁸ erg s⁻¹, respectively; Strong et al. 2010). The γ-ray luminosity of the Milky Way at VHE turns out to be one to two orders of magnitude lower than in those two lower energy bands. This demonstrates that the presented models are compatible within the available energy budget constraints and indicates a drop in luminosity between the HE and VHE ranges.

5. Conclusions

We present models of the VHE γ-ray source population of the Milky Way, based on different assumptions of the spatial source distributions. Power-law indices of luminosity and radius functions of the population are derived from a subsample of the HGPS source catalogue (namely, extended sources with known distances) and its sensitivity. We pay special attention to correction of the observation bias. The validation of this bias correction is done with simulated toy models and demonstrates very good reconstruction capabilities. Furthermore, the simulations demonstrate that relying on the detected set of sources with no bias correction gives more or less arbitrary results. In this context it has to be noted that the limitation to the range of completeness does not completely avoid this problem because the completeness relates to point-like sources and does not apply to sources of larger extension.

A comparison of the source models with HGPS observations demonstrates a reasonable agreement. Despite a lack of asymmetry in the latitude distribution seen in all models as compared to the observed distribution (see Fig. 7), simulations approximately reproduce the spatial HGPS source distribution as well as the distribution of source fluxes and extents. Only the model mSp2B is disfavoured owing to its distinctively different longitude distribution in the HGPS sensitivity range and the clear under-prediction of the total flux from VHE γ-ray sources. All models under-predict the fraction of extended sources in the detectable sample, which can be attributed to either effects in the construction of the HGPS catalog (e.g. source confusion) or shortcomings of the model (e.g. invalidness of underlying assumptions) or a combination of both.

Despite the rather limited statistics in the sample of HGPS sources, our data-driven approach, which minimises the degrees of freedom, allows for relatively good predictions to be made. The derived models can be used to study possibilities and limitations of VHE observations. Examples are expectations for future instruments, an assessment of the amount of yet unresolved sources in a large-scale diffuse emission measurement depending on the sensitivity threshold, and the study of observational challenges like source confusion. We derived the total number of VHE sources in the Milky Way and the total luminosity and flux of the Milky Way in VHE γ-rays. Disregarding the one disfavoured model, of the four remaining models, mPWN, mSNR, and mSp4 do not deviate substantially from one another. While model mFE also yields similar results regarding the distribution of source properties, the predicted number of sources within the Milky Way exceeds the predictions of the other models by a factor > 3. Thus, the predicted range of VHE sources in our Galaxy is 800 to 7000. The scatter in the total luminosity and the total flux lies with a factor ∼4 in the range of (1.6 − 6.3) × 10³⁶ ph s⁻¹ and (3 − 15) × 10⁻¹⁰ ph cm⁻² s⁻¹, respectively. A significant fraction, (13 − 32)%, of the γ-ray emission of sources within the HGPS region is attributed to yet unresolved sources and contributes to the measured diffuse emission. With a foreseen sensitivity of 2 mCrab in the central Galactic plane, CTA should be able to increase the known Galactic VHE γ-ray source sample by a factor between 5 and 9.

It should be noted that the derived properties depend not only on the chosen spatial distribution model but also on the luminosity range covered by the model. Especially, the low-luminosity limit can have a significant impact on the expected number of sources and the energy budget. The choice of this limit was made to encompass the HGPS data, and it is possible that weaker sources (or a subdominant class of less luminous sources) can be on the verge of detection and not being identified with the present, inhomogeneous sensitivity. Expanding the dynamical range towards lower luminosities in the models increases the fraction of sources that lie outside the sensitivity limit and thus increases the expected number of VHE sources in the Galaxy and simultaneously increases (although to a lesser degree) predictions for fluxes and luminosities. Another factor that affects the estimation of the size of the population is source confusion. Given the predicted number of sources in the Galaxy, confusion of sources in the FOV and their ascribed fluxes seems to be inevitable. If source confusion already underlies the HGPS detected data set, the true size of the VHE source population can be expected to be even larger and luminosity and radius functions are subject to an additional bias.

Consequently, improvements of these models planned for the future include the correct treatment of source confusion in the comparison between simulations and data. Additionally, inclusion and proper treatment of sources with incomplete information (point-like sources/no known distance estimation) can increase the statistics of the sample. The incorporation of spectral information will allow for comparisons with neighbouring energy bands, making use of the wealth of information at Fermi-LAT energies and the complementarity of wide FOV HAWC observations.

In order to study the physics of VHE γ-ray sources, differentiation is required between the various source classes, such as SNRs and PWNe as the likely dominating contributors, and smaller contributions from other classes. This goes in conjunction with a widening of the parameter space that is needed to include physical source properties rather than the phenomenological observables used in this work. Although weakly constrained with the currently available data, this inclusion of physical source modelling for the characterisation of the population properties of VHE source classes is the logical next step after the in-depth studies of single objects and the description of the detected samples of sources so far performed with HGPS data (see e.g. Cristofari et al. 2013, 2017). The presented models can serve as a basis for such studies. With future surveys of more sensitive instruments such as CTA the amount of available data will increase, allowing for more accurate characterisations of the SNR and PWN populations.

Acknowledgments

We gratefully acknowledge fruitful discussions and feedback by Olaf Reimer, Ralf Kissmann, and Julia Thaler. For a very useful and constructive report which helped to improve the paper we would like to thank the anonymous referee. This work was conducted in the context of the CTA Consortium. This paper has gone through internal review by the CTA Consortium for which we like to thank Matthieu Renaud and Ullrich Schwanke.

References

All Tables

All Figures

	Fig. 1. Distribution of VHE sources in the Milky Way following the four-arm spiral model by Steiman-Cameron et al. (2010) as described in Sect. 2 (grey contours). Additionally, the HGPS detection horizon for point-like sources with a luminosity of 1033 ph s−1 is plotted on top (blue contours).
In the text

	Fig. 2. Estimated correction function for mSp4 showing the fraction of sources in the FOV of the HGPS that for a given luminosity and radius are detectable according to Eq. (9). The low values for small sources result from the rejection of point-like sources in the analysis. Additionally those sources of the HGPS that fulfil the selection criterion are shown.
In the text

	Fig. 3. Distribution of reconstructed αL (top) and αR (bottom) for 600 toy models using different spatial models and true values of αL, true = −2 and αR, true = −1. True values are denoted by the vertical line.
In the text

	Fig. 4. One-dimensional luminosity distribution (top) and radius distribution (bottom) for simulated populations. Original distributions are given in blue; distribution of detected sources in orange and the bias-corrected reconstruction of the original distribution are shown in green. Details are given in the text.
In the text

	Fig. 5. Kernel density estimation for the flux-extent distribution of the HGPS sample (grey) and mSp4 (orange). The contour lines indicate the 1σ, 2σ, and 3σ containment fraction. Flux is given in “crab units” (c.u.).
In the text

	Fig. 6. Comparison of the longitudinal source distribution between HGPS sources (grey bars) and predictions of different models (coloured lines). The dotted lines represent the modelled source densities within the HGPS region, arbitrarily scaled for the sake of visibility. In addition, the dynamic range of the HGPS sensitivity is shown on the top right panel. Details are given in the text.
In the text

	Fig. 7. Comparison of the latitudinal source distribution as in Fig. 6.
In the text

Fig. 8. Cumulative source distribution over flux (log(N)−log(F)). The HGPS source distribution is given by grey points with horizontal error bars depicting the bin width and vertical error bars showing the Poissonian error. Coloured lines represent the distribution of the whole Galactic source population averaged over the different realisations of the respective model. As the whole population also includes sources outside the FoV or too extended to be detectable, most models also overshoot the HGPS data in the range of completeness for point-like sources. Only mSp2B is in clear conflict with HGPS data.

In the text

	Fig. 9. Fraction of detectable sources in the FOV of the HGPS cumulatively over source flux. The shaded regions represent the interquartile range; the lines depict the median.
In the text

	Fig. 10. Example realisation of a synthetic population of VHE γ-ray sources for model mSp4. Source luminosities are given by the colour scale; source radii are proportional to circle radii. Detectable sources within the HGPS sensitivity are indicated with orange circles.
In the text

	Fig. 11. Sky maps of the VHE source population in the HGPS region. Top: from the HGPS catalogue. This map has been constructed by assuming a spherical source shape for better comparability with the simulations. Middle: detectable sources from the synthetic population shown in Fig. 10, based on the mSp4 model. Bottom: all sources from this simulated realisation.
In the text