© The Authors 2026

1. Introduction

Pulsars, the rapidly rotating and highly magnetized remnants of massive stellar explosions, lose their rotational energy predominantly through the emission of a relativistic wind of particles and magnetic field. When this wind interacts with the surrounding medium, it forms a pulsar wind nebula (PWN) – a bubble of relativistic plasma bounded by a termination shock where particles, primarily leptons, are accelerated to extreme energies. The PWNe radiate broadband non-thermal emission predominantly via synchrotron and inverse Compton (IC) processes, extending from radio to peta-electronvolt γ rays (see Gaensler & Slane (2006), Slane (2017), Olmi & Bucciantini (2023) for comprehensive reviews). Because they probe both the pulsar central engine and the interaction with the host supernova remnant (SNR) and interstellar medium (ISM), PWNe are unique laboratories for high-energy astrophysics, relativistic plasma physics, and particle acceleration under extreme conditions.

The PWNe constitute the most numerous class of firmly identified very high-energy (VHE, E > 100 GeV) Galactic γ-ray sources detected by current instruments such as High Energy Stereoscopic System (H.E.S.S.) (H. E. S. S. Collaboration 2018a,b), Major Atmospheric Gamma-ray Imaging Cherenkov (MAGIC) (Anderhub et al. 2010; Rico 2016), and Very Energetic Radiation Imaging Telescope Array System (VERITAS) (Konopelko 2008; Aliu et al. 2013, 2014; Mukherjee 2016). They are also key targets for wide-field detectors such as High-Altitude Water Cherenkov (HAWC) and the Large High Altitude Air Shower Observatory (LHAASO) (Albert et al. 2020; Cao et al. 2024). Understanding their global properties, demographics, evolution, and radiative output is therefore essential both for interpreting current surveys and for guiding the expectations of, for example, the upcoming Cherenkov Telescope Array Observatory (CTAO; see, e.g., (de Oña-Wilhelmi et al. 2013)) and Southern Wide-field Gamma-ray Observatory (SWGO; see, e.g., Conceição (2023)).

The evolution of PWNe is intimately tied to the dynamical stage of its host SNR, itself a function of the past stellar wind history of its progenitor and of the local conditions of the ISM. In the early free-expansion phase, the nebula expands within the cold, homologously expanding ejecta, and its size and energetics closely track the pulsar’s spin-down luminosity. However, the majority of Galactic PWNe are not young systems: they are middle-aged, mostly beyond ∼10⁴ yr. Early on in this stage, the SNR reverse shock travels back toward the explosion center and collides with the nebula. This event initiates the reverberation phase, in which the PWN is compressed, heated, and later re-expands within the surrounding ejecta (Reynolds & Chevalier 1984; Gelfand et al. 2009; Torres & Lin 2018; Bandiera et al. 2023b). Its effects are profound: the nebular magnetic field can be amplified by orders of magnitude, particle populations can be reenergized, and radiative losses are drastically enhanced. These processes strongly reshape the PWN broadband spectrum and temporal evolution.

Despite its importance, reverberation has remained the most challenging phase of PWN evolution to model. On the one hand, one-zone thin-shell models have been widely employed in population studies and spectral fittings, owing to their computational simplicity and speed (e.g. Pacini & Salvati 1973; Reynolds & Chevalier 1984; Becker & Huang 2007; Zhang et al. 2008; Qiao et al. 2009; Gelfand et al. 2009; Fang & Zhang 2010; Tanaka & Takahara 2010, 2011; Bucciantini et al. 2011; Martín et al. 2012; Vorster et al. 2013; Torres et al. 2013, 2014; Martín et al. 2016; Martin & Torres 2022; van Rensburg et al. 2018; Zhu et al. 2018; Torres & Lin 2018; Torres et al. 2019; Fiori et al. 2020; De Sarkar et al. 2022; Kolb et al. 2017; Temim et al. 2015). In the best of cases, these models approximate the swept-up shell as infinitely thin and impose analytic prescriptions for the confining SNR. While robust during the free-expansion phase − when the SNR ejecta are cold and self-similar − the reliability of these assumptions rapidly decreases during reverberation. At this stage, the interaction between the PWN and the SNR’s reverse shock becomes highly nonlinear, involving complex processes such as shell thickening, shock reflections, and dynamical feedback that lie beyond the reach of simple analytical prescriptions. On the other hand, multidimensional hydrodynamic (HD) and magnetohydrodynamic (MHD) simulations can accurately follow the complex interplay of shocks, turbulence, and magnetic fields (Blondin et al. 2001; van der Swaluw et al. 2001; Bucciantini et al. 2003; Komissarov & Lyubarsky 2004; Del Zanna et al. 2006; Porth et al. 2014; Olmi et al. 2016; Olmi & Torres 2020) and have recently included detailed descriptions of the circumstellar medium (CSM; Meyer & Meliani 2022; Meyer et al. 2024; Meyer & Torres 2025; Meyer et al. 2025). However, relativistic MHD simulations are prohibitively expensive for a long PWN evolution: they typically consume millions of CPU hours to model only a few hundred years and have constraints in volume too. As a result, they cannot be readily applied to the long-term evolution of hundreds or thousands of PWNe, as required for statistical studies or for comparison with survey data.

This modeling gap and large computational cost has motivated the development of the hybrid TIDE+L framework (Bandiera et al. 2020, 2021, 2023a,b). TIDE+L combines the efficiency of the thin-shell approach in the free-expansion phase with a Lagrangian treatment of the SNR structure during reverberation. This strategy allows the PWN–SNR system to be evolved self-consistently across all stages, including the complex compression episodes of reverberation, while also following the evolution of the particle spectrum subject to radiative and adiabatic losses. The code is fast enough to simulate thousands of systems, yet accurate enough to capture the essential physics of reverberation – at least when compared with 1D HD simulations. This means that TIDE+L represents a step forward for population synthesis studies of PWNe connecting pulsar birth properties to the observable γ-ray sky. The need for such an approach is particularly acute when addressing detection prospects. Because middle-aged PWNe dominate the Galactic population, and because their spectra and luminosities are heavily shaped by reverberation, predictions that neglect this phase risk achieving biased conclusions. Accurate forecasts are necessary to interpret present source catalogs, assess selection effects, and anticipate the yield of surveys with CTAO, SWGO, LHAASO, and other observatories.

In this work, we present our efforts to achieve population synthesis of Galactic middle-aged PWNe, including a more realistic treatment of reverberation. We constructed a synthetic population of PWNe sampled from distributions of the initial spin period, magnetic field, ejecta mass, and other SNR and ambient properties, and evolved each system with TIDE+L up to a maximum age of 10⁵ yr. After computing the radiative and dynamical evolution of each PWN, we compared the results with the flux sensitivity thresholds, while taking into account the sky coverage of current and upcoming facilities (H.E.S.S., HAWC, CTAO, LHAASO, and SWGO) as well as the sensitivity degradation effect due to source extension, thereby quantifying, in a statistically robust way, the fraction of the Galactic PWNe population that would be detectable. This paper is the first of a two-part series. Here, in Paper I, we examine the detectability of Galactic PWNe and establish the statistical baseline for their visibility in γ-ray surveys. Paper II will focus on the most extreme outcome of reverberation – superefficiency – where the nebular luminosity can temporarily exceed the pulsar’s spin-down power. Together, we aim for these papers to provide a unified framework for understanding the Galactic PWN population, its evolutionary diversity, and its observational signatures in the era of next-generation γ-ray astronomy.

This paper is organized as follows. In Sect. 2, we describe the underlying distributions that define the initial population and briefly outline the method used to evolve it. Sect. 3 presents and discusses the main results of the population synthesis, while Sect. 4 summarizes our conclusions.

2. The population and its evolution

2.1. Number of sources in the population

The synthetic population of Galactic PWNe was generated to reproduce the statistical properties of middle-aged systems, while remaining consistent with pulsar birth distributions and core-collapse supernova (SN) progenitors. In total, ∼2400 sources were simulated. Out of the total number of sources, we randomly considered 1600 sources, corresponding to the expected number of PWNe formed in the Galaxy, assuming a core-collapse SN rate of about 1.6 per century, reported as the best-fit value in Rozwadowska et al. (2021), implying the formation of about 1600 PWNe over the past 10⁵ yr. We neglected systems older than 10⁵ yr, which are expected to be too diluted to contribute significantly. This SN explosion rate is a factor of ∼3–4 larger than recent progenitor-based estimates of the core-collapse SN rate derived from near-complete censuses of OB stars in the solar neighborhood (Quintana et al. 2025). However, such lower values would be difficult to reconcile with the observed population of young PWNe, unless as-yet-unknown systematic biases in pulsar or PWN population modeling, such as tera-electronvolt detectability times or radiative efficiencies, significantly affect current estimations.

To account for the uncertainty resulting from individual population realizations, we performed 1000 random realizations of samples containing 1600 sources each, drawn from our full population. The ensemble of realizations provides the mean expected number of detectable PWNe, as well as the associated uncertainties arising purely from population stochasticity and sampling variance.

2.2. Spatial distribution

For the spatial distribution of the synthetic PWNe population, we used the same approach as is presented in Cristofari et al. (2017), whereby the core collapse SNRs are positioned in the Galaxy according to the spatial distribution of the Galactic pulsar population as modeled by Faucher-Giguère & Kaspi (2006). A similar approach was adopted by Fiori et al. (2022), as the method was optimized to reproduce the observed population of Galactic γ-ray SNRs. Distances to each of the PWNe were assigned to mimic the large-scale pulsar distribution (Fiori et al. 2022), ensuring a spread across longitude and latitude consistent with observed populations.

2.3. Age distribution

The age distribution of the population was randomly sampled from the SN birth rate, with a cut excluding systems younger than 100 years. This was assumed to maintain observational consistency, as no confirmed nearby SN explosion has been observed in the last 100 years, and any such nearby event would likely be detected unless severely obscured by dust.

2.4. SN explosion energy distribution

The explosion energy is typically considered to be fixed at E_SN = 10⁵¹ erg; however, in reality, the core collapse SN energetics show variability (Hamuy 2003; Nadyozhin 2003; Zampieri et al. 2003; Müller et al. 2017). To account for this variation, the distribution of SN explosion energy was calculated using a truncated Gaussian distribution. After imposing lower and upper truncation limits at 10⁵⁰ and 10⁵² erg, respectively, the considered truncated distribution has a mean ⟨log₁₀(E_SNR/erg)⟩ = 51 and a standard deviation of 0.54 (Batzofin et al. 2024, 2025; Martinez et al. 2022).

2.5. Interstellar medium characteristics

The ISM density n_ISM(=ρ_ISM/1.27 m_p) follows a uniform distribution in the range of 0.01 ≤ n_ISM ≤ 10  cm⁻³. Note that the factor of 1.27 accounts for a standard chemical composition including hydrogen, helium, and traces of heavier elements.

Our range of ISM densities is different to the values assumed by Cristofari et al. (2017), in which the authors consider a wider range (10⁻⁵ ≤ n_ISM ≤ 10  cm⁻³). A reasonable minimum cutoff was implemented in the ISM distribution according to expectations of typical Galactic environments, and to maintain the characteristic radius of the underlying SNR within a physically consistent range. As an example, for a choice of n_ISM = 10⁻⁵ cm⁻³, the SNR characteristic radius would be around 200-400 pc, depending on the value of the ejecta mass, which is inconsistent with the properties of currently observed Galactic SNRs. Such low-density environments may, in principle, host physically large, extremely radio-faint, diluted remnants that fall below the sensitivity of current surveys, thereby easing the tension regarding the total number of SNRs expected based on standard ccSN rates. Nevertheless, in this paper, we adopt a lower cutoff in n_ISM at 0.01 cm⁻³ to probe tenuous environments, while avoiding the production of overly diluted SNRs with unrealistically large characteristic radii, for which there is currently no clear observational support in Galactic SNR catalogs.

2.6. Estimation of ejecta mass from progenitors

To determine the ejecta mass, M_ej, we relied on the relation between the zero-age main-sequence (ZAMS) mass and ejecta mass for nonrotating core-collapse SN progenitors at solar metallicity, provided by tabulated models generated by the GENEVA stellar evolution code (Eggenberger et al. 2008; Ekström et al. 2012). GENEVA¹ simulates stellar structures in a one dimensional fashion, taking into account the microphysical processes of convection and diffusion in the stellar interior, using the Schwarzschild mixing length criterion as well as the Zahn viscosity diffusion coefficient for angular momentum transport. The physics of the mass-loss rate is described in Georgy (2010). This relation, illustrated in Fig. 1 of Batzofin et al. (2024), was used to map each progenitor ZAMS mass to its corresponding M_ej. To do this, the ZAMS masses were sampled from an initial mass function (IMF), which represents the distribution of stellar masses at birth. For high-mass stars, the IMF is well approximated by a power law of the form p(M)∝M^−α, where α = 2.35 corresponds to the Salpeter slope. We considered stellar masses in the range M ∈ [8,50] M_⊙. The higher end is unlikely to end in neutron stars (instead of black holes) but it is used here so that there is the possibility of finding a resulting maximum ejecta mass (from GENEVA tracks) at higher values. With a maximum ZAMS mass of 25 M_⊙, the ejecta results in around 12 M_⊙. Since a higher ejecta mass has been estimated from already-established PWNe, we have increased it to 50 M_⊙ to include slightly higher ejecta masses. Nevertheless, we recall that the number for the higher-mass ZAMS that we consider is much reduced due to the IMF.

Fig. 1. Spatial distribution of the synthetic PWN population in the Milky Way Galaxy for a random realization of 1600 sources. The corresponding ejecta masses (in units of M⊙) of each PWN are color-coded as indicated in the color bar. The right panel of the figure shows the Fermi-LAT γ-ray skymap and the distribution of the same PWNe population with respect to it. The background top-down (face-on) view of the Milky Way (Credit: NASA/JPL-Caltech/R. Hurt (SSC/Caltech)) and the all-sky γ-ray (Aitoff) skymap (Credit: NASA/DOE/Fermi LAT Collaboration) were generated using the mw-plot code (Hurt 2025).

We then used the inverse transform sampling technique (Devroye 1986) that converts uniform random numbers between 0 and 1 into ZAMS masses that are distributed according to the desired IMF. See Appendix A for details. We generated ZAMS values within the progenitor mass range, which shows a greater number of lower-mass stars due to the steep negative exponent of the IMF, as expected. Once the ZAMS mass was drawn, the ejecta mass, M_ej, was computed via interpolation from the GENEVA model outputs and randomly assigned to each of the synthetic PWNe. As can be seen from Fig. 2, the ejecta mass in this work ranges from 7 to 15 M_⊙ resulting from our adopted ZAMS mass range.

Fig. 2. Histograms of the parameter distributions employed in this work. Each histogram shows the mean number of sources per bin, obtained by averaging the bin counts over 1000 random realizations, each built by drawing 1600 PWNe without replacement from the full synthetic parent population of 2400 systems. The bin-by-bin uncertainty corresponds to the 1σ scatter among realizations, and typically ranges from ≈7% to ≈15% for the parameters shown. Note that the final magnetic field distribution, BPWN, shown in the figure was not assumed a priori, but instead emerges from the simulation. The histogram was constructed using the PWN magnetic field evaluated at the current age of each source.

Note that previous efforts (see, e.g., Abe et al. 2024; Fiori et al. 2022) have considered the ejecta mass to be distributed according to a Gaussian distribution centered at 13 M_⊙ with σ = 3 M_⊙, ranging from a minimum of 5 M_⊙ to 20 M_⊙. The latter high-end of the ejecta mass would in turn imply a very large ZAMS mass, with an unclear evolutionary path to a neutron star.

In reality, the distribution of ejecta mass is barely constrained. Different models (or the same model with different parameters) result in a huge variation of the final masses. For example, from Sukhbold et al. (2016) or Renzo et al. (2017), it can be seen that the variation in the final mass considering different wind models and/or different efficiencies is extremely wide. However, the ZAMS masses are better constrained, justifying our approach.

The impact of the distribution of the ejecta mass on the overall population, for example, in the characteristic luminosity–time (L₀/L_ch–τ₀/t_ch) plane, which will be discussed later, is modest. For example, a uniform ejecta mass distribution ranging from 5 M_⊙ to 20 M_⊙ would not change the population distribution on the L₀/L_ch–τ₀/t_ch plane drastically from what is reported in this paper. The spatial distribution of a single random realization of synthetic PWNe population of 1600 sources on the backdrop of the Milky Way Galaxy is shown in Fig. 1, with each point color-coded with the corresponding ejecta mass, M_ej.

2.7. Pulsar characteristics

Each simulated PWN is powered by a pulsar whose initial properties are drawn from observationally motivated distributions following Bandiera et al. (2023a). The initial spin period was sampled from a Gaussian distribution with mean ⟨P₀⟩ = 100 ms and standard deviation σ_P0 = 80 ms, truncated at a minimum value of 10 ms. This choice lies intermediate between the distributions proposed by Watters & Romani (2011), Johnston et al. (2020), and Faucher-Giguère & Kaspi (2006), thereby sampling both the γ-ray and radio pulsar populations without overemphasizing either extreme, and remaining within the parameter space where associated PWNe are expected to form. The pulsar initial magnetic field follows a lognormal distribution with mean ⟨log₁₀(B₀/G)⟩ = 12.3 and dispersion σ = 0.25, similar to Faucher-Giguère & Kaspi (2006) and Gullón et al. (2015), but with a cut at a minimum value of 10¹² G.

To calculate the initial period derivative, spin-down luminosity, and characteristic spin-down timescale, we employed the following relations (e.g., see Gaensler & Slane (2006)). The rotational evolution of a pulsar is primarily governed by the loss of its rotational kinetic energy due to electromagnetic torques. Assuming that the pulsar behaves as an isolated rotating magnetic dipole, the spin-down luminosity is given by

$$\begin{array}{c}\hfill {\dot{E}}_{\mathrm{rot}}=-I\mathrm{\Omega }\dot{\mathrm{\Omega }},\end{array}$$(1)

where Ω = 2π/P is the angular velocity, $\dot{\mathrm{\Omega }}$ is its time derivative, and I is the moment of inertia of the neutron star. The evolution of Ω follows a braking law of the form

$$\begin{array}{c}\hfill \dot{\mathrm{\Omega }}=-K{\mathrm{\Omega }}^n,\end{array}$$(2)

where n is the braking index and K is a constant that encapsulates the dependence on the magnetic field strength, radius, inclination angle, and moment of inertia of the neutron star. For a pure magnetic dipole radiation in vacuum, the braking index is n = 3. In our simulations, we adopted a constant braking index of n = 2.33, which corresponds to the secular braking index of Crab (Lyne et al. 2015; Horvath 2019), to account for observed deviations from the ideal dipole model. Rewriting the spin-down equation in terms of the period, P, and integrating under the assumption that K and n remain constant yields the period evolution, we obtained

$$\begin{array}{c}\hfill P(t)=P_0{(1+\frac{t}{{\tau }_0})}^{\frac{1}{n-1}},\end{array}$$(3)

where P₀ is the initial spin period. Differentiating the above equation gives the time-dependent period derivative

$$\begin{array}{c}\hfill \dot{P}(t)=\frac{P_0}{{\tau }_0(n-1)}{(1+\frac{t}{{\tau }_0})}^{\frac{2-n}{n-1}},\end{array}$$(4)

from which one can get the initial spin-down timescale, τ₀:

$$\begin{array}{c}\hfill {\tau }_0=\frac{P_0}{(n-1){\dot{P}}_0}.\end{array}$$(5)

The relationship between the initial period, period derivative, and magnetic field can be obtained using the following expression:

$$\begin{array}{c}\hfill B_0\approx \sqrt{\frac{3I{c}^n}{2{(2\pi )}^{(n-1)}{R}^{(n+3)}}}\times \sqrt{\dot{P_0}{P}_0^{(n-2)}}\mathrm{G},\\ \hfill \end{array}$$(6)

which is a general formula with the dimension of the magnetic field, valid for any arbitrary braking index, n. For example, for the magnetic dipole radiation in vacuum case (n = 3), with I = 10⁴⁵ g cm² and R = 10⁶ cm, Eq. (6) boils down to the well-known relation $B_0\approx 3.2\times {10}^{19}\sqrt{P_0\dot{P_0}}\mathrm{G}$.

The age of a pulsar with present-day period P(t) and period derivative Ṗ(t) can also be computed by integrating the spin-down law, leading to the expression

$$\begin{array}{c}\hfill t=\frac{P(t)}{(n-1)\dot{P}(t)}[1-{\left(\frac{P_0}{P(t)}\right)}^{n-1}].\end{array}$$(7)

The current age, period, and period derivative at the current age of the sources in the population are compatible with the formulation presented in this section.

The initial spin-down luminosity of a pulsar can be calculated from its initial spin period, P₀, and period derivative, Ṗ₀, using the relation

$$\begin{array}{c}\hfill L_0=4{\pi }^{2}I\frac{{\dot{P}}_0}{P_0^3}\mathrm{erg}{\mathrm{s}}^{-1}.\end{array}$$(8)

Once L₀ is determined, the time evolution of the spin-down luminosity can be described by

$$\begin{array}{cc}\hfill L(t)& =4{\pi }^{2}I\frac{\dot{P}(t)}{P{(t)}^3}\mathrm{erg}{\mathrm{s}}^{-1}\hfill \\ \hfill & =L_0{(1+\frac{t}{{\tau }_0})}^{-\frac{n+1}{n-1}}\mathrm{erg}{\mathrm{s}}^{-1},\hfill \end{array}$$(9)

where τ₀ is the initial spin-down timescale and n is the braking index. It is worth noting that the true value of the braking index likely reflects a more complex and poorly constrained spin-down history. Given the limited number of measured braking indices available, we adopt this representative value throughout the paper; this also follows Bandiera et al. (2023a,b).

The histograms highlighting the distributions of the aforementioned parameters are shown in Fig. 2. The sampling procedure employed in the parameter distribution histograms is discussed in Appendix B.

2.8. Assumed model parameters

Additional parameters relevant to the nebular emission were sampled from distributions motivated by previous studies (Fiori et al. 2022; Torres et al. 2014). These include the energy break in the injected particle spectrum (sampled from a log-normal distribution, with log₁₀(γ_b) having mean 5.74 and standard deviation 0.37, which aptly covers the plausible γ_b range from 5 × 10⁴ to 1 × 10⁷), low- and high-energy spectral indices (uniform distribution in 1.0 < α₁ < 1.7 and 2.0 < α₂ < 2.7, respectively), the containment factor (uniform distribution in 0.01 < R_L/R_TS < 0.9, where R_L and R_TS are the Larmor radius and termination shock radius, respectively), and the magnetic fraction (log-normal distribution with log₁₀(η) = − 1.80 and standard deviation 0.35, with the additional condition η < 1, to cover the plausible range from 10⁻³ to 0.2. This is also consistent with the typical magnetic fraction posited in Torres et al. (2014). For each system, the local interstellar radiation field (ISRF) was estimated by interpolating tabulated near- and far-infrared (NIR and FIR) energy densities according to the Galactic position (Porter et al. 2006). To account for the higher-energy densities inferred from spectral modeling of observed PWNe (Torres et al. 2014), multiplicative correction factors, randomly selected from an uniform distribution ranging from 2 to 10, were applied to both NIR and FIR components, with fixed temperatures of T_NIR = 3000 K and T_FIR = 70 K. Such enhancements are physically plausible, since PWNe are commonly located in complex environments within the Galactic plane, near star-forming regions, heated dust, and molecular clouds, where the radiation field can be significantly stronger than the large-scale, azimuthally averaged ISRF implemented in Porter et al. (2006). Altogether, this setup provides a statistically robust and physically motivated ensemble of PWNe, with initial conditions covering the plausible range of pulsar and its environmental properties. These synthetic systems form the basis for the subsequent dynamical and radiative evolution calculations, which are briefly discussed below.

2.9. Evolving the population

The evolution of each synthetic PWN was computed using the hybrid framework TIDE+L, which combines the time-dependent one-zone code TIDE (Martín et al. 2012; Torres et al. 2014; Martín et al. 2016; Martin & Torres 2022) with the Lagrangian reverberation modules introduced by Bandiera et al. (2020, 2021, 2023a,b). This approach addresses the limitations of thin-shell models, with analytical formalism of the host SNR structure assumed a priori, during the interaction between the PWN and SNR reverse shock, where nonlinear feedback dominates the dynamics and radiative outcome. During the early free-expansion phase, the nebula evolves under the thin-shell approximation. However, during the reverberation stage, when the SNR reverse shock approaches the nebular radius, earlier models have overestimated the degree of compression. In particular, the treatment of the swept-up shell and the host SNR structure becomes especially important during the reverberation phase, since the reverse-shock-driven compression of the PWN is transmitted through the surrounding SNR structure. Unlike earlier thin-shell approaches, which prescribe the surrounding SNR structure a priori, TIDE+L evolves the HD structure of the SNR self-consistently together with the PWN, including the dynamical feedback of the shell itself. This provides a more realistic treatment of the pressure transfer during reverberation and avoids the overly abrupt or artificially strong compression that may arise in more simplified thin-shell prescriptions (see, e.g., Figs. 12 and 13 of Bandiera et al. 2023a). Later prescriptions yielded smoother evolution but could not follow the shell thickness or internal shocks. The Lagrangian scheme, as implemented in Bandiera et al. (2023a,b), overcomes these issues by tracking discrete ejecta shells, capturing interface thickening, delayed compression, and multiple shock reflections. This enables realistic modeling of the PWN-SNR dynamics and their observational signatures.

The radiative evolution is also coupled to the dynamics through the time-dependent particle distribution, accounting for synchrotron, IC, bremsstrahlung, synchrotron-self-Compton (SSC), and adiabatic losses. The injected spectrum, powered by the pulsar’s spin-down luminosity, partitions the available energy between relativistic particles and the magnetic field according to a chosen magnetization parameter, from which the resulting non-thermal emission is subsequently computed.

3. Results

3.1. PWN evolutionary stages

In this work, we have divided the evolution of the PWN-SNR system into three phases:

1. The free expansion phase, in which the PWN is freely expanding within the unshocked SNR ejecta,

2. The reverberation compressing phase, in which the PWN is getting compressed by the reverse shock, but the minimum radius has not yet been reached, and

3. The reverberation post-compression phase, in which the first minimum has been reached, and the PWN is evolving further.

Figure 3 presents a pie chart illustrating the mean fraction of the simulated PWNe currently in each evolutionary stage. The sampling procedure employed to estimate the mean count of sources in each evolutionary stage and to produce the pie chart is discussed in Appendix B. Note that for PWNe in post-compression, the evolution of PWNe radii is complex, and to have them all in a single category is a simplification, as subsequent compression or expansion events can, in principle, happen. However, since these will likely be damped by instabilities or anisotropies in the medium, this global measure is enough for our aims.

Fig. 3. Pie chart distinguishing sources in free expansion (green), reverberation in the compressing phase (orange), and reverberation in the post-compression phase (red). The plotted counts correspond to the mean number of sources in each evolutionary stage, obtained by averaging over 1000 realizations of 1600 randomly drawn PWNe at the final simulation time. The 1σ uncertainties for each phase are also provided in the legend.

In order to simplify the characterization of the PWN-SNR evolution and the reverberation phase properties, and following Bandiera et al. (2023a,b), we introduced the SNR characteristic radius (R_ch), time (t_ch), and luminosity (L_ch). These characteristic scales are dependent on the values of M_ej, E_SN, and ρ_ISM, and can be written

$$\begin{array}{cc}\hfill R_{\mathrm{ch}}& =M_{\mathrm{ej}}^{1/3}{\rho }_{\mathit{ISM}}^{-1/3},\hfill \end{array}$$(10)

$$\begin{array}{cc}\hfill t_{\mathrm{ch}}& =E_{\mathrm{SN}}^{-1/2}{M}_{\mathrm{ej}}^{5/6}{\rho }_{\mathit{ISM}}^{-1/3},\hfill \end{array}$$(11)

$$\begin{array}{cc}\hfill L_{\mathrm{ch}}& =\frac{E_{\mathrm{SN}}}{t_{\mathrm{ch}}}\hfill \end{array}$$(12)

(see, e.g., Truelove & McKee 1999).

The combined evolution of the PWN–SNR system can be characterized by the ratios τ₀/t_ch and L₀/L_ch throughout its evolution (Bandiera et al. 2023a), so that each source can be represented as a point in the τ₀/t_ch–L₀/L_ch plane. Although this parametrization is strictly valid only in the non-radiative case and when the dimensionless ejecta structure is held fixed across all systems (Bandiera et al. 2023a,b), it remains a useful approximate framework more generally. In particular, it captures not only the pulsar birth properties, but also, through the adopted normalization, the influence of the SN explosion and the ambient medium, while reducing the number of explicit parameters needed to describe the evolution. Within this plane, the upper region is populated by energetic, long-lived pulsars, whose nebulae expand rapidly and undergo only mild reverberation, resembling Crab-like systems with comparatively large radii and weak compression. By contrast, the lower region corresponds to weak, rapidly fading pulsars, whose nebulae are more strongly affected by the SNR reverse shock, leading to stronger compression and more pronounced reverberation. Intermediate regions describe transitional cases, in which the PWN grows substantially before reverberation sets in, but still experiences a significant compression phase, resulting in moderate compression factors (CFs) and a delayed spectral evolution.

3.2. Population on the L₀/L_ch–τ₀/t_ch plane

Figure 4 shows a single random realization of 1600 synthetic PWNe on the log₁₀(L₀/L_ch)–log₁₀(τ₀/t_ch) plane. We further used the Gaussian kernel density estimation (KDE) to produce a smooth, continuous estimate of the density by using Gaussian (bell-shaped) kernels centered on each point. This let us identify regions where sources are more densely clustered on the plane, and thus visualize the spatial density distribution of the population. We also show ellipses encompassing the 10th, 50th, and 98th percentiles of the entire population, along with constant ratios of PWN-SNR energetics L₀τ₀/E_SN.

Fig. 4. The distribution of the population on the characteristic log10(L0/Lch)–log10(τ0/tch) plane. The left panel shows a single random realization of the synthetic PWN population containing 1600 sources on the plane, with ellipses marking the isolevels enclosing 10, 50, and 98 percent of the population. Each point is color-coded by density, derived from Gaussian KDE, with colors ranging from red (denser regions) to blue (sparser ones). The right panel displays the same plane along with the isolevels, but color-coded by the CF (CF = Rmax/Rmin) at the first compression, as indicated in the color bar. In this plot, only PWNe in the post-compression phase are shown. Colored triangles in this plot highlight two example PWNe within a 50% enclosing ellipse with different CFs used for individual simulations. Straight lines in both plots denote loci of four constant pulsar–SNR energetics.

We further investigated the dynamical compression experienced by PWNe by computing the CF, defined as the ratio between the PWN radius at the first local maximum and the subsequent local minimum in the evolution of the PWN radius, i.e., R_max/R_min always computed at the first compression event. For each simulated source, we identified these turning points in the radius-time profile and calculated the corresponding CF, capturing the degree of reverberation-induced contraction. By mapping the CF as a color gradient for all points over the characteristic pulsar initial spin-down luminosity and timescale ( log₁₀(L₀/L_ch) - log₁₀(τ₀/τ_ch) ) plane, we visualize how reverberation strength varies across different pulsar configurations. This enabled the identification of regimes where strong compression is more likely, and offered a framework for comparing model predictions with observed PWNe properties. The CF distribution on log₁₀(L₀/L_ch) − log₁₀(τ₀/τ_ch) plane is plotted in the right panel of Fig. 4 along with 10th, 50th, and 98th percentiles enclosing ellipses and constant ratios of PWN-SNR energetics.

3.3. Individual examples of PWN evolution

To show the time evolution of individual sources from our population, we selected two representative PWNe within the 50% enclosing ellipse, each characterized by a different CF, extracted from the distribution shown in Fig. 4 (their respective locations on the L₀/L_ch–τ₀/t_ch plane are indicated in the figure). The corresponding input parameters for these sources are listed in Table 1. In each case, we examined the multiwavelength (MWL) spectral energy distribution (SED) evolution over a series of snapshot ages to explore the effects of reverberation. Additionally, we tracked the evolution of the SNR forward shock (FS), reverse shock (RS), and PWN radius to illustrate the onset and impact of the reverberation phase on the PWN structure. As the system enters reverberation and undergoes compression, the magnetic field strengthens, leading to an increase in synchrotron luminosity and enhanced cooling. This results in a steepening of the SED, particularly in the optical to X-ray regime. These individual PWN results are shown in Fig. 5. To capture this behavior, we also present the time evolution of the photon index in the optical (4 × 10¹⁴ − 7.5 × 10¹⁴ Hz) and soft X-ray band (0.1–1 keV) for both sources, computed at the same snapshot ages as above.

Fig. 5. Results of the individual source simulations based on the CF are shown in this figure. The upper (lower) row corresponds to example 1 (example 2). The left panel in each row shows the evolution of the SED with time. The solid (dashed) lines correspond to the SED before (after) the reverberation, whereas the dotted line corresponds to the SED at the reverberation age (tbeg, rev). The shaded red (blue) region signifies the energy range in which the optical (X-ray) photon index is calculated. The middle panel shows the time evolution of the radius of the PWN, SNR forward and reverse shocks. In the right panel, we show the time evolution of the optical (X-ray) photon index in red (blue) for both sources, while also providing 10% error bars for all photon index values.

Following the onset of reverberation, the MWL SED of the sources exhibits complex, time-dependent evolution. In the high-compression case (Example 1; CF ≈ 50.3, upper panel of Fig. 5), at the onset of reverberation, the PWN radius evolution clearly shows a sharp turnover and pronounced contraction of the PWN radius, marking a brief but powerful compression phase before gradual recovery. This leads to a sharp increase in synchrotron luminosity, peaking in the optical to X-ray regime, resulting from an amplified magnetic field. This is followed by a rapid steepening of the spectrum at higher energies as the amplified magnetic field drives strong radiative cooling. The optical photon index initially hardens slightly at the onset of reverberation, as the particle population is rapidly reenergized during compression. This is followed by a brief but pronounced softening, caused by enhanced radiative losses in the amplified magnetic field during the early stages of reverberation. At later times, the optical spectrum hardens again, as continued particle injection and the cooling of higher-energy particles replenish the electron population radiating at optical energies. The X-ray photon index shows a similar initial mild hardening at the beginning of compression, but subsequently softens more gradually as the evolution proceeds and the highest-energy particles continue to lose energy in the strengthened magnetic field. In the high-CF case, the imprint of reverberation on the SED evolution is fast.

In contrast, the low-compression case (Example 2; CF ≈ 5.8, lower panel of Fig. 5) exhibits a milder and smoother evolution: the PWN compression due to reverberation is delayed, and further evolution is comparatively gradual. Consequently, the magnetic field amplification, and thus the variation in the synchrotron component and cooling at the higher energies, is also comparatively delayed and gradual across time. The overall SED evolution does not show a significant effect of reverberation long after the actual reverberation epoch (t_beg, rev), owing to the slower response of the PWN to the crushing of the reverse shock. The optical photon index shows a mild hardening shortly after the onset of reverberation, followed by a delayed softening, and then hardens again as higher-energy particles cool and progressively repopulate the lower-energy, optical-emitting range. By contrast, the X-ray photon index also exhibits an initial mild hardening, but then continues to soften gradually with time. This reflects the progressive depletion of the particles responsible for the X-ray emission, whose spectrum shifts to lower energies as the increasing magnetic field enhances synchrotron losses. Note that in this case, the softening of the X-ray photon index sets in later, but is more pronounced and long-lived than in the high-CF case, reaching higher values and evolving over a more extended time interval. Altogether, these examples illustrate that the magnitude of the CF primarily regulates the strength and duration of the luminosity amplification and spectral hardening or softening episodes during reverberation, highlighting how the diversity in CF across the population leads to markedly different observable signatures among middle-aged PWNe, which can be captured by X-ray instruments such as Chandra and XMM-Newton or detectors observing in the optical wavelengths.

3.4. Source extension as a function of distance

Figure 6 shows the distribution of the current age PWNe radii against the source distances for a single random realization of 1600 synthetic PWNe, along with the maximum and minimum extension found in the H.E.S.S. Galactic Plane Survey (HGPS). We also indicate the sources that would be detectable by CTAO and H.E.S.S. from our population in blue and red (purple if detectable by both), respectively. See Subsect. 3.5 for details on the applied flux thresholds and observable sky coverage for each instrument.

Fig. 6. Current PWN radius versus distance plot for a single random realization of the synthetic population containing 1600 sources. The red circles (if any) denote the sources that will be detectable by H.E.S.S., the blue circles represent the detectable sources by CTAO, and the purple circles denote sources that will be detectable by both. The gray circles signify the rest of the population. The dashed and solid lines represent the maximum and minimum angular extension, respectively, as posited by H. E. S. S. Collaboration (2018b) from HGPS-detected PWNe. In the plot, firmly identified PWNe in the HGPS (firebrick squares), candidate HGPS PWNe (gold circles), and PWNe outside the HGPS (orange diamonds) are also shown.

In Fig. 6, we also show the firmly detected and candidate PWNe by H.E.S.S., as well as PWNe outside HGPS. The extensions of the H.E.S.S.-detectable sources in our synthetic population remain broadly consistent with the range of extensions inferred from the HGPS data. Furthermore, our results indicate that CTAO is expected to detect a larger number of sources within the extension range (between 0.03 deg and 0.6 deg) posited by HGPS. Notably, most of the simulated sources do not exceed the maximum angular extension observed in the HGPS, while a significantly larger fraction are predicted to have radii smaller than the minimum range reported in the HGPS. This is most probably the consequence of the treatment adopted here, which neglects mixing and enforces spherical symmetry in the coupled PWN–SNR evolution. On the one hand, this suppresses internal SNR anisotropies and mixing (see, e.g., Ferreira & de Jager 2008); on the other, it limits distortions at the PWN contact discontinuity, thereby naturally favoring more compact and regular nebular shapes.

Additionally, we adopt a braking index of n = 2.33 rather than n = 3, implying a faster decline of the pulsar spin-down luminosity than in the pure dipole case. Furthermore, note that in Fiori et al. (2022) the evolution of each PWN was artificially modified at the typical escape time from the SNR (see the insets of Fig. 3 in Fiori et al. 2022), which tends to produce larger late-time nebulae. More generally, our present setup does not include additional mechanisms that could increase the apparent extension of the tera-electronvolt-emitting region, such as a pulsar proper motion, which may generate displaced or elongated structures through the coexistence of relic and freshly injected particles, or the formation of extended pulsar halos beyond the nebula proper.

In addition, the simulations are followed only up to a finite maximum age, which may further limit the development of the largest systems. We therefore regard the present treatment as a reasonable compromise between physical realism and computational feasibility, while noting that it may still retain a residual bias toward underestimating the largest observed extensions. Even without accounting for additional effects mentioned above that could further increase the sizes of PWNe, our simulations indicate that, given the angular resolution of CTAO, roughly half of the detectable sources should be spatially resolvable.

3.5. Flux predictions and detectability

In this section, we estimate the detectability of the simulated PWN population on a source-by-source basis, using the present-day gamma-ray output of each nebula together with the observational characteristics of each instrument. For every simulated source, the model provides a spectrum at its current evolutionary age, and this present-age spectrum was used to evaluate whether the source would be observable once the relevant sky coverage, flux sensitivity, and angular-extension effects were taken into account. The resampling procedure used to derive the predicted number of detections and the associated uncertainty is described in Appendix B. We discuss below the detectability criteria for different instruments considered in this paper.

CTAO: We assumed the alpha array configuration for CTAO. We estimated detectability using the present-day integrated photon flux above 125 GeV. Only the simulated sources lying inside the adopted CTAO sky coverage region were considered further. In the implementation used here, the latitude coverage was taken to be −6° ≤b ≤ 6° (Cherenkov Telescope Array Consortium 2019; Abe et al. 2024), while the longitude dependence of the survey sensitivity was treated in five segments; namely, [ − 60° ,60° ], [60° ,150° ], [150° , − 150° ], [ − 150° , − 120° ], and [ − 120° , − 60° ]². To each of these regions, we assigned different point-source sensitivities above 125 GeV, i.e., 1.8, 2.7, 3.8, 3.1, and 2.6 mCrab, respectively, according to Cherenkov Telescope Array Consortium (2019)³. Note that the effective exposure times will depend on the details of CTAO survey implementation. An example can be seen in Tables 6.3 and 6.4 of Cherenkov Telescope Array Consortium (2019). Thus, the CTAO sensitivity depends explicitly on the location of the source within the survey footprint.

A second key ingredient is the treatment of source extension. For each PWN, we derived an angular size from its final modeled radius and distance, using the geometrical relation θ = tan⁻¹(R_PWN/d), and we used this directly as the characteristic source size (σ_src) in the detectability calculation. Here, note that the adopted σ_src is a lower limit to the real size of the PWNe, which neglects effects such as mixing and the pulsar proper motion, which will lead to more extended PWNe as discussed in Subsect. 3.4. The point-source threshold in the appropriate longitude segment was then degraded according to source extension following the prescription given in H. E. S. S. Collaboration (2018b),

$$\begin{array}{c}\hfill F_{min}({\sigma }_{\mathrm{src}})=F_{min,0}\frac{\sqrt{{\sigma }_{\mathrm{src}}^2+{\sigma }_{\mathrm{PSF}}^2}}{{\sigma }_{\mathrm{PSF}}},\end{array}$$(14)

where σ_PSF = 0.07° is the adopted CTA point-spread-function scale and F_min, 0 is the point-source sensitivity for the relevant longitude interval. In addition, we imposed an upper angular-size cut and excluded sources with σ_src > 2° (Martin et al. 2022), since such very extended systems would fall outside the regime for which the adopted survey sensitivity can be meaningfully applied. Additionally, similar to Martin et al. (2022), we also considered an upgraded version of CTAO, dubbed CTAO+ in the paper, where the point-source sensitivity was improved by a factor of 2 compared to that presented in Cherenkov Telescope Array Consortium (2019).

H.E.S.S.: For H.E.S.S., we followed the same source-by-source approach, using the present-day integrated photon flux above 1 TeV for each simulated nebula. We restricted the source population to the adopted HGPS sky coverage; namely, |b|< 3° and l < 65° or l > 250° (H. E. S. S. Collaboration 2018b). The point-source sensitivity was taken to be 10 mCrab above 1 TeV (Martin et al. 2022), corresponding to 50 h of observation time. In this case, the mCrab conversion was obtained from the same adopted Crab spectrum (Eq. (13)), but now integrated from 1 TeV upward.

Similarly to CTAO, the source extension was explicitly taken into account. The angular size of each PWN was derived from its modeled final radius and distance, and the effective H.E.S.S. threshold was degraded according to Eq. (14), where here we used σ_PSF = 0.08° (Martin et al. 2022). We also imposed an upper angular-size cut of σ_src ≤ 0.7°, such that more extended systems were excluded from the detectable sample.

HAWC: For HAWC, we again applied the same source-by-source procedure using the present-day integrated photon flux above 1 TeV. To evaluate the sky visibility, the simulated Galactic coordinates were converted to equatorial declination, and only sources lying within the adopted HAWC field of view, −20° ≤δ ≤ 60°, were retained (Albert et al. 2020). The point-source sensitivity was taken to be 50 mCrab above 1 TeV (Abeysekara et al. 2017), corresponding to 5 years of observation time. As above, the mCrab conversion was obtained from the adopted Crab spectrum Eq. (13), here integrated from 1 TeV upward. Source extension was treated with the same prescription introduced above, with σ_PSF = 0.3° (Abeysekara et al. 2013) and an upper angular size cut at 2° (arbitrary) being adopted in this case.

LHAASO: For LHAASO, we treated WCDA and KM2A within the same general framework, but unlike CTAO, H.E.S.S., and HAWC, the detectability was not evaluated through a single (or a series of) mCrab thresholds with an explicit PSF-based extension correction. Instead, for each source we computed the flux directly from the model spectrum at a fixed reference energy; namely, F₃ at 3 TeV for WCDA and F₅₀ at 50 TeV for KM2A. The source positions were converted from Galactic coordinates to equatorial declination, and only objects within the adopted LHAASO field of view, −20° ≤δ ≤ 80° (Cao et al. 2024), were considered. The source angular size was again inferred from the modeled final radius and distance and expressed as r₃₉. The observation time for LHAASO was considered to be the same as in Cao et al. (2024), i.e., 508 days for WCDA and 933 days for KM2A.

For WCDA, detectability was determined by comparing F₃ with a size- and declination-dependent threshold surface reconstructed from the sensitivity curves given in Cao et al. (2024), using the E^−2.5 case. These curves provide the dependence of the threshold on declination for r₃₉ = 0° and 0.5°, and on r₃₉ for fixed declinations of −15° and 30°. We performed a two-dimensional interpolation in the (δ, r₃₉) plane using the tabulated sensitivity grid. A source was then counted as detectable by WCDA if its modeled F₃ exceeds the corresponding interpolated threshold. KM2A was treated analogously, but using the threshold surface for F₅₀ at 50 TeV, again based on the E^−2.5 sensitivity curves in declination and angular size. In this case, a source was counted as detectable when its modeled F₅₀ is larger than the corresponding interpolated KM2A threshold.

SWGO: For SWGO, we followed an approach analogous to that adopted for LHAASO; namely, evaluating the source flux at a fixed reference energy rather than an integrated flux above threshold. In this case, we computed F₁ at 1 TeV directly from the resulting model spectra. The simulated Galactic coordinates were converted to equatorial declination, and only sources within the adopted SWGO sky coverage, −70° ≤δ ≤ 20° (Scharrer et al. 2025; SWGO Collaboration 2025), were considered. The point-source sensitivity was taken to be 5 mCrab at 1 TeV (Angüner & Ergin 2024), corresponding to 5 years of observation time, where the mCrab conversion, in this case, was defined from the differential Crab spectrum at 1 TeV (see Eq. (13)). Source extension was then accounted for using the same prescription introduced above, adopting here σ_PSF = 0.1°, and we additionally imposed an upper size cut of σ_src ≤ 1° (Scharrer et al. 2025).

Using the detectability criteria for different instruments discussed above and after applying the resampling procedure across 1000 realizations, we present the detectability statistics resulting from our synthetic population in Table 2.

3.5.1. Reverberation impact on detectability

To further investigate how the detectability of PWNe depends on the treatment of reverberation, we compared the results obtained with the hybrid TIDE+L framework to those derived with TIDE (Martín et al. 2012, 2016; Martin & Torres 2022), which adopts a simplified, one-zone, purely thin-shell description. For this comparison, we applied the same instrument-dependent detectability criteria and the same sampling procedure discussed in Appendix B, so that the differences can be directly attributed to the underlying dynamical treatment. The resulting predictions are summarized in Table 2. We find that the detectability is systematically reduced in the TIDE case across all instruments. This behavior is consistent with the fact that the simplified treatment does not capture the reverberation phase with the same realism as TIDE+L, and may therefore overestimate the compression and the resulting magnetic field amplification, leading to excessively strong synchrotron losses, a reduced surviving high-energy particle population, resulting in a lower detectability fraction. The comparison thus reinforces the importance of modeling the reverberation phase appropriately when predicting the Galactic population of tera-electronvolt-emitting PWNe.

To illustrate this effect more directly, we compare the broadband SEDs of individual sources between the two approaches. Specifically, we first identify the PWNe that are predicted to be detectable by CTAO in a single random realization of TIDE+L simulated population results, and then compare their broadband SEDs with those obtained for the corresponding sources in the TIDE setup. This allows the impact of the reverberation treatment to be assessed at the level of the emitted spectrum itself, rather than only through the total number of detectable objects. These comparisons are shown in Fig. 7 together with the relevant instrumental sensitivity curves. As can be seen there, the same sources that are detectable by CTAO in the TIDE+L simulation generally lie closer to, or even below, the sensitivity curves when modeled with TIDE. This demonstrates that a more realistic treatment of reverberation has a substantial effect on the predicted tera-electronvolt emission, and consequently on the inferred size of the Galactic tera-electronvolt PWN population.

Fig. 7. Comparison of the SEDs of sources predicted to be detectable by CTAO for a single random realization of the simulated population taking the reverberation effects into consideration. For each source, the SED obtained with TIDE+L (left) is compared with that of the corresponding source evolved with the simplified TIDE framework (right). The relevant instrumental sensitivity curves are also shown in both panels.

3.5.2. Detectability as a function of PWN evolutionary stages

We also performed the sampling method discussed in Appendix B for each evolutionary phase, and for each realization we applied the instrument-specific flux thresholds and sky-coverage cuts as discussed above. Table 3 reports the mean detections averaged over all realizations, and the quoted uncertainties correspond to the sample 1σ standard deviation of the sampled ensemble. A more realistic treatment of reverberation increases the expected number of detections because it avoids the overly strong compression that may arise in simplified treatments. In such cases, the magnetic field can be overamplified during the compression phase, leading to excessively severe synchrotron losses that deplete the high-energy particle population and suppress the γ-ray emission. By evolving the reverberation phase more self-consistently, TIDE+L yields a milder and more realistic compression, so that synchrotron cooling is not artificially overestimated and more sources remain detectable. By contrast, sources in the post-compression stage are typically older, larger, and intrinsically dimmer, and therefore their detectability is generally lower. In addition, the post-compression evolution is likely more sensitive to effects such as anisotropy, which may become important at late times but are not included in the present model. For these reasons, only a small fraction of post-compression systems are predicted to be detectable.

3.6. Cumulative distribution of the population

To further assess the detectability of PWNe as a function of flux, we constructed cumulative distributions of sources, normalized to the Crab Nebula flux, above a certain energy threshold. For two chosen energy thresholds of > 100 GeV and > 1 TeV, we computed the expected integrated flux with respect to the Crab spectrum as given in Eq. (13). Integrating this expression yielded the corresponding flux of 1 crab unit (CU) in photons cm⁻² s⁻¹, with values of 1 CU (> 100 GeV) ≈ 6.5 × 10⁻¹⁰ photons cm⁻² s⁻¹ and 1 CU (> 1 TeV) ≈ 2.26 × 10⁻¹¹ photons cm⁻² s⁻¹, which were then used to normalize the integrated PWNe fluxes above a corresponding energy threshold.

The resulting normalized fluxes were sorted and used to construct the cumulative distribution curves. Here as well, instead of relying on a single realization, we employed the sampling method discussed in Appendix B to calculate the model cumulative distribution curves. As discussed in Appendix B, we used the mean cumulative distribution across all realizations of all the sources in the visible range of the respective instruments as the central model prediction in both cases, also including an error band computed as $\sqrt{N}$, where N is the mean cumulative count for each flux threshold. In both cases, we compared the model curves with HGPS results by overlaying cumulative distributions for different subsets: firmly identified PWNe, composite candidates, and unidentified (UNID) sources. This enabled a direct comparison between the simulated and observed source populations. For comparison, we selected all the sources from our population that are within H.E.S.S. and CTAO footprints, separately. The cumulative distribution plots are given in Fig. 8.

Fig. 8. Mean cumulative flux distribution out of 1000 realizations of synthetic population containing 1600 sources emitting above two specific energies, i.e., > 100 GeV (left panel), and > 1 TeV (right panel). The flux has been normalized with respect to 1 Crab unit. The lighter colored areas represent the errors of each curve. In both plots, the mean distribution of sources across the realizations of the synthetic population within the CTAO (in blue) and H.E.S.S. (in red) footprints are directly compared with the firmly identified PWNe (in orange), the same added with PWNe in a composite SNR (in green), and the sum of these two with the known HGPS unidentified sources (in black).

By comparing both > 100 Gev and > 1 TeV distributions with the observables, it can be seen that both the simulated curves corresponding to CTAO and H.E.S.S. are consistent, within errors, with observational PWNe and PWNe+composite curves until 10⁻² CU. Below that value, the instrumental sensitivity decreases, resulting in the flattening of the observational curves. The CTAO cumulative curve is consistent with the PWNe+composite+UNID curves at the highest flux (≳0.2 CU), indicating that these unidentified sources at the highest fluxes have a strong possibility of being PWNe. The simulated curves also deviate slightly from the PWNe+composite+UNID curve between around 0.02 CU to 0.1 CU in both cases (comparatively more pronounced in the > 1 TeV case), which may indicate that the SNRs, pulsar halos, or γ-ray binaries could contribute in that energy range. The relative contributions of the middle-aged PWNe to the observable tera-electronvolt γ-ray sky will be directly testable with future CTAO observations.

3.7. Tera-electronvolt luminosity distribution as a function of distance

Figure 9 shows the integrated 1-10 TeV γ-ray luminosity for one random realization of 1600 sources against distance (left panel) and spin-down power (right panel). In the luminosity-distance plot, the red curve represents the H.E.S.S. detection threshold, whereas the blue curve corresponds to that of the CTAO. As expected, a vast amount of sources are outside the H.E.S.S. detection threshold, whereas a much larger fraction will be detectable by CTAO in the future. We further infer that about 80% of the sources lying above the H.E.S.S. and CTAO threshold have ages < 40 kyr. Furthermore, about 95% of the sources above the H.E.S.S. threshold and 85% of those above the CTAO threshold are located within 15 kpc. This trend is consistent with the presently firmly identified PWNe, which can all be found within ∼15 kpc and have an age range within ≲40 kyr. We also find that only six sources above the H.E.S.S. threshold and nine above the CTAO threshold fall in the age range 80 < t_age < 100 kyr. This suggests that comparatively old PWNe are likely to lie below the CTAO detection threshold. In turn, this may indicate that, if older systems are detected by CTAO, they are more likely to correspond to halo-like objects than to classical PWNe.

Fig. 9. Distribution of integrated γ-ray luminosities in the 1-10 TeV band with source distance (left panel) and spin-down power (right panel) for a single random realization of synthetic population containing 1600 sources. In the left panel, the dashed red and blue curves signify the flux detection thresholds for H.E.S.S. (H. E. S. S. Collaboration 2018a) and CTAO (Abe et al. 2024), respectively. In the right panel, the solid line represents the slope corresponding to the L1−10 TeV − Ė relation, as obtained from H. E. S. S. Collaboration (2018b), whereas the dashed gray line represents that calculated from our population. The population distributions in all panels have been color-coded based on the current age of the corresponding sources. In all of these plots, firmly identified PWNe in HGPS (firebrick squares), candidate HGPS PWNe (gold circles), and PWNe outside HGPS (orange diamonds) are also shown.

Note that in previous literature, such as de Oña-Wilhelmi et al. (2013), it was posited that ∼300–600 PWNe would be detectable by CTA, especially in the CTA-I and CTA-D configurations. However, the authors considered three PWNe as representative of the entire TeV PWNe population, and did not account for the time evolution of the PWNe sizes and luminosities, which is very complex, especially when reverberation is considered. Hence, previous estimations naturally lead to an overestimation, and at best, can be considered as an upper limit. Our estimates, while still in the hundreds of PWNe, are more conservative as well as more solidly based on PWN physics.

3.8. Tera-electronvolt luminosity as a function of spin-down luminosity

A correlation between X-ray luminosity and spin-down power was suggested by Kargaltsev et al. (2013). In the 0.1–100 Gev Fermi band, a possible relation L_0.1−100 GeV with Ė can be inferred from Acero et al. (2013), Acharyya (2025), but data are too scarce to come to a meaningful conclusion. In the tera-electronvolt energy range, the dependence was found to be L₁₋₁₀ TeV ∝ Ė^0.59±0.21 (H. E. S. S. Collaboration 2018b), also with significant dispersion.

Here, the right panel of Fig. 9 shows an apparent correlation between our integrated tera-electronvolt luminosity and the spin-down power. When integrated γ-ray luminosity within 1-10 TeV is fit with a power law in logarithmic space, our population yields L₁₋₁₀ TeV ∝ Ė^1.15. The HGPS slope can be retrieved if we consider only the sources with spin-down power Ė > 5.9 × 10³⁵ erg s⁻¹, or sources with an age of fewer than 5.5 kyr. This indicates that instrumental sensitivity limitations can introduce biases in the inferred slope. The luminosity threshold required to recover this correlation reflects the absence of survey-specific biases (e.g., flux- or extension-limited detectability) in our synthetic sample, while the younger age range suggests that the model predicts a more rapid decline or spectral softening of the tera-electronvolt emission with time – likely driven by the evolution following reverberation – than what is inferred from the HGPS data. Note that we do not report any uncertainties on both the estimation of slopes from our population, or after applying selection cuts to our population, where the inferred slope matches that posited by HGPS. Given the large dispersion in the sample data points in which the slope is being fit in both cases, the uncertainties in the slopes do not bear much statistical meaning. Nevertheless, this exercise demonstrates that the HGPS results can be plausibly reproduced, provided that appropriate selection cuts are applied. For comparison, as before, we also show the firmly detected and candidate PWNe by H.E.S.S., as well as PWNe outside HGPS, in all plots of Fig. 9.

4. Concluding remarks

Given that PWNe are expected to dominate the tera-electronvolt γ-ray sky revealed by future observatories, accurately incorporating the physics governing their detectability is crucial. Our results are based on a more realistic treatment of reverberation compared to similar former approaches in the literature, implemented through the hybrid TIDE+L framework (Bandiera et al. 2023b). We find that neglecting reverberation dynamics leads to an underestimation of the γ-ray output, and thus to a lower fraction of detectable sources. Magnetic amplification, particle reheating, and subsequent losses induced by reverberation emerge as key ingredients shaping both luminosity distributions and SED diversity. Even though we caveat that reverberation is not a 1D process, and despite the improvement in our treatment, we are not addressing inhomogeneous compressions and the influence of the CSM in the evolution. For timescales of tens of thousands of years, these can only be addressed by 3D simulations such as the ones in Meyer et al. (2025) using MHD, or Olmi et al. (2016), Porth et al. (2014), and Kolb et al. (2017) provided future advances in computing power allow for larger and longer simulations in relativistic magnetohydrodynamics (RMHD), which is currently unrealistic.

This study is not intending to reproduce the P−Ṗ diagram, as the latter is mostly populated by radio sources, much older than the ones we are interested in here. In addition, we are using rotational losses for the pulsars, with a fixed braking index. Only a few values of the latter are known, and its range is vast. Moreover, it can fluctuate and secularly vary across the pulsar evolution. Although its influence on the evolution of the PWN is minor, it does significantly affect the pulsar period at sufficiently late times. At such ages, other effects, such as magneto-thermal cooling, may be relevant to determine the P−Ṗ diagram. These effects, never included in prior PWN studies, are also neglected in the present study as they are not expect to influence the population results, given the relatively young pulsar ages. We shall explore in more detail if there is any possible impact of magneto-thermal cooling on the evolution of PWNe in a future work.

Qualitatively, the predicted number of H.E.S.S. detectable PWNe in our population is consistent with HGPS statistics (H. E. S. S. Collaboration 2018b) if the majority of the HGPS UNID sources are actually associated with pulsar-powered objects, while CTAO is expected to increase detections by an order of magnitude. The updated CTAO+ configuration is expected to detect more than 200 PWNe in the Galaxy, in agreement with previous estimates. We also predict that future facilities such as SWGO will play an important role in expanding and confirming the Galactic tera-electronvolt PWN population. Our predicted LHAASO detectability is consistent with Cao et al. (2024), which reports 35 sources associated with pulsars. Of these, 24 are linked to pulsars younger than 10⁵ yr. Moreover, 22 out of the 35 sources are classified as ultrahigh-energy emitters, with emission reaching or exceeding 50 TeV, again in line with our predictions. The cumulative flux distributions suggest that many unidentified HGPS sources at the highest fluxes may correspond to evolved PWNe, while also leaving room for pulsar halos – targets that future CTAO and SWGO surveys will directly test.

The derived correlation between integrated 1–10 TeV luminosity and pulsar spin-down power (L₁₋₁₀ TeV ∝ Ė^1.15) indicates that the well-known L_γ–Ė relation arises intrinsically from the coupled pulsar–nebula evolution rather than from survey biases. The spin-down luminosity threshold (> 5.9 × 10³⁵ erg s⁻¹) needed to recover the HGPS survey trend underscores the absence of flux- or extension-limited biases in our synthetic sample, while younger systems (t_age < 5.5 kyr) reproduce the steeper relation consistent with observations.

We predict that most γ-ray-bright PWNe that are to be detected are likely younger than ∼40 kyr, while any older systems (80 < t_age < 100 kyr) detected are likely to be associated with halo-like morphologies. This work hints at a natural evolutionary link between middle-aged PWNe and the emergence of extended tera-electronvolt halos around mature pulsars. Individual case studies (Fig. 5) further demonstrate how compression governs spectral hardening or softening and transient brightening, explaining the observed diversity of Galactic PWN SEDs.

Overall, our study provides a statistically robust picture of the Galactic PWNe population and its detectability, while incorporating a realistic formalism for reverberation. This approach bridges the connection between pulsar birth properties and survey observables, offering theoretical guidance for CTAO, LHAASO, and SWGO. However, although comprehensive, our approach in this paper is not without caveats. The escape of accelerated particles from the nebula into the host SNR, and subsequently into the surrounding medium, can modify the radiative signatures and extension of the system (Martin et al. 2024, 2025). Moreover, a proper treatment of the CSM surrounding the PWN–SNR complex can significantly influence its dynamical evolution (Meyer et al. 2024, 2025). In addition, the total energy budget available from the pulsar spin-down includes a non-negligible pulsed component, which should be accounted for when evaluating the effective energy injection into the nebula. We also did not take into account the possible escape of the pulsars from their host SNRs. Another point that remains to be evaluated is the role of mixing, especially in the reverberation phase, and if there is a possibility of including it in one-zone models using an ad hoc recipe.

Acknowledgments

We thank the anonymous reviewer for helpful suggestions and constructive criticism, which vastly improved the quality of this work. This work has been supported by the grant PID2024-155316NB-I00 funded by MICIU /AEI /10.13039/501100011033 / FEDER, UE, and CSIC PIE 202350E189. This work is also partly supported by the Spanish program Unidad de Excelencia María de Maeztu CEX2020-001058-M, financed by MCIN/AEI/10.13039/501100011033, and by the MaX-CSIC Excellence Award MaX4-SOMMA-ICE, and also supported by MCIN with funding from European Union NextGeneration EU (PRTR-C17.I1). ADS is supported by the grant Juan de la Cierva JDC2023-052168-I, funded by MCIU/AEI/10.13039/501100011033 and by the ESF+. D.F.T. acknowledges the T. D. Lee Institute, where part of this research was done, for hospitality. BO and NB acknowledge partial support by the INAF Mini-grant HYPNOTIC87A: Hidden Young Pulsar Nebula Occupying The Inner Core of 87A and by the European Union – NextGenerationEU RRF M4C2 1.1 grant PRIN-MUR 2022TJW4EJ.

References

Appendix A: Sampling stellar masses from a power-law initial mass function

To generate stellar masses following a prescribed IMF, we use inverse transform sampling (Devroye 1986). This approach is well suited to power-law IMFs, since it allows direct sampling once the cumulative distribution function (CDF) is inverted.

We assume a power law probability density function of the form

$$\begin{array}{c}\hfill p(M)=k{M}^{-\alpha },\end{array}$$(A.1)

defined over the finite interval [M_min, M_max], where the normalization constant k satisfies

$$\begin{array}{c}\hfill {\displaystyle {\int }_{M_{\min }}^{M_{\max }}}p(M)dM=1.\end{array}$$(A.2)

For α ≠ 1, this gives

$$\begin{array}{c}\hfill k=\frac{1-\alpha }{M_{\max }^{1-\alpha }-M_{\min }^{1-\alpha }}.\end{array}$$(A.3)

The corresponding cumulative distribution function is

$$\begin{array}{c}\hfill F(M)=\frac{M^{1-\alpha }-M_{\min }^{1-\alpha }}{M_{\max }^{1-\alpha }-M_{\min }^{1-\alpha }}.\end{array}$$(A.4)

Setting F(M) = U, where U ∈ [0, 1] is a uniform random variable, and solving for M yields

$$\begin{array}{c}\hfill M(U)={[U(M_{\max }^{1-\alpha }-M_{\min }^{1-\alpha })+M_{\min }^{1-\alpha }]}^{1/(1-\alpha )}.\end{array}$$(A.5)

This expression is then used to draw ZAMS masses distributed according to the desired power law IMF. In Fig. A.1, we generated 1600 samples within the mass range [8, 50] M_⊙, and the resulting distribution is compared with the analytical IMF. As expected, the sample reproduces the power-law slope, with a greater number of lower-mass stars due to the steep negative exponent. Once the ZAMS mass is drawn, the ejecta mass M_ej is computed via interpolation from the GENEVA model outputs, which is also shown in Fig. A.1.

Fig. A.1. The left panel shows the analytical IMF and 1600 randomly sampled ZAMS masses between 8-50 M⊙ to ensure that we are sampling the IMF correctly. The right panel shows the relation between ZAMS mass and ejecta mass of nonrotating core-collapse SN progenitors at solar metallicity in the GENEVA library (Ekström et al. 2012).

Appendix B: Sampling procedure and statistical uncertainties

To quantify the statistical uncertainty arising from the finite size of the synthetic Galactic PWN population, we adopt a common sampling framework that is applied consistently to all quantities discussed in this work. The procedure is designed to capture the intrinsic variance expected when drawing finite realizations from a larger synthetic parent population, and it is used to estimate uncertainties on parameter distributions, evolutionary phase fractions, detectability statistics, and cumulative flux distributions.

All sampling procedures described below are based on a parent population of N_avail simulated PWNe (2400 in this work). From this parent set, we repeatedly draw mock Galactic realizations by randomly selecting a fixed number N_samp = 1600 sources without replacement, ensuring that no source appears more than once in a given realization. Each realization therefore represents an equally plausible finite Galactic PWN population. This sub-sampling is repeated N_real = 1000 times, yielding an ensemble of independent realizations drawn from the same underlying synthetic population.

B.1. Sampling of parameter distributions

For each parameter of interest, histogram bin edges are defined prior to the sampling, and the same bins are used for all realizations. In a given realization k, the sampled parameter values are binned to obtain a count c_k, j in bin j. Repeating this process for all realizations produces an ensemble of counts {c_k, j} that reflects fluctuations arising solely from finite-sample statistics.

The expected distribution is taken to be the mean count per bin,

$${\overline{c}}_j=\frac{1}{N_{\mathrm{real}}}{\displaystyle \sum _{k=1}^{N_{\mathrm{real}}}}{c}_{k,j},$$

while the associated statistical uncertainty is estimated from the sample standard deviation,

$${\sigma }_j=\sqrt{\frac{1}{N_{\mathrm{real}}-1}{\displaystyle \sum _{k=1}^{N_{\mathrm{real}}}}{(c_{k,j}-{\overline{c}}_j)}^2}.$$

In well-populated bins, the distribution of c_k, j approaches a Gaussian by virtue of repeated random sampling, making σ_j a meaningful 1σ estimate of the uncertainty. This procedure therefore yields a mean predicted distribution together with an intrinsic uncertainty arising from the finite size of the population.

B.2. Sampling of evolutionary phase fractions

Each simulated PWN is classified into one of three evolutionary phases: free expansion, reverberation–compressing, or reverberation–post–compression. The sampling procedure described above is then applied to estimate the expected fraction of Galactic PWNe in each phase.

For each realization k, the sampled population is reduced to a set of counts (c_k, 0, c_k, 1, c_k, 2), corresponding to the number of sources in each evolutionary phase. Across all realizations, this yields an ensemble {c_k, i} for each phase i. The mean expected number of PWNe in phase i is given by

$${\overline{c}}_i=\frac{1}{N_{\mathrm{real}}}{\displaystyle \sum _{k=1}^{N_{\mathrm{real}}}}{c}_{k,i},$$

with an associated uncertainty

$${\sigma }_i=\sqrt{\frac{1}{N_{\mathrm{real}}-1}{\displaystyle \sum _{k=1}^{N_{\mathrm{real}}}}{(c_{k,i}-{\overline{c}}_i)}^2}.$$

Because each count results from many independent random draws from a large population, the central limit theorem ensures that the distributions of c_k, i are approximately Gaussian, allowing σ_i to be interpreted as a 1σ uncertainty.

B.3. Sampling for detectability

For a given instrument considered, every simulated PWN is assigned an instrument-specific detectability flag based on its present-day model spectrum, sky position, and angular size. The sampling procedure is then applied to estimate the expected number of detectable sources and the associated sampling variance.

In each realization, we draw a population of size N_samp without replacement from the full eligible parent pool and evaluate the detectability criteria of the corresponding instrument for all sampled sources. Because the same realization is used to assess all instruments, correlated fluctuations between different observatories are naturally preserved. For a given instrument, the total number of detected sources in the kth realization is denoted by n_k, inst. From the ensemble of N_real realizations, we compute the mean predicted number of detections as

$${\overline{n}}_{\mathrm{inst}}=\frac{1}{N_{\mathrm{real}}}{\displaystyle \sum _{k=1}^{N_{\mathrm{real}}}}{n}_{k,\mathrm{inst}},$$

and the corresponding standard deviation as

$${\sigma }_{\mathrm{inst}}=\sqrt{\frac{1}{N_{\mathrm{real}}-1}{\displaystyle \sum _{k=1}^{N_{\mathrm{real}}}}{(n_{k,\mathrm{inst}}-{\overline{n}}_{\mathrm{inst}})}^2}.$$

These values are quoted as the expected detectability and its 1σ uncertainty. In practice, the distributions of n_k, inst are generally close to Gaussian, and their fit widths are consistent with the directly computed standard deviations, supporting the interpretation of σ_inst as the sampling uncertainty on the predicted number of detectable sources. In this paper, the resampling procedure is applied directly to the detectability flags, since the detectability analysis is carried out on a source-by-source basis using instrument-specific flux thresholds, sky-coverage constraints, and extension corrections.

B.4. Detectability across evolutionary phases

The detectability of PWNe as a function of evolutionary phase is estimated using the same sampling framework. In each realization, the sampled sources are first grouped according to their evolutionary stage. For each phase i and instrument “inst”, we compute the number of detectable systems $c_{i,\mathrm{inst}}^{(k)}$.

The mean number of detectable PWNe in phase i for a given instrument is

$${\overline{c}}_{i,\mathrm{inst}}=\frac{1}{N_{\mathrm{real}}}{\displaystyle \sum _{k=1}^{N_{\mathrm{real}}}}{c}_{i,\mathrm{inst}}^{(k)},$$

with uncertainty

$${\sigma }_{i,\mathrm{inst}}=\sqrt{\frac{1}{N_{\mathrm{real}}-1}{\displaystyle \sum _{k=1}^{N_{\mathrm{real}}}}{(c_{i,\mathrm{inst}}^{(k)}-{\overline{c}}_{i,\mathrm{inst}})}^2}.$$

As in the previous cases, the large number of realizations ensures that the resulting distributions are approximately Gaussian, and σ_i, inst provides a meaningful 1σ uncertainty.

B.5. Sampling of cumulative flux distributions

To compare model predictions with the observed HGPS cumulative source distributions, we construct ensembles of cumulative flux curves for CTAO and H.E.S.S. For each realization, instrument-specific sky coverage is applied to the sampled sources, and their integral fluxes above each flux threshold, for two different energy threshold cases (E > 1 TeV or E > 100 GeV) are computed and expressed in Crab Units.

For each realization k, sources are sorted by flux thresholds and used to construct a cumulative distribution N_k(> F). Repeating this procedure for all realizations yields an ensemble of cumulative curves. The expected cumulative distribution is taken as the mean,

$$\overline{N}(>F)=\frac{1}{N_{\mathrm{real}}}{\displaystyle \sum _{k=1}^{N_{\mathrm{real}}}}{N}_k(>F).$$

The associated uncertainty is approximated as $\sqrt{\overline{N}(>F)}$, reflecting the typical Poisson-like variation and enabling a direct comparison with observational results.

All Tables

All Figures

Fig. 1. Spatial distribution of the synthetic PWN population in the Milky Way Galaxy for a random realization of 1600 sources. The corresponding ejecta masses (in units of M⊙) of each PWN are color-coded as indicated in the color bar. The right panel of the figure shows the Fermi-LAT γ-ray skymap and the distribution of the same PWNe population with respect to it. The background top-down (face-on) view of the Milky Way (Credit: NASA/JPL-Caltech/R. Hurt (SSC/Caltech)) and the all-sky γ-ray (Aitoff) skymap (Credit: NASA/DOE/Fermi LAT Collaboration) were generated using the mw-plot code (Hurt 2025).

In the text

Fig. 2. Histograms of the parameter distributions employed in this work. Each histogram shows the mean number of sources per bin, obtained by averaging the bin counts over 1000 random realizations, each built by drawing 1600 PWNe without replacement from the full synthetic parent population of 2400 systems. The bin-by-bin uncertainty corresponds to the 1σ scatter among realizations, and typically ranges from ≈7% to ≈15% for the parameters shown. Note that the final magnetic field distribution, BPWN, shown in the figure was not assumed a priori, but instead emerges from the simulation. The histogram was constructed using the PWN magnetic field evaluated at the current age of each source.

In the text

Fig. 3. Pie chart distinguishing sources in free expansion (green), reverberation in the compressing phase (orange), and reverberation in the post-compression phase (red). The plotted counts correspond to the mean number of sources in each evolutionary stage, obtained by averaging over 1000 realizations of 1600 randomly drawn PWNe at the final simulation time. The 1σ uncertainties for each phase are also provided in the legend.

In the text

Fig. 4. The distribution of the population on the characteristic log10(L0/Lch)–log10(τ0/tch) plane. The left panel shows a single random realization of the synthetic PWN population containing 1600 sources on the plane, with ellipses marking the isolevels enclosing 10, 50, and 98 percent of the population. Each point is color-coded by density, derived from Gaussian KDE, with colors ranging from red (denser regions) to blue (sparser ones). The right panel displays the same plane along with the isolevels, but color-coded by the CF (CF = Rmax/Rmin) at the first compression, as indicated in the color bar. In this plot, only PWNe in the post-compression phase are shown. Colored triangles in this plot highlight two example PWNe within a 50% enclosing ellipse with different CFs used for individual simulations. Straight lines in both plots denote loci of four constant pulsar–SNR energetics.

In the text

Fig. 5. Results of the individual source simulations based on the CF are shown in this figure. The upper (lower) row corresponds to example 1 (example 2). The left panel in each row shows the evolution of the SED with time. The solid (dashed) lines correspond to the SED before (after) the reverberation, whereas the dotted line corresponds to the SED at the reverberation age (tbeg, rev). The shaded red (blue) region signifies the energy range in which the optical (X-ray) photon index is calculated. The middle panel shows the time evolution of the radius of the PWN, SNR forward and reverse shocks. In the right panel, we show the time evolution of the optical (X-ray) photon index in red (blue) for both sources, while also providing 10% error bars for all photon index values.

In the text

Fig. 6. Current PWN radius versus distance plot for a single random realization of the synthetic population containing 1600 sources. The red circles (if any) denote the sources that will be detectable by H.E.S.S., the blue circles represent the detectable sources by CTAO, and the purple circles denote sources that will be detectable by both. The gray circles signify the rest of the population. The dashed and solid lines represent the maximum and minimum angular extension, respectively, as posited by H. E. S. S. Collaboration (2018b) from HGPS-detected PWNe. In the plot, firmly identified PWNe in the HGPS (firebrick squares), candidate HGPS PWNe (gold circles), and PWNe outside the HGPS (orange diamonds) are also shown.

In the text

Fig. 7. Comparison of the SEDs of sources predicted to be detectable by CTAO for a single random realization of the simulated population taking the reverberation effects into consideration. For each source, the SED obtained with TIDE+L (left) is compared with that of the corresponding source evolved with the simplified TIDE framework (right). The relevant instrumental sensitivity curves are also shown in both panels.

In the text

Fig. 8. Mean cumulative flux distribution out of 1000 realizations of synthetic population containing 1600 sources emitting above two specific energies, i.e., > 100 GeV (left panel), and > 1 TeV (right panel). The flux has been normalized with respect to 1 Crab unit. The lighter colored areas represent the errors of each curve. In both plots, the mean distribution of sources across the realizations of the synthetic population within the CTAO (in blue) and H.E.S.S. (in red) footprints are directly compared with the firmly identified PWNe (in orange), the same added with PWNe in a composite SNR (in green), and the sum of these two with the known HGPS unidentified sources (in black).

In the text

Fig. 9. Distribution of integrated γ-ray luminosities in the 1-10 TeV band with source distance (left panel) and spin-down power (right panel) for a single random realization of synthetic population containing 1600 sources. In the left panel, the dashed red and blue curves signify the flux detection thresholds for H.E.S.S. (H. E. S. S. Collaboration 2018a) and CTAO (Abe et al. 2024), respectively. In the right panel, the solid line represents the slope corresponding to the L1−10 TeV − Ė relation, as obtained from H. E. S. S. Collaboration (2018b), whereas the dashed gray line represents that calculated from our population. The population distributions in all panels have been color-coded based on the current age of the corresponding sources. In all of these plots, firmly identified PWNe in HGPS (firebrick squares), candidate HGPS PWNe (gold circles), and PWNe outside HGPS (orange diamonds) are also shown.

In the text

	Fig. A.1. The left panel shows the analytical IMF and 1600 randomly sampled ZAMS masses between 8-50 M⊙ to ensure that we are sampling the IMF correctly. The right panel shows the relation between ZAMS mass and ejecta mass of nonrotating core-collapse SN progenitors at solar metallicity in the GENEVA library (Ekström et al. 2012).
In the text