EDP Sciences logo
Subscriber Authentication Point
Search
Menu
Home All issues Volume 712 (August 2026) A&A, 712 (2026) A56 Full HTML
Open Access
Issue
A&A
Volume 712, August 2026
Article Number A56
Number of page(s) 7
Section Interstellar and circumstellar matter
DOI https://doi.org/10.1051/0004-6361/202659784
Published online 03 August 2026

© The Authors 2026

Licence Creative CommonsOpen Access article, published by EDP Sciences, under the terms of the Creative Commons Attribution License (https://creativecommons.org/licenses/by/4.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

This article is published in open access under the Subscribe to Open model. This email address is being protected from spambots. You need JavaScript enabled to view it. to support open access publication.

1 Introduction

Pulsar wind nebulae (PWNe) are bubbles of relativistic pair plasma that originate from the confinement of the pulsar (PSR) wind by the parent supernova remnant (SNR) of its stellar progenitor (Gaensler & Slane 2006). While young systems are mostly observed as powerful nonthermal synchrotron sources, typically in radio and X-rays (Kargaltsev et al. 2015), older systems are mostly revealed in the TeV gamma-ray domain. This is where the relic low-energy leptons accumulated over the life of the system (Olmi 2023) produce radiation through inverse-Compton scattering on seed photons, most frequently on photons from the cosmic microwave background (CMB) or thermal infrared photons in the far- and near-infrared (FIR and NIR). Based on the estimated rate of core-collapse supernovae (Rozwadowska et al. 2021), we expect roughly 2000 potential gamma-ray emitting systems in the Galaxy (Fiori et al. 2022; De Sarkar et al. 2026), of which we currently detect and identify about 30 (Kargaltsev et al. 2013). Moreover, the population of ultrahigh-energy gamma-ray sources (with emission above 100 TeV up to PeVs) recently detected by the Large High Altitude Air Shower Observatory (LHAASO) approaches 100 to date, with the indication that many of them are directly connected with pulsars (Cao et al. 2024). These numbers are bound to increase substantially in the future with the operation of new TeV observatories such as the Cherenkov Telescope Array Observatory (CTAO Cherenkov Telescope Array Consortium 2019) or the ASTRI Mini-Array (Vercellone et al. 2022), as well as with the continuation of LHAASO operations. A correct modeling of the expected emission from these systems is then mandatory to assess the possibility of future detection and to properly characterize any observable target.

Current strategies for the modeling of old systems mostly rely on the so-called one-zone thin-shell approach (Ostriker & Gunn 1971; Reynolds & Chevalier 1984; Gelfand et al. 2009; Bucciantini et al. 2011; Martín et al. 2012). While this approach has proved quite robust for young systems, establishing its reliability for older systems has always been problematic, with different groups resorting to different recipes that often led to contrasting results (Gelfand et al. 2009; Bucciantini et al. 2011; Martín et al. 2012). Much of this stems from the fact that while for young systems the PWN and the SNR are causally disconnected and can be evolved separately, their evolution at later times becomes coupled, and their mutual feedback is highly non-linear. Even when the pulsar proper motion or gradients in the density of the ambient medium are neglected, there is no simple analytical way of parameterizing this interaction.

Simple 2D and even 3D hydrodynamical simulations of old PWNe have been presented in the past (Blondin et al. 2001; van der Swaluw et al. 2004; Kolb et al. 2017; Meyer et al. 2024; Meyer & Torres 2025), but no attempt to use them to derive their expected spectral properties has been made, mostly because in order to make these simulations numerically manageable, several strong simplifications have to be made (adiabaticity, no magnetic field, nonrelativistic fluids, and no pulsar spin-down), to the point that they can only be relied on at a qualitative level, but cannot be used to make robust quantitative predictions.

The standard one-zone thin-shell models that are used to describe the evolutions of PWNe inside SNRs are based on the assumption that during the initial free-expansion phase, all of the SNR ejecta swept up by the PWN accumulate in a thin shell at the very boundary of the nebula itself (Ostriker & Gunn 1971; Reynolds & Chevalier 1984; Gelfand et al. 2009; Bucciantini et al. 2011; Martín et al. 2012). In a series of papers (Bandiera et al. 2020, 2023a,b) it has recently been shown that in the following reverberation phase, this shell roughly maintains its integrity, at least as long as the dynamics is 1D. These studies also show that the shell does not accumulate matter during the reverberation phase and that it can be treated as a separate structure with respect to the PWN on the inside and the SNR shell on the outside. These findings inspired the development of a new strategy for modeling the dynamical and spectral evolution of PWNe inside SNRs: the PWN is treated using the standard one-zone approach, with or without radiation losses, as a piston that acts on a thin shell of swept-up ejecta, and the SNR is instead evolved using the equation of ideal hydrodynamics and is subject to inner boundary conditions dictated by the dynamics of the thin shell (this might be easily done adopting a Lagrangian formalism). The thin shell itself evolves according to the mass- and momentum-conservation laws under the action of the PWN on its inner and the SNR on its outer side. The set of equations forms a closed system that can be evolved in time (Bandiera et al. 2023b; Bucciantini & Olmi 2026; De Sarkar et al. 2026).

However, in a recent paper (Bucciantini & Olmi 2026, Paper I hereafter), this entire picture has been questioned in relation to the early free-expansion phase and possibly also to the later reverberation phase. It has been known, for a long time indeed, that the swept-up ejecta in the Crab nebula (and 3C 58) are not confined in a thin shell at the boundary of the relativistic PWN, but form a layer of thermally emitting filaments/knots that penetrate deeply inside the PWN (Trimble 1968; Hester et al. 1996; Fesen et al. 2008; Owen & Barlow 2015; Martin et al. 2021, 2025). Moreover, recent mass estimates seem to suggest that most of the swept-up mass is in these deeply embedded filaments and not in the outer skin (Owen & Barlow 2015). Driven by these and other considerations, we developed a new formalism that allows us to salvage the one-zone approach, with its versatility in terms of spectral evolution, and still model the development of a thick layer of mixed ejecta (Hester et al. 1996; Jun 1998; Bucciantini et al. 2004; Porth et al. 2014). Application of the new formalism to the Crab nebula suggested that when a layer of filaments was accounted for, the evolution of a PWN in the free-expansion could be quantitatively different from the evolution predicted using the standard thin-shell approach. That work, however, left the question unanswered of how relevant a thick layer, in place of a thin shell, might be for the following reverberation phase.

Unfortunately, while there is some evidence of mixing also in evolved system (Ma et al. 2016), foremost the large extension observed in a few TeV PWNe (H. E. S. S. Collaboration 2019a; Aharonian et al. 2006; H. E. S. S. Collaboration 2019b), there is really no system that can offer us some guidance on how to modify our theoretical approach, on par with the Crab nebula. Moreover, due to the large numerical cost of simulating old PWNe, even numerical results (e.g. Blondin et al. 2001; Kolb et al. 2017; Meyer & Torres 2025) do not provide a robust and reliable picture of their evolution. In the light of this uncertainty, we present here a modified version of the standard one-zone thinshell approach to the reverberation phase that accounts for the fact that when the PWN reaches the SNR reverse shock, there is no massive shell, but an extended layer of filaments.

The paper is organized as follows: in Sect. 2 the equations for the model are introduced, in Sect. 3 we present and discuss our results for a selected sample of models spanning the expected PWN-SNR population and compare them with the canonical thin-shell results. Finally, in Sect. 4 we summarize our conclusions.

Thumbnail: Fig. 1 Refer to the following caption and surrounding text. Fig. 1

Schematic evolution of a PWN inside an SNR in the presence of an extended layer of high-density filaments that penetrate the PWN itself. Top left panel: Initial structure in the free expansion phase. An SNR shell is present, bounded on the inside by a reverse shock (RS), while a high-density PWN with extended high-density fingers expandings outward. Top right panel: moment of beginning reverberation, at which the PWN also reaches the reverse shock. Bottom right panel: compression of the PWN by the SNR following reverberation. The reverse shock has disappeared, and matter in the filaments now accumulates in a thin shell at the boundary between PWN and SNR. Bottom left panel: compression of the PWN at maximum, and almost all matter in the filaments has now been accumulated in the thin shell, even though some filaments might still be present in the central region (not shown in figure).

2 The model

We model the PWN-SNR system with three fluids (see Fig. 1 for a schematic representation of the structure of this system during its evolution), with different levels of coupling among them.

  • The PWN is treated as a uniform bubble containing a relativistically hot massless plasma whose energy varies due to adiabatic losses (and radiation) and injection by the PSR. During the initial free-expansion phase, the PWN is not well coupled to the surrounding cold ejecta, which, instead of accumulating in a thin shell at the PWN boundary, finally forms a layer of filaments that penetrate it, and over which it exerts a reduced drag. When the swept-up mass forms fingers that penetrate inside and expand at a lower speed than the PWN edge (Martin et al. 2025), then the kinetic energy of this material is lower than in the canonical model, where it entirely accumulates in a thin shell at the edge. This means that the energy transferred by the PWN is also reduced, implying a lower coupling between the PWN and the cold SNR ejecta. During the later reverberation phase, on the other hand, the PWN is coupled with ejecta that have been shock-heated by the SNR reverse shock. The reason for this difference can be attributed to the different magnetization of the latter to the former. While it is expected that the cold ejecta that give rise to the thermal filaments are almost completely unmagnetized (and this holds for many of the high density thermal knots (Hester et al. 1996; Owen & Barlow 2015)), there is some evidence that the plasma downstream of the reverse shock is magnetized (Gotthelf et al. 2001; DeLaney et al. 2002; Rho et al. 2002), which will ensure a much better coupling with the nonthermal PWN. This means that a well-defined boundary will from between the PWN and the reverse-shock-heated ejecta, and in the late reverberation phase, this boundary can be taken as the radius of the PWN, Rpwn.

  • The SNR shell is treated as an ideal classical fluid, obeying the equations of 1D radial fluid dynamics. The matter in the SNR shell, especially the matter downstream of the reverse shock, is coupled to the PWN and the filaments. The level of coupling of the two is unclear, and it is possible that the highest-density knots, if produced by Rayleigh–Taylor instability from the interaction of the PWN and SNR, might freely penetrate the SNR in the form of shrapnel (Park et al. 2004; DeLaney et al. 2010; Miceli et al. 2013). For simplicity, we assume that having similar densities, the filaments do not efficiently penetrate the SNR shell and that they accumulate at its inner boundary Rin at most, which is also the edge of the PWN.

  • The filaments are treated as a cold pressure-less dense plasma, characterized by a high compressibility, which, while uncoupled to the PWN, interact efficiently with the SNR shell accumulating in a thin shell at a radius Rsh, coincident with the SNR inner boundary. We assume that the filaments penetrate all the way to the center. In reality, the fingers might only partially penetrate the PWN, as in the case of the Crab nebula, where they are seen to extend from 0.3–0.4Rpwn outward (Martin et al. 2021). However, as shown in Paper I, the presence of such a small central region, free of filaments, does not change the overall dynamics significantly. In this sense, our assumption amounts to maximizing the mixing, and it can be seen as the opposite limit with respect to the canonical thin-shell model, where the swept-up matter entirely accumulates at the edge.

The tree components are then all coupled together at Rpwn = Rin = Rsh. There, a shell will form. The equations describing the evolution of this shell are the same as are used in the standard thin-shell formalism. The mass in the shell grows by accumulating the matter from the filaments, ddtMsh(t)=4πRsh(t)2vrel(t)ρfil(Rsh,t)Θ[vrel(t)],Mathematical equation: $\[\frac{\mathrm{d}}{\mathrm{dt}} M_{\mathrm{sh}}(t)=4 \pi R_{\mathrm{sh}}(t)^2 v_{\mathrm{rel}}(t) \rho_{\mathrm{fil}}\left(R_{\mathrm{sh}}, t\right) \Theta\left[v_{\mathrm{rel}}(t)\right],\]$(1)

where vrel(t)=vfil(Rsh,t)−R˙sh(t)Mathematical equation: $\[v_{\text {rel}}(t)=v_{\text {fil}}(R_{\text {sh}}, t)-\dot{R}_{\text {sh}}(t)\]$ is the relative velocity between the shell and the filaments, whose velocity profile is given by vfil(r, t), and whose density profile is ρfil(r, t). On the other hand, the shell momentum evolves through momentum transfer from the accumulating filaments and under the action of the PWN pressure on the inside and the SNR pressure on the outside, ddt[Msh(t)R˙sh(t)]=4π[Ppwn(t)−Psnr(t)]Rsh(t)2+dMsh(t)dtvfil(Rsh,t).Mathematical equation: $\[\begin{aligned}\frac{\mathrm{d}}{\mathrm{dt}}\left[M_{\mathrm{sh}}(t) \dot{R}_{\mathrm{sh}}(t)\right]= & 4 \pi\left[P_{\mathrm{pwn}}(t)-P_{\mathrm{snr}}(t)\right] R_{\mathrm{sh}}(t)^2 \\& +\frac{\mathrm{d} M_{\mathrm{sh}}(t)}{\mathrm{dt}} v_{\mathrm{fil}}\left(R_{\mathrm{sh}}, t\right).\end{aligned}\]$(2)

When a fully developed homogeneously expanding filamentary layer of uniform density is assumed, as the one considered in Paper I, then we can introduce the auxiliary radius, R~(t)=Rpwn(tbeg,rev)+R˙pwn(tbeg,rev)(t−tbeg,rev),Mathematical equation: $\[\tilde{R}(t)=R_{\mathrm{pwn}}(t_{\mathrm{beg}, \mathrm{rev}})+\dot{R}_{\mathrm{pwn}}(t_{\mathrm{beg}, \mathrm{rev}})(t-t_{\mathrm{beg}, \mathrm{rev}}),\]$(3)

where tbeg,rev is the time at which reverberation begins, defined by the condition Rpwn(t) = Rrs(t), with Rrs(t) the reverse-shock radius. R~Mathematical equation: $\[\tilde{R}\]$ is the maximum radius that the filaments (equal to the radius of the PWN) would have had if the PWN had not collided with the reverse shock, but had instead moved outward at the velocity it had at the moment of collision with the reverse shock. As in Paper I, we assume that during the free-expansion phase (t < tbeg,rev), the outer edge of the filaments corresponds to the outer edge of the PWN and that it moves at the same speed. Then vfil(Rsh,t)=R˙pwn(tbeg,rev)Rsh(t)R~(t),Mathematical equation: $\[v_{\mathrm{fil}}(R_{\mathrm{sh}}, t)=\dot{R}_{\mathrm{pwn}}(t_{\mathrm{beg}, \mathrm{rev}}) \frac{R_{\mathrm{sh}}(t)}{\tilde{R}(t)},\]$(4)

and for a standard core-envelope profile for the SNR ejecta (again, as in Paper I), Eq. (1) can be solved, and we find Msh(t)=M~ej(Rpwn(tbeg,rev)vctbeg,rev)3−δ[1−(RshR~(t))3],Mathematical equation: $\[M_{\mathrm{sh}}(t)=\tilde{M}_{\mathrm{ej}}\left(\frac{R_{\mathrm{pwn}}\left(t_{\mathrm{beg}, \mathrm{rev}}\right)}{v_{\mathrm{c}} t_{\mathrm{beg}, \mathrm{rev}}}\right)^{3-\delta}\left[1-\left(\frac{R_{\mathrm{sh}}}{\tilde{R}(t)}\right)^3\right],\]$(5)

where vc is the velocity at the boundary of the core region of the ejecta, enclosing a mass M~ejMathematical equation: $\[\tilde{M}_{\mathrm{ej}}\]$. These equations can be generalized to more complex distributions of the filaments, including the case of partial penetration.

To close this system, an equation for the PWN pressure is required. In the nonradiative case, it can be computed according to 4πddt[Rsh3(t)Ppwn(t)]=Lpsr(t)−4πPpwn(t)Rsh(t)2R˙sh(t).Mathematical equation: $\[4 \pi \frac{\mathrm{~d}}{\mathrm{dt}}\left[R_{\mathrm{sh}}^3(t) P_{\mathrm{pwn}}(t)\right]=L_{\mathrm{psr}}(t)-4 \pi P_{\mathrm{pwn}}(t) R_{\mathrm{sh}}(t)^2 \dot{R}_{\mathrm{sh}}(t).\]$(6)

The SNR pressure just outside the PWN, Psnr(t) is derived by solving the full 1D evolution of the SNR, for which we employ the Lagrangian scheme described in Bandiera et al. (2023a).

3 Comparison with the canonical thin-shell approach

3.1 Population synthesis

In order to make a meaningful comparison between the canonical approach to the reverberation and our new strategy, we first computed a representative population of PWN-SNR systems. This was done following recipes similar to those present in the literature (Fiori et al. 2022; Bandiera et al. 2023a; De Sarkar et al. 2026), but without accounting for their distribution within our Galaxy, their magnetization, or any other input quantity related to their emissivity because we were not interested in computing their observability and operated in the nonradiative regime. In this case, the only parameters of interest are the pulsar initial spin-down power L0, the spin-down time τsd, the SNR energy Esn, the ejecta mass Mej, and the ambient medium density ρism (which can be used to define a characteristic time Tch and radius Rch, Truelove & McKee 1999). The pulsar spin-down properties are given by its initial period and magnetic field, whose distributions were taken from Watters & Romani (2011). For the initial period, we set a lower limit of 10 ms, as in the original work by Watters & Romani (2011), while we set an upper limit of 500 ms, above which there are mostly magnetars (Mereghetti et al. 2015). For the magnetic field, the lower limit was set at 1011 G, compatible with the minimum values of observed PSRs within PWNe (lower values are more characteristic of thermal neutron stars that show no evidence of pair-creation activity, Pavlov et al. 2004; Ho 2013), and an upper limit at 4.4 × 1013 G given that pair creation above the QED critical field is suppressed by photon splitting (Timokhin & Harding 2015, 2019). For the spin-down law, we assumed a standard dipole braking index, n = 3. Moreover, we adopted a standard PSR radius of 10 km (Özel & Freire 2016), a typical PSR mass Mpsr = 1.40 M⊙ (Biswas & Rosswog 2025), and a moment of inertia I = 1045 g cm2 (Breu & Rezzolla 2016). The SNR energy was taken from a uniform distribution in the range (0.2–2.0) × 1051 erg s−1 (Sukhbold et al. 2016). The ejecta mass was equal to Mfin − Mpsr, where in order to account for mass losses, the final mass of the progenitor star Mfin was a function of the zero-age main-sequence mass Mzams according to Sukhbold et al. (2016), and the latter was sampled from a Salpeter initial mass function (Salpeter 1955) in the range 9–25 M⊙ (Sukhbold et al. 2016). Given the large uncertainty in pre-supernova mass loss rates and more recent estimates that suggested reduced losses due to wind clumpiness (Renzo et al. 2017), we also considered the case of no mass loss at all (Mfin = Mzams) and found that the distribution of the final population was not altered in any significant way. The ejecta density followed a core-envelope profile, as in Paper I, with the inner core slope δ = 0 and an envelope slope ω = 12 (Matzner & McKee 1999). For the ISM density, we instead used the numerical results by Kim & Ostriker (2017) for the local Galactic medium. We assumed that one-third of the supernovae are, in the nomenclature of Kim & Ostriker (2017), runaway (i.e., they sample the ISM in an unbiased way), and the rest were classified as cluster (Chu & Gruendl 2008; Carretero-Castrillo et al. 2024). Again, we compared this choice with the simple analytical prescription by Elwood et al. (2019), based on SNR modeling, and found no significant difference in the distribution of the PWN-SNR population. Our synthetic population allowed us to define a region of interest in the characteristic plane L0τsd/Esn–τsd/Tch (Bandiera et al. 2023a) that contains 98% of the expected Galactic PWN-SNR systems, as shown in Fig. 2.

The upper left part of this region represents very energetic PSRs that also release all of their energy on timescales much shorter than the SNR characteristic time. They behave in an explosive fashion, and their energetics is strong enough to affect the global SNR dynamics. However, we expect no more than a few systems like these to have formed in the past 100 kyr in the Galaxy. This means that it is unlikely that we will be able to see any of them early enough to appreciate this explosive interaction. On the other hand, systems in the bottom right corner of the population have PSRs that not only inject little energy, but do so on very long timescales. For these, we expect the SNR evolution to be almost unaffected by the PWN.

3.2 Conditions at the beginning of the reverberation

In Fig. 2, we plot in this region of interest the difference in the time it takes the system to reach reverberation, tbeg,rev, and the difference in the pressure of the PWN at the beginning of reverberation between the canonical thin-shell case and the case of a fully developed homologously expanding filamentary layer (corresponding to X = 0.6 in the equations for the thick-shell of Paper I). This gives an overview of the differences already present at the beginning of the reverberation for systems characterized by the exact same parameters, due just to different recipes for the behavior of the swept-up ejecta.

In the presence of a thick filamentary layer, the PWN reaches the reverse shock earlier. In the case of very energetic systems, this can be about 50% earlier than in the thin-shell case, and it typically occurs when the reverse shock is still moving outward, while for low-energy systems, the difference is 10% at most and reverberation begins when the reverse shock is moving inward. The difference in tbeg,rev implies a difference in the size of the PWN at reverberation, Rbeg,rev. It is also evident that there is a clear trend for the time difference along the major axis of our elliptical region. However, this difference is not strongly pronounced because the dynamical trajectory of the reverse shock, whose radius reaches a maximum and then begins to recede, changes by less than 20% between 0.5 and 1.75 Tch. For the vast majority of systems, the maximum PWN size is unaffected by the presence or absence of a thick filamentary layer. Table 1 lists a comparison for equivalent models in the two regimes for the two sequences of Fig. 2, tracing systems along the major and minor axis of the region of interest. The size of the PWN at tbeg,rev does not change significantly in the region of the minor axis, while some stronger variations are seen along the major axis, with models at the beginning and end of the sequence reverberating with a smaller PWN radius because they encounter the reverse shock very early or very late. In Fig. 2 we also show the difference in the energy content of the PWN. Although in the presence of a thick layer, the PWN reaches reverberation earlier, its energy content can be almost twice higher for energetic systems. For low-energy systems, the trend is instead reversed. However, the differences are smaller than 15% for the vast majority of systems, and they are almost zero at the peak of the distribution (the center of the region).

Thumbnail: Fig. 2 Refer to the following caption and surrounding text. Fig. 2

Comparison of canonical thin-shell and thick-shell models (X = 0.6) for the PWN-SNR population. The black ellipses represent the region in the characteristic plane containing 98% of the expected Galactic PWN-SNR population. The points represent the models of Table 1: the red dot at the center is C0, yellow triangles represent models of the sequence S#, and magenta squares show those of L#. The upper panel represents the relative difference in the time at which reverberation begins Δtbeg,rev=tbeg,revthin−tbeg,revthickMathematical equation: $\[\Delta t_{\text {beg,rev}}=t_{\text {beg,rev}}^{\text {thin}}-t_{\text {beg,rev}}^{\text {thick}}\]$. The lower panel represents the relative difference in the energy content of the PWN at the time at which reverberation begins ΔEpwn=Epwnthick−EpwnthinMathematical equation: $\[\Delta E_{\mathrm{pwn}}=E_{\mathrm{pwn}}^{\text {thick}}-E_{\mathrm{pwn}}^{\text {thin}}\]$. The green line separates systems with a positive difference from systems with a negative difference.
Table 1

Models.

3.3 Evolution during reverberation

In Fig. 3, we show the evolution of the radius of the PWN for models C0, S1, S4, along the minor axis, which span a wide range of compressions from just 30% to about a factor 30, as shown in Table 2. A similar trend also holds along the major axis from model L4 to model L1, and similar differences between the canonical thin-shell and our thick-shell formalism are present. Model C0 shows that in the thick case, reverberation is not only anticipated, but the PWN tends to expand to a slightly larger radius, and this is reached at a later time than the maximum of the thin case. Interestingly, the evolution of the radius in the thick case appears to be a delayed and stretched version of the thin case. The first compression takes place at ~5Tch and is followed by a re-expansion, and then, a second compression sets in. The series of compressions and re-expansions is very similar to what is typically found for the equivalent thin-shell model, and the radii of the first minima are roughly the same. There is some indication, however, that the amplitude of these oscillations is reduced. The reason might be that unlike in the thin-shell case, where the mass of the shell is saturated at the beginning of the reverberation, the mass of the shell at the PWN outer surface in the thick-layer case grows in time, probably leading a more relaxed behavior. Interestingly, for a thick-layer case, it takes no more than one characteristic time for the mass accumulated at the outer PWN boundary to almost saturate, and this mass reaches values that are about twice what is found in the thin-shell case.

This overall qualitative trend also holds for less and more compressive systems. For model S1, where compression during reverberation is marginal, the dynamics of the PWN radius in the presence of a thick-layer does not show the even minor repeated compressions and re-expansions of the thin-shell case, and it seems to relax quite smoothly to the long-term trend. For the strongly compressive S4 model, the radial evolution follows the thin-shell case very closely, but in a somewhat delayed fashion. Interestingly, the mass accumulated at the PWN edge in all cases is roughly twice the mass of the thin-shell case. That the maximum and minimum radius of the PWN are not much affected by the way the swept-up mass behaves implies that the evolution of the PWN pressure (and hence, the magnetic energy) is also quite similar. We found indeed that the maximum pressure during compression changes by less than a factor of 2.

Thumbnail: Fig. 3 Refer to the following caption and surrounding text. Fig. 3

Evolution of the PWN radius (blue) and mass in the outer PWN shell (red) for the standard one-zone thin-shell formalism (dashed) and in the case of a thick homogeneous fully developed filamentary layer (solid). The yellow dots represent the moment when reverberation begins tbeg,rev. From top to bottom: models C0, S1, and S4. The dotted (green) curve describes the trajectory of the reverse shock for an SNR with δ = 0 and ω = 12. The characteristic mass is by definition the mass Mch of the SNR ejecta, while the characteristic radius Rch can be derived from the SNR energy, ejecta mass, and ISM density (Bandiera et al. 2021).
Table 2

Models and comparison of the properties of compression during reverberation.

4 Conclusions

Following the work presented in Paper I, where we introduced a new formalism for investigating the dynamics of PWNe inside SNRs in the presence of a thick filamentary layer that develops as the PWN expands inside the cold SNR ejecta, we have extended this analysis here to the following reverberation phase in a three-structure model, where the filaments accumulate in a shell at the boundary between the PWN and SNR. We implemented this formalism within a thin-shell Lagrangian code that followed the approach first introduced with the new TIDE-L code.

The main findings of our work are listed below:

  • In the presence of a thick filamentary layer of cold ejecta, the duration of the free-expansion phase is shorter than in the standard thin-shell formalism by up to about 50% for very energetic systems. This does not translate into major differences in the typical size of PWNe at reverberation, however, which is typically 0.65 Rch for a large portion of the population;

  • The qualitative behavior during reverberation is unaltered: after reaching the reverse shock, the PWN expands slightly farther up to a maximum radius that differs by ~10% at most from the thin-shell case; then, a compression sets in until a minimum radius is reached, which again does not differ substantially from the value in the thin-shell approximation. Afterward, the PWN undergoes a series of minor compressions and re-expansions;

  • Interestingly, the main difference with respect to the thinshell case appears to be in the time it takes from the beginning of the reverberation to reach the first minimum (duration of the compression phase). In the presence of a thick layer of filaments, the duration of the first major compression is lengthened and the minimum is typically reached ~ 50% later in time. This smoother behavior is mainly a consequence of the difference in the way the mass accumulates at the PWN boundary;

  • After reaching the first minimum, the PWN extension typically oscillates around values that are the same as in the thin-shell case. Most of the difference between the two cases is limited to the first compression;

  • There is no appreciable difference in the measured compression factors over a wide range of compressibilities. Because this is the main quantity that affects the spectral evolution of medium age and old PWNe, potentially driving them to super-efficiency (Torres et al. 2019) and fast cooling, we do not expect the presence of a thick layer of massive filaments to affect the long-term spectral evolution of PWNe strongly.

We recall that this formalism and the approach we introduced here, while a step forward in understanding and modeling the role played by a filamentary layer in the dynamics and evolution of PWNe, is still limited to a pure 1D radial evolution. In this sense, it cannot answer questions related to the multidimensional stability of the compression phase. Simulations and observations (Blondin et al. 2001; Kolb et al. 2017; Ma et al. 2016; Meyer & Torres 2025) clearly indicate that mixing of the PWN and SNR is likely to be important during reverberation. Complex physics like this can only handled using multidimensional numerical codes for fluid dynamics. However, from a numerical perspective, it is still not feasible from a numerical point of view to follow the entire evolution of a PWN to very late times (from a few years to several tens of kiloyears) in the appropriate regimes (combining a relativistic PWN with a nonrelativistic SNR) with the correct energetics (including the PSR spin-down) and with a resolution sufficient to resolve the formation of a filamentary layer, not to mention the inclusion of magnetic field. It is too demanding in terms of time and resources. It still remains to be assessed how effective and how strong the role of any possible instability would be. This also holds for its effect on the spectral properties of PWNe.

Acknowledgements

The authors thank D. F. Torres for fruitful comments and suggestion on this work. N.B. and B.O. 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

  1. Aharonian, F., Akhperjanian, A. G., Bazer-Bachi, A. R., et al. 2006, A&A, 456, 245 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  2. Bandiera, R., Bucciantini, N., Martín, J., Olmi, B., & Torres, D. F. 2020, MNRAS, 499, 2051 [NASA ADS] [CrossRef] [Google Scholar]
  3. Bandiera, R., Bucciantini, N., Martín, J., Olmi, B., & Torres, D. F. 2021, MNRAS, 508, 3194 [NASA ADS] [CrossRef] [Google Scholar]
  4. Bandiera, R., Bucciantini, N., Martín, J., Olmi, B., & Torres, D. F. 2023a, MNRAS, 520, 2451 [NASA ADS] [CrossRef] [Google Scholar]
  5. Bandiera, R., Bucciantini, N., Olmi, B., & Torres, D. F. 2023b, MNRAS, 525, 2839 [NASA ADS] [CrossRef] [Google Scholar]
  6. Biswas, B., & Rosswog, S. 2025, Phys. Rev. D, 112, 023045 [Google Scholar]
  7. Blondin, J. M., Chevalier, R. A., & Frierson, D. M. 2001, ApJ, 563, 806 [NASA ADS] [CrossRef] [Google Scholar]
  8. Breu, C., & Rezzolla, L. 2016, MNRAS, 459, 646 [NASA ADS] [CrossRef] [Google Scholar]
  9. Bucciantini, N., Amato, E., Bandiera, R., Blondin, J. M., & Del Zanna, L. 2004, A&A, 423, 253 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  10. Bucciantini, N., Arons, J., & Amato, E. 2011, MNRAS, 410, 381 [NASA ADS] [CrossRef] [Google Scholar]
  11. Bucciantini, N., & Olmi, B. 2026, A&A, submitted [arXiv:2606.18912] [Google Scholar]
  12. Cao, Z., Aharonian, F., An, Q., et al. 2024, ApJS, 271, 25 [NASA ADS] [CrossRef] [Google Scholar]
  13. Carretero-Castrillo, M., Ribó, M., Paredes, J. M., & Benaglia, P. 2024, Bull. Soc. Roy. Sci. Liege, 93, 103 [Google Scholar]
  14. Cherenkov Telescope Array Consortium (Acharya, B. S., et al.) 2019, Science with the Cherenkov Telescope Array [Google Scholar]
  15. Chu, Y.-H., & Gruendl, R. A. 2008, in Astronomical Society of the Pacific Conference Series, 387, Massive Star Formation: Observations Confront Theory, eds. H. Beuther, H. Linz, & T. Henning, 415 [Google Scholar]
  16. De Sarkar, A., Torres, D. F., Olmi, B., Bucciantini, N., & Meyer, D. M.-A. 2026, A&A, 711, A296 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  17. DeLaney, T., Koralesky, B., Rudnick, L., & Dickel, J. R. 2002, ApJ, 580, 914 [NASA ADS] [CrossRef] [Google Scholar]
  18. DeLaney, T., Rudnick, L., Stage, M. D., et al. 2010, ApJ, 725, 2038 [Google Scholar]
  19. Elwood, B. D., Murphy, J. W., & Díaz-Rodríguez, M. 2019, MNRAS, 483, 4551 [NASA ADS] [CrossRef] [Google Scholar]
  20. Fesen, R., Rudie, G., Hurford, A., & Soto, A. 2008, ApJS, 174, 379 [Google Scholar]
  21. Fiori, M., Olmi, B., Amato, E., et al. 2022, MNRAS, 511, 1439 [NASA ADS] [CrossRef] [Google Scholar]
  22. Gaensler, B. M., & Slane, P. O. 2006, ARA&A, 44, 17 [Google Scholar]
  23. Gelfand, J. D., Slane, P. O., & Zhang, W. 2009, ApJ, 703, 2051 [NASA ADS] [CrossRef] [Google Scholar]
  24. Gotthelf, E. V., Koralesky, B., Rudnick, L., et al. 2001, ApJ, 552, L39 [Google Scholar]
  25. H. E. S. S. Collaboration (Abdalla, H., et al.) 2019a, A&A, 621, A116 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  26. H. E. S. S. Collaboration (Abdalla, H., et al.) 2019b, A&A, 627, A100 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  27. Hester, J. J., Stone, J. M., Scowen, P. A., et al. 1996, ApJ, 456, 225 [Google Scholar]
  28. Ho, W. C. G. 2013, in IAU Symposium, 291, Neutron Stars and Pulsars: Challenges and Opportunities after 80 years, ed. J. van Leeuwen, 101 [Google Scholar]
  29. Jun, B.-I. 1998, ApJ, 499, 282 [Google Scholar]
  30. Kargaltsev, O., Rangelov, B., & Pavlov, G. 2013, in The Universe Evolution: Astrophysical and Nuclear Aspects, eds. I. Strakovsky, & L. Blokhintsev (Nova Science Publishers), 359 [Google Scholar]
  31. Kargaltsev, O., Cerutti, B., Lyubarsky, Y., & Striani, E. 2015, Space Sci. Rev., 191, 391 [NASA ADS] [CrossRef] [Google Scholar]
  32. Kim, C.-G., & Ostriker, E. C. 2017, ApJ, 846, 133 [NASA ADS] [CrossRef] [Google Scholar]
  33. Kolb, C., Blondin, J., Slane, P., & Temim, T. 2017, ApJ, 844, 1 [CrossRef] [Google Scholar]
  34. Ma, Y. K., Ng, C.-Y., Bucciantini, N., et al. 2016, ApJ, 820, 100 [Google Scholar]
  35. Martín, J., Torres, D. F., & Rea, N. 2012, MNRAS, 427, 415 [NASA ADS] [CrossRef] [Google Scholar]
  36. Martin, T., Milisavljevic, D., & Drissen, L. 2021, MNRAS, 502, 1864 [NASA ADS] [CrossRef] [Google Scholar]
  37. Martin, T., Milisavljevic, D., Temim, T., et al. 2025, MNRAS, 544, 1 [Google Scholar]
  38. Matzner, C. D., & McKee, C. F. 1999, ApJ, 510, 379 [Google Scholar]
  39. Mereghetti, S., Pons, J. A., & Melatos, A. 2015, Space Sci. Rev., 191, 315 [NASA ADS] [CrossRef] [Google Scholar]
  40. Meyer, D. M.-A., Meliani, Z., & Torres, D. F. 2024, A&A, 692, A207 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  41. Meyer, D. M.-A., & Torres, D. F. 2025, MNRAS, 537, 186 [Google Scholar]
  42. Miceli, M., Orlando, S., Reale, F., Bocchino, F., & Peres, G. 2013, MNRAS, 430, 2864 [Google Scholar]
  43. Olmi, B. 2023, arXiv e-prints [arXiv:2303.12019] [Google Scholar]
  44. Ostriker, J. P., & Gunn, J. E. 1971, ApJ, 164, L95 [NASA ADS] [CrossRef] [Google Scholar]
  45. Owen, P. J., & Barlow, M. J. 2015, ApJ, 801, 141 [NASA ADS] [CrossRef] [Google Scholar]
  46. Özel, F., & Freire, P. 2016, ARA&A, 54, 401 [Google Scholar]
  47. Park, S., Hughes, J. P., Slane, P. O., et al. 2004, ApJ, 602, L33 [Google Scholar]
  48. Pavlov, G. G., Sanwal, D., & Teter, M. A. 2004, in IAU Symposium, 218, Young Neutron Stars and Their Environments, eds. F. Camilo, & B. M. Gaensler, 239 [Google Scholar]
  49. Porth, O., Komissarov, S. S., & Keppens, R. 2014, MNRAS, 443, 547 [CrossRef] [Google Scholar]
  50. Renzo, M., Ott, C. D., Shore, S. N., & de Mink, S. E. 2017, A&A, 603, A118 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  51. Reynolds, S. P., & Chevalier, R. A. 1984, ApJ, 278, 630 [CrossRef] [Google Scholar]
  52. Rho, J., Dyer, K. K., Borkowski, K. J., & Reynolds, S. P. 2002, ApJ, 581, 1116 [Google Scholar]
  53. Rozwadowska, K., Vissani, F., & Cappellaro, E. 2021, New Astron., 83, 101498 [NASA ADS] [CrossRef] [Google Scholar]
  54. Salpeter, E. E. 1955, ApJ, 121, 161 [Google Scholar]
  55. Sukhbold, T., Ertl, T., Woosley, S. E., Brown, J. M., & Janka, H.-T. 2016, ApJ, 821, 38 [Google Scholar]
  56. Timokhin, A. N., & Harding, A. K. 2015, ApJ, 810, 144 [NASA ADS] [CrossRef] [Google Scholar]
  57. Timokhin, A. N., & Harding, A. K. 2019, ApJ, 871, 12 [NASA ADS] [CrossRef] [Google Scholar]
  58. Torres, D. F., Lin, T., & Coti Zelati, F. 2019, MNRAS, 486, 1019 [Google Scholar]
  59. Trimble, V. 1968, AJ, 73, 535 [Google Scholar]
  60. Truelove, J. K., & McKee, C. F. 1999, ApJS, 120, 299 [Google Scholar]
  61. van der Swaluw, E., Downes, T. P., & Keegan, R. 2004, A&A, 420, 937 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  62. Vercellone, S., Bigongiari, C., Burtovoi, A., et al. 2022, J. High Energy Astrophys., 35, 1 [NASA ADS] [CrossRef] [Google Scholar]
  63. Watters, K. P., & Romani, R. W. 2011, ApJ, 727, 123 [NASA ADS] [CrossRef] [Google Scholar]

All Tables

Table 1

Models.

In the text
Table 2

Models and comparison of the properties of compression during reverberation.

In the text

All Figures

Thumbnail: Fig. 1 Refer to the following caption and surrounding text. Fig. 1

Schematic evolution of a PWN inside an SNR in the presence of an extended layer of high-density filaments that penetrate the PWN itself. Top left panel: Initial structure in the free expansion phase. An SNR shell is present, bounded on the inside by a reverse shock (RS), while a high-density PWN with extended high-density fingers expandings outward. Top right panel: moment of beginning reverberation, at which the PWN also reaches the reverse shock. Bottom right panel: compression of the PWN by the SNR following reverberation. The reverse shock has disappeared, and matter in the filaments now accumulates in a thin shell at the boundary between PWN and SNR. Bottom left panel: compression of the PWN at maximum, and almost all matter in the filaments has now been accumulated in the thin shell, even though some filaments might still be present in the central region (not shown in figure).
In the text
Thumbnail: Fig. 2 Refer to the following caption and surrounding text. Fig. 2

Comparison of canonical thin-shell and thick-shell models (X = 0.6) for the PWN-SNR population. The black ellipses represent the region in the characteristic plane containing 98% of the expected Galactic PWN-SNR population. The points represent the models of Table 1: the red dot at the center is C0, yellow triangles represent models of the sequence S#, and magenta squares show those of L#. The upper panel represents the relative difference in the time at which reverberation begins Δtbeg,rev=tbeg,revthin−tbeg,revthickMathematical equation: $\[\Delta t_{\text {beg,rev}}=t_{\text {beg,rev}}^{\text {thin}}-t_{\text {beg,rev}}^{\text {thick}}\]$. The lower panel represents the relative difference in the energy content of the PWN at the time at which reverberation begins ΔEpwn=Epwnthick−EpwnthinMathematical equation: $\[\Delta E_{\mathrm{pwn}}=E_{\mathrm{pwn}}^{\text {thick}}-E_{\mathrm{pwn}}^{\text {thin}}\]$. The green line separates systems with a positive difference from systems with a negative difference.
In the text
Thumbnail: Fig. 3 Refer to the following caption and surrounding text. Fig. 3

Evolution of the PWN radius (blue) and mass in the outer PWN shell (red) for the standard one-zone thin-shell formalism (dashed) and in the case of a thick homogeneous fully developed filamentary layer (solid). The yellow dots represent the moment when reverberation begins tbeg,rev. From top to bottom: models C0, S1, and S4. The dotted (green) curve describes the trajectory of the reverse shock for an SNR with δ = 0 and ω = 12. The characteristic mass is by definition the mass Mch of the SNR ejecta, while the characteristic radius Rch can be derived from the SNR energy, ejecta mass, and ISM density (Bandiera et al. 2021).
In the text

Current usage metrics show cumulative count of Article Views (full-text article views including HTML views, PDF and ePub downloads, according to the available data) and Abstracts Views on Vision4Press platform.

Data correspond to usage on the plateform after 2015. The current usage metrics is available 48-96 hours after online publication and is updated daily on week days.

Initial download of the metrics may take a while.