A&A 380, 309-317 (2001)
DOI: 10.1051/0004-6361:20011437

Pulsar wind nebulae in supernova remnants

Spherically symmetric hydrodynamical simulations

E. van der Swaluw ^1,2 - A. Achterberg ¹ - Y. A. Gallant^1,3,4 - G. Tóth ⁵

1 - Astronomical Institute, Utrecht University, PO Box 80000, 3508 TA Utrecht, The Netherlands
2 - Dublin Institute for Advanced Studies, 5 Merrion Square, Dublin 2, Republic of Ireland
3 - Osservatorio Astrofisico di Arcetri, Largo E. Fermi 5, 50125 Firenze, Italy
4 - Service d'Astrophysique, CEA Saclay, 91191 Gif-sur-Yvette Cedex, France
5 - Department of Atomic Physics, Pázmány Péter sétány 1, 1117 Budapest, Hungary

Received 19 December 2000 / Accepted 27 September 2001

Abstract
A spherically symmetric model is presented for the interaction of a pulsar wind with the associated supernova remnant. This results in a pulsar wind nebula whose evolution is coupled to the evolution of the surrounding supernova remnant. This evolution can be divided in three stages. The first stage is characterised by a supersonic expansion of the pulsar wind nebula into the freely expanding ejecta of the progenitor star. In the next stage the pulsar wind nebula is not steady; the pulsar wind nebula oscillates between contraction and expansion due to interaction with the reverse shock of the supernova remnant: reverberations which propagate forward and backward in the remnant. After the reverberations of the reverse shock have almost completely vanished and the supernova remnant has relaxed to a Sedov solution, the expansion of the pulsar wind nebula proceeds subsonically. In this paper we present results from hydrodynamical simulations of a pulsar wind nebula through all these stages in its evolution. The simulations were carried out with the Versatile Advection Code.

Key words: ISM: supernova remnants - stars: pulsars: general - hydrodynamics - shock waves

1 Introduction

The explosion of a massive star at the end of its life as a supernova releases an amount of energy roughly equal to 10⁵³ erg. Roughly 99% of this energy is radiated away in the form of neutrinos as a result of the deleptonization of the $\sim$1 $M_{\odot}$ stellar core as it collapses into a neutron star. The remaining (mechanical) energy of about 10⁵¹ erg is contained in the strong shock propagating through the stellar mantle, and ultimately drives the expansion of the supernova remnant (SNR). In those cases where a rapidly rotating neutron star (pulsar) remains as a "fossil'' of the exploded star, a pulsar wind, driven by the spindown luminosity of the pulsar, can be formed. The precise magnetospheric physics leading to such a pulsar wind is not fully understood, but it is believed that a major fraction of the spin-down luminosity of the pulsar is converted into the mechanical luminosity of such a wind.

The total rotational energy released by a Crab-like pulsar over its lifetime is of order 10⁴⁹-10⁵⁰ erg. This is much less than the mechanical energy of $\sim$10⁵¹ erg driving the expansion of the SNR. Therefore, the dynamical influence of the pulsar wind on the global evolution of the supernova remnant itself will be small. From an observational point of view, however, the presence of a pulsar wind can lead to a plerionic supernova remnant, where the emission at radio wavelengths shows an extended, flat-spectrum central source associated with a Pulsar Wind Nebula (PWN). The best-known example of such a system is the Crab Nebula, and about a half-dozen other PWNs are known unambiguously around pulsars from radio surveys (e.g. Frail & Scharringhausen 1997). These surveys suggest that only young pulsars with a high spindown luminosity produce observable PWNs at radio wavelengths. At other than radio wavelengths, in particular X-rays, there are about ten detections of PWN around pulsars both in our own galaxy and in the large Magellanic cloud (LMC) (e.g. Helfand 1998; Table 1 of Chevalier 2000). The expansion of an isolated SNR into the general interstellar medium (ISM) can be divided in four different stages (Woltjer 1972; see also Cioffi 1990 for a review): the free expansion stage, the Sedov-Taylor stage, the pressure-driven snowplow stage and the momentum-conserving stage. In the models presented here we will only focus on the evolution of a pulsar wind nebula (PWN) during the first two stages of the SNR: the free expansion stage and the Sedov-Taylor stage. We will assume that the pulsar is stationary at the center of the remnant, excluding such cases as CTB80 (e.g. Strom 1987; Hester & Kulkarni 1988), PSR1643-43 and PSR1706-44 (Frail et al. 1994), where the pulsar position is significantly excentric with respect to the SNR, presumably due to a large kick velocity of the pulsar incurred at its birth, assuming of course that SNR and pulsar are associated and we are not dealing with a chance alignment of unrelated sources. The case of a pulsar moving through the remnant with a significant velocity will be treated in a later paper.

In this paper we compare (approximate) analytical expressions for the expansion of a PWN in a supernova remnant with hydrodynamical simulations carried out with the Versatile Advection Code (VAC). We confirm earlier analytical results (Reynolds & Chevalier 1984; Chevalier & Fransson 1992) which state that the PWN is expanding supersonically when it is moving through the freely expanding ejecta of the SNR. Due to deceleration of the expanding SNR ejecta by the interstellar medium (ISM), a reverse shock propagates back to the center of the SNR (e.g. McKee 1974; Cioffi et al. 1988). Due to the presence of reverberations of the reverse shock in the SNR, the expansion of the PWN goes through an unsteady phase when this reverse shock hits the edge of the PWN. After these reverberations have decayed, the expansion of the PWN through the ejecta of the SNR progenitor star continues subsonically with the PWN almost in pressure equilibrium with the interior of the SNR.

This paper is organised as follows. In Sects. 2 and 3 we discuss the aforementioned two stages of the PWN/SNR system. In Sect. 4 the hydrodynamical simulations will be presented and compared with the analytical expressions from Sects. 2 and 3.

2 Pulsar Wind Nebula in a freely expanding Supernova Remnant

In the early stage of the evolution of a PWN, the SNR consists mostly of the stellar ejecta expanding freely into the interstellar medium. The PWN expands into these ejecta. The sound velocity in the interior of the SNR is much smaller than the expansion velocity of the PWN. The supersonic expansion of the PWN results in a shock propagating into the ejecta (see Fig. 1).

An analytical equation for the radius of this shock can be derived for a constant spindown luminosity. Using this solution, the assumption of supersonic expansion will be checked a posteriori. For simplicity we assume that the ejecta have a uniform density,

$$% \rho_{\rm ej}(t) = \frac{3M_{\rm ej}}{4 \pi R_{\rm ej}^{3}},$$ (1)

and a linear velocity profile as a function of radius,

$$% V_{\rm ej}(r) = \frac{r}{t} = V_{0} \: \left( \frac{r}{R_{\rm ej}} \right),$$ (2)

with $R_{\rm ej}=V_0t$ the radius of the front of the ejecta. The value of V₀ is determined by the requirement that the kinetic energy of the ejecta equal the total mechanical energy E₀ of the SNR:

$$% E_0 = \mbox{{$\displaystyle{\frac{1}{ 2}}$ }} \: \rho_{\rm ... ...ox{{$\displaystyle{\frac{3}{10}}$ }} \: M_{\rm ej} V_{0}^{2} .$$ (3)

This yields:

$$% V_0=\sqrt{\frac{10}{3} \: \frac{E_0}{M_{\rm ej}}}\cdot$$ (4)

We assume that the stellar ejecta swept up by the strong shock which bounds the PWN collect in a thin shell, and that this material moves with the post-shock velocity. Neglecting the contribution of the thermal energy we can write the total (kinetic) energy of this shell, $E_{\rm shell}$, as:

$$% E_{\rm shell}(t)= \mbox{{$\displaystyle{\frac{1}{2}}$ }} \:... ...(t) + \mbox{$\frac{1}{4}$ } \frac{R_{\rm pwn}}{t} \right)^2 ,$$ (5)

where

$$% M_{\rm sw}(t) \equiv M_{\rm ej} \: \left({R_{\rm pwn}\over V_0 t}\right)^3$$ (6)

is the ejecta mass swept up by the pulsar wind nebula. In deriving the post-shock velocity, we assumed that the ejecta behave as an ideal non-relativistic gas with adiabatic heat ratio $\gamma_{\rm ej} = 5/3$ and used the Rankine-Hugoniot jump conditions for a strong shock.

The interior of the PWN is dominated by thermal energy. The sound speed in a realistic PWN is close to the speed of light c, while the expansion velocity is much less than c. Perturbations in the pressure will be smoothed out rapidly, on a sound crossing time $t_{\rm s} \sim R_{\rm pwn}/c$, much less than the expansion time scale $t_{\rm exp} \sim R_{\rm pwn}/\dot{R}_{\rm pwn}$. Therefore, we can assume a nearly uniform pressure $P_{\rm pwn}$ in the PWN. The internal energy of the PWN then equals

$$% E_{\rm pwn}={4\pi\over 3(\gamma_{\rm pwn} -1)} \:R_{\rm pwn}^3 \: P_{\rm pwn}.$$ (7)

Here we take $\gamma_{\rm pwn} =4/3$ because the pulsar wind nebula material is relativistically hot. The pressure of the interior of the PWN is a few times ($\sim$4 when $\gamma_{\rm ej} = 5/3$) the pressure in the shocked ejecta just downstream of the outer shock of the pulsar wind nebula (Chevalier 1984, Fig. 4). In the Appendix, we give a derivation of this result using the thin shell approximation for the motion of the shocked ejecta, assuming a power-law expansion: $R_{\rm pwn} \propto t^{\alpha}$. There we show that the ratio of the shock radius $R_{\rm s}$ and the radius $R_{\rm c}$ of the contact discontinuity separating the pulsar wind material from the shocked ejecta satisfies (Eq. (A.2)) $R_{\rm s}/R_{\rm c} = (2 \alpha + \gamma_{\rm ej} - 1)/(\alpha (\gamma_{\rm ej} + 1))$(=24/23 when $\gamma_{\rm ej} = 5/3$ and $\alpha = 6/5$, see below), and that the pressure of the pulsar wind bubble equals (Eq. (A.6)):

$$% P_{\rm pwn}(t) = \frac{(\alpha - 1)(4 \alpha - 3)}{3} \: \rho_{\rm e}(t) \: \left( \frac{R_{\rm pwn}}{t} \right)^{2} .$$ (8)

Combining these relations yields:

$$% E_{\rm pwn} = \frac{(4 \alpha - 3)(\alpha - 1)}{3(\gamma_{\... ...1)} \: M_{\rm sw} \: \left( \frac{R_{\rm pwn}}{t} \right)^{2}.$$ (9)

Energy conservation for the PWN system reads:

$$% E_{\rm shell}(t) + E_{\rm pwn}(t)= E_{\rm init}(t) + L_0 t .$$ (10)

Here $E_{\rm init}(t)$ is the kinetic energy which the swept-up ejecta would have if they were freely expanding. This quantity can be obtained by integrating the kinetic energy density of ejecta in a sphere with radius $r < R_{\rm pwn}$ if there was no PWN. This yields:

$$% E_{\rm init}(t)=E_0 \; \left(\frac{R_{\rm pwn}}{R_{\rm ej}} \right)^5.$$ (11)

After some algebra using the Eqs. (1)-(11) one can obtain a power-law solution for the radius of the pulsar wind bubble:

$$% R_{\rm pwn}(t)=C\left({L_0t\over E_0}\right)^{1/5}V_0t\propto t^{6/5},$$ (12)

where C is a numerical constant of order unity:

$$% C =\left({6 \over{\displaystyle 15 (\gamma_{\rm pwn}-1)}}+\frac{289}{240} \right)^{-1/5} \simeq 0.839$$ (13)

assuming $\gamma_{\rm ej} = \frac{5}{3}$ and $\gamma_{\rm pwn} = \frac{4}{3}$. Chevalier (1977) and Reynolds & Chevalier (1984) already obtained this $R(t)\propto t^{6/5}$ expansion law. It can easily be checked that the expansion velocity obtained in this manner is indeed much larger than the sound velocity in the freely expanding supernova remnant.

Figure 1: Schematic representation of PWN in a freely expanding SNR. There are a total of four shocks and two contact discontinuities. From left to right one can see: the pulsar wind termination shock $R_{\rm ts}$ (dashed line), the first contact discontinuity $R_{\rm cd1}$ (dotted line) separating shocked pulsar wind material from shocked ejecta, the PWN shock $R_{\rm pwn}$ (solid line) bounding the PWN. For the SNR we have a reverse shock $R_{\rm rs}$ (dashed line), the second contact discontinuity $R_{\rm cd2}$ (dotted line) separating shocked ejecta from shocked ISM, and the SNR shock $R_{\rm snr}$ (solid line) which is the outer boundary of the PWN/SNR system.

Open with DEXTER

3 Pulsar Wind Nebula in a Sedov-Taylor remnant

3.1 Pulsar Wind Nebula expansion for a constant wind luminosity

Towards the end of the free expansion stage a reverse shock is driven deep into the interior of the SNR. This reverse shock reheats the stellar ejecta, and as a result the sound velocity increases by a large factor. When the reverberations due to reflections of the reverse shock have almost completely dissipated, one can approximate the interior of the SNR by using the analytical Sedov solution (Sedov 1959). The interaction with the reverse shock influences the evolution of the pulsar wind nebula quite dramatically. Cioffi et al. (1988) have already shown in their 1D simulation of a pure shell SNR that the reverse shock gives rise to all kinds of sound waves and weak shocks traveling back and forth through the ejecta before the interior relaxes towards a Sedov solution. We will show that during the process of relaxation the radius of the pulsar wind nebula contracts and expands due to reverberations of the reverse shock. Compression waves are partly reflected and partly transmitted at the edge of the PWN. We will come back to this point when we discuss results from hydrodynamics simulations in Sect. 4, which allow a more detailed picture of this process. In this section we consider a fully relaxed Sedov SNR. The PWN expands subsonically into the remnant because the interior of the SNR has been re-heated by the reverse shock. For the case of a constant (mechanical) luminosity driving the pulsar wind an analytical expression for the radius of the PWN can be easily obtained. In this stage of the PWN evolution, we associate its radius $R_{\rm pwn}$ with the contact discontinuity separating pulsar wind material from the ejecta of the progenitor star (see Fig. 2).

Figure 2: Schematic representation of a PWN in a Sedov SNR. There are a total of 2 shocks and 2 contact discontinuities. From left to right one can see: the pulsar wind termination shock $R_{\rm ts}$ (dashed line), the first contact discontinuity $R_{\rm pwn}$ (dotted line) separating shocked pulsar wind material from shocked ejecta, bounding the PWN. Furthermore there is another contact discontinuity $R_{\rm cd}$ (dotted line) separating shocked ejecta from shocked ISM, and the SNR shock $R_{\rm snr}$ (solid line) which bounds the PWN/SNR system.

Open with DEXTER

We first present an order-of-magnitude calculation which leads to the correct power-law solution for the radius of the PWN. The assumption of subsonic expansion implies approximate pressure equilibrium between the wind material and the stellar ejecta at the edge the PWN. In the interior of the SNR the pressure scales as

$$% P_{\rm snr}\propto E_0/R_{\rm snr}^3.$$ (14)

On the other hand, the pressure in the interior of the PWN scales as

$$% P_{\rm pwn}\propto L_0t/R_{\rm pwn}^3,$$ (15)

with L₀ the mechanical luminosity driving the wind. Pressure equilibrium at the contact discontinuity at $R_{\rm pwn}$ implies the following relation for the radius of the PWN as a function of time:

$$% R_{\rm pwn}(t) = {\bar C} \left({L_0t\over E_0}\right)^{1/3}R_{\rm snr}(t)\propto t^{11/15},$$ (16)

with the constant of proportionality $\bar{C}$ to be determined below.

A more detailed derivation uses the first law of thermodynamics, assuming once again a constant energy input L₀ into the PWN by the pulsar-driven wind:

$$% {\rm d}E_{\rm th}=L_{0} \: {\rm d}t - P_{\rm i} \: {\rm d}{\cal V}_{\rm pwn}.$$ (17)

Here $E_{\rm th}$ is the thermal energy of the PWN, $P_{\rm i}$ its internal pressure, and ${\cal V}_{\rm pwn}$ its volume. This yields the following equation describing the energy balance of a slowly expanding PWN:

$$% {{\rm d}\over {\rm d}t}\left({4\pi\over 3}{P_{\rm i}R_{\rm ... ...}\: \left( \frac {{\rm d} {R}_{\rm pwn}}{{\rm d} t} \right) ,$$ (18)

or equivalently

$$% {{\rm d}\over {\rm d}t}\left({4\pi\over 3} {\gamma_{\rm pwn... ...n}^3 \: \left( \frac{{\rm d} {P}_{\rm i}}{{\rm d} t} \right) .$$ (19)

This equation has a power-law solution for $R_{\rm pwn}(t)$provided the internal pressure $P_{\rm i}(t)$ in the SNR behaves as a power-law in time so that the relation

$$% R_{\rm pwn}^3 \: \left( \frac{{\rm d} P_{\rm i}}{{\rm d} t} \right) \; = \; \mbox{constant}$$ (20)

can be satisfied. For a Sedov SNR expanding into a uniform ISM one has $P_{\rm i} \propto t^{-6/5}$ and one finds:

$$% R_{\rm pwn}(t)= D \left({L_{0} t\over P_{\rm i}(t)}\right)^{1/3} \propto t^{11/15} ,$$ (21)

where

$$% D = \left[ \frac{4\pi}{3} \: \left({\gamma_{\rm pwn}\over\gamma_{\rm pwn} -1}+\frac{6}{5} \: \right) \right]^{- 1/3} .$$ (22)

If $R_{\rm pwn} \ll R_{\rm snr}$ we can use the central pressure from the Sedov solution with $\gamma_{\rm ism} =5/3$ for the interior pressure in the SNR which confines the PWN (e.g. Shu 1992):

$$% P_{\rm i}(t)=P_{\rm snr}(t) \; \simeq \; 0.074\: \left( {E_0\over R_{\rm snr}^3} \right)\propto t^{-6/5}.$$ (23)

We find the same result for $R_{\rm pwn}(t)$ as in the order-of-magnitude calculation (Eq. (16)), determining the constant in that expression as $\bar C\simeq 0.954$ for a non-relativistic fluid ( $\gamma _{\rm pwn} =5/3$) and $\bar C\simeq 0.851$ for a relativistic fluid ( $\gamma_{\rm pwn} =4/3$). By comparing the sound speed with the expansion velocity at the edge of the PWN, we confirm that the expansion remains subsonic.

Figure 3: Comparison between results from numerical simulations and analytical result for the radius of the PWN, i.e. Eq. (12). The dashed line indicates the radius for the PWN obtained from numerical simulations. The solid line corresponds to Eq. (12) with $C \simeq 0.889$, as appropriate for $\gamma _{\rm pwn} =5/3$. The different physical parameters are as indicated in Table 1 (Simulation 1). The injected mass of the pulsar wind has been chosen in such a way that the termination velocity of the pulsar wind equals the speed of light. One can see that in the simulation the radius $R_{\rm pwn}$ is about 5% smaller than predicted by the analytical result, but the power-law behaviour $R_{\rm pwn} \propto t^{6/5}$ is correctly reproduced.

Open with DEXTER

3.2 Pulsar Wind Nebula expansion for varying wind luminosity

The constant wind luminosity assumption is not very realistic by the time the effects of the reverse shock and its associated reverberations have vanished. The spin-down luminosity of the pulsar is more realistically described by the luminosity evolution from a rotating magnetic dipole model:

$$% L(t)=\frac{L_0}{\displaystyle \left(1+{t\over\tau}\right)^2}\cdot$$ (24)

Therefore we now consider the more realistic case of a time-dependent luminosity given by (24). The energy balance equation for the PWN reads:

$${{\rm d}\over {\rm d}t}\left({4\pi\over 3} {P_{\rm i}R_{\rm p... ...a_{\rm pwn} -1)}\right) = {L_0\over\left( 1+{t/\tau}\right)^2} % -4\pi R_{\rm pwn}^2 P_{\rm i} \left( \frac{{\rm d} {R}_{\rm pwn}}{{\rm d} t} \right).$$ (25)

We solve this equation numerically using a fourth-order Runge-Kutta method (e.g. Press et al. 1992). As an initial condition we take the radius of the PWN equal to zero at the start of the evolution, neglecting the initial stage when the PWN is expanding supersonically. For the pressure $P_{\rm i}$, we use the pressure at the center of the Sedov SNR (23). We find that the solution for $R_{\rm pwn}$ converges to $R_{\rm pwn}\propto t^{3/10}$. The same result was obtained by Reynolds & Chevalier (1984) who noted that this behaviour follows from the adiabatic expansion of the relativistic pulsar wind bubble after the energy input of the pulsar has almost ceased, where the pressure scales as $P_{\rm pwn} \propto R^{-3 \gamma_{\rm pwn}} = R^{-4}$, and the behaviour of the pressure $P_{\rm i} \propto t^{-6/5}$ in the Sedov remnant. Pressure equilibrium then yields $R_{\rm pwn} \propto t^{-3/10}$. Figure 9 shows this semi-analytical result together with results from hydrodynamical simulations. For the semi-analytical equation we use $\gamma _{\rm pwn} =5/3$, because the hydrodynamics code also uses this value (see Sect. 4.1 below).

4 Numerical simulations

4.1 Method

Our simulations were performed using the Versatile Advection Code (VAC, Tóth 1996) which can integrate the equations of gas dynamics in a conservative form in 1, 2 or 3 dimensions. We used the TVD-MUSCL scheme with a Roe-type approximate Riemann solver from the numerical algorithms available in VAC (Tóth & Odstrcil 1996); a discussion of this and other schemes for numerical hydrodynamics can be found in Leveque (1998). In this paper our calculations are limited to spherically symmetric flows.

We use a uniform grid with a grid spacing chosen sufficiently fine to resolve both the shocks inside the PWN and the larger-scale shocks associated with the SNR. Table 1 gives the physical scale associated with the grid size for the simulations presented here. An expanding SNR is created by impulsively releasing the mechanical energy of the SN explosion in the first few grid cells. The thermal energy and mass deposited there lead to freely expanding ejecta with a nearly uniform density, and a linear velocity profile as a function of radius.

Figure 4: Pressure profile of the PWN/SNR system as a function of radius, at time t=1000 years after the SN explosion. Physical parameters are as indicated in Table 1 (Simulation 2). Moving outwards in radius one can see the wind termination shock, the shock bounding the PWN, the reverse shock of the SNR and the shock bounding the SNR. The interior of the PWN is nearly isobaric. There is a sudden increase in pressure of the ejecta behind the SNR reverse shock.

Open with DEXTER

	Figure 5: Density profile for the same PWN/SNR system as in Fig. 4.
Open with DEXTER

A realistic shocked pulsar wind is presumably highly relativistic with an adiabatic heat ratio $\gamma_{\rm pwn} =4/3$. The (shocked) stellar ejecta on the other hand are non-relativistic with $\gamma_{\rm ej} = 5/3$.

The VAC code does not currently include relativistic hydrodynamics. Therefore, the best approach available to us is to keep $\gamma _{\rm pwn} =5/3$, but to take a luminosity for the pulsar wind, L(t), and an associated mass injection, $\dot M_{\rm pw}(t)$, such that the terminal velocity obtained from these two parameters,

$$% v_\infty =\sqrt{2L(t)/\dot M_{\rm ej}(t)} ,$$ (26)

roughly equals the speed of light. Since the pulsar wind material downstream of the termination shock moves with only a mildly relativistic bulk speed we expect our results to be qualitatively correct. Thermal energy and mass are deposited continuously in a small volume as a source for the wind. The hydrodynamics code then develops a steady wind reaching the terminal velocity $v_\infty$ well before the (cold) wind is terminated by the standing termination shock.

We trace the total mass injected into the PWN by the pulsar wind in order to determine the radius of the contact discontinuity which separates the pulsar wind material from the SN ejecta ( $R_{\rm cd1}$ in Fig. 1 and $R_{\rm pwn}$ in Fig. 2). We also determine the position of the shock bounding the PWN during the stage of supersonic expansion. This enables us to compare the numerical results with the analytical expressions derived in Sects. 2 and 3 for the PWN radius.

  

Table 1: Simulation parameters.

	Figure 6: Sound velocity profile as a function of radius, for the same case as in Figs. 4 and 5. Because of the Rankine-Hugoniot jump conditions at the wind termination shock, the sound velocity of the shocked wind material in the PWN bubble is close to the speed of light. Behind the contact discontinuity, where the bubble consists of swept-up ejecta, the sound speed has a smaller value.
Open with DEXTER

	Figure 7: Velocity profile for the PWN/SNR system with the same parameters. The terminal wind velocity, $v_\infty$, is close to the speed of light. The large jump in the velocity at a radius $\sim$4 pc is the reverse shock which is still propagating forwards in the laboratory frame. The velocity jump of the PWN shock at radius $\sim$2.5 pc is much smaller.
Open with DEXTER

As a test of the code we have calculated a pulsar wind driven by a constant luminosity L₀ (Simulation 1 in Table 1). We let the PWN evolve until the reverse shock propagating in the SNR hits its outer edge. Figure 3 shows the radius of the shock of the PWN together with the analytical Eq. (12). We take $\gamma _{\rm pwn} =5/3$ in the analytical expressions for comparison with the numerical results. Although the analytical result of Eq. (12) is not reproduced exactly (the radius is about 5% smaller), the power-law expansion law $R_{\rm pwn} \propto t^{6/5}$ is reproduced.

4.2 Evolution of the PWN-SNR system into the Sedov phase

Our simulations of the evolution of a pulsar wind nebula inside a supernova remnant employ the parameters listed in Table 1.

In the early stage of its evolution the PWN is bounded by a strong shock propagating through the ejecta of the progenitor star. In Figs. 4-7 one can clearly identify the four shocks indicated schematically in Fig. 1. Moving outward in radius one first encounters the pulsar wind termination shock; this termination shock is followed by the PWN shock. In the sound velocity profile of Fig. 6 one can see a large jump between these two shocks: the contact discontinuity separating shocked pulsar wind material from shocked ejecta. Further outward one encounters the SNR reverse shock, which at this stage of the SNR evolution is still moving outwards from the point of view of a stationary outside observer. The whole PWN-SNR system is bounded by the SNR blast wave. Figure 9 shows the evolution of the contact discontinuity radius $R_{\rm cd}$, which can be identified with the radius of the PWN in the subsonic expansion stage. One can clearly see the moment at $t\simeq 1.75$kyr when the reverse shock hits the edge of the PWN: the expansion becomes unsteady with the PWN contracting and expanding due to the interaction with the pressure pulses associated with the reverberations of the reverse shock. When these reverberations have almost dissipated the expansion of the PWN relaxes to a steady subsonic expansion. In this stage, we can fit the radius of the PWN obtained from the simulations with the (semi-)analytical solution obtained from a numerical integration of Eq. (18), as shown in this figure.

The interaction of the PWN with the reverse shock and the associated reverberations is quite complicated. We will therefore describe this process in more detail.

	Figure 8: The radius of the PWN contact discontinuity as a function of time (solid line). We compare with the semi-analytical solution from Eq. (25) (dashed line). Here one can see that the expansion of the PWN is unsteady, due to the reverberations of the reverse shock. This simulation was done with the parameters listed in Table 1 (Simulation 3).
Open with DEXTER

Figure 9: The radius of the PWN as a function of time, together with the ratio of the total energy in the PWN bubble with respect to the total energy input by the pulsar wind. The solid line represents the radius of the contact discontinuity of the PWN (in arbitrary units), the open squares represent the aforementioned ratio of energy. This simulation was done with the parameters as listed in Table 1 (Simulation 3).

Open with DEXTER

4.3 The influence of reverse-shock reverberations

The reverse shock initially encounters the PWN in its supersonic expansion stage. After the collision between the reverse shock and the PWN shock a reflected shock propagates back towards the outer (Sedov-Taylor) blast wave of the SNR. A transmitted shock propagates into the shocked ejecta inside the PWN. When this shock hits the contact discontinuity bounding the pulsar wind material a similar reflection/transmission event occurs: a shock moves radially outwards, and a compression wave moves into the pulsar wind material. The latter wave is rapidly dissipated in the pulsar wind bubble because of the high sound speed in the shocked pulsar wind. After a few sound crossing times the pulsar wind bubble contracts adiabatically in response to the pressure increase inside the SNR. After this contraction it regains pressure equilibrium with the surrounding SNR and the PWN expands subsonically henceforth. This chain of events can be clearly seen in Fig. 8 where we plot the radius of the PWN. The whole process takes a time comparable with the duration of the initial supersonic expansion stage.

4.4 Subsonic expansion stage

When the PWN has more or less relaxed to a steady subsonic expansion the PWN has gained energy as a result of the interaction with the reverse shock. Consequently, the radius of the PWN is roughly 20% larger than the value predicted by the semi-analytical solution obtained from Eq. (18) in Sect. 2. In Fig. 9 we show the ratio between the (mostly thermal) energy of the pulsar wind bubble, i.e. the part of the PWN that consists of shocked pulsar wind material, and the total mechanical energy deposited by the pulsar. One can clearly see the increase in the energy content of the pulsar wind bubble. A large fraction of the energy deposited by the pulsar wind in the stage when the expansion is supersonic is contained in the kinetic energy of the shocked stellar ejecta in the PWN shell. When the reverse SNR shock is interacting with the PWN bubble, energy is apparently transferred from this thin shell to the interior of the bubble by the process of adiabatic compression. One can see the effect of adiabatic compression in Fig. 9 where the radius of the PWN and the energy content of the PWN are anti-correlated.

5 Conclusions and discussion

We have considered a spherically symmetric PWN/SNR system in the early and middle stages of its evolution, well before cooling of the SNR shell becomes dynamically important and before a significant disruption of spherical symmetry due to a possible (large) kick velocity of the pulsar can take place. The expansion of the PWN is coupled with the dynamics of the expanding SNR, leading to two distinct evolutionary stages separated by an unsteady transition phase:

• When the PWN is surrounded by the freely expanding ejecta of the SNR, the expansion of the PWN is supersonic. In this stage the pressure in the interior of the PWN bubble is slightly larger than one would expect from ram pressure of the surrounding ejecta alone, using the Rankine-Hugoniot relations at the PWN shock. This is due to the thin shell of shocked, swept-up ejecta which needs to be accelerated by the outward force due to the interior pressure of the PWN.
• This stage of supersonic expansion is ultimately followed by a subsonic expansion of the PWN. This happens after the reverse shock has encountered the shock of the PWN.
• The transition between these two stages is unsteady due to the interaction of the PWN with the reverse shock and its associated reverberations. From the hydrodynamical simulations we see that the time scale for adjustment to the pressure of the surrounding SNR is determined by the sound speed of the ejecta shocked by the PWN in the first stage.

Two of the prototypes of plerionic SNRs are the Crab Nebula and 3C58. In the Crab, there is a decrease in the radio flux of $0.167 \pm 0.015\%$ yr^-1 (Aller & Reynolds 1985). By contrast, 3C58 shows an increase in its flux density at radio frequencies between 1967 and 1986 (Green 1987). This increase might be the result of the reverse shock which has encountered the PWN shock around 3C58; the PWN is being compressed and therefore the flux density is going up.

Our numerical simulations are different from the results presented by Jun (1998). This author concentrates on the details of the PWN in the supersonic expansion stage, and in particular on the formation of Rayleigh-Taylor fingers in his two-dimensional simulations. Our simulations include the whole supernova remnant, but can not address the development of Rayleigh-Taylor instabilities due to our assumption of spherical symmetry.

In future work we will discuss how these results change when the influence of a significant kick velocity of the pulsar is taken into account. If this is taken into account, the model presented here will lose its validity at a certain time: one can calculate when the motion of a pulsar will become supersonic in a Sedov stage. One can show that this will happen when the pulsar is about halfway from the explosion center to the edge of the SNR: a bow shock is expected to result from this and clearly the model presented here will break down. Observationally there is evidence that this is the case for the pulsar associated with the SNR W44.

Appendix A: Thin shell approximation

We briefly derive an expression for the pressure in the PWN, and for the thickness of the shell of shocked ejecta material surrounding the PWN in the first (supersonic) expansion stage of PWN evolution.

The starting point of our derivation is the assumption that velocity of the material in the shell moves at a constant velocity. This assumption is justified by detailed calculations by Chevalier (1984) who calculates the velocity, pressure and density using similarity variables. We denote the radius of the contact discontinuity separating pulsar wind material from ejecta material by $R_{\rm c}$, and the radius of the shock wave propagating into the ejecta by $R_{\rm S}$. The constant velocity assumption implies that the velocity of the contact discontinuity must equal the fluid velocity directly behind the shock. Assuming the shock to be strong with compression ratio $(\gamma_{\rm ej}-1)/(\gamma_{\rm ej} +1)$ this implies:

$$% \frac{{\rm d} R_{\rm c}}{{\rm d} t} = \frac{2 }{\gamma_{\rm... ..._{\rm ej} -1}{\gamma_{\rm ej} +1} \: \frac{R_{\rm s}}{t} \cdot$$ (A.1)

Here we used the fact that upstream of the shock the velocity of the uniformly expanding ejecta equals $V_{\rm ej}(R_{\rm s}) = R_{\rm s}/t$. Assuming a self-similar solution where $R_{\rm c} \: , \: R_{\rm s} \propto t^{\alpha}$ the above condition implies

$$% R_{\rm c} = \frac{2 \alpha + \gamma_{\rm ej} - 1}{\alpha (\gamma_{\rm ej} +1)} \: R_{\rm s}.$$ (A.2)

If we assume $\gamma_{\rm ej} = 5/3$ and use that $\alpha = 6/5$ if the mechanical power L of the pulsar is constant we find

$$% R_{\rm c} = \frac{23}{24} \: R_{\rm s} ,$$ (A.3)

and the thickness of the layer of shocked ejecta material is $\Delta r = R_{\rm s}/24$. The thin layer approximation is therefore a good one.

We now neglect the difference between $R_{\rm c}$ and $R_{\rm s}$, putting

$$R_{\rm c} \sim R_{\rm s} \equiv R_{\rm pwn}$$

as is done in the main paper. The equation of motion for the thin layer of shocked ejecta can be approximated by (e.g. Ostriker & Gunn 1971)

$$% M_{\rm sw}(t) \: \frac{{\rm d}^{2} R_{\rm pwn}}{{\rm d} t^{... ...eft[ P_{\rm pwn} - \rho_{\rm ej}(t) \: V_{\rm rel}^{2}\right],$$ (A.4)

where

$$% V_{\rm rel} = \frac{{\rm d} R_{\rm pwn}}{{\rm d} t} - \frac{R_{\rm pwn}}{t}$$ (A.5)

is the relative velocity between the shell and the expanding ejecta. The outward motion of the shell is driven by the difference between the presure in the PWN and the ram-pressure $\rho_{\rm e} V_{\rm rel}^{2}$. Here $M_{\rm sw}(t)$ is the mass of the swept-up ejecta. Using Eqs. (1) and (6) of the main paper we can write

$$% M_{\rm sw}(t) = \frac{4 \pi}{3} \: \rho_{\rm ej}(t) \: R_{\rm pwn}^{3}.$$ (A.6)

Assuming that the PWN expands as $R_{\rm pwn} \propto t^{\alpha}$, the above equation can be solved for the pressure $P_{\rm pwn}$ of the pulsar wind material at the contact discontinuity:

$$% P_{\rm pwn}(t) = \frac{(\alpha - 1)(4 \alpha - 3)}{3} \: \rho_{\rm ej}(t) \left( \frac{R_{\rm pwn}}{t} \right)^{2} .$$ (A.7)

The pressure just behind the strong shock at the outer edge of the layer is given by

$$% P_{\rm sh} = \frac{2 }{\gamma_{\rm ej} + 1} \: \rho_{\rm ej... ... R_{\rm pwn}}{{\rm d} t} - \frac{R_{\rm pwn}}{t} \right)^{2}.$$ (A.8)

Using once again $R_{\rm pwn} \propto t^{\alpha}$, this implies that the ratio of the pressure in the PWN and the post-shock pressure equals

$$% \frac{P_{\rm pwn}}{P_{\rm sh}} = \frac{(\gamma_{\rm ej} +1... ... \mbox{$\frac{3}{2}$ } \: \left( \gamma_{\rm ej} + 1 \right),$$ (A.9)

where we assumed $\alpha = 6/5$ in the second equality. For $\gamma_{\rm ej} = 5/3$ one has $P_{\rm pwn} = 4 P_{\rm sh}$, in good agreement with the calculations of Chevalier (1984) (see his Fig. 4). The pressure rises considerably in the layer of shocked ejecta material, where the sound speed is low. This result implies that the energy of the hot material in the PWN, defined in Eq. (7) of the main paper, equals

$$E_{\rm pwn} \frac{(4 \alpha - 3)(\alpha - 1)}{3(\gamma_{\rm pwn} - 1)} \: M_{\rm sw} \: \left( \frac{R_{\rm pwn}}{t} \right)^{2}$$ = (A.10)

$$\frac{3 M_{\rm sw} }{25 (\gamma_{\rm pwn}-1)} \: \left( \frac{R_{\rm pwn}}{t} \right)^{2} ,$$ =

where we used $\alpha = 6/5$ in the second equality.

Acknowledgements

The authors would like to thank L. Drury for discussions on the subject, T. Downes for discussing several numerical techniques, D. Frail for clarifying the nature of W44, and R. Keppens for his assistance with the code VAC. We thank the referee, Dr. R. A. Chevalier, for his helpful comments. The Versatile Advection Code has been developed as part of the project on "Parallel Computational Magneto-Fluid Dynamics'', funded by the Dutch Science Foundation (NWO) from 1994 to 1997. EvdS is currently supported by the European Commission under the TMR programme, contract number ERB-FMRX-CT98-0168. Y.A.G. acknowledges support from The Netherlands Foundation for Research in Astronomy (GBE/MPR grant 614-21-008) and the Italian Ministry of Universities and Research (grant Cofin-99-02-02). G.T. is currently supported by the Hungarian Science Foundation (OTKA, grant No. D 25519) and the Education Ministry of Hungary (grant No. FKFP-0242-2000).

References

 

• Aller, H. D., & Reynolds, S. P. 1985, ApJ, 293, L73 In the text NASA ADS  
• Chevalier, R. A. 1977, in Supernovae, ed. D. N. Schramm (D. Reidel Publ. Co., Dordrecht, Holland), 53 In the text  
• Chevalier, R. A. 1984, ApJ, 280, 797 In the text NASA ADS  
• Chevalier, R. A. 2000, ApJ, 539, L45 In the text NASA ADS  
• Chevalier, R. A., & Fransson, C. 1992, ApJ, 395, 540 In the text NASA ADS  
• Cioffi, D. F., McKee, C. F., & Bertschinger, E. 1988, ApJ, 334, 252 In the text NASA ADS  
• Cioffi, D. F. 1990, in Physical Processes in Hot Cosmic Plasmas, ed. W. Brinkmann, A. C. Fabian, & F. Giovanelli (Kluwer Academic Publ., Dordrecht, Holland), 1 In the text  
• Frail, D. A., Goss, W. M., & Whiteoak, J. B. Z. 1994, ApJ, 437, 781 In the text NASA ADS  
• Frail, D. A., & Scharringhausen, B. R. 1997, ApJ, 480, 364 In the text NASA ADS  
• Green, D. A. 1987, MNRAS, 225, 11 In the text  
• Helfand, D. J. 1998, MSAI, 69, 791 In the text  
• Hester, J. J., & Kulkarni, S. R. 1988, ApJ, 331, L121 In the text NASA ADS  
• Jun, B.-I. 1998, ApJ, 499, 282 In the text NASA ADS  
• Leveque, R. J. 1998, in Computational Methods for Astrophysical Fluid Flow, Saas Fee Adv. Course 27, ed. R. J. Leveque, D. Mihalas, E. A. Dorfi, & E. Müller (Springer Verlag, Heidelberg), 1 In the text  
• McKee, C. F. 1974, ApJ, 188, 335 In the text NASA ADS  
• Ostriker, J. P., & Gunn, J. E. 1971, ApJL, 164, L95 In the text NASA ADS  
• Press, W. H., Teukolsky, S. A., Vertterling, W. T., & Flanery, N. P. 1992, Numerical recipes in C (Cambridge University Press, New York) In the text  
• Reynolds, S. P., & Chevalier, R. A. 1984, ApJ, 278, 630 In the text NASA ADS  
• Sedov, L. I. 1959, Similarity and Dimensional Methods in Mechanics (Academic, New York) In the text  
• Shu, F. H. 1992, The Physics of Astrophysics: vol. II Gas Dynamics (University Science Books, Mill Valley) In the text  
• Strom, R. G. 1987, ApJ, 319, L103 In the text NASA ADS  
• Tóth, G. 1996, Astrophys. Lett. & Comm., 34, 245 In the text  
• Tóth, G., & Odstrcil, D. 1996, J. Comp. Phys., 128, 82 In the text  
• Woltjer, L. 1972, ARA&A, 10, 129 In the text NASA ADS