2.1 PSR evolution

The engine driving the dynamics of the PWN is the central pulsar. It loses energy, injecting it in the form of a relativistic plasma, following a typical magnetic dipole spin-down rate (or luminosity, see Pacini & Salvati Reference Pacini and Salvati1973):

(1) \begin{equation} \dot{E} = L = 4 \pi^2 I \frac{\dot{P}}{P^3}\,,\end{equation}

where $I\simeq 10^{45}$ g cm $^2$ is the pulsar moment of inertia, while P is its rotational period and $\dot{P}$ its time derivative. For canonical pulsars, the period varies typically between $0.01$ and 1 s, while the period derivative is in the range $10^{-15}$ – $10^{-11}$ ss $^{-1}$ . As one might expect, X-ray bright synchrotron nebulae, likely younger, tend to be associated with energetic PSRs ( $\dot{E} \gtrsim 10^{36}$ erg s $^{-1}$ ). Older objects instead tend to be associated with less powerful engines (see e.g. Figure 8 in the following Section 5).

To model the evolution of the pulsar energy injection with time, it is customary to assume that its spin slows down from an initial value $P_0$ according to: $\dot{P} \propto P^{2-n}$ , where the coefficient n is called the braking index. This is usually assumed to be in the range $2\leq n\leq3$ , with $n=3$ the case of pure dipole spin-down. However, recently larger values have also been suggested (as in the case of H. E. S. S. J1640-465, reported in Archibald et al. Reference Archibald2016), even if their physical meaning is still not fully understood (Parthasarathy et al. 2020). If n remains constant during the pulsar life, the variation of the pulsar energy input is (Pacini & Salvati Reference Pacini and Salvati1973):

(2) \begin{equation} L(t) = L_0 \left( 1 + \frac{t}{\tau_0} \right)^{-(n+1)/(n-1)}\,,\end{equation}

where $\tau_0=P_0/[(n-1)\dot{P_0}]$ is the initial spin-down time of the pulsar. From this equation, it is easy to see that the pulsar input can be considered as constant only for $t\ll \tau_0$ , while properly accounting for its temporal evolution becomes important for the correct modelling of the nebula beyond $\tau_0$ .

2.2 SNR evolution

A very important point that we want to make before discussing the phenomenology of the various evolutionary phases, is that the age of the PWN (or the pulsar) is not a good indicator of the evolutionary stage of the PWN. A much better indication is instead given by the SNR characteristic time, first introduced by Truelove & McKee (Reference Truelove and McKee1999), providing the typical time scale for the evolution of the SNR:

(3) \begin{equation} t_{\rm{ch}} =E_{\rm{sn}}^{-1/2} M_{\rm{ej}}^{5/6} \rho_{\rm{ism}}^{-1/3}\,,\end{equation}

where $E_{\rm{sn}}$ is the SN explosion energy (typically of order of $10^{51}$ erg), $M_{\rm{ej}}$ is the mass of the ejecta (in the range $\sim$ 6–20 $\,{\rm M}_\odot$ for SNRs hosting a PWN, see e.g. Smartt Reference Smartt2009) and $\rho_{\rm{ism}}$ the density of the interstellar medium (ISM).

We wish to remark here that, contrary to the widely used Truelove & McKee (Reference Truelove and McKee1999) model, which predicts large variations, depending on the progenitor structure, in the evolution of the SNR reverse shock (RS, an important factor in setting the various PWN evolutionary stages), and in particular on the time it takes to reach the center of the remnant (an event that we name implosion, happening at time $t_{\mathrm{implo}}$ ), a much smoother trend has been recently found by Bandiera et al. (Reference Bandiera, Bucciantini, Martín, Olmi and Torres2021), using fully Lagrangian simulations. In particular, large variations of the implosion time only appear when drastically chancing the structure of the ejecta core (see e.g. Figure 4 of that work), while for the largely used case in PWN+SNR models (ejecta with a flat core plus a steep envelope), implosion happens at $t_{\mathrm{implo}}\simeq 2.4 \, t_{\mathrm{ch}}$ . This value clearly represents the maximum time the PWN can remain in free-expansion (see next subsection), and can be used as a sort of theoretical upper limit for the duration of this first phase see Section 2.4. In reality, for typical PSR energy injection, the duration of the free-expansion phase is no longer than about half of the implosion time. The implosion time is more representative of the time at which the compression due to the interaction of the PWN with the SNR reaches the maximum. It should be clear that, depending on the properties of the SNR, this time might be very different for diverse systems, and this is the reason why the age of the PWN (or the pulsar) is not a good indicator of the stage of its evolution. Indeed, a change by a factor of 2 in the mass of the ejecta leads to a 1.8 factor of difference in the implosion time: if we assume, for example, a fixed SN energy of $E_{\rm{sn}}=10^{51}$ erg, and an ISM with number density 1 particle cm $^{-3}$ , a PWN in a $6\,{\rm M}_\odot$ SNR could stay at most for $\sim$ 5,000 yr in the first stage, while a PWN powered by the same pulsar in a $10\,{\rm M}_\odot$ remnant could stay in free-expansion for $\sim$ 9,000 yr.

Bandiera et al. (Reference Bandiera, Bucciantini, Martín, Olmi and Torres2023) show that the entire evolution of each PWN-SNR system is fully determined—in the absence of significative radiative losses—by the two quantities: $[\tau_0/t_{\mathrm{ch}}\,,\; L_0t_{\mathrm{ch}}/E_{\mathrm{sn}}]$ . These in practice weight the pulsar time and energetics with respect to the SNR ones.

In the following of this section, we limit ourselves to a phenomenological discussion of the properties of the various phases.

2.3 Free expansion

In the first phase the PWN expands, with a mild acceleration, in the cold freely expanding ejecta of the SNR core, hence its name. The typical PWN expansion speed is of the order of few thousands km s $^{-1}$ , much higher than the average velocity the PSR can acquire during the SN event (the kick velocity, ranging in 100–500 km s $^{-1}$ , see e.g. Faucher-Giguere & Kaspi 2006), so that one can safely consider the PSR to be stationary at this stage. During this phase, the evolution of the PWN is independent from that of the SNR shell, since no interaction has been established yet between the two (Reynolds & Chevalier Reference Reynolds and Chevalier1984). Based on the results by Jun (1998), showing that the PWN collects material in a thin-shell at its boundary during its initial expansion ( $t\ll \tau_0$ ), a simplified description of the PWN evolution can be obtained using the following formulas:

(4) \begin{align}\frac{d}{dt}\left(4\pi\,P(t)R(t)^4\right) {}= {} L(t)\,R(t)\,, \end{align}

(5) \begin{align}\frac{d}{dt}\left(M(t)\frac{dR(t)}{dt}\right)\;\; {}= 4\pi\,P(t)R(t)^2+\frac{dM(t)}{dt}\frac{R(t)}{t},\end{align}

known as the ‘thin-shell approximation’. Here P(t) is the PWN pressure, M(t) the mass of the thin-shell and R(t) the shell radius, that within the thin-shell approximation, can be taken as that of the PWN.

The freely expanding ejecta have a density profile characterised by a flat core surrounded by a steep envelope. It is customary to use power laws to describe them: the steep envelope as $r^{-\omega}$ (with $\omega > 5$ ) and the shallow core as $r^{-\delta}$ (with $\delta<3$ , see Bandiera et al. Reference Bandiera, Bucciantini, Martín, Olmi and Torres2021 and references therein):

(6) \begin{equation}\rho_{\mathrm{ej}}(r,t)=\begin{cases}A\,(v_{\mathrm{t}}/r)^\delta/t^{3-\delta}, {}& \text{if } r < v_{\mathrm{t}} t\,, \\ \\[-8pt] A\,(v_{\mathrm{t}}/r)^\omega t^{\omega-3}, {} &\text{if } r \geq v_{\mathrm{t}} t \,,\end{cases}\end{equation}

with the parameters A and $v_{\mathrm{t}}$ (the expansion velocity of the ejecta core) that depend on the SN energy and mass of the ejecta as:

(7) \begin{align}A {}= {}\frac{(5-\delta)(\omega-5)}{2 \pi (\omega-\delta)}\frac{E_{\mathrm{sn}}}{v_{\mathrm{t}}^5},\end{align}

(8) \begin{align}v_{\mathrm{t}} {}= {} \sqrt{\frac{2(5-\delta)(\omega-5)}{(3-\delta)(\omega-3)}\frac{E_{\mathrm{sn}}}{M_{\mathrm{ej}}}}\,\,.\end{align}

For the simplified case of the flat density profile ( $\delta=0$ ) plus a steep envelope ( $\omega=\infty$ ), actually the most commonly used, the mass of the shell can be expressed as: $M(t)=4\pi R^3(t) A/(3t^3)$ . At early enough times, when the pulsar input can be still considered as constant (namely $t \ll \tau_0$ , so that $L(t)\simeq L_0$ ) an analytical solution for Equations (4)–(5) can be found (Bandiera et al. Reference Bandiera, Bucciantini, Martín, Olmi and Torres2023):

(9) \begin{equation} R_{\textrm{PWN}}(t)\Bigg\vert_{t\ll \tau_0}=\left[\frac{(3-\delta)(5-\delta)^3}{(9-2\delta)(11-2\delta)}\frac{L_0}{4\pi}\frac{t^{6-\delta}}{A v_t^\delta}\right]^{1/(5-\delta)}\,.\end{equation}

This reduces to the well-known trend found for the specific case of a flat core ( $\delta=0$ ) by Reynolds & Chevalier (Reference Reynolds and Chevalier1984) and later by van der Swaluw et al. (Reference van der Swaluw, Achterberg, Gallant and Tóth2001): $R_{\textrm{PWN}}(t)\propto t^{6/5}$ . On the contrary an analytic solution valid at all times up to reverberation has not been found, while a tentative solution based on the series expansion of the spin-down law was investigated in Bucciantini et al. (Reference Bucciantini, Blondin, Del Zanna and Amato2003). For the flat core case, Bandiera et al. (Reference Bandiera, Bucciantini, Martín, Olmi and Torres2023) have recently shown that a simple solution can be determined based on the fitting of numerical results:

(10) \begin{equation} R_{\textrm{PWN}}(t)\Bigg\vert_{\delta=0} \simeq \mathcal{V}_0\tau_0\frac{\left[ 1+(c_l t/\tau_0 )^{b}\right]^{1/b}}{\left[1+\left(0.82 \,t/\tau_0 \right)^{-a}\right]^{6/(5a)}}\,,\end{equation}

where the parameters $c_l,\, a$ and b depend on the value of the braking index n, while $\mathcal{V}_0\tau_0\simeq 1.91 (L_0\tau_0/E_{\mathrm{sn}})^{1/5} (\tau_0/t_{\mathrm{ch}})R_{\mathrm{ch}}$ , and $R_{\mathrm{ch}}=M_{\mathrm{ej}}^{1/3}\rho_{\mathrm{ism}}^{-1/3}$ is the characteristic radius. A very interesting result is that the evolution of the PWN in the free-expansion phase is poorly affected by the value of the braking index, at least in the range $1.8 \leq n \leq 4.1$ . This can be seen in Figure 3, comparing the radial evolution computed with Equation (10) for $n=1.8$ ( $c_l=0.86310,\, a=0.78492,\, b=0.76490$ ), $n=2.33$ ( $c_l=0.9517,\, a=0.73014,\, b=0.71979$ ), $n=3$ ( $c_l=1.0329,\, a=0.66355,\, b=0.65937$ ) and $n=4.1$ ( $c_l=1.1334,\, a=0.59129,\, b=0.59210$ ), up to a maximum time of $5-10 \tau_0$ , much longer than the duration of this initial phase for typical systems.

Figure 3. Left panel: time evolution of the PWN radius (in units of $V_0\tau_0$ and $\tau_0$ ) for four values of the braking index in the range $n=\{1.8;\,4.1\}$ . In magenta colour we show the pure dipole braking ( $n=3$ ), while in orange the value measured for Crab ( $n=2.33$ , Lyne et al. Reference Lyne, Jordan, Graham-Smith, Espinoza, Stappers and Weltevrede2015; Horvath 2019). Right panel: ratio of the previous two.

As we will better see in Section 3, being characteristic of young objects, and hence of bright systems, the free-expansion phase has received a lot of attention in the past in terms of modelling, with a variety of approaches based on different strategies (Reynolds & Chevalier Reference Reynolds and Chevalier1984; van der Swaluw et al. Reference van der Swaluw, Achterberg, Gallant and Tóth2001; Bucciantini et al. Reference Bucciantini, Blondin, Del Zanna and Amato2003; Gelfand, Slane, & Zhang Reference Gelfand, Slane and Zhang2009; Martín, Torres, & Rea Reference Martín, Torres and Rea2012; Fiori et al. Reference Fiori, Olmi, Amato, Bandiera, Bucciantini, Zampieri and Burtovoi2022; Komissarov & Lyubarsky Reference Komissarov and Lyubarsky2004; Del Zanna et al. Reference Del Zanna, Volpi, Amato and Bucciantini2006; Porth, Komissarov, & Keppens Reference Philippov and Spitkovsky2014; Olmi et al. Reference Olmi, Del Zanna, Amato, Bandiera and Bucciantini2014, Reference Olmi, Del Zanna, Amato, Bucciantini and Mignone2016).

2.4 Reverberation

When the PWN outer boundary hits the SNR RS, the free-expansion phase ends and reverberation starts. In reality, due to Rayleigh-Taylor like instabilities (R-T), that develop at the PWN boundary, as well as asymmetries of different nature (e.g. high proper motion of the pulsar, strong ISM density gradients, asymmetries induced from the stellar explosion), this transition is not as sharp as simplified 1D models would suggest, neither is the following reverberation phase. This phase and its transition have received much less attention due to several reasons: their extreme complexity, the fact that the SNR evolution and properties become more relevant, the PSR kick velocity begins to play a role. Moreover, there are few well characterised systems from an observational point of view that can be taken as benchmark to evaluate the accuracy of a model. Only recently a detailed study of the transition between free-expansion and reverberation have been presented (Bandiera et al. Reference Bandiera, Bucciantini, Martín, Olmi and Torres2023), still limited to a simplified 1D evolution. The time at which this transition happens (that we name $t_{\mathrm{begrev}}$ ) is technically obtained as the intercept between the RS radius and the PWN one. It depends on the relation between the PWN energetics ( $L_{0}\tau_0$ ) and $E_{\mathrm{sn}}$ , and can be computed with some accuracy in the limit $\tau_0 \ll t_{\mathrm{ch}}$ using the following formula by Bandiera et al. (Reference Bandiera, Bucciantini, Martín, Olmi and Torres2023), expressed for simplicity in terms of the variable $\lambda_E=\log_{10}[{(L_0\tau_0)/E_{\mathrm{sn}}}]$ :

(11) \begin{equation} \!\!\!\!\!t_{\mathrm{begrev}}(\lambda_E)\simeq 2.4102\frac{1-\exp(-0.1494+1.1606\,\lambda_E)}{1+\exp(1.6831 + 0.6805\,\lambda_E)}\,t_{\mathrm{ch}}.\end{equation}

An extension to larger $\tau_0$ is possible substituting to $L_0\tau_0$ the energetic of a massive shell interacting with the SNR in place of the PWN, as discussed in Section 5.1 of Bandiera et al. (Reference Bandiera, Bucciantini, Martín, Olmi and Torres2023). The time $t_{\mathrm{begrev}}$ typically ranges between 1 and 2 characteristic times for most of the systems (see Figure 2 of Bandiera et al. Reference Bandiera, Bucciantini, Martín, Olmi and Torres2023).

When reverberation starts, the PWN begins to interact directly with the shocked ejecta in the SNR. The pressure exerted by the ejecta generally induces a strong deceleration of the shell accumulated at the PWN boundary van der Swaluw et al. Reference van der Swaluw, Achterberg, Gallant and Tóth2001) which, in many cases, turns into a compression of the nebula. During this compression the internal energy of the PWN, and its magnetic field, increase. Ultimately, depending on the efficiency of radiation losses and magnetic energy dissipation, the total internal pressure rises to become comparable to the external one, and the compression is suddenly reverted into a new expansion.

One-zone and 1D models predict a number of successive oscillations and expansions before the system reaches a new equilibrium in the Sedov-Taylor phase, and the PWN keeps expanding without bouncing back anymore. This oscillatory behaviour is what gives the name reverberation to this phase. As already pointed out in Bucciantini, Arons, & Amato (Reference Bucciantini, Arons and Amato2011), this extended oscillatory behaviour is an artifact of the 1D approximation. Higher dimensional models in fact show that during reverberation the PWN boundary is largely subjected to R-T like instabilities at the boundary (Blondin, Chevalier, & Frierson Reference Blondin, Chevalier and Frierson2001; Kolb et al. Reference Kolb, Blondin, Slane and Temim2017), producing effective mixing between the PWN material and the outer one. This acts as a sort of viscous term that might halt oscillations already after the first compression. This is why in one-zone models reverberation is usually stopped artificially after the first bounce (Martín et al. Reference Martín, Torres and Rea2012; Torres et al. Reference Torres, Cillis, Martín and de Oña Wilhelmi2014), or when the PWN pressure reaches equilibration with the outer one (Bucciantini et al. Reference Bucciantini, Arons and Amato2011). In simplified 1D models, the effect of the reverberation phase on the dynamics of the PWN can be measured in terms of the compression factor (CF), namely the ratio between its maximum radius (close to the begin of reverberation) and the radius at the maximum of compression, i.e the minimum one. More generally the parameter of interest is the change in the total volume taken by the relativistic non-thermal plasma of the PWN, since it translates in a variation of the nebular magnetic field ( $B\propto R^{-2}$ ). It is in fact this quantity that directly impacts the chance an old PWN might become visible again in X-rays.

Roughly speaking, the effects of the reverberation can be divided into two extreme cases: (i) the PWN is sufficiently energetic, then the compression is relatively small, or even not appreciable, and the PWN continues its expansion almost undisturbed; (ii) the PWN is weak and then overwhelmed by the SNR pressure, contracting down to very small radii (volumes). Low energetic systems are the most critical in terms of modelling, since they might undergo violent compression with important modifications of their multi-wavelength spectral properties. In fact compression enhances the magnetic field and energises particles, leading to an increase in radiative losses. For very high compression efficiency ( $CF\gtrsim 1\,000$ ), the PWN can enter a fast cooling regime, where a large fraction of the particle energy is lost, and the particle energy distribution is strongly modified (Torres, Lin, & Coti Zelati Reference Torres, Lin and Coti Zelati2019). Bandiera et al. (Reference Bandiera, Bucciantini, Martín, Olmi and Torres2023) show that these extreme behaviour is not expected to be common within the present PWNe population, with only a relatively limited number of objects having possible CF larger than few hundreds. On the contrary, the vast majority of the systems undergo a rather small compression, with CF ranging from a few to tens. However, Bandiera et al. (Reference Bandiera, Bucciantini, Martín, Olmi and Torres2023) do not consider the effect of radiative losses and thus, this estimate might change depending on the relevance of losses in the first evolutionary phase.

2.5 Post-reverberation

The reverberation phase, together with the pulsar proper motion, is what shapes the PWN in the later evolutionary phases. Gradients in the ambient medium density are known to impact the evolution of the SNR reverse shock (Ferreira & de Jager Reference Ferreira and de Jager2008; Kolb et al. Reference Kolb, Blondin, Slane and Temim2017) and, in general, one expects that as time passes the level of asymmetries in the system grows. This implies that simplified one-zone or 1D models became progressively less accurate in the description of these systems, and attempts to include these extra effects (Gelfand et al. Reference Gelfand, Slane and Zhang2009) less and less reliable. On top of this, the mixing due to R-T instability can be so strong as to disrupt completely the PWN as a coherent object. It is then clear that a PWN in this stage would likely be very far from being spherically symmetric, rather having a complex distorted shape. As a consequence, modelling this phase can be very demanding: the need to properly account for asymmetries requires the use of multi-D models. Being the shape and properties of the evolved PWN strongly dependent on its previous history, one has to follow its entire evolution, which implies a large dynamical range in term of temporal and spatial scales. Moreover the morphology of the system, especially in the presence of mixing and instabilities, is very dependent on the model dimensionality and a comprehensive description can only be done in 3D. A very beautiful example of a 3D model of a largely asymmetric evolved system can be found in Kolb et al. (Reference Kolb, Blondin, Slane and Temim2017).

Unfortunately there are not that many systems, with a detailed and robust characterisation, that can be used to benchmark our theoretical models for this evolutionary phase. Some of them like G327.1-1.1 (Temim et al. Reference Temim, Slane, Kolb, Blondin, Hughes and Bucciantini2015; Ma et al. Reference Lyutikov, Komissarov and Porth2016), IC443 (Swartz et al. Reference Swartz2015), W44 (Frail et al. Reference Frail, Giacani, Goss and Dubner1996; Petre, Kuntz, & Shelton Reference Petre, Kuntz and Shelton2002), have been the targets of previous analyses (van der Swaluw, Downes, & Keegan Reference van der Swaluw, Downes and Keegan2004; Bucciantini et al. Reference Bucciantini, Arons and Amato2011), but they are too few to provide a significative sample.

From the one-zone description, one expects that when reverberation ends, with the oscillatory behaviour almost completely dissipated, the PWN relaxes to a steady subsonic expansion (van der Swaluw et al. Reference van der Swaluw, Achterberg, Gallant and Tóth2001). This was found to happen on a long time-scale of $\sim$ $20t_{\mathrm{ch}}$ . This time is possibly even longer, since the relaxation of the contact discontinuity of the SNR has been recently show to happens on times of $\sim$ $35t_{\mathrm{ch}}$ , while between $20t_{\mathrm{ch}}$ and $30t_{\mathrm{ch}}$ small oscillations are still present in the SNR radius (Bandiera et al. Reference Bandiera, Bucciantini, Martín, Olmi and Torres2021).

2.6 Late phase—bow shocks

Ultimately PWNe are likely to end their life as bow shock nebulae (BSPWN), due to the fact that a large fraction of pulsars is born with high kick velocity (100–500 km s $^{-1}$ , Cordes & Chernoff Reference Cordes and Chernoff1998; Arzoumanian, Chernoff, & Cordes Reference Arzoumanian, Chernoff and Cordes2002; Faucher-Giguère & Kaspi Reference Faucher-Giguère and Kaspi2006; Sartore et al. Reference Sartore, Ripamonti, Treves and Turolla2010;Verbunt, Igoshev, & Cator Reference Verbunt, Igoshev and Cator2017), and as such it is bound to emerge out of the parent SNR before the pulsar spin-down luminosity becomes so weak as to make particle acceleration and non-thermal synchrotron emission negligible. Considering the typical SNR decelerated expansion of the Sedov-Taylor phase, one can easily estimate in a few tens of thousand of years the time the pulsar takes to emerge the SNR bubble, to be compared with the longer typical age of pulsars in the Galaxy ( $\sim$ $10^6$ yr). From that moment on, the pulsar interacts directly with the ambient medium. Given the high speed of the pulsar, larger than the typical sound speed of the surrounding medium, its motion generally becomes supersonic already inside the SNR (Gaensler & Slane Reference Gaensler and Slane2006). IC443 and W44 are exemplary cases of BSPWNe still confined within their parent SNR (Swartz et al. Reference Swartz2015; Frail et al. Reference Frail, Giacani, Goss and Dubner1996; van der Swaluw et al. Reference van der Swaluw, Downes and Keegan2004). The supersonic motion induces the formation of a bow shock around the pulsar and its nebula, reshaping drastically the PWN. Now it does not appear as a bubble anymore but, on the contrary, it assumes an elongated cometary-like shape.

In BSPWNe the pulsar is located at the bright head of a very elongated tail, extending in the direction opposite to the pulsar motion. The new morphology of the system depends on the balance of the ram pressure of the pulsar wind with respect to that of the incoming (in the frame of the PSR) ambient medium. This balance actually determined the thickness of the PWN at the head front, representing the characteristic dimension of a BSPWN, known as stand off distance:

(12) \begin{equation} d_0=\sqrt{\frac{L}{4\pi c \,\rho\, v^2_{\mathrm{psr}}}}\,,\end{equation}

where $\rho$ is the density of the ambient medium. Gradients in the ambient medium reflect into asymmetries in the bow shock shape (Vigelius et al. Reference Vigelius, Melatos, Chatterjee, Gaensler and Ghavamian2007). On the contrary, even larger asymmetries in the pulsar wind energy distribution tend to have minor effects on the overall bow-shock shape (Olmi & Bucciantini Reference Olmi and Bucciantini2019a), remaining mostly concentrated in the head, that is typically not resolved even in the deepest Chandra’s observations.

The PSR wind, shocked in the head, is now diverted along the tail, where the plasma becomes less magnetised and more turbulent with increasing distance from the pulsar (Wilkin Reference Wilkin1996; Bucciantini & Bandiera Reference Bucciantini and Bandiera2001; Bucciantini Reference Bucciantini2002). How turbulent and magnetised the tail is depends not only on the PSR wind magnetisation, but also on the level of asymmetry in its energy distribution, on the inclination of the PSR spin axis with respect to the direction of the kick velocity, and on the relative inclination and strength of the ambient medium magnetic field (Olmi & Bucciantini Reference Olmi and Bucciantini2019a).

2.7 Other environments

PWNe represent proto-typical systems where one can witness the interaction of a relativistic magnetised outflow with a confining surrounding environment. As such they form a paradigm for a variety of different classes of high-energy astrophysical objects. Moreover, PWNe can be found more or less anywhere a NS is active, even if the activity is just a transient one. Indeed, from an historical point of view, they have often been ‘the first ones’ where high-energy astrophysical processes have been discovered/studied/understood, and much of the theory developed for their study has found application elsewhere. Here, we briefly review some of these other environments.

• PWNe, from the point of view of fluid dynamics, acceleration properties, and emission mechanisms, are representative of the very broad class of relativistic wind bubbles. The physics that have been thoroughly investigated to explain the observed presence of non-thermal particles or the properties of relativistic outflows at large distances from their engine (Lyubarskij Reference Lyubarskij1992; Li Reference Li1993; Barkov & Komissarov Reference Barkov and Komissarov2008; Chen & Zhang Reference Chen and Zhang2021; Zhu, Zhang, & Fang Reference Zhu, Zhang and Fang2019; Kagan et al. 2015; Sironi, Keshet, & Lemoine Reference Sironi, Keshet and Lemoine2015) have proven to be quite general and with a much larger applicability, in terms of astrophysical systems ranging from gamma ray bursts (GRBs, Thompson Reference Thompson1994) to active galactic nuclei jets and radio lobes (Uzdensky Reference Uzdensky, Gonzalez and Parker2016, Reference Uzdensky2018).

• Bright persistent X-ray emitting PWNe require a steady engine, providing the necessary particle and energy injection. It was held for a long time that only canonical PSRs, with their active magnetospheres, capable of supporting pair creation cascades with high multiplicity (Timokhin & Harding Reference Timokhin and Harding2019) could be surrounded by such nebulae. However, recently has emerged a more dynamical picture, where transient PWNe can form and shine, even around neutron stars with inefficient pair creation, as in the case of magnetars and/or rotating radio transients (RRATs, Camero-Arranz et al. Reference Camero-Arranz2013; Vink & Bamba Reference Vink and Bamba2009; Tanaka Reference Tanaka2016; Torres Reference Torres2017; Granot et al. 2017; Blumer, Safi-Harb, & McLaughlin Reference Blumer, Safi-Harb and McLaughlin2017; Torres Reference Torres, Perera, Preston and Sanidas2018; Margalit & Metzger Reference Margalit and Metzger2018). These transient nebulae, typically associated with bursting/active phases, are still electromagnetic powered, and shine in synchrotron, even if it is still unclear if the neutron star activity leads to freshly injected particles, or simply re-energizes particles already present and injected at earlier times. Even the driving energy reservoir is not well constrained in general, though magnetic energy has been suggested for the case of magnetars. Moreover, violent events, as giant flares, can produce transient nebulae, whose evolution resembles, even if at a faster pace, that of regular PWNe (Gaensler et al. Reference Gaensler2005; Gelfand et al. Reference Gelfand2005).

• PWNe can also form in binary systems, where the pulsar wind is confined by the ram pressure of the companion star. In case of massive stars, with strong equatorial flows, this interaction leads to the formation of a bow-shock, whose dynamics, and by consequence emission, can be both highly variable due to the possible large orbital eccentricity, and to the fact that it is also strongly affected by centrifugal and Coriolis forces. As in normal bow-shocks, particles can be accelerated to high energy, but unlike regular systems, the orbital dynamics can lead to a substantial mixing with the stellar material, a high level of turbulence, and the development of multiple shocks, with distributed acceleration. Moreover the presence of a bright source makes IC a relatively strong cooling process, to the point that it is unclear if the observed X-ray emission is synchrotron or IC (Neronov & Chernyakova Reference Nakamori2007; Kargaltsev et al. 2014; Bednarek & Sitarek Reference Bednarek and Sitarek2013; Zdziarski, Neronov, & Chernyakova Reference Wilkin2010; Paredes-Fortuny et al. Reference Oort and Walraven2015; Bosch-Ramon et al. Reference Bosch-Ramon, Barkov and Khangulyan2012) These systems offer us a unique opportunity to investigate the PSR wind and its acceleration properties, at distances much smaller than in regular PWNe, and can have important consequences on a large set of dissipative wind models (Kirk & Skjæraasen Reference Kirk, Lyubarsky and Petri2003; Kirk Reference Khangulyan, Arakawa and Aharonian2006; Lyubarsky Reference Lyubarsky2010; Takamoto, Inoue, & Inutsuka Reference Sudoh, Linden and Beacom2012; Cerutti & Philippov Reference Cerutti and Philippov2017; Cerutti, Philippov, & Dubus Reference Cerutti, Philippov and Dubus2020; Huber et al. 2021a; Huber, Kissmann, & Reimer 2021b).

• Proto-magnetar wind nebulae, and in general relativistic outflows from millisecond rotating newly born magnetars, have been invoked to explain both long (Usov Reference Torres, Lin and Coti Zelati1992; Thompson Reference Temim, Slane, Kolb, Blondin, Hughes and Bucciantini2007; Bucciantini et al. Reference Bucciantini, Quataert, Arons, Metzger and Thompson2007, Reference Bucciantini, Quataert, Arons, Metzger and Thompson2008; Metzger et al. Reference Mestre, Torres, de Oña Wilhelmi and Martí2011b) and short GRBs (Bucciantini et al. Reference Bucciantini, Metzger, Thompson and Quataert2012; Wang Reference Volpi, Del Zanna, Amato and Bucciantini2017; Metzger & Piro Reference Metzger, Giannios, Thompson, Bucciantini and Quataert2014; Rezzolla & Kumar 2015), mostly because the detection of the so called ‘late activity’ (Gompertz, O’Brien, & Wynn 2014; Rowlinson et al. Reference Romani and Ng2013) has made almost mandatory to assume long lived engines, hardly compatible with the timescale for disk accretion in a stellar mass black hole. The key idea is that the relativistic outflow that we think is at the origin of the prompt gamma-ray emission, is nothing else than the magneto-centrifugally driven neutrino-wind, coming from the cooling proto-NS. The dynamics of such wind reaches rapidly the force-free limits, as the baryon loading rapidly drops in time (Metzger et al. Reference Melatos2011a), and its interaction with the surrounding layers of the progenitor star or the ejecta of the binary merger (for SGRBs), leads to the formation of a hot relativistic wind bubble, that is both a reservoir of energy, to be released at later time, as well as a reservoir of high energy photons that can lead to the appearance of a so called kilo-nova (Ciolfi & Siegel Reference Ciolfi and Siegel2015; Siegel & Ciolfi Reference Scuderi2016a,b). Lower energy system have been considered also to justify a possible continuum of explosion phenomenology down to super-luminous and broad line Ib/c supernovae (Dessart et al. Reference Dessart, Hillier, Waldman, Livne and Blondin2012; Moriya, Chen, & Langer Reference Moran, Shearer, Mignani, Słowikowska, De Luca, Gouiffès and Laurent2017; Soker & Gilkis Reference Slane, Helfand, van der Swaluw and Murray2017; Milisavljevic et al. Reference Metzger and Piro2018; Wang,Wang, & Dai Reference Wang2019; Vurm & Metzger Reference Vink and Bamba2021; Shankar et al. Reference Schweizer, Bucciantini, Idec, Nilsson, Tennant, Weisskopf and Zanin2021; Dessart Reference Dessart2019; Margalit et al. Reference Margalit and Metzger2018).

3.1 One-zone models

In one-zone (or equivalently 0+1) models, the PWN is described as a uniform bubble in interaction with the surrounding SNR, being subjected to adiabatic, radiation and possibly particles losses. Energy and magnetic field are injected into the bubble with constant relative efficiencies by the PSR, following its spin-down luminosity (Equation (2). The magnetic field energy within the bubble can be evolved either assuming magnetic flux conservation (Pacini & Salvati Reference Pacini and Salvati1973; Gelfand et al. Reference Gelfand, Slane and Zhang2009; Torres et al. Reference Torres, Cillis, Martín and de Oña Wilhelmi2014) or a constant ratio between magnetic and plasma pressures (Bucciantini et al. Reference Bucciantini, Arons and Amato2011). In the early free-expansion phase, if radiation losses are negligible, the two approaches are equivalent (however they might have different implications in terms of wind magnetisation), but once radiation losses become important they might lead to sizeable different results (Bucciantini et al. Reference Bucciantini, Arons and Amato2011).

One-zone models are all based on the thin-shell approximation: the shell is evolved considering mass and momentum conservation (Equations (4)–(5) in Section 2.3) and its radius is equated with that of the PWN. To the time evolution of the radius, one then couples a description of the evolution of the internal particles distribution function subject to injection and losses, from which one finally derives the observed PWN multi-wavelength spectrum.

Spectral evolution models can be traced back to the original work by Pacini & Salvati (Reference Pacini and Salvati1973) that, assuming an approximated description of the temporal evolution of the PWN, limited to its initial free-expanding phase, derived the evolution of the spectral energy distribution function showing its relation with the injected particle distribution. This work was later extended with a model also describing the interaction with the SNR, but still limited to the interaction with the unshocked ejecta (Reynolds & Chevalier Reference Reynolds and Chevalier1984), and considering possible variations to the particles injection (Atoyan Reference Atoyan1999).

Despite their evident over-simplification, one-zone models have been—and still are—widely used. Results from one-zone models have in fact proved to give a good description of the global properties of young PWNe, allowing one to model the evolution of these systems and their spectral properties, and to sample a large parameter space, in ways that are not possible with current multi-dimensional models. They can be easy implemented and are not much demanding in terms of numerical resources (Gelfand et al. Reference Gelfand, Slane and Zhang2009; Tanaka & Takahara Reference Takamoto, Inoue and Inutsuka2010; Bucciantini et al. Reference Bucciantini, Arons and Amato2011; Martín et al. Reference Martín, Torres and Rea2012; Tanaka & Takahara Reference Tanaka and Takahara2013; Torres et al. Reference Torres, Cillis, Martín and de Oña Wilhelmi2014; Gelfand Reference Gelfand and Torres2017; Fiori et al. Reference Fiori, Zampieri, Burtovoi, Caraveo and Tibaldo2020). These characteristics make them very much appealing for large populations studies (e.g. Torres et al. Reference Torres, Cillis, Martín and de Oña Wilhelmi2014; Fiori et al. Reference Fiori, Olmi, Amato, Bandiera, Bucciantini, Zampieri and Burtovoi2022), but on the other hand they must be used with care to follow the evolution beyond free-expansion. In (Bandiera et al. Reference Bandiera, Bucciantini, Martín, Olmi and Torres2020, Reference Bandiera, Bucciantini, Martín, Olmi and Torres2023) it has been in fact shown that one-zone models predict excessive compressions of the nebula during reverberation, possibly leading to the burn-off of a huge amount of the emitting particles, changing dramatically the spectral evolution. Torres et al. (Reference Torres, Lin and Coti Zelati2019) show that in cases of extreme compressions ( $CF\gtrsim1\,000$ ) a super-efficient phase can appear, when the PWN luminosity is so enhanced to make it visible again at X-rays at later times. One-zone models tend to overestimate the PWN compression mainly because the simplified description of the pressure in the SNR, in general assumed to be equal to the central pressure from the Sedov solution, or the pressure at the SNR FS scaled with some arbitrary factor (for a complete description see Bandiera et al. Reference Bandiera, Bucciantini, Martín, Olmi and Torres2023). This kind of approximation indeed introduces large errors when the PWN starts to interact with the SNR. Recently Bandiera et al. (Reference Bandiera, Bucciantini, Martín, Olmi and Torres2023) have in particular shown that the pressure in the SNR is very far from the Sedov solution during the entire first compression, and that in the aforementioned assumptions lead to a consistent overestimation of the outer pressure. This means that the number of systems undergoing a super-efficient phase, might be much smaller than expected, if a correct description of the SNR pressure is considered during reverberation.

A preliminary attempt to extend one-zone models beyond free-expansion has been recently made by Fiori et al. (Reference Fiori, Olmi, Amato, Bandiera, Bucciantini, Zampieri and Burtovoi2022), where the SNR pressure is shaped on the results of 1D hydrodynamic simulations of the PWN-SNR interaction. Conversely to Bandiera et al. (Reference Bandiera, Bucciantini, Martín, Olmi and Torres2023), in this case there was no specific focus in reproducing correctly the dynamical effects of the interaction between the PWN and the surrounding SNR.

The extension of one-zone models is generally limited to the first series of compressions and re-expansions of the PWN, while in few cases they have been extended to longer times halting by hands the oscillations to match the SNR pressure from the Sedov solution (Torres et al. Reference Torres, Lin and Coti Zelati2019).

3.2 1D models

In 1D models the evolution and properties of the PWN and the SNR are given as a function of the radial coordinate. The first 1D model of a PWN was put forward by Rees & Gunn (Reference Porth, Komissarov and Keppens1974) to describe, as many others later on, the Crab nebula. Despite the very simplified description of the system, that model was able to predict a number of features and to give them a physical interpretation: the appearance of the ultra-relativistic un-shocked wind as an under-luminous area in the inner nebula, later observed by Weisskopf et al. (Reference Wang, Wang and Dai2000) at X-rays; the nebula shrinkage with increasing frequency as sign of the synchrotron emission and central particle injection; the nebular magnetic field close to the equipartition value; the Lorentz factor of the wind, its magnetisation and the injection rate of particles from the synchrotron luminosity. This preliminary model was then elaborated and extended by Kennel & Coroniti (Reference Kargaltsev, Rangelov, Pavlov, Strakovsky and Blokhintsev1984a, b), that provide estimates of many characteristic quantities of the Crab nebula based on the full solution of the relativistic MHD equations within the PWN itself. Few years later, Emmering & Chevalier (Reference Emmering and Chevalier1987) found a time dependent analytic solution of the same problem.

Many other 1D models have been then developed in later years, all based on the numerical implementation of the equations describing the PWN-SNR interaction and evolution, both in the classic and relativistic regimes and in the hydrodynamic (HD) or MHD frameworks, extending to a longer evolution than the free-expansion phase (van der Swaluw et al. Reference van der Swaluw, Achterberg, Gallant and Tóth2001, Reference van der Swaluw, Downes and Keegan2004; Bucciantini et al. Reference Bucciantini, Blondin, Del Zanna and Amato2003; de Jager, Ferreira, & Djannati-Ataï Reference de Jager, Ferreira, Djannati-Ataï, Aharonian, Hofmann and Rieger2008; Bandiera et al. Reference Bandiera, Bucciantini, Martín, Olmi and Torres2023).

One of the longer standing problems in PWNe physics, the so called sigma-paradox (Melatos Reference Martín, Torres and Rea1998), arose as consequence of these 1D models, and puzzled the community for more than thirty years: in order to explain the existence of the TS, the average magnetisation $\sigma$ in the nebula must be quite low: $\sigma\simeq \rm{few}\times 10^{-3}$ (Kennel & Coroniti Reference Kennel and Coroniti1984b). But theoretical models of pulsar magnetospheres predict a much larger magnetisation at the light cylinder of the star: $\sigma\sim 10^{4}$ (Kirk et al. Reference Kirk2009; Arons Reference Arons2012). To make compatible these two opposite predictions, a huge amount of magnetic dissipation must be considered along the path separating the light cylinder and the pulsar wind termination shock, converting the pulsar wind from a Poynting dominated outflow to a particle dominated one. Magnetic dissipation is actually expected to occur in the current sheet of the alternating pulsar wind (Lyubarsky Reference Lyubarskij2003; Sironi & Spitkovsky Reference Siegel and Ciolfi2011), but $\sim$ 7 orders of magnitude are simply too many to be accounted for with any known process. We want to remark here that the magnetisation at the wind termination shock not only affects the PWN dynamics, but it is also an important parameter that controls the efficiency of the various acceleration mechanisms (Amato Reference Amato2014). As we will discuss in the following paragraphs, going multi-D is the way to reduce—if not solve—this long standing problem, since the augmented dimensionality allows increasing the amount of magnetic dissipation.

Results of one-zone and 1D models are in excellent agreement in the free-expansion phase, while the accordance disappears as soon as the PWN enters reverberation. As for one-zone models, the reliability of 1D models is in any case limited to at most the beginning of reverberation. Unlike one-zone, however, 1D Lagrangian models (Bandiera et al. Reference Bandiera, Bucciantini, Martín, Olmi and Torres2023) can be extended to the following evolution (the first compression and the subsequent sequence of multiple compressions and re-expansions). Recall that, in any case, these are just an artifact of the reduced dimensionality. 1D simulations show a sort of damped oscillations, with compressions generally becoming less severe as time passes by. These behaviour is not found in a multi-D description of the PWN/SNR system, where the mixing produced by boundary instabilities helps in mitigating oscillations. However, one should be careful not to confuse the observed angular size/dimension of the PWN with its effective volume (the one that matters for compression), once mixing becomes important.

3.3 2D models

Overcoming the limitations of 1D models became urgent with the first detailed images of the inner Crab nebula at optical (Hubble Space Telescope, Hester et al. 1995) and X-rays (Chandra, Weisskopf et al. Reference Wang, Wang and Dai2000). The Crab nebula was in fact the first source where a beautiful bright jet-torus structure, with an equatorial torus and two opposite polar jets, was identified. Later on other systems showed up the same properties (Gaensler et al. Reference Gaensler, Arons, Kaspi, Pivovaroff, Kawai and Tamura2002; Lu et al. Reference Leung and Ng2002; Romani & Ng Reference Roberts, Romani, Johnston and Green2003; Camilo et al. Reference Camilo, Gaensler, Gotthelf, Halpern and Manchester2004; Slane et al. Reference Sironi and Spitkovsky2004; Romani et al. Reference Roberts, Romani and Johnston2005), now believed to be a common feature of young PWNe. It is clear that such a complex structure cannot be reproduced with the simplified geometry of 1D models, neither by slighting modifying the 1D approach. Indeed the first 2D analytic model of the Crab nebula by Begelman & Li (Reference Begelman and Li1992), considering the toroidal structure of the magnetic field, was only capable to explain its observed prolate shape.

A critical point for the theoretical description of the inner nebula was the appearance of the polar jets so close to the pulsar. They in fact seem to form in the un-shocked relativistic wind, where magnetic collimation is known to be poorly efficient. Theoretical and numerical models of the relativistic wind emanating from the pulsar (Begelman & Li Reference Begelman and Li1992; Contopoulos, Kazanas, & Fendt Reference Contopoulos, Kazanas and Fendt1999; Bogovalov Reference Bogovalov1999, and later Komissarov Reference Kolb, Blondin, Slane and Temim2006; Spitkovsky Reference Spitkovsky2006; Timokhin Reference Thompson2006) already predicted a non-uniform distribution of the wind energy flux at the termination shock, with most of the energy concentrated in the equatorial plane (the so called split-monopole models), but there was no evidence for a collimation of part of the flux in the polar regions. The solution to the jet formation was found only few years later: modelling the dynamics of the plasma with an anisotropic distribution of the energy flux in the wind, the terminations shock was found to become oblate, with larger extension in the equatorial region than in the polar one (Lyubarsky Reference Lyne, Jordan, Graham-Smith, Espinoza, Stappers and Weltevrede2002; Bogovalov & Khangoulian Reference Bogovalov and Khangoulian2002a; Bogovalov & Khangoulyan Reference Bogovalov and Khangoulyan2002b; Khangoulian & Bogovalov Reference Kennel and Coroniti2003). Polar jets can then form due to magnetic hoop stresses in the post-shock plasma, immediately beyond the polar front of the shock, and appear closer to the pulsar than the torus thanks to the oblate morphology of the TS itself.

This theoretical prediction was verified later on with the use of relativistic 2D MHD numerical simulations (Del Zanna, Amato, & Bucciantini Reference Del Zanna, Amato and Bucciantini2004; Komissarov & Lyubarsky Reference Komissarov and Lyubarsky2004; Bogovalov et al. Reference Bogovalov, Chechetkin, Koldoba and Ustyugova2005; Del Zanna et al. Reference Del Zanna, Volpi, Amato and Bucciantini2006). Moreover, 2D models show that the formation of jets is only possible starting from a minimum wind magnetisation in the pulsar wind of $\sigma\sim 10^{-2}$ , one order of magnitude larger than that originally set by 1D models. The increase of dimensionality then appears as the first possible way to alleviate the sigma-paradox.

Thanks to its luminosity and vicinity, the largest part of 2D models were made to investigate the properties of the Crab nebula that was—and still is—the perfect source to look at for a detailed comparison. They were extremely successful at reproducing many of its features down to very fine details, especially in the inner nebula: maps of the synchrotron emission (Del Zanna et al. Reference Del Zanna, Volpi, Amato and Bucciantini2006; Olmi et al. Reference Olmi, Del Zanna, Amato, Bandiera and Bucciantini2014); polarisation properties (Bucciantini et al. Reference Bucciantini, del Zanna, Amato and Volpi2005); the complete spectrum, from radio to gamma-rays (Volpi et al. Reference Volpi, Del Zanna, Amato and Bucciantini2008); the variability at small scales in the inner nebula (Camus et al. Reference Camus, Komissarov, Bucciantini and Hughes2009; Olmi et al. Reference Olmi, Bucciantini and Morlino2015; Lyutikov, Komissarov, & Porth 2016). In particular the multi-wavelength appearance of variable arc-like bright structures (the so called wisps), forming close to the TS and moving outwards at a consistent fraction of c, was shown to be a perfect tool to trace the properties of the underlying plasma and to constrain the location and mechanism of particle acceleration (Olmi et al. Reference Olmi, Bucciantini and Morlino2015).

On the contrary, results of 2D models at larger scales are limited by the imposed axisymmetry, with only the global properties of the nebula, such as the extension and the shrinkage with the increasing energy, reproduced reasonably well. Axisymmetry is in fact the major limitation of 2D models: it induces the artificial accumulation of magnetic loops along the polar axis of the nebula, producing enhanced compression of the magnetic field. This imposes an upper limit to the wind magnetisation of $\sigma \lesssim 0.1$ , otherwise the polar compression of the field becomes so strong that the collimated magnetic flux punches the nebular shell at the polar boundaries. The consequence of this non-redistribution of magnetic field is that polar regions have an (excessively) intense magnetic field, while the rest of the nebula is under-magnetized, with an average magnetic field, for the Crab, well below the expected value of $\sim$ $150\,\mu$ G. To reproduce the PWN spectrum then one is forced on one side to inject an excessive number of particles, on the other to steepen by hand the injection spectrum of the X-ray emitting component to alleviate the average lower energy losses, and match the synchrotron emission.

The problem introduced by the wrong geometry of the magnetic field remains very evident if looking at the morphology of the nebula at large scales: the radio emission, expected to be rather uniform, indeed resembles the shape of the magnetic field (it is concentrated around the polar axis); the IC emission is overestimated, approximately by the same factor by which the averaged magnetic field is underestimated (Olmi et al. Reference Olmi, Del Zanna, Amato, Bandiera and Bucciantini2014).

Despite these limitations, 2D models remain appealing in many cases: their results are in good agreement with more complex 3D models if limited to the inner nebula and small scales; they allow for a longer evolution of the PWN and for a larger investigation in terms of number of sources and physical parameters with respect to their equivalent in 3D, due to the lower numerical cost. A possibility to use 2D models to infer the large scale properties of PWNe has been investigated in Olmi & Torres (Reference Olmi, Del Zanna, Amato and Bucciantini2020), where the HD scheme has been preferred to the MHD one to avoid the aforementioned problems with the field geometry, affecting the large scale emission. In that work the magnetic field is excluded from the dynamics but traced numerically (with a recipe linking it to the thermal pressure), approximately accounting for the particles losses during the evolution.

Finally, 2D HD and MHD models have also been successfully used to investigate the formation of bow shock PWNe produced by fast moving pulsars escaped from their SNRs, and their variations depending on the properties of the ambient medium (Toropina et al. Reference Timokhin and Harding2001; Bucciantini Reference Bucciantini2002; Bucciantini et al. Reference Bucciantini, del Zanna, Amato and Volpi2005; Olmi, Bucciantini, & Morlino Reference Olmi and Bucciantini2018; Toropina, Romanova, & Lovelace Reference Timokhin2019).

3.4 3D models

With the first 3D MHD simulations of the Crab nebula (Porth, Komissarov, & Keppens Reference Petre, Kuntz and Shelton2013) it became apparent that the solution to the sigma-paradox was finally in reach. We have already seen how moving from 1D to 2D permitted to gain more than one order of magnitude in wind magnetisation. Thanks to the development of efficient processes of magnetic dissipation in the nebula, a magnetisation $\sigma\geq 1$ becomes finally accessible with 3D simulations, drastically reducing the impact of the sigma-paradox. Configurations with a toroidal magnetic field are subject to current driven instabilities (kink-like, Begelman Reference Begelman1998; Nalewajko & Begelman Reference Morlino, Lyutikov and Vorster2012), and signs of such a process have been detected in different PWNe (Mori et al. Reference Milisavljevic, Patnaude, Chevalier, Raymond, Fesen, Margutti, Conner and Banovetz2004; Pavlov et al. Reference Pavan2003). Numerical simulations of the Crab nebula jet proved not only that the kink instability efficiently develops in the nebula, but that it is also responsible for the variations seen in the jets at different epochs (Mignone et al. Reference Metzger, Giannios, Thompson, Bucciantini and Quataert2013). Mixing of the magnetic field induced by kink-like instabilities is so efficient that, even if the initial field configuration is fully toroidal, a poloidal component raises rapidly in 3D, becoming even dominant immediately outside the inner nebula (see e.g. Porth et al. Reference Philippov and Spitkovsky2014; Olmi et al. Reference Olmi, Del Zanna, Amato, Bucciantini and Mignone2016). The mixing causes on one side the magnetic field geometry to become much more complex than in 2D, but also results in a more uniform distribution of the magnetic field in the nebular volume (see Figure 4).

Figure 4. 3D plot of the magnetic field lines (left side) and slice of the magnetic field intensity (right side) from the 3D simulation of the Crab nebula presented in Olmi et al. (Reference Olmi, Del Zanna, Amato, Bucciantini and Mignone2016), at a simulated age of $\sim$ 250 yr. Stream lines are coloured with the intensity of the field as from the right panel, but using a different colour scale to highlight the different structures. Both figures have been produced using VisIt (Childs et al. Reference Childs2012).

Despite being the best way to account for the properties of a PWN, the use of 3D models is limited by the huge amount of resources (time/numerical) they need. The spatial scales that must be reproduced are extremely different, from the injection region of the pulsar wind (smaller than the TS) to the PWN contact discontinuity, and the surrounding SNR bubble in case of young systems. Moreover to correctly reproduce the pulsar wind, the injection site must be solved with a minimum number of grid cells ( $\gtrsim$ 10–20), and the injection region must always remain detached from the radius of the termination shock, otherwise the correct jump conditions at the shock cannot be ensured. This translates in the necessity of a very high resolution at the center of the numerical domain, mapping a zone which represents only $\sim$ 1/100 of the global size of the system. The grid is then usually optimized with the use of an adaptive mesh technique, able to increase or decrease the resolution as needed, allowing to save a large amount of time, or with expanding grids. Nevertheless the resources needed to run such models remain large and still prohibitive for long term evolution. To date the longest 3D MHD simulation of the Crab nebula reproduces $\sim$ $1/4$ of the age of the source (Olmi et al. Reference Olmi, Del Zanna, Amato, Bucciantini and Mignone2016) and it required few millions of core/hours of computational time and several months to be run. A longer evolution ( $\sim$ 7 500 yr) has been reached with HD simulations and expanding grids in the case of a high speed PWN ( $v_{\mathrm{kick}}\simeq 300$ km s $^{-1}$ ) interacting with a composite SNR (Kolb et al. Reference Kolb, Blondin, Slane and Temim2017), producing a strongly asymmetric system out of the reverberation phase.

A different approach has been applied to the modelling of bow shock pulsar wind nebulae. For these systems, going 3D is particularly important to correctly capture the magnetic field topology and the development of turbulence in the tail, which are strongly affected by geometric constraints and relevant for the observational properties of the source. A first attempt to 3D modelling of bow shocks, limited to the classical HD regime, was presented by Vigelius et al. (Reference Vigelius, Melatos, Chatterjee, Gaensler and Ghavamian2007). Then, in recent years, the extension to the MHD relativistic regime was investigated by different groups (Barkov, Lyutikov, & Khangulyan Reference Barkov, Lyutikov and Khangulyan2019; Olmi & Bucciantini Reference Olmi and Bucciantini2019a). In all these cases, the bow shocks have been modelled directly assuming the pulsar in interaction with the ISM, not considering then the transitional phase when the star is emerging from the SNR bubble. That transitional phase was indeed investigated in van der Swaluw et al. (Reference Uzdensky, Gonzalez and Parker2003).

The PSR reference frame is generally used, with the star being positioned at a specific point of the computational domain. The star then sees the ambient medium as an incoming, cold flow, that might be magnetized or un-magnetized, depending in the specific case. A large sample of different configurations has been investigated, varying the inclination of the magnetic field and pulsar-spin axis, the direction of the pulsar motion compared to the first two, the magnetisation and level of anisotropy of the pulsar wind. The evolution is then followed for enough time to ensure the system dynamics has reached a relaxed quasi-steady regime. Emitting properties can be computed using the same approach generally used for young systems. Particles responsible for the emission are injected at the wind termination shock with a broken power law distribution and their emissivity is computed as discussed in the following Section 4, with different possible assumptions on their density in the PWN (e.g. they can be considered as uniformly distributed or with a distribution shaped on the thermal pressure, see Olmi & Bucciantini Reference Ng and Romani2019b).

In a recent work (Olmi & Bucciantini Reference Olmi and Bucciantini2019c) we have addressed the possibility for particles to escape the bow shock, modelling the evolution of the particle trajectories in the electric and magnetic fields of the MHD simulations.

3.5 Highlights from 3D MHD models of young PWNe

Here we summarize the most important findings of 3D numerical simulations of young—Crab like—PWNe, according to Porth et al. (Reference Petre, Kuntz and Shelton2013, Reference Philippov and Spitkovsky2014) and Olmi et al. (Reference Olmi, Del Zanna, Amato, Bucciantini and Mignone2016).

• Solution of the sigma-paradox: values of the wind magnetisation at injection larger than unity become possible thanks to the efficient magnetic dissipation produced in 3D. Actually the average value of the magnetic field when the system has reached a self-similar evolution (at ages $\gtrsim$ 150 yr) is still a factor of $\sim$ 1.5–2 lower than the expected one (see Figure 6 in Olmi et al. (Reference Olmi, Del Zanna, Amato, Bucciantini and Mignone2016).

• Development of poloidal magnetic field: despite the magnetic field at injection is purely toroidal, turbulence and high-speed polar flows rapidly modify its topology. A polar component easily develops immediately beyond the inner nebula, and becomes even dominant in the outer nebula (see e.g. Figure 11 in Porth et al. Reference Philippov and Spitkovsky2014 or Figure 8 in Olmi et al. Reference Olmi, Del Zanna, Amato, Bucciantini and Mignone2016). The complex structure of the magnetic field topology and its intensity can be seen in Figure 4.

• Despite the complex topology of the magnetic field, the magnetic pressure in the nebula is rather uniform on large scales (see Figure 9 in Olmi et al. Reference Olmi, Del Zanna, Amato, Bucciantini and Mignone2016), as well as the total pressure (magnetic and thermal, see Figure 3 in Porth et al. Reference Philippov and Spitkovsky2014). This is a net difference with 2D MHD axisymmetric models, and explains why an HD approach seems more suitable to reproduce bulk and macroscopic properties of PWNe (Olmi & Torres Reference Olmi, Del Zanna, Amato and Bucciantini2020).

• The inner nebula is reasonably well described by 2D MHD axisymmetric models, including the Crab wisps and knot (Komissarov & Lyubarsky Reference Komissarov and Lyubarsky2004; Lyutikov et al. 2016). A direct comparison shows that, if one sticks to the inner—toroidal—nebula, the description obtained with 2D models does not differ much from the structure simulated in 3D (see Figure 5). This of course is a very important result, supporting the reliability of previous models in 2D and meaning that, if only interested in the properties of the inner nebula, one can safely use less demanding 2D simulations.

• Time variability of the inner nebula: the wisp-like variability close to the shock position is reproduced correctly and is comparable with what previously found with 2D models (Camus et al. Reference Camus, Komissarov, Bucciantini and Hughes2009; Olmi et al. Reference Olmi, Bucciantini and Morlino2015). In 3D a number of variable non-axisymmetric structures is also found in the outer nebula, very similar to what observed at large scales in radio (Bietenholz et al. Reference Bietenholz, Hester, Frail and Bartel2004).

• Lacking of correct spectral description: despite the huge improvements in the description of the fluid structure, and in general in the global dynamics of the nebula, the emitting properties are not fully reproduced, both in terms of morphology and integrated spectrum. The cause may be both in the—still—too low magnetic field on average, and/or in the approximated reconstruction—a posteriori—of the properties and history of the emitting particles.

Figure 5. A selection of images for a PWN in the free-expansion phase, with maps from MHD simulations specific for the Crab nebula, elaborated from the simulations described in Del Zanna et al. Reference Del Zanna, Volpi, Amato and Bucciantini2006 (middle panel—2D) and Olmi et al. Reference Olmi, Del Zanna, Amato, Bucciantini and Mignone2016 (panel on the right—3D). In particular, from left to right: toy model of the inner structure of the nebula, showing the oblate shape of the wind termination shock (the yellow spiral marks the striped wind) and the formation of jets due to hoop stresses at the polar front of the shock; 1 keV X-ray synthetic map of the Crab nebula, reproducing its inner morphology (map in logarithmic scale and intensity scaled to the maximum); intensity map of the 3D velocity field from a 3D simulation of the Crab nebula, with the evident formation of kinking jets.

3.6 Highlights from 3D MHD models of bow shocks

In this subsection we summarize the results of 3D relativistic MHD simulations of bow shock nebulae. The discussion is mostly based on Barkov et al. (Reference Barkov, Lyutikov and Khangulyan2019) and Olmi & Bucciantini (Reference Olmi and Bucciantini2019a) for the dynamics, on Olmi & Bucciantini (Reference Ng and Romani2019b) for the emission and polarisation properties and on Olmi & Bucciantini (Reference Olmi and Bucciantini2019c) for the properties of the particle escape.

• The large scale structure of the bow shock is rather independent on the variation of the geometry and wind properties, with major deviations caused by different inclinations of the pulsar spin-axis and the direction of motion (maximum for $45^\circ$ ). The differences appear as small extrusions and blobs, in some cases resulting in periodic perturbation of the forward shock.

• The dynamics in case of an isotropic wind is very similar to the one found in 2D HD for intermediate values of the magnetisation ( $\sigma\simeq 0.1$ ).

• The dynamics in the tail is largely dominated by the level of magnetisation and by the wind anisotropy: unlike anisotropy, the lower is the magnetisation, the higher is the turbulence. Anisotropic models then are more turbulent than isotropic ones, showing strong mixing also for large magnetisation (see Figure 3 of Olmi & Bucciantini Reference Olmi and Bucciantini2019a). In general, anisotropic and low magnetized systems are fully dominated by turbulence even close to the termination shock, resulting in the complete loss of information from the injection region along the tail. On the other hand, isotropic and highly magnetized systems show a coherent structure of the magnetic field, with the injection properties still affecting the dynamics far from the shock along the tail, in very good agreement with simplified semi-analytic models (Bucciantini Reference Bucciantini2018).

• No efficient magnetic amplification from turbulence. Even in presence of strong turbulence, the magnetic field tends to reach equipartition with the turbulent kinetic energy, in general smaller then the thermal energy in the tail. The only sign of field enhancement is found close to the contact discontinuity, possibly resulting from efficient shear instability amplification.

• Low turbulence cases—emission and polarisation properties: a strong correlation between conditions at injection and surface brightness is found. The variety of observational morphologies is very wide: from bright heads to bright tails or, in some cases, bright wings. The polarisation fraction is higher in the tail for higher magnetisation.

• High turbulence cases—emission and polarisation properties: once magnetisation drops (and anisotropy grows), turbulence starts to dominate also in the appearance of surface brightness, and distinguishing between the various cases becomes hard. When turbulence increases, the polarisation fraction drops.

• The escape of particles from the bow shock is found to be an energy dependent process: the threshold for escape is set by the condition that the particle Larmor radius ( $r_L$ ) in the equipartition magnetic field ( $B_{\mathrm{eq}}\sim \mathrm{few}\,\mu$ G) is equal to the typical size of the bow-shock in the head (the stand off distance $d_0$ ). This translates in a Lorentz factor of the particles of: $\gamma\sim e B_{\mathrm{eq}} \, d_0/(m_e c^2)$ , that for typical systems^{Footnote d} corresponds to $\gamma\gtrsim \mathrm{few}\times 10^7$ , or in a particle energy $\gtrsim$ 10 TeV.

• There is a transition in the escape process: particles manage to escape more easily if injected at the frontal polar region of the pulsar wind, while the others tend to remain confined in the tail. At lower energies, particles escape only in the presence of reconnection points at the magnetopause between the shocked pulsar wind and the ISM, and this might give rise to the appearance of one-sided jet-like features. Particles show an increasingly more diffusive escape with energy: the outflow becomes more uniform (see Figure 6), but charge separation increases.

Figure 6. Comparison of two different regimes of particle escaping from a bow shock nebula from a 3D simulation: from massive and almost diffusive escape for high Lorentz factor particles ( $\gamma=10^8$ , left-blue coloured side) to directional escape along external field lines of lower energy particles (with $\gamma=10^7$ , right-orange coloured side). Data come from the simulation presented in Olmi & Bucciantini (Reference Olmi and Bucciantini2019c) and have been displayed using VisIt (Childs et al. Reference Childs2012). The inclination of the image is shown with the bottom triad, while the position of the bow shock is marked by the dashed white line.

4.1 Injection at the shock

The particle injection spectrum is thus typically modeled as a broken power-law in the particle energy E:

(15) \begin{equation} Q(E,t)=Q_0(t)\left(\frac{E}{E_b}\right)^{-p_i} \,,\end{equation}

consisting of two distinct families, characterised by different injection indices $p_i$ , at energies below or above the injection break ( $E_b$ ). The normalisation function $Q_0(t)$ is determined by the requirement that the power injected in particles is a fixed fraction of the spin-down luminosity, namely: $(1-\eta) L(t) = \int_{E_{\mathrm{min}}}^{E_{\mathrm{max}}} Q(E,t) E \,dE$ , where $\eta$ is the magnetic fraction (i.e. how much of the injection goes into magnetic energy), linked to the magnetisation (Equation (13) as: $\eta=\sigma/(\sigma+1)$ .

A rather flat injection index ( $p_{\mathrm{low}}\sim 1$ – $1.5$ ) is characteristic of the lower energy component, leading to a synchrotron spectrum with flux density $S_\nu(\nu)\propto \nu^{-\alpha_{\mathrm{low}}}$ and spectral index $\alpha_{\mathrm{low}}\sim 0$ – $0.3$ . The higher energy part indeed shows steeper spectra at injection ( $p_{\mathrm{low}}\sim 2$ – $2.7$ ), leading to a synchrotron spectrum with $\alpha_{\mathrm{high}}\sim 1$ – $1.2$ or, as more commonly reported, a photon index $\Gamma=1+\alpha \sim 2$ – $2.2$ . The recent PeV observations by LHAASO (Cao et al. Reference Cao2021a, b) confirm, as originally suggested by Bucciantini et al. (Reference Bucciantini, Arons and Amato2011) and recently investigated by de Oña Wilhelmi et al. (Reference de Oña Wilhelmi, López-Coto, Amato and Aharonian2022), that high energy particles can be accelerated up to the pulsar voltage (see also Khangulyan, Arakawa, & Aharonian Reference Khangoulian and Bogovalov2020).

The possibility that particles responsible for the radio, optical and X-ray emissions belong to separate families, accelerated via different mechanisms, is supported also by the multi-wavelength variability observed in the Crab nebula. The arc-like bright structures named wisps, that appear very close to the TS location and move outwards with mild relativistic velocity, have been observed at multi-wavelengths (Hester et al. 2002; Bietenholz, Frail, & Hester Reference Bietenholz, Frail and Hester2001; Bietenholz et al. Reference Bietenholz, Hester, Frail and Bartel2004), and shown to neither be spatially coincident nor characterized by the same velocity (Schweizer et al. Reference Sartore, Ripamonti, Treves and Turolla2013). In the MHD framework, which has provided a very good description of the nebular dynamics, wisps trace the structure of the underlying plasma—the magnetic field in particular—and as such they can only be non-coincident if particles with different energies are produced at different locations of the shock. This likely means that acceleration processes act differently in different sectors of the shock (Olmi et al. Reference Olmi, Bucciantini and Morlino2015).

The acceleration mechanisms proposed so far are mainly three: (i) diffusive shock acceleration, or Fermi-I like processes; (ii) diffuse acceleration due to stochastic magnetic reconnection, or Fermi-II, in MHD turbulence; (iii) acceleration conveyed by driven magnetic reconnection. Actually a fourth mechanism has been invoked: resonant absorption of ion cyclotron waves, which however requires the presence of ions in the wind, still not confirmed (nor excluded).

Diffusive shock acceleration requires a very low magnetisation to be effective ( $\sigma \lesssim 10^{-3}$ , Sironi et al. Reference Sironi, Keshet and Lemoine2015), a condition that can only be sustained at the equatorial sector of the shock, where the striped wind ensures a huge dissipation of the field, or very close to the polar axis, where the field naturally vanishes. The power law index at injection of optical/X-ray emitting particles is compatible with what predicted by Fermi-I acceleration, and MHD models for the X-ray wisps are also in agreement with a scenario in which that particles are injected mainly at the equatorial front of the oblique termination shock. On the other hand, driven magnetic reconnection requires a much higher magnetisation ( $\sigma \gtrsim 30$ ) and a large pair multiplicity ( $\kappa\gtrsim 10^8$ ), difficult to account within present models of pulsars magnetospheres (Sironi & Spitkovsky Reference Siegel and Ciolfi2011; Timokhin & Harding Reference Timokhin and Harding2019). The power law index of the radio emitting particles at injection is on the other hand compatible with both reconnection and Fermi-II acceleration, while no particular information arises from wisps models in this case. A possible conclusion is that radio particles are accelerated at higher latitudes along the shock, out of the equatorial sector, or via a more distributed acceleration in the nebula (Olmi et al. Reference Olmi, Bucciantini and Morlino2015; Lyutikov et al. Reference Lyubarsky2019).

4.2 Radiation mechanisms

Once injected in the nebula, the particle energy E evolves according to the following equation:

(16) \begin{equation} \frac{\partial E}{\partial t} = - \frac{E}{R}\frac{\partial R}{\partial t} - c_2 E^2 \left( \frac{B^2}{8\pi} + U_{\mathrm{rad}}\right)\,,\end{equation}

where $c_2=4/3 \,\sigma_{\mathrm{th}}/(m_e^2 c^3)$ , with $\sigma_{\mathrm{th}}$ the Thomson cross section, and $U_{\mathrm{rad}}$ is the energy density in radiation. The first term on the right side of the equation represents the energy variation due to adiabatic expansion or contraction (depending on the evolutionary phase), while the second one describes radiation losses (synchrotron and IC). From this equation one can easily find that the cooling energy $E_c$ is given in general by

(17) \begin{equation} E_c = \frac{1}{R}\frac{\partial R}{\partial t}\frac{1}{c_2 \left( B^2/8\pi + U_{\mathrm{rad}}\right)}\,.\end{equation}

Particles with $E<E_c$ are then dominated by adiabatic processes, losses or gains depending on the expansion or compression of the PWN. All particles with energy above $E_c$ are instead dominated by radiation losses. In the free-expansion phase $E_c$ determines the cooling spectral break, that it is commonly found to move to higher energies with time e.g. Pacini & Salvati Reference Pacini and Salvati1973). During reverberation instead $E_c$ separates those particles that gain energy due to compression ( $E<E_c$ ) from those that are still loosing energy due to radiation losses ( $E>E_c$ ).

As already mentioned, the spectral energy distribution of a PWN is fully non-thermal, with the primary emitting mechanism, responsible for emission from radio up to few hundreds of MeV, being synchrotron radiation. Higher energies are produced with the second emitting process characteristic of those systems: IC scattering between local photons and the same leptons responsible for the synchrotron emission. An example of a full spectral energy distribution coming from a one-zone modelling can be found, in the specific case of the Crab nebula, in Figure 7.

Figure 7. Non-thermal spectral energy distribution of the Crab nebula, elaborated from the original plot in Bucciantini et al. (Reference Bucciantini, Arons and Amato2011), computed with a one-zone radiative model. Details of the assumed parameters and origin of data points can be found in the reference paper (see Figure 1 and its caption). The different coloured areas highlight the emission at various energy bands. From left to right: light-red for radio and IR; light green for optical and UV; light blue for X-rays (considering 0.1–100 keV); light yellow for the low energy gamma rays (up to the synchrotron limit $\sim$ 250 MeV) and darker yellow for the high energy gamma-rays (fully due to IC). The blue line is for the synchrotron component, the red one for the IC component, to which major contributions come from self-synchrotron Compton (in cyan) and scattering with CMB photons (in yellow).

The main contribution to IC in general comes from the interaction with the photons of the cosmic microwave background (CMB), with minor contributions from the interstellar radiation field and the synchrotron photons from the PWN itself. An exception in this respect is the Crab nebula: due to its young age and intense magnetic field ( $\sim$ 150–200 $\mu$ G), its IC spectrum has in fact a consistent contribution from self-synchrotron radiation. The Crab nebula is moreover the only known case where the maximum energy of the accelerated particles is limited by radiation losses in the Galaxy, and is smaller than the maximum energy inferred by the available potential from the pulsar.

Once the injection spectrum has been defined (Equation (15), the PWN luminosity can be computed as:

(18) \begin{equation} L_{\textrm{PWN}}(\nu)=4\pi \int_{V_{\textrm{PWN}}} \left\{j_\nu^{\mathrm{SYNC}}(\nu) + j_\nu^{\mathrm{IC}}(\nu) \right\} dV\,,\end{equation}

where $V_{\textrm{PWN}}$ is the PWN volume, while $j_\nu^{\,(i)}$ is the emissivity, obtained integrating over the particle distribution function either for synchrotron or IC radiation as $j_\nu^{\,(i)}=\int_{E_{\mathrm{min}}}^{E_{\mathrm{max}}} \tilde{Q}(E,t) \, \mathcal{P}_\nu^{\,(i)}(E,\nu) \,dE$ , where $\tilde{Q}(E,t)$ is the evolved particle spectrum in the nebula, from the injection one defined in Equation (15) (see e.g. Bucciantini et al. Reference Bucciantini, Arons and Amato2011). When computing the emissivity, in the 3D case, one should of course take also into account the spatial dependence due to orientation of the line of sight. General expressions for the synchrotron and IC power $\mathcal{P}_\nu(E,\nu)$ can be found in many textbooks, e.g. Rybicki & Lightman (1979).

5.1 Young systems

To date we have identified less than 20 sources as PWNe in their free-expansion phase. They constitute a large part of the catalog of the PWNe detected mainly at X-rays and not associated with a fast moving pulsar ( $\sim$ 40 sources), reported in Table A.1. From the morphological point of view, synchrotron dominated systems are characterised by a larger extension at lower energies than at higher ones (see e.g. Figure 9), with sort of a spherical/elliptical shape. X-rays highlight the inner nebula, revealing the presence of a jet-torus structure (see the zoom in of Figure 9), believed to be a rather common feature in young PWNe, first discovered in the Crab and successively identified in another bunch of systems.

Figure 9. Left panel: Composite optical image of the Crab nebula, Credits: ESO. Right panel: Combined optical (in red—from Hubble) and X-ray (in blue—from Chandra) images of the Crab nebula, Credits: Optical—NASA/HST/ASU/J. Hester et al. ; X-Ray—NASA/CXC/ASU/J. Hester et al.

Another common feature of the observed torii is the appearance of enhanced brightness at one side (effect of the relativistic Doppler boosting of the emitting particles moving towards the observer) and the presence of variable, both in brightness and position, arc-like (the wisps, Scargle Reference Rybicki, Lightman, Matheson and van Leeuwen1969; Bietenholz, Frail, & Hankins Reference Bietenholz, Frail and Hankins1991; Bietenholz et al. Reference Bietenholz, Frail and Hester2001; Helfand, Gotthelf, & Halpern 2001; Pavlov et al. 2001; Bietenholz et al. Reference Bietenholz, Hester, Frail and Bartel2004) or point like (knots, Lou 1998) structures marking the high variability of the inner nebula, where (most of?) the particles are accelerated. As mentioned in section 3.3, the discovery of the complex inner structure of the Crab nebula was what prompted the move from 1D models to 2D MHD simulations.

As discussed previously (Section 4), young PWNe are characterized by extremely broad band spectra, extending from radio to gamma-rays; with the advent of LHAASO, now the Crab spectrum has been further extended above PeV energies (Cao et al. Reference Cao2021a).

5.2 Middle aged—reverberating systems

The reverberation phase, characterized by the interaction with the SNR reverse shock, is hardly identifiable. At present we only know a handful of systems showing clear evidence of being in that stage, among which Vela X (Blondin et al. Reference Blondin, Chevalier and Frierson2001), the Boomerang nebula (Kothes, Reich, & Uyanıker 2006) and the Snail in G327.1-1.1 (Temim et al. Reference Tanaka and Takahara2009, Reference Temim, Slane, Kolb, Blondin, Hughes and Bucciantini2015), shown in Figure 10.

Figure 10. Composite IR (for the stellar field), radio (in red colour) and X-ray (in blue) image of the PWN G327.1-1.1, one of the very few systems in clear interaction with the SNR reverse shock. Credits: X-Ray—NASA/CXC/SAO/T. Temim et al. and ESA/XMM-Newton; Radio: SIFA/MOST and CSIRO/ATNF/ATCA; IR: UMass/IPAC-Caltech/NASA/NSF/2MASS.

Independently of the pulsar speed, middle aged systems are expected to show large asymmetries. The interaction with the SNR reverse shock is not expected to happen spherically, as simplified one-zone models are forced to assume (see Section 3.1). The compactness of the PWN contact discontinuity itself is partially destroyed by R-T like instabilities (Section 2.5) well before the onset of reverberation. Of course, in case of a high proper motion of the star, the asymmetry is even larger. The asymmetry introduced in the PWN morphology in this stage is then expected to be a common feature of almost all middle-aged systems, as well as of old ones, if they have not become bow shock nebulae. The original PWN bubble might be fragmented, with radio (and gamma-ray) separated bubbles surrounding the remaining nebula and expanding under adiabatic forces. The background of the PWN can be then very noisy, making it difficult to be identified at radio or gamma-rays. On the other hand, X-ray emission might simply be too faint.

Estimating the level of fragmentation and mixing of the PWN after reverberation can be quite complex: HD simulations typically do not converge, because in the HD regime the fastest growing scale of R-T like instabilities is set purely by numerical viscosity at the grid scale (Blondin et al. Reference Blondin, Chevalier and Frierson2001; Bucciantini et al. Reference Bucciantini, del Zanna, Amato and Volpi2005), while MHD simulations have not been performed for old systems; observationally the mixing has only been estimated for the Snail (Ma et al. Reference Lyutikov, Komissarov and Porth2016), suggesting a pulsar wind filling factor of order of 50%. For a discussion of the possible impact of mixing on thin-shell modelling see Bucciantini et al. (Reference Bucciantini, Arons and Amato2011).

5.3 Bow shock nebulae

Through the combination of radio, H $_\alpha$ (in case of a partially ionized ambient medium) and, especially, X-ray observations, nowadays we have identified 25 fast moving pulsars with an associated bow shock nebula. They are listed in Table A.2 and are marked with blue coloured circles in Figure 8.

Few bow shock nebulae show an elongated X-ray tail, in some cases associated with an even more extended radio tail. Only a part of them shows a spectral variation along the X-ray tail, in particular a softening indicating synchrotron cooling (e.g. the Mouse and the Lighthouse, Kargaltsev et al. Reference Hare2017). To date no TeV emission has been detected from bow shocks directly, while UV has only been detected in correspondence with H $_\alpha$ emission, probably coming from the heated shocked ISM. The bow shock head appears not to have a standard morphology, with even drastic variations from one object to another. As shown by 3D MHD models (Barkov et al. Reference Barkov, Lyutikov and Khangulyan2019; Olmi & Bucciantini Reference Olmi and Bucciantini2019a,b), these variations can be ascribed to intrinsic differences in the geometry of the pulsar magnetosphere, in the orientation of the pulsar spin-axis with respect to the pulsar direction of motion, and the orientation with respect to the observer’s line of sight.

In some cases, the structure of the bow shock appears to be modified in the so called head-and-shoulder shape: the bow shock shows an evident widening with distance from the pulsar, with possibly a periodic structure, as the famous example of the Guitar nebula (Chatterjee & Cordes Reference Chatterjee and Cordes2004; van Kerkwijk & Ingle Reference van der Swaluw, Achterberg, Gallant and Tóth2008, see e.g. Figure 11). This is believed to be the sign of the mass loading of ambient neutral atoms into the bow shock through the shocked ISM (Morlino et al. Reference Mori, Burrows, Hester, Pavlov, Shibata and Tsunemi2015); those atoms then interact with the pulsar wind and modify its dynamics (Bucciantini & Bandiera Reference Bucciantini and Bandiera2001; Bucciantini Reference Bucciantini2002). This effect has been proved through numerical simulations by Olmi et al. (Reference Olmi and Bucciantini2018), showing that the lateral expansion of the bow shock tail is a function of the pulsar Mach number only, namely it increases with the Mach number as the effect of the augmented ram pressure exerted by the ISM on the bow shock nebula contact discontinuity.

Figure 11. Composite optical (in blue colour) and X-ray (in magenta) image of the Guitar BSPWN, generated by the fast moving pulsar J2225+6535, and its extended X-ray misaligned tail. The Guitar nebula is also one of the systems characterized by the modification of the tail structure in the so called head-and-shoulder shape, possibly indicating mass loading from the ambient medium into the tail. Credits: X-Ray—NASA/CXC/UMass/S.Johnson et al; Optical: NASA/STScI & Palomar Observatory 5-m Hale Telescope.

In recent years bow shock nebulae have gained renewed interest thanks to the detection of collimated, extended and generally highly misaligned (with respect to the direction of motion), jet-like features, only visible at X-rays (observed in the Chandra band: 0.5–8 keV), usually referred to as misaligned tails (De Luca et al. 2011; Pavan et al. Reference Paredes-Fortuny, Bosch-Ramon, Perucho and Ribó2014; Klingler et al. Reference Klepser2016; de Vries & Romani Reference de Vries and Romani2020; Wang Reference Vurm and Metzger2021; de Vries et al. Reference de Vries2022; de Vries & Romani Reference de Vries and Romani2022). These structures were already observed few years ago surrounding a couple of systems (e.g. the Lighthouse nebula), but with the recently increased number of detection they now seem a rather common feature of these evolved nebulae. At present a generally accepted interpretation (Bandiera Reference Bandiera2008) is that they are produced by high energy particles (close to the maximum limit of the potential drop) leaking from the bow shock nebula, and then producing emission via synchrotron radiation in the local magnetic field. The observed asymmetry appears to be related to the mutual inclination of the spin-axis of the pulsar, its magnetic field, the pulsar speed and the direction of the magnetic field lines in the ambient medium (Olmi & Bucciantini Reference Olmi and Bucciantini2019c). Indeed an open question is how to amplify the magnetic field from the ISM value to the order of few tens of $\mu$ G required to produce the observed emission.

Evolved pulsars are also associated with the formation of TeV halos (Abeysekara et al. Reference Abeysekara2017), that have been again interpreted as high-energy particles escaping from the bow shock nebula, and then diffusing (with some suppression) in the ambient medium to form bright gamma-ray bubbles, much more extended than the original bow shock (see e.g. a sketch for the case of Geminga in Figure 12).

Figure 12. Sketch roughly comparing the size of the TeV halo around the Geminga pulsar wind nebula (from HAWC measures at 100 TeV) and that of the X-ray pulsar wind nebula (image adapted from the original composite picture at X-rays—Chandra—and IR—Spitzer). Credits for the PWN map: X-Ray—NASA/CXC/PSU/B. Posselt et al; IR: NASA/JPL-Caltech.

6.1 Observational prospects

The actual population of PWNe is expected to increase largely in the next future, especially thanks to very high energy ( $>$ 100 GeV) observations. The actual Galactic plane survey from the H.E.S.S. telescope (Abdalla et al. Reference Abdalla2018a) found that more than half of the 24 extended sources detected are identifiable with PWNe (14, mainly thanks to their multi-wavelength counterpart). In the Fermi-LAT 3FGL catalog (Abdo et al. Reference Abdo2013) the unidentified sources are $\sim$ 20% of the total and it is plausible that most (all?) of them are actually PWNe with no direct association with a known pulsar. As we already said, thanks to their longer lifetime as gamma-ray emitters, middle aged PWNe will likely dominate the very high energy sky, possibly representing up to $\sim$ 60% of the Galactic sources, and a huge number of new detection ( $\sim$ 200) may be expected with the upcoming Cherenkov Telescope Array, thanks to its unprecedented angular and energy resolution (Remy et al. Reference Remy, Tibaldo, Acero, Fiori, Knödlseder, Olmi and Sharma2022; Fiori et al. Reference Fiori, Olmi, Amato, Bandiera, Bucciantini, Zampieri and Burtovoi2022). A very challenging problem will be that of the source confusion in the crowded Galactic plane, that will reduce the number of identified sources with respect to theoretical estimate (Mestre et al. Reference Martin, Tibaldo, Marcowith and Abdollahi2022). Moreover an additional source of confusion might be that of TeV halos associated with evolved pulsars. Given their extended and weak emission, they are difficult to identify, and likely constitute a background noise for other sources. However the number of expected halos in the Galaxy is still matter of debate, ranging from few hundreds to a few, depending on their physical interpretation (Sudoh et al. Reference Soker and Gilkis2019; Giacinti et al. Reference Giacinti, Mitchell, López-Coto, Joshi, Parsons and Hinton2020; Martin et al. Reference Margalit, Metzger, Thompson, Nicholl and Sukhbold2022).

A very exciting recent result is the detection of PeV photons coming from 12 sources in the Galaxy, plus the Crab nebula, by (Cao et al. Reference Cao2021a,b). Unfortunately, the limited angular resolution of LHAASO does not permit to identify the exact location of the source (and origin of these energetic photons), except for the case of the Crab. One or more pulsars can be found in the same region covered by the PSF of the instrument for all the 12 sources. Out of these, 11 result to be theoretically compatible with being powered by a pulsar (de O.a Wilhelmi et al. 2022), meaning that pulsars and their nebulae might also be the most numerous class of extremely high energy emitters in the Galaxy.

PeV data are also fundamental to either confirm or exclude the presence of an hadronic component in the pulsar wind (Atoyan & Aharonian Reference Atoyan and Aharonian1996; Bednarek & Protheroe Reference Bednarek and Protheroe1997; Bednarek & Bartosik Reference Bednarek and Bartosik2003; Amato et al. Reference Amato, Guetta and Blasi2003). If present, hadrons might show up only at the very high energies, where the leptonic emission from IC drastically falls due to Klein-Nishina suppression (Amato & Olmi Reference Amato and Olmi2021), thus the next generation of imaging atmospheric Cherenkov telescopes will have a crucial role in answering this question.

The recently launched IXPE satellite (Weisskopf et al. Reference Weisskopf, Silver, Kestenbaum, Long and Novick2022), operating in the range 2–8 keV, would enable us to sample and image, for the first time, the magnetic field structure (not just its strength) and the level of turbulence in a handful of PWNe, resolving the magnetic field pattern in the central region of the torus-arcs characteristic of young objects, and possible in the jets. The three main targets of the mission for space resolved polarimetry are: the Crab nebula, Vela and MSH 15-52. Other PWNe will be observed in the second and third year of operations. Prior to IXPE the Crab nebula was the only object to have been observed with a detected average X-ray polarisation of 19%. IXPE will open a new observational window into the way we understand and characterize these objects. Preliminary results are coming in these same days (e.g Bucciantini et al. Reference Bucciantini2022). By the end of the mission we expect to have a much better understanding of the dynamical conditions in these relativistic accelerators.

With the approaching end of operations of the Chandra telescope, a very important instrument for future observations of PWNe will be the LYNX X-ray observatory,^{Footnote f} that thanks to its improved sensitivity and field of view promises to open new windows of opportunity to investigate the structure and details of X-ray sources (an example is the possibility to finally detect the compact object in SNR 1987A, see (e.g. Greco et al. 2021).

6.2 Modelling prospects

Despite the astonishing progresses made in the last two decades, there are still many important open questions in our understanding of the physical processes operating in PWNe. If a part of the historical open problems seem nowadays to be solved, as the case of the long standing sigma-paradox, some are still not properly answered, among which:

• What are the physical mechanisms responsible for particle acceleration at the different energies? All the proposed mechanisms have strengths and weaknesses and require precise—and very diverse—physical properties of the nebular plasma to be viable.

• What is the origin of the gamma-ray flares in the Crab nebula? And is the Crab the unique source producing flares? Many possibilities had been proposed and investigated, most of them requiring a mG magnetic field at the emitting region (this is the case for powerful reconnection events in $\sigma\gg 1$ regions), much larger than that estimated from PeV emission. A very detailed discussion of this point, and of the proposed models, has been recently reported in Amato & Olmi (Reference Amato and Olmi2021).

• Are hadrons present in the pulsar wind? Despite being certainly a minority by number, hadrons could even be—if present—energetically dominant in the wind, changing completely our understanding of its properties.

In the latest years new questions have been added, especially thanks to observations of evolved systems that have revealed a number of unexpected features:

• How and with what efficiency particles can escape from evolved PWN?

• How escaped particles produce misaligned tails or extended halos?

• Do pulsars accelerate particles at the theoretical limit of the maximum potential drop?

• How does the interaction with the reverse shock of the SNR modify the emitting properties and morphology of evolved PWNe?

• Which ingredients of current modelling need to be modified/improved to interpret the upcoming large amount of gamma-ray observations and to manage source confusion? Is there a way to theoretically model the possible different morphologies and spectra of evolved PWNe that can help in their identification, especially in the lack of a multi-wavelength counterpart?

Recent modelling efforts go in the direction of trying to answer these questions. How particles escape from evolved systems seems now to be clarified, as discussed in Olmi & Bucciantini (Reference Olmi and Bucciantini2019c). The efficiency of the escape process is strictly linked to the particle energy, with only the most energetic particles, close to the maximum energy achievable, able to escape in large fractions. The possibility for them to be revealed as extended and diffuse, or asymmetric and thin, structures is also partially accounted for with numerical models, that predict a variety of escaping processes depending on the properties of the system and those of the surrounding medium.

We still lack an understanding of is what happens once particles have been injected in the ambient medium. What process causes the reduction of the Galactic diffusion length in the vicinity of pulsars (if this is the case), to produce TeV halos with the observed size through diffusion? And also, what amplifies the magnetic field to the value needed for the observed synchrotron emission from misaligned X-ray tails, a factor of 2–10 larger than that expected in the ISM? To answer these questions new dedicated modelling of particle propagation must be investigated, considering diffusion properties, development of self-turbulence and the onset of instabilities able to modify the properties of the ISM, possibly using a hybrid approach between pure MHD and Particle In Cell (PIC) techniques, that correctly accounts for the evolution of particles in the plasma.

A better understanding of the composition, and of the physical properties, of the pulsar wind passes through a correct interpretation of high energy gamma-ray data. In particular, the firm identification of the observed PeVatrons will shed light on the possibility for pulsars to efficiently accelerate particles very close to the theoretical limit. A refined modelling of the spectral properties through evolutionary phases will also be extremely important.

A very challenging point will be to develop models that help in disentangling the PWNe contribution to the gamma-ray Galactic emission from that of other sources. Part of this requires to understand how to model evolved systems, both in terms of emission and morphology, with fast and light enough approaches to reproduce a large sample of the expected population, but enough refined to be reliable. In Bandiera et al. (Reference Bandiera, Bucciantini, Martín, Olmi and Torres2023) a first step to better include the reverberation phase in the description of the PWNe evolution has been made, based on the modification of the standard one-zone description. Nevertheless these results, despite being much more accurate than previous models, still are far for being definitive. A 3D MHD modelling of reverberation is necessary to understand which role the third spatial dimension, and the structure of the magnetic field, play in shaping the PWN during its interaction with the SNR.

Moreover one-zone models by construction cannot account for the formation of asymmetric systems, that we expect will constitute the largest part of middle aged to evolved PWNe. Then, models must be somehow generalised to account for different geometries, but how without running expensive 3D—or even 2D—models is absolutely unclear, with only few preliminary studies presented for the moment (e.g. Olmi & Torres Reference Olmi, Del Zanna, Amato and Bucciantini2020). This will also be the flip of the coin for the interpretation of gamma-ray data.