Diffusion in pulsar wind nebulae: an investigation using magnetohydrodynamic and particle transport models

Abstract

We study the transport of high-energy particles in pulsar wind nebulae (PWN) using three-dimensional magnetohydrodynamic (MHD) and test-particle simulations, as well as a Fokker–Planck particle transport model. The latter includes radiative and adiabatic losses, diffusion, and advection on the background flow of the simulated MHD nebula. By combining the models, the spatial evolution of flux and photon index of the X-ray synchrotron emission is modelled for the three nebulae G21.5−0.9, the inner regions of Vela, and 3C 58, thereby allowing us to derive governing parameters: the magnetic field strength, average flow velocity, and spatial diffusion coefficient. For comparison, the nebulae are also modelled with the semi-analytic Kennel & Coroniti model but the Porth et al. model generally yields better fits to the observational data. We find that high velocity fluctuations in the turbulent nebula (downstream of the termination shock) give rise to efficient diffusive transport of particles, with average Péclet number close to unity, indicating that both advection and diffusion play an important role in particle transport. We find that the diffusive transport coefficient of the order of ∼ 2 × 10²⁷(L_s/0.42 Ly) cm² s^{− 1} (L_s is the size of the termination shock) is independent of energy up to extreme particle Lorentz factors of γ_p ∼ 10¹⁰.

diffusion, MHD, turbulence, pulsars: individual: 3C 58, pulsars: individual: G21.5−0.9, pulsars: individual: Vela

1 INTRODUCTION

Pulsars produce highly relativistic winds that consist of electrons and positrons, and frozen-in magnetic field (e.g. Kirk, Lyubarsky & Petri 2009). Where the ram pressure of this magnetized wind equals the confining pressure of the ambient medium, a termination shock is formed (Rees & Gunn 1974; Kennel & Coroniti 1984a). Charged particles are accelerated at the termination shock (e.g. Kennel & Coroniti 1984b; Reynolds & Chevalier 1984), presumably via the Fermi acceleration mechanism but also magnetic reconnection of the striped pulsar wind could play an important role (see e.g. Sironi & Spitkovsky 2011). Downstream of the termination shock, relativistic leptons produce synchrotron radiation from radio to X-ray wavelengths. Particles are continuously transported away from the termination shock by the downstream flow, inflating a luminous nebula around the pulsar commonly known as a pulsar wind nebula (PWN).

Overall, the synchrotron energy spectra can be described using a combination of power laws (see e.g. De Jager & Djannati-Ataï 2009). Spatially resolved X-ray observations of PWNe show that the photon index Γ evolves as a function of distance from the pulsar (see e.g. Slane et al. 2000; Bocchino et al. 2001; Mangano et al. 2005; Schöck et al. 2010). As shown by e.g. Vorster & Moraal (2013), the evolution of the lepton spectrum is not only determined by the energy loss rate, but also by the rate at which particles are transported towards the outer parts of the nebula. While particle evolution models of PWNe typically only take convection into account, it has been suggested that diffusion may also play an important role in the transport process (see e.g. Gaensler & Slane 2006). This idea is supported by the results of Tang & Chevalier (2012), who found that they were able to better model the evolution of Γ observed for G21.5−0.9 and 3C 58 by including both convection and diffusion in their model, in contrast to only including convection.

The exact nature of diffusion is still an open question, with both kinetic and turbulent processes at its core. Adopting a kinetic view, one expects diffusion to be related to the magnetic field geometry, while there exists some uncertainty as to the correct topology to use for PWNe. If one assumes the purely toroidal magnetic field from the classic Kennel & Coroniti (1984a) model, then all outward diffusion is perpendicular to the field lines which render kinetic diffusive transport inefficient. In contrast to this simple picture, observations of filamentary structures in the outer regions of Crab nebula (e.g. Hester 2008, and references therein) support the presence of a poloidal field component that wraps around and along the filaments. This impact on the transport of particles, which can then slide outwardly along the field lines. In fact, deep Chandra X-ray observations by Seward, Tucker & Fesen (2006) find indications that spectral index variations along filaments are less pronounced than in the cross-field direction, highlighting the importance of the magnetic topology in particle transport.

There is now little doubt that the prominent outer filaments of Crab nebula are due to the Rayleigh–Taylor instability (RTI) of the interface between the PWN and the supernova remnant (Chevalier & Gull 1975; Jun 1998; Porth, Komissarov & Keppens 2014b; Bietenholz & Nugent 2015). In general, the RTI can be suppressed for strong and ordered field (e.g. Bucciantini et al. 2004) or when the interface has stopped accelerating in old PWNe with age of ∼ 10⁴ yr (see also Chevalier & Reynolds 2011). Thus it is uncertain whether the RTI on its own can account for the destruction of ordered field and seed turbulence throughout the nebula.

However, fluid instabilities in the bulk of the nebula, triggered by the causally connected PWN jet and equatorial shear-flow contribute to a destruction of the ordered toroidal field that is injected with the relativistic wind (e.g. Begelman 1998; Mizuno et al. 2011; Porth et al. 2014a). Perhaps more than by changing the geometry, fluid instabilities can directly enhance particle transport through driving a turbulent cascade in which particles are transported diffusively as we will show in Section 2.

The present paper aims to address two specific points: the first is to use the magnetohydrodynamic (MHD) results of Porth et al. (2014a) to study particle transport and derive diffusion coefficients for PWNe. The second is to further investigate the role of diffusion in the evolution of the leptons spectra by extending the results of Tang & Chevalier (2012). These authors focused on modelling the spatial evolution of Γ, while neglecting the corresponding evolution of the X-ray synchrotron flux. Naively one might expect that a good fit to Γ would automatically produce a good fit to the flux. However, Vorster (2013) found that the evolution of Γ is, to a degree, determined by the strength of the magnetic field, the convection velocity, and the diffusion coefficient relative to each other. It is therefore possible to find a similar evolution of Γ by suitably varying the above-mentioned parameters. The paper therefore addresses this issue by taking into account the spatial evolution of the flux in the modelling, where possible.

While the magnetic field and flow velocity calculated by Porth et al. (2014a, hereafter PKK14) are the primary choice for investigating the role of diffusion, it is useful to compare the results with those found from using the classic model developed by Kennel & Coroniti (1984a, hereafter KC84), where the spatial dependence of the magnetic field and flow velocity are determined by the ratio σ of magnetic to particle energy in the PWN. The data fitting will therefore allow one to derive the value of this important parameter within the context of the KC84 model. Thus far this value has only been rigorously derived for the Crab nebula (Kennel & Coroniti 1984b).

The rest of the paper is structured as follows. Section 2 presents the diffusion coefficients derived from the PKK14 model and characterizes the turbulence obtained in the simulations. Section 3 introduces the particle transport model that is used to fit the X-ray data, with the modelling results presented in Section 4. Lastly, a discussion of the results together with the main conclusion drawn from the modelling can be found in Section 5.

2 DIFFUSION IN PULSAR WIND NEBULAE

The coefficients of ordinary diffusion in the r-direction can be defined as

\begin{eqnarray} D_{rr}(\Delta t) = \frac{\langle (\Delta r)^2 \rangle }{2 \Delta t}, \end{eqnarray}

(1)

where the average is taken over the differences in positions Δr = r(t + Δt) − r(t). For the diffusion approximation to hold, many small angle deflections within Δt are required, and result in a stochastic behaviour of the particles. It is therefore required that Δt ≫ τ, where τ is the collision time-scale. Under these circumstances D_rr becomes independent of Δt. On the other hand, for anomalous diffusion the mean squared displacement follows a more general relation 〈Δr²〉 ∝ Δt^β, with β ≠ 1 such that the running diffusion coefficient, as given by equation (1), will show an inherent time dependence. Whether particle transport in simulated PWN obeys relation (1) or a more general law will be investigated in Section 2.3. We first discuss two diffusion processes likely to occur in PWN and then move to the numerical results.

2.1 Ordered field: Bohm diffusion

Analytical (KC84) as well as two-dimensional MHD models of PWNe adopt a purely toroidal field geometry. Radial transport can then occur either due to (turbulent) advection or kinetic cross-field diffusion. A useful, upper limit for cross-field diffusion is provided by the empirical Bohm law:

\begin{eqnarray} D^{\rm B} = \frac{1}{3}r_{\rm g}^2 \omega _{\rm g} = 1.7 \times 10^{26} \left(\frac{\Gamma _{\rm p}}{10^9}\right) \left(\frac{B}{100{\rm \,\mu G}}\right)^{-1}\, \rm cm^2\, s^{-1}, \end{eqnarray}

(2)

with diffusion time

\begin{eqnarray} t_{\rm B} = \frac{L^2}{2D^{\rm B}} = 3 \times 10^{3} \left(\frac{\Gamma _{\rm p}}{10^9}\right)^{-1} \left(\frac{B}{100\,{\rm \mu G}}\right) \left(\frac{L}{2\, \rm pc}\right)^2\, \rm yr. \end{eqnarray}

(3)

Even for PeV particles this is larger than the age of the Crab nebula.

On the other hand, the advective time-scale t_a can be estimated from the flow velocity in the nebula. In the laminar model of KC84, the flow velocity monotonously connects the expansion speed v_n with the fast flow near the termination shock |$v_{\rm f}\sim c/\sqrt{3}$|⁠, hence

\begin{eqnarray} 4\times 10^{8} \left(\frac{L}{2\,\rm pc}\right) \left(\frac{v_{\rm f}}{0.5 c}\right)^{-1}\, {\rm s}\le t_{\rm a} \end{eqnarray}

(4)

\begin{eqnarray} \le 4\times 10^{10} \left(\frac{L}{2\,\rm pc}\right) \left(\frac{v_{\rm n}}{1500\, \rm km\,s^{-1}}\right)^{-1}\, \rm s, \end{eqnarray}

(5)

which is shorter than t_B for particles up to 1 PeV. Thus, kinetic particle diffusion with D^B is negligible in PWNe.

2.2 Diffusion in a turbulent flow

Next to the field geometry, an important factor in particle diffusion in PWNe is the presence of vigorous turbulence. In this case even particles tied to a given field line experience special diffusion due to the effect of field line wandering (e.g. Salu & Montgomery 1977; Kirk, Melrose & Priest 1994; Matthaeus et al. 1995). This ‘frozen-in’ diffusion should not be mistaken with energy-dependent kinetic diffusion due to wave–particle interactions as it is often encountered in theories for cosmic ray diffusion (e.g. Shalchi 2009).

Let us now provide a back of the envelope estimate for the effective diffusivity that is connected to the turbulent nebula flow. The role of the velocity field can be easily seen at the simplified axisymmetric system where |$\boldsymbol {B}=B \boldsymbol {\hat{e}}_{\phi }$|⁠, |$\boldsymbol {v}=v_{r}\boldsymbol {\hat{e}}_{r}+v_{\theta }\boldsymbol {\hat{e}}_{\theta }$|⁠. This reflects our plasma injection conditions in good approximation to the pulsar wind far away from the source. Here, the drift velocity |$\boldsymbol {v}_{{\rm D}}=\boldsymbol {E}\times \boldsymbol {B}/B^{2}$| reduces to the fluid velocity |$\boldsymbol {v}_{{\rm D}}=\boldsymbol {v}$| and to first order, particles are simply advected with the flow. In case the flow exhibits turbulent motion, particles will follow the complex flow field which in our application leads to a diffusive transport of the particles (see Section 2.3).

For the largest eddies on the scale of the termination shock L_s (the driving scale of the turbulence), we obtain a typical velocity of v_f ∼ 0.5c close to the sound speed in an ultrarelativistic gas |$c_{\rm s}=c/\sqrt{3}$|⁠. The corresponding eddy turnover time is

\begin{eqnarray} t_{\rm Ls} = \frac{L_{\rm s}}{v_{\rm f}} = 2.8\times 10^7 \left(\frac{L_{\rm s}}{0.42\, \rm Ly}\right)\left(\frac{v_{\rm f}}{0.5c}\right)^{-1}\, \rm s, \end{eqnarray}

(6)

which is much shorter than the nebula age. Thus interpreting t_Ls as the collision time, we have Δt ≫ τ and expect a diffusive transport of particles tied to the eddies. If the nebula expands self-similarly, the latter will remain true during all of its evolution. The effective diffusion coefficient due to eddies on the scale L_s is hence estimated with

\begin{eqnarray} D^{\rm E}_{\rm Ls} = \frac{1}{3} v_{\rm f} L_{\rm s} = 2.1 \times 10^{27}\left(\frac{v_{\rm f}}{0.5c}\right)\left(\frac{L_{\rm s}}{0.42\, \rm Ly}\right)\, \rm cm^2\, s^{-1}, \end{eqnarray}

(7)

where we have adopted a mean free path of L_s and mean velocity v_f. We again estimate the diffusion time-scale for typical parameters of a young PWN:

\begin{eqnarray} t_{\rm E} = \frac{L^2}{2D^{\rm E}} = 8.57\times 10^{9} \left(\frac{v_{\rm f}}{0.5c}\right)^{-1}\left(\frac{L_{\rm s}}{0.42\, \rm Ly}\right)^{-1}\left(\frac{L}{2\, \rm pc}\right)^2\, \rm s, \end{eqnarray}

(8)

which is less than 300 yr for the Crab nebula. This corresponds to the ‘escape time’ of particles injected at the termination shock.

Adopting a Kolmogorov (1941) scaling for the velocities on the two scales λL and L: v_λL ∼ λ^1/3v_L, we obtain the scalings for the eddy turnover time t_λL ∼ λ^2/3t_L and for the diffusion coefficient |$D^{\rm E}_{\lambda L}\sim \lambda ^{4/3} D^{\rm E}_{L}$|⁠. One can see that the largest scales dominate the diffusion process which is why we could neglect the contribution of smaller scales in our estimate (7). Indeed, if we define a scale-averaged diffusion coefficient,

\begin{eqnarray} \langle D^{\rm E} \rangle _{\lambda 0} \equiv \frac{1}{1-\lambda _0}\int _{\lambda 0}^1 \lambda ^{4/3} D_{\rm Ls}^{\rm E}\, {\rm d}\lambda , \end{eqnarray}

(9)

we see that the effect of allowing smaller scales is rather small: in the limit λ₀ → 0, we merely obtain |$\langle D^{\rm E} \rangle _{\lambda }\simeq 0.43 D_{{\rm Ls}}^{\rm E}$|⁠.

Note that the diffusion time in terms of the eddy turnover time is just t_E/t_Ls = 3/2(L/L_s)² which is roughly conserved in the self-similar expansion of the nebula. In Crab, we have L/L_s ∼ 14, thus particles witness several hundred turnover times before escaping from the system.

The coefficient given by equation (7) is quite universal as its governing dependence is just the termination shock size – it could thus easily be adapted to other PWNe. Since D^E > D^B, for energies up to ∼ 6 PeV, we do not expect an energy dependence of the diffusion coefficient for the large majority of particles. In particular, no energy dependence for X-ray synchrotron emitting particles is expected.

2.3 Test-particle diffusion in 3D simulations of PWN

To estimate the diffusion coefficient in the MHD simulations of PKK14, we follow orbits of test particles in the MHD flow, similar to the approach by Wisniewski, Spanier & Kissmann (2012) in application to heliospheric plasma. Monoenergetic electrons are injected uniformly and isotropically into a simulation snapshot and the particles are advanced according to the Lorentz-force resulting from the MHD fields |$\boldsymbol {E},\boldsymbol {B}$| (e.g. Landau & Lifshitz 1960):

\begin{eqnarray} \frac{{\rm d} \boldsymbol {u}}{{\rm d}t} = \frac{q}{mc} \left( \boldsymbol {E}+\frac{\mathbf {u}\times \boldsymbol {B}}{c\gamma _{\rm p}}\right) + \boldsymbol {g}, \end{eqnarray}

(10)

where |$\boldsymbol {u}=\Gamma _{\rm p}\boldsymbol {v}/c$| is the particles four-velocity and q/m signifies the electron charge to mass ratio.

Although implemented in the code, here we omit the radiation reaction force which is responsible for particle cooling and set |$\mathbf {g}=\boldsymbol {0}$|⁠. Considering a typical average field strength of 100 μm G and particles with the highest relevant Lorentz factor Γ_p = 10⁹, the synchrotron lifetime given by equation (17) becomes t_syn = 2.5 yr and thus approaches the adopted test-particle integration time of ∼ 1 yr. However, radiative losses are not relevant for the measurement of the transport coefficients for given particle energy which is the main goal of this section. Later in the spectral modelling discussed in Section 3 the synchrotron cooling is consistently taken into account.

Equation (10) is advanced with a second-order symplectic Boris scheme (e.g. Birdsall & Langdon 1991) where the fields |$\boldsymbol {E},\boldsymbol {B}$| at the particles position are obtained via linear interpolations in space and time between the CFL limited fluid steps. The electric field is taken consistent with the ideal MHD condition |$\boldsymbol {E}=1/c\, (\mathbf {B}\times \boldsymbol {v})$|⁠. In fact, when interpolating the fields to the particles position, we ensure to machine precision that the ideal MHD condition is satisfied by interpolating both |$\boldsymbol {B}$| and |$\boldsymbol {v}$| separately and then calculating the corresponding electric field. This way, particle acceleration due to a spurious parallel field component is excluded.

For numerical stability, each particle is advanced with individual time step ensuring that one gyration is resolved by at least 60 steps. This particle treatment is implemented in MPI-AMRVAC¹ (Keppens et al. 2012; Porth et al. 2014c) allowing for massively parallel calculations and adaptive mesh refinement. To obtain sufficient statistics, we follow the orbits of 5 × 10⁵ particles. The resulting diffusion coefficients can then be measured for different time intervals Δt according to equation (1). Here we restrict ourselves to the models with parameters as ‘B3D’ of PKK14. Thus we adopt an obliqueness angle of α = 45°, and an initial wind magnetization σ₀ = 1.

To check if normal diffusion is an accurate description for the transport process, we first investigate the convergence of the running diffusion coefficient. Fig. 1 (top panel) thus shows the mean radial diffusion coefficient obtained within the PWN volume as a function of the time interval. Here D_rr(Δt) is calculated by simply replacing ΔxΔy in equation (1) with Δr² for the radial transport. For the relatively short averaging times employed, transport due to mean-flow advection (see also Section 3.2) is negligible at least in the outer parts of the nebula and has not been subtracted in the determination of D_rr.

Figure 1.

Top: convergence in time of the mean radial diffusion coefficient 〈D_rr(t)〉 for particles with γ_p = 10⁸ for increasing resolutions of the MHD flow, where Δx = 2 × 10¹⁶ cm. Convergence in time is obtained after ≈1 yr. After 1.5 yr, the estimated diffusion coefficient decreases again, most likely due to particle escape. Bottom: radial profiles of the radial diffusion coefficient for γ_p = 10⁸ for various time intervals t.

After an initial ballistic regime where D_rr ∝ t, convergence is obtained at approximately 1 yr (roughly three t_Ls) at a value of 〈D_rr〉 ≈ 3.5 × 10²⁶ cm² s^{− 1}. Taking into account the termination shock size of L_s ≃ 0.2 Ly in the simulation, this agrees reasonably well with the expectation given by equation (7). As the numerical resolution of the MHD simulation is increased by a factor of 4, the temporally converged diffusion coefficient decreases by 14 per cent from a value of 3.5 × 10²⁶ to 3 × 10²⁶ cm² s^{− 1}. We suspect that this moderate dependency reflects the fact that our simulations are not yet able to reach perfect convergence of the velocity power spectrum.

The bottom panel of Fig. 1 illustrates the radial profile of D_rr. In the inner (laminar) region of the PWN, particle diffusion is suppressed with coefficient rising to a plateau reached at ≈1/3 of the nebula radius. The systematic radial advection present in the inner nebula leads to an outward shifting offset of the profiles as longer time intervals are considered. At the outer edge, particle escape leads to a systematic decrease of the measured coefficient as the time interval increases. This cautions us to restrict the averaging times to ≈1 yr maximum.

Finally, we show flow characteristics and a map of the corresponding 〈D_rr〉_ϕ in Fig. 2. The out-of-plane magnetic field indicates that the slowly varying striped wind (with obliqueness angle α = 45°) forms a current-sheet in the torus region of the nebula. Magnetic polarity mixes in the bulk of the nebula as can be seen at the example of the blue strand of positive flux in the north-western quadrant. None the less, the net polarity in each hemisphere is given by the injected one. Equatorial shear flow and the violently unstable polar jet generate turbulence and small-scale flow within one to two termination shock radii. This way, the instantaneous flow velocity differs substantially from the radial expansion of the average flow. Only the latter has to match with the slow expansion velocity of the nebula. This is a major difference to the laminar model of KC84.

Figure 2.

Flow characteristics in highest resolution run with Δx = 5 × 10¹⁵ cm, shown in the y = 0 plane. Left: out-of-plane component of the magnetic field illustrating the formation of an equatorial current-sheet j_r in the torus region of the nebula. The net polarity in the two hemispheres is still given by the injected field, indicating a toroidal guide field. Middle: velocity magnitude illustrating development and decay of small-scale turbulent flow within 1–2 termination shock radii. Right: map of the radial diffusion coefficient as obtained from the test particles and averaged over azimuthal direction. High transport coefficients are realized especially in the polar plume. The drop at the outer radius is a systematic introduced by particle escape.

Particle transport is highest in the region of the polar plume and near the equator. This is expected, as in these regions, violent mixing is triggered via kink instability and fast shear-flow. We obtain a variation of the transport coefficients by almost two orders of magnitude in the nebula where the minimum is found at intermediated polar angles of 45°.

The variation of the local diffusion coefficients with particle energy are shown in Fig. 3.

Figure 3.

Transport coefficient for various particle energies. Diffusion remains energy independent until the gyroradius becomes larger than the scale L_S at Γ_p > 10⁹. Results for Γ_p = 10⁷ and 10⁸ overlap almost identically. The x-axis extends all the way to the nebula boundary r_n ≃ 1.2 × 10¹⁸ cm.

As expected, the spatial diffusion coefficient increases once the gyroradius,

\begin{eqnarray} r_{\rm g}=\frac{p_\perp c}{e B} = 1.7\times 10^{16} \left(\frac{\Gamma _{\rm p}}{10^9}\right)\left(\frac{B}{100\,\mu \mathrm{G}}\right)^{-1}\, \rm cm, \end{eqnarray}

(11)

becomes comparable to the scale L_S. At a typical magnetic field strength of 100 μG, this is the case for ultrarelativistic particles with Γ_p ≈ 10¹⁰. These unrealistically high energy particles surpass even the Crab-flare particles of Γ_p ≈ 1 × 10⁹–3 × 10⁹ (Abdo et al. 2011; Tavani et al. 2011) which is why we will pay no further mind to them. However, for X-ray emitting particles with energies of Γ_p ≲ 10⁷, we do not observe an energy dependence of the diffusion coefficient. It is clear that the numerical resolution of Δx = 2 × 10¹⁶ cm does not suffice to resolve variations in the flow on the order of the gyroradius for Γ_p ≲ 10⁹ and the particle transport reduces to the motion of the guiding centres. In this sense, the lacking energy dependence comes as no surprise. However, as the relevant transport process for X-ray emitting particles is turbulent eddy diffusion which in turn is dominated by the well-resolved largest scale of the cascade (Section 2.2), we expect the results obtained for D_rr to be in the physically relevant regime.

2.4 Characterization of the turbulent flow

To derive quantities of interest for the turbulent transport of particles, we recover an averaged and a fluctuating part of the MHD vector fields:

\begin{eqnarray} \langle B\rangle &=& \sqrt{\langle B_\phi \rangle ^2+\langle B_r\rangle ^2+\langle B_\theta \rangle ^2}, \end{eqnarray}

(12)

\begin{eqnarray} \delta B &=& \sqrt{\delta B_r^2+\delta B_\theta ^2+\delta B_\phi ^2}, \end{eqnarray}

(13)

\begin{eqnarray} \delta B_i &=& \sqrt{\langle B_i^2\rangle - \langle B_i\rangle ^2}, \end{eqnarray}

(14)

where the averages are performed over several correlation time-scales and the azimuthal direction. Since the initial and boundary conditions are independent of the ϕ-coordinate, we can interpret the ϕ-average as an ensemble average which improves the statistics of the data.

Maps of the average and fluctuating magnetic field are shown in Fig. 4. As discussed already in PKK14, the magnetic field in the nebula varies by at least a factor of 10 from the area around the termination shock to the outer parts. The magnetic field strength is lowest in the equatorial and polar regions where changes of polarity occur. These are also the regions of the strongest turbulence as characterized by δB/〈B〉. The fluctuating component can become more than 10 times the average field strength. In the more moderate mid-latitude regions of the outer nebula, we still obtain a range of δB/〈B〉 ∈ [0.3–1]. Close to the termination shock however the field is regular and δB/〈B〉 ∼ 0.1. This is true in particular for the flow in the arch-shock where the knot emission is suggested to originate (Yuan & Blandford 2015; Lyutikov, Komissarov & Porth 2016).

Figure 4.

Characterization of the turbulence in magnetic (top) and velocity fields (bottom). Left: logarithm of the average fields (magnetic field given in Gauss). We also show stream-lines of the poloidal velocity vector fields as white lines. Two large-scale counter-streaming vortices form on each hemisphere. Middle: magnitude of the fluctuating field component. Right: fluctuating components in terms of the average. We average over the azimuthal angle and 21 snapshots, corresponding to 20 yr.

Comparing Figs 4 with 3, one can see that fluctuations of the magnetic field do not correlate necessarily with regions with strong diffusion of particles as would be suggested by quasi-linear theory (QLT) and field line random walk (FLRW) of particle transport (e.g. Shalchi et al. 2012). The high values of δB/B ≳ 10 in the outer equatorial region are not seen as enhancements in the diffusion map.

The mean velocity in the nebula highlights the fast jet and equatorial flow and we recover two large-scale vortices per hemisphere that re-cycle this material. Comparing the fluctuating component with the average velocity shows that the fluctuations in the recycling region are typically stronger than the mean flow and hence the large-scale vortices cannot be a stationary feature. In regions of fast systematic flow velocity, fluctuations are relatively weak, however, velocity fluctuations are largest with δv/v ≳ 10 in the corresponding shear layers. It is interesting to point out that right at the base of the jet, both velocity and magnetic fluctuations indicate strong turbulence. As suggested by e.g. PKK14, this turbulent region should be identified with the highly variable ‘anvil’ observed in the Crab nebula (Hester et al. 1995, 2002; Tavani et al. 2011).

To test the hypothesis outlined in Section 2.2 that particles are transported mainly via turbulent flow in the nebula, we apply the Taylor–Green–Kubo (TGK) formulation (e.g. Shalchi & Döring 2007, and references therein) directly to the MHD velocity field:

\begin{eqnarray} \hat{D}_{xy}(\Delta t) \equiv \int _{0}^{\Delta t} {\rm d}\tau \langle (v_{x}(\tau )-\langle v_{x} \rangle ) (v_{y}(0) - \langle v_{y} \rangle ) \rangle , \end{eqnarray}

(15)

where we subtract the background flow velocities 〈v_x,y〉. In the following, we will only study the radial transport, thus x, y = r. Assuming homogeneity in time and purely diffusive behaviour, this definition would be identical to the general form equation (1) if the velocities were taken as particle velocities. On the other hand, as long as the magnetic field remains toroidally dominated, the components of particle drift velocity |$(\boldsymbol {E}\times \boldsymbol {B})\, \boldsymbol {\hat{e}}_{r}$| and flow velocity v_r can be interchanged. It is hence instructive to see how equation (15) compares to the direct measurement from test particles using equation (1).

To obtain convergence of the integral |$\hat{D}_{rr}$|⁠, we need to ensure that the upper bound Δt is chosen larger than the correlation time following from:

\begin{eqnarray} R(\Delta t) = \frac{\langle \left(v_{r}(\Delta t)-\langle v_{r} \rangle \right) \left(v_{r}(0) - \langle v_{r} \rangle \right) \rangle }{\langle v_{r}^{2}\rangle - \langle v_{r} \rangle ^{2}}. \end{eqnarray}

(16)

The right-hand panel of Fig. 5 shows the resulting coefficient for radial diffusion |$\hat{D}_{rr}$| after an integration over 5 yr.² We overplot white contours of R = 1/2 to visually separate correlated from uncorrelated regions. At this time, the flow is uncorrelated in most regions which indicates convergence of |$\hat{D}_{rr}$|⁠. Qualitatively, we obtain good agreement with the measurement from test particles. The region of strongest diffusion is near the poles, diffusion is suppressed at polar angles of 45° and it becomes stronger again in the equatorial region. In the polar jet, we also obtain reasonably good quantitative agreement with the test particle simulation and diffusivities reach the maximum of D_rr ≈ 3 × 10²⁷ cm² s^{− 1}.

Figure 5.

Diffusion resulting from fluctuations of the velocity field according to equation (15) for Δt = 5. White contours indicate the regions where the radial velocity field is still correlated (R = 1/2). Regions where |$\hat{D}_{rr}\le 0$| are masked out. We average over the azimuthal angle and 21 snapshots, corresponding to 20 yr.

However, there are also clear differences in the results obtained via formulations (1) and (15). The total average value |$\langle \hat{D}_{rr}\rangle _{r,\theta }\simeq 1.2\times 10^{26}\, \rm cm^{2}\, s^{-1}$| is almost three times lower than the corresponding test particle result. In particular for the outer regions the flow velocity field tends to underpredict D_rr by an order of magnitude. We suggest that this is related to the generation of a poloidal field component such that the radial particle transport decouples from the flow velocities. Particles can then be transported outwardly also along the field lines, subject only to pitch-angle scattering.

Finally, we investigate the power spectrum of the flow in the nebula proper. Fig. 6 shows the temporally averaged quantity |$P_{\rm K}(k) = \langle \boldsymbol {\hat{u}}\cdot \boldsymbol {\hat{u}}^{{\ast }}\rangle _{t}$| where |$\boldsymbol {\hat{u}} = \boldsymbol {\hat{u}}(k)$| denotes the shell average of the Fourier-transformed four-velocity. Note that the we omit the factor 2π in the wavenumber k ≡ 1/λ and thus the termination shock corresponds to a wavenumber of k_S ≈ 5 × 10^{− 18} cm^{− 1}. In principle we would expect a peak in the power spectrum at the driving scale k_S. Since we did not subtract the mean velocity, it is difficult to see this in the data as the spectrum is contaminated by mean flows on larger scales, e.g. in the jet. It is reassuring to see that the large scales are well converged between the two resolutions. Both resolutions show a notable bump at k = 10^{− 17} cm^{− 1} which is approximately the transverse size of the fast equatorial flow. Apart from this bump, the spectrum on scales smaller than L_S can be approximated by a Kolmogorov 5/3 law. (Relativistic turbulence does show the classical Kolmogorov scaling; Radice & Rezzolla 2013.)

Figure 6.

Velocity power spectrum in the 3D domain for a numerical resolution of Δx = 2 × 10¹⁶ cm (blue) and Δx = 1 × 10¹⁶ cm (red). The variable termination shock radius corresponds to a wavenumber of k_S ≈ 5 × 10^{− 18} cm^{− 1}. To guide the eye, we also show a Kolmogorov 5/3 power law which matches the overall cascade quite well.

3 THE TRANSPORT MODEL

Having derived diffusion coefficients, we now turn our attention to implementing these coefficients in a particle transport model, with the ultimate goal of modelling the X-ray emission observed from three PWNe.

The electrons and positrons responsible for emitting the observed X-ray emission have a relatively short synchrotron lifetime compared to the age of the nebula, with an estimate for the synchrotron lifetime is given by (see e.g. De Jager & Djannati-Ataï 2009)

\begin{equation} t_{\rm {syn}} = 1.3 \times 10^{3} \left(\frac{B}{100 \,\mu \rm G}\right)^{-2} \left(\frac{E}{\rm TeV}\right)^{-1}\, \rm yr. \end{equation}

(17)

For an average magnetic field of B = 300 μG and a 1 TeV electron, the synchrotron lifetime is t_syn ∼ 100 yr. This implies that new particles must continually be injected at the termination shock of the PWN. Consequently one may assume that this particular particle population has reached a steady state. Furthermore, as the evolution of the PWN occurs on time-scales that are larger than the above-mentioned synchrotron lifetime, one may also assume a steady state for the magnetic field, convection velocity, and diffusion coefficient. Lastly, the X-ray emission will be modelled under the assumption that the PWNe have a spherically symmetric morphology. The validity of this last assumption varies between PWNe, and will therefore be discussed when the results for the individual sources are presented.

In a spherically symmetric, steady-state system where convection, diffusion, adiabatic losses, and synchrotron emission are important, the evolution of the lepton spectrum is described by the Fokker–Planck transport equation (see e.g. Ginzburg & Syrovatskii 1964; Parker 1965):

\begin{eqnarray} &&{\kappa \frac{\mathrm{\partial} ^2 f}{\mathrm{\partial} r^2} + \left[\frac{1}{r^2}\frac{\mathrm{\partial} }{\mathrm{\partial} r}\left(r^2\kappa \right)-V\right]\frac{\mathrm{\partial} f}{\mathrm{\partial} r}} \nonumber\\ &&{\quad+ \left[\frac{1}{3r^2}\frac{\mathrm{\partial} }{\mathrm{\partial} r}\left(r^2V\right)+z_{\rm p} p\right]\frac{\mathrm{\partial} f}{\mathrm{\partial} \ln p}+4z_{\rm p}pf = 0,} \end{eqnarray}

(18)

where r is the radial position, p is the particle's momentum, and f(r, p, t) is the omnidirectional distribution function. Furthermore, κ(r) is the diffusion coefficient, V(r) the convection velocity, and

\begin{equation} z_{\rm p}(r)=\frac{4}{3}\sigma _{\rm {T}}\frac{1}{\left(m_{\rm e}c\right)^2}\frac{B^2}{8\pi } \end{equation}

(19)

a coefficient related to the synchrotron loss rate. In the equation above σ_T represents the Thomson cross-section, m_e the mass of the electron, and B(r) the magnetic field. The spatial dependence for B(r), V(r), and κ(r) that are used for the modelling will be discussed in more detail in Sections 3.1 and 3.2.

The electrons and positrons injected at the termination shock are included by employing a suitable boundary condition. The number of particles per momentum interval that flow through the inner radial boundary must be equal to the total number of particles produced per time and momentum interval by the pulsar, Q_s = Q*(p/p*)^{− (α + 2)}, with Q* a normalization constant. This implies that

\begin{equation} \oint \boldsymbol {S}_{\rm s}\cdot {\rm d}\boldsymbol {A}_{\rm s} = Q_{\rm s}, \end{equation}

(20)

where |${\rm d}\boldsymbol {A}_{\rm s}=r_{\rm s}^2\sin \theta \, {\rm d}\theta \,{\rm d}\phi \, \boldsymbol {e}_r$| is the surface element on the inner boundary, and

\begin{equation} \boldsymbol {S} = 4\pi p^2(Vf-\kappa \cdot \nabla f). \end{equation}

(21)

Integrating equation (20) leads to the inner boundary condition,

\begin{equation} CV_{\rm s} f(r_{\rm s},p) - \kappa _{\rm s}\frac{\mathrm{\partial} f(r_{\rm s},p)}{\mathrm{\partial} r} = Q^{\ast }\frac{1}{4\pi r_{\rm s}^2}\left(\frac{p_{\rm s}}{p^{\ast }}\right)^{-(\alpha +2)}, \end{equation}

(22)

that is solved simultaneously with equation (18).

The normalization constant Q* is determined by the requirement that the total momentum (energy) in the particle spectrum injected into the nebula must be some fraction η of the spin-down luminosity L of the pulsar:

\begin{equation} \int _{p_{\min }}^{p_{\max }} Q^{\ast }\left(p/p^{\ast }\right)^{-(\alpha +2)} p^3 {\rm d}p = \frac{\eta L}{c}. \end{equation}

(23)

The value of η is left as a free parameter, but it should be noted that this quantity only controls the normalization of the particle spectrum, and therefore has no influence on the spatial evolution of either the synchrotron photon index Γ or the X-ray flux. Consequently this parameter is only of secondary importance.

Observations indicate that PWNe consist of two particle populations (see e.g. De Jager & Djannati-Ataï 2009, and references therein): the first is responsible for producing the radio synchrotron emission, whereas the second population is responsible for producing the X-ray emission. The present modelling is concerned with only the latter population, and p_min in equation (23) therefore represents the minimum particle momentum of the second population. As the value of p_min is not constrained by observations, p_minc = 0.1 TeV is used, comparable to the values p_minc = 0.08–0.15 TeV derived by Zhang, Chen & Fang (2008) for three PWNe.

To calculate the maximum particle momentum p_max, or alternatively the maximum energy E_max, two possible constraints can be used. The first is that the rate of acceleration of a particle at a shock will eventually be balanced by synchrotron losses, with De Jager et al. (1996) deriving the limit

\begin{equation} E_{\max } = \gamma _{\max } m_{\rm e}c^2 = 44 B_{\rm s}^{-1/2}\, {\rm erg}, \end{equation}

(24)

where γ_max is the maximum Lorentz factor of the accelerated particles and B_s is the magnetic field at the termination shock. Additionally, particles can only be accelerated as long as they are confined within the shock, leading to the gyroradius limit (see e.g. De Jager & Djannati-Ataï 2009)

\begin{equation} E_{\max }=q_{\rm e} B_{\rm s} r_{\rm s}, \end{equation}

(25)

where r_s is the radius of the termination shock. For the modelling the value of E_max is chosen as the smaller of the two limits.

These limits are based on the assumption that the X-ray emitting electrons and positrons are accelerated through the Fermi mechanism. One well-known characteristic of this process is that the maximum index of the particle spectrum is α = 2, with this value chosen for the modelling.

For the outer boundary it is assumed that the particles can freely escape from the PWN. A more realistic model would allow for the confinement of particles, but also requires the use of a time-dependent model. This was investigated by Vorster (2013), who found that both boundary conditions lead to identical results. This makes sense if it is kept in mind that the X-ray emitting particles have a very short lifetime, and one would therefore not expect the number of these particles to increase with time near the outer boundary of the PWN.

The transport equation (18) is solved numerically using the second-order accurate, finite-difference, Crank–Nicolson scheme. For more details on this model, Vorster & Moraal (2013) and Vorster (2013) can be consulted, where it is also shown how the various energy loss and transport processes influence the spatial evolution of the particle spectrum.

As will be discussed in Sections 3.1 and 3.2, the model has a minimum of four free parameters (B_s, V_s, K_s, and η), and it would thus be too computationally intensive to investigate the whole possible parameter space. Rather, the search method developed by Nelder & Mead (1965) is used to find a set of parameters that result in a minimum χ² value. In order to ensure that the minimum found is not a local minimum, the search method is started at three different positions in the possible parameter space.

Lastly, for completeness it should be noted that the X-ray data were extracted from either circular or annular regions of increasing size, and that the emission observed from the inner circular regions of the nebula will contain some fraction of the emission emitted from regions further out. This is taken into account when modelling the data by performing a line-of-sight integral on the emission predicted by the model.

3.1 The KC84 model

This model is based on the well-known results of KC84. In this model the magnetic field B_s and the advection velocity V_s at the termination shock are determined by the ratio of electromagnetic to particle energy σ. This ratio also determines the spatial dependence of both B and V. However, this model does not provide a direct prediction for the diffusion coefficient. Furthermore, as a toroidal magnetic field is assumed, diffusion will necessarily occur perpendicular to the magnetic field. Based on this consideration, the choice of diffusion coefficient is motivated as discussed below.

It has long been known that the perpendicular diffusion coefficient for charged particles in a collisionless plasma must in some way depend on their Larmor radius, and thus implicitly on the magnitude of the background magnetic field (see e.g. Chen 1974). The nature of this magnetic field dependence is, however, not perfectly clear. Many studies assume Bohm diffusion, where κ ∝ B⁻¹ (see equation 2). This relation has been found, through numerical simulations of perpendicular diffusion coefficients, to be a reasonable choice in the presence of strong turbulence (Hussein & Shalchi 2014). The field line random walk (FLRW) diffusion coefficient represents another limiting case. Here it is assumed that particles exclusively follow magnetic field lines, which themselves wander due to turbulence. The perpendicular diffusion coefficient in this case scales as B⁻², and is energy independent (see e.g. Shalchi 2009, and references therein). Much theoretical work on perpendicular diffusion has been done in recent years, taking into account more thoroughly the effects of turbulence on perpendicular scattering. These recent theories, being for the most part refinements of the non-linear guiding centre (NLGC) theory proposed by Matthaeus et al. (2003), tend to yield results that are relatively intractable functions of, amongst other turbulence quantities such as the magnetic variance, the parallel mean free path (see e.g. Shalchi, Li & Zank 2010; Engelbrecht 2013; Engelbrecht & Burger 2015). Shalchi, Bieber & Matthaeus (2004), however, derive approximate, analytical expressions for the NLGC perpendicular mean free path for various limiting scenarios, including the FLRW limit described above, which have been used in heliospheric cosmic ray modulation studies (see e.g. Burger et al. 2008). These expressions require as input the parallel mean free path. Assuming the quasi-linear result for this quantity (see e.g. Burger et al. 2008), this perpendicular mean free path will scale as B^−4/3 at the high energies of interest to this study, and remain energy independent.

Based on the above considerations, two scenarios will be investigated. The first assumes a B⁻¹ dependence for the diffusion coefficient [hereafter referred to as the KC84(a) scenario], and the second a B⁻² dependence [hereafter referred to as the KC84(b) scenario]. These values represent limiting cases, with a more realistic case expected to fall somewhere in between. Both scenarios additionally assume energy-independent coefficients. The value of the diffusion coefficient at the termination shock κ_s is left as a free parameter, while κ is limited by equation (7).

3.2 The PKK14 model

This model is based on the MHD simulations of PKK14. To implement these simulations into the spherically symmetric transport model, B and V (predicted by the MHD model) are averaged over time for a number of spherical shells, with the resulting averages shown in Fig. 7. The averaged data are fitted with the Gaussian function:

\begin{equation} G(r) = \mathcal {A} \exp \left(-\frac{r^2}{2\mathcal {B}^2}\right)+\mathcal {C}, \end{equation}

(26)

with this function having the advantage that it can easily be incorporated into the transport model. The fits are performed from the termination shock at r_s ≈ 0.15r_n out to the nebula radius r_n. It yields the parameters |$\mathcal {A} = 1.536c$|⁠, |$\mathcal {B} = 1.56727 \times 10^{-2}c$|⁠, and |$\mathcal {C} = 0.1005$| for the velocity, and |$\mathcal {A} = 5.239 \times 10^{-4}\,{\rm G}$|⁠, |$\mathcal {B} = 7.241 \times 10^{-5}\,{\rm G}$|⁠, and |$\mathcal {C} = 0.2059$| for the magnetic field. Fits using power-law profiles were also investigated, but it was found that the Gaussian profiles lead to an improvement in the fit.

Figure 7.

Time-averaged profiles of the magnetic field magnitude and radial velocity using the model ‘B3D’ of PKK14. A Gaussian was fitted to profiles averaged from four snapshots at t = 60, 61, 62, 63 yr after the start of the simulation.

As the PWNe that are modelled have varying ratios of r_s/r_n, the following characteristics are kept constant for the different PWNe: (1) the value of B_s is chosen as a free parameter; (2) the magnetic field decreases by a factor of ∼5 within approximately three termination shock radii, before reaching a constant value; (3) the velocity at the inner boundary is V_s = 0.51c; (4) the velocity decreases rapidly within about two termination shock radii, before reaching a constant value; and (5) the constant velocity after two termination shock radii is chosen to match the observed expansion speed of the PWNe.

One aspect of the PKK14 model that has to be investigated is the effect of the rapid decrease in V in the inner regions of the PWN, as this will lead to adiabatic heating. This requires the use of a time-dependent model that is able to take into account both energy losses and gains. Using the time-dependent transport model developed by Vorster & Moraal (2013) a large number of scenarios with varying parameters were investigated, and it was found that adding the effect of adiabatic heating only changes the normalization of the spectra. However, the same spatial evolution is predicted, regardless of whether the adiabatic heating in the inner regions of the PWN is included or neglected. Ideally one would do the fitting with this time-dependent model. However, compared to the steady-state model it is not only significantly more computationally intensive, but it is also difficult to obtain the same numerical accuracy as the steady-state model. The fitting will therefore be done using the steady-state model, and neglecting adiabatic heating in the inner regions.

4 THE X-RAY PULSAR WIND NEBULAE

The three PWNe that are chosen for the modelling are the two young sources G21.5−0.9 and 3C 58, as well as the inner regions of the Vela PWN. The parameters derived for the nebulae are listed in Table 1. In order to better compare the results between the models and the selected PWNe, spatially averaged values have also been calculated for the various parameters, and are listed in the second part of the table. The last line of the table lists the average Péclet number (see e.g. Prandtl 1953), defined as

\begin{equation} \xi = \frac{Vr}{\kappa }. \end{equation}

(27)

This dimensionless number is an indication of the importance of advection relative to diffusion. When ξ ≫ 1 the system is advection dominated, while being diffusion dominated for ξ ≪ 1.

Table 1.

Values derived for the free parameters: KC84 represent values found using the model of Kennel & Coroniti (1984a) – Scenario (a) corresponds to the case κ ∝ B⁻¹ and Scenario (b) to the case κ ∝ B⁻² – and PKK14 the values found using the model of PKK14. The first part of the table lists the values of the various parameters at the termination shock, and the second part values averaged over the PWN.

	G21.5−0.9			Vela			3C 58
Parameter	KC84(a)	KC84(b)	PKK14	KC84(a)	KC84(b)	PKK14	KC84(a)	KC84(b)	PKK14
Bs (μG)	33	33	283	264	264	38	8	8	300
Vs (units of c)	0.36	0.36	0.51	0.52	0.52	0.51	0.35	0.35	0.51
κs (1026 cm2 s−1)	8.5	6.0	5.7	0.5	0.5	1.4	27.8	17.1	13.3
σ (10−3)	1.3	1.3	–	142	143	–	0.6	0.6	–
η (10−2)	8.8	22.5	4.5	0.3	0.3	2.1	–	–	–
|$\bar{B}$| (μG)	158	158	43	30	30	5.8	63	63	46
|$\bar{V}$| (10−3, units of c)	4.2	4.2	3.1	159	159	3.3	4.5	4.5	2.6
|$\bar{\kappa }$| (1026 cm2 s−1)	1.2	0.3	5.7	1.1	0.9	1.4	1.8	0.3	13.3
|$\bar{\xi }$|	2.0	9.7	0.3	129	186	2.1	2.0	15	0.2

4.1 G21.5−0.9

With a spin-down luminosity of L = 3.3 × 10³⁷ erg s⁻¹ (Camilo et al. 2006), the pulsar in the supernova remnant G21.5−0.9 is one of the most energetic pulsars in the Galaxy. The PWN is located at a distance of 4.8 kpc (Tian & Leahy 2008), with an estimated age of 870 yr (Bietenholz & Bartel 2008). X-ray observations (Slane et al. 2000; De Rosa et al. 2009; Tsujimoto et al. 2011) reveal a bright, highly spherical nebula with a radius of r_pwn ∼ 40 arcsec, while Slane et al. (2000) estimated a value of r_s ≳ 1.5 arcsec for the radius of the termination shock. In order to keep the number of free parameters to a minimum, the value r_s = 1.5 arcsec is used for fitting the data.

Bietenholz & Bartel (2008) measured that the PWN is expanding at a velocity of V_pwn = 910 ± 160 km s⁻¹. With the assumed value of r_s, the KC84 model predicts that the above-mentioned expansion velocity in the outer regions of the PWN is compatible with the range of values 1.0 × 10⁻³ ≤ σ ≤ 1.6 × 10⁻³. For the modelling the intermediate value σ = 1.3 × 10⁻³ is chosen, leading to magnetic field strength of B_s = 33 μG at the termination shock. This range of σ values is also comparable to the value σ = 3 × 10⁻³ derived for the Crab nebula (Kennel & Coroniti 1984b). For the PKK14 model, the value V_pwn = 910 km s⁻¹ is used.

The data that are modelled is the 2–8 keV observations reported by Tsujimoto et al. (2011). For these observations data were extracted from circular regions of increasing size, with the centre of each region placed on the position of the pulsar. The observations therefore represent cumulative data. A 3–45 keV observation for the region r ≤ 30 arcsec has also recently been reported by Nynka et al. (2014). The authors found statistically significant evidence for a spectral break at ∼9 keV, deriving Γ = 1.852 ± 0.0011 for the 3–9 keV energy range, and |$\Gamma =2.099^{+0.019}_{-0.017}$| for the 9–45 keV energy range. However, the 3–9 keV measurement is not entirely in agreement with the results of Tsujimoto et al. (2011), who found Γ = 1.78 ± 0.02 for the same region. Despite this discrepancy, the single 9–45 keV observation of Nynka et al. (2014) was taken into account for the modelling, as it was found that this does help to constrain the possible free parameters. As such the single observation should have only a limited influence on the χ² value calculated for the fit.

The best-fitting model prediction using the KC84 model is shown in Fig. 8, with the best-fitting values listed in Table 1. The figure shows that the KC84 model is not able to reproduce the evolution of Γ. Additionally, there is also little variation in the results when a diffusion coefficient is used that scales as either κ ∝ B⁻¹ or κ ∝ B⁻² [henceforth referred to as KC84(a) and KC84(b)]. The KC84 model is however capable of modelling the evolution of the flux to some degree, with the largest deviation between the model prediction and the data occurring in the inner parts of the nebula. The average magnetic field derived from the fit is |$\bar{B}=158\,\mu {\rm G}$|⁠, comparable to the equipartition value B = 180 μG estimated by Safi-Harb et al. (2001), and used by Tang & Chevalier (2012) in their modelling of G21.5−0.9. However, using both X-ray and very high energy gamma-ray observations, de Jager, Ferreira & Djannati-Ataï (2008a) derived the lower value |$\bar{B}=25\,\mu {\rm G}$|⁠, comparable to the values |$\bar{B}=47$| and 64 μG derived by Tanaka & Takahara (2011) using a spatially independent transport model.

Figure 8.

G21.5−0.9: fit to the data using the magnetic field and velocity profiles from the model of Porth et al. (2014a, PKK14), and using the profiles from the model of Kennel & Coroniti (1984a, KC84). The X-ray data are taken from Tsujimoto et al. (2011) and corresponds to the 2–8 keV energy range.

The average diffusion coefficients |$\bar{\kappa }=1.2\times 10^{26}$| and 2.5 × 10²⁵ cm² s⁻¹, corresponding to the KC84(a) and KC84(b) scenarios, respectively, are smaller than the value |$\bar{\kappa }=2\times 10^{27}\,{\rm cm}^2\,{\rm s}^{-1}$| derived by Tang & Chevalier (2012). The average Péclet numbers |$\bar{\xi }=2.0$| and 9.7, for the KC84(a) and KC84(b) scenarios, respectively, show that convection is the more dominant particle transport process. As such one would expect advection to be more dominant for the KC84(b) scenario compared to the KC84(a) scenario: B increases as a function of radius for the largest part of the nebula, and therefore the B⁻² scaling implies that κ should decrease faster than for the KC84(b) scaling. However, Fig. 8 shows that the different scalings used for κ have only a small influence on the results, as one would expect in a scenario where advection is more dominant.

Fig. 8 also shows the fit obtained using the PKK14 model, with the best-fitting parameters listed in Table 1. The figure shows that the PKK14 model can reproduce the evolution of Γ significantly better than the KC84 model, but has a problem in predicting the correct flux in the inner regions of the nebula. This is not necessarily surprising as the flow and magnetic field have a complexity (see Fig. 2) that is not taken into account in the spherically symmetric transport equation. The average magnetic field derived from the fit is |$\bar{B}=43\,\mu {\rm G}$|⁠, comparable to the value derived from the KC84 model. However, it may be argued that the value B_s = 283 μG predicted by the PPK14 model is more realistic given the large luminosity of the pulsar powering the G21.5−0.9. The average diffusion coefficient |$\bar{\kappa }=5.7\times 10^{26}\,{\rm cm}^2\,{\rm s}^{-1}$| is once again limited by equation (7). As is the case for the KC84 model, the PKK14 model predicts that diffusion should be the dominant transport process, with the value |$\bar{\xi }=0.34$| calculated from the fit. Lastly, the PKK14 model predicts that Γ = 1.89 in the 9–45 keV energy range.

To illustrate the role of diffusion, it is useful to fit the data with diffusion neglected. The results of this comparison are shown in Fig. 9. As such it is difficult to make any conclusions based on the results from the KC84 model. Not only does the model fail to fit the evolution of Γ, but the values ξ > 1 in Table 1 predict that convection should be more important that diffusion.

Figure 9.

The same as Fig. 8, but with diffusion neglected.

Within the context of the PKK14 model, comparison of Figs 8 and 9 shows that diffusion leads to a moderately better fit to Γ, but a significantly better fit to the evolution of the flux. However, there are two important aspects to keep in mind. The first is that, in order to improve the fit to the data, the spectral index injected at the shock was changed to α = 2.2, i.e. neglecting diffusion can to some degree be compensated for by changing certain parameters. Secondly, and more importantly, the average magnetic field in the absence of diffusion is |$\bar{B}=6.9\,\mu {\rm G}$|⁠, lower than the prediction of any previous model. The consequence is that the conversion efficiency of the pulsar's spin-down luminosity to particle energy must exceed 100 per cent, i.e. η > 1.

The discussion presented in the previous paragraph also implicitly illustrates an important consequence of diffusion. In order for electrons and positrons to reach the outer regions of G21.5−0.9, the transport time-scale has to be shorter than the synchrotron lifetime of the particles. Neglecting diffusion increases the transport time-scale, and subsequently the synchrotron lifetime of the particles has to decrease, implying that |$\bar{B}$| has to become smaller. Thus, with diffusion included it becomes possible for the particles to propagate to distances from the pulsar that would otherwise not have been possible due to synchrotron losses.

4.2 The inner regions of the Vela PWN

The Vela PWN, located at a distance of 290 pc (Dodson et al. 2003), is a highly asymmetric PWN, as shown in X-ray (Markwardt & Ögelman 1995) and very high energy gamma-ray observations (Aharonian et al. 2006), with the age of the nebula estimated to be 11 400 yr (Manzali, De Luca & Caraveo 2007). Although the Vela PWN is a morphologically complex source with an extension of ∼1°, 3–10 keV observations (Mangano et al. 2005) of the inner region show a bright X-ray nebula that is approximately spherically symmetric. X-ray observations also show a double torus and jet-like structure at the centre of the PWN (Helfand, Gotthelf & Halpern 2001), commonly associated with the radius of the termination shock. Using a geometrical model to fit the double torus, Ng & Romani (2004) derived a value of r_s ∼ 21 arcsec. The aim is not to model the entire PWN, but only the above-mentioned inner region. The 3–10 keV data (Mangano et al. 2005) that is modelled was again extracted from circular regions of increasing size up to a radius of r = 15 arcmin.

As the velocity V is unknown, the value of σ, and consequently B_s, is now an additional free parameter in the KC84 model. For the PKK14 model, it is assumed the velocity reaches the constant value V = 1000 km s⁻¹ after two termination shock radii (see also discussion in Section 3.2), comparable to the value measured for G21.5−0.9.

Fig. 10 shows the fits to the data using the KC84 model. It can be seen that the model can fit the Γ data moderately well, while providing a good fit to the flux data for the largest part of the nebula. However, there is a discrepancy between the model prediction and flux data in the outer regions of the nebula. Both κ scenarios (κ ∝ B⁻¹ and κ ∝ B⁻²) predict an average magnetic field of |$\bar{B}=30\,\mu {\rm G}$| for the inner region of the nebula, larger than the average magnetic field of |$\bar{B}=5\,\mu {\rm G}$| estimated by De Jager, Slane & LaMassa (2008b) for the whole PWN. The best-fitting value σ = 0.14 × 10⁻¹ is also in agreement with the range σ = 0.05–0.5 estimated by Sefako & de Jager (2003), and comparable to the range σ = 0.01–0.1 estimated by Bogovalov et al. (2005). Average diffusion coefficients of |$\bar{\kappa }=1.1\times 10^{26}$| and 8.7 × 10²⁵ cm² s⁻¹ are predicted for the κ ∝ B⁻¹ and κ ∝ B⁻² scenarios, respectively. These values correspond to Péclet numbers of |$\bar{\xi }=130$| and 186, indicating that advection is the dominant process.

Figure 10.

Inner regions of Vela: fit to the data using the PKK14 and KC84 models. The X-ray data are taken from Mangano et al. (2005) and correspond to the 3–10 keV energy range.

A salient feature of Fig. 10 is the similarity between the two κ scenarios in the KC84 model. The large value of σ implies that B decreases in the inner regions of the PWN, before approaching a constant value. The constant value therefore implies that the radial dependence of κ will also be constant, regardless of the scaling used.

Apart from the 3–10 keV data, Mangano et al. (2005) also extracted a single data point from a r = 15 arcmin circular region in the 20–100 keV energy range. The photon index and flux was measured as Γ = 2.00 ± 0.05 and F_X = 15.9 ± 0.4 erg cm⁻² s⁻¹, respectively. This data were also taken into account for the modelling, with the KC84 model fit in Fig. 10 predicting Γ = 1.78 and a flux of F_X = 10.6 erg cm⁻² s⁻¹ for both diffusion coefficients.

Fig. 10 also shows that while the PKK14 model can produce a better fit to the Γ data, compared to the KC84 model, the former model has trouble in predicting the correct flux in the inner regions of the nebula. A smaller average magnetic field of |$\bar{B}=5.8\,\mu {\rm G}$| and Péclet number |$\bar{\xi }=2.1$| is found for the PKK14 model, although a diffusion coefficient |$\bar{\kappa }=1.4\times 10^{26}\,{\rm cm}^2\,{\rm s}^{-1}$| comparable to those of the KC84 model is predicted. Lastly, the PKK14 model predicts Γ = 2.01 and a flux of F_X = 17.2 erg cm⁻² s⁻¹ in the 20–100 keV energy range. Note that this is an improvement over the prediction made by the KC84 model (cf. previous paragraph).

4.3 3C 58

The PWN 3C 58, located at a distance of 3.2 kpc (Roberts et al. 1993), is generally associated with the supernova SN 1181, implying an age of 833 yr (see e.g. Slane et al. 2004). X-ray observations in the 2.2–8 keV energy range not only reveal a nebula with an elliptical morphology similar to that of the Crab nebula, but also a show jet-like and a possible torus structures in the inner regions of the PWN (Slane et al. 2004). The torus structure is elongated, and extends between 2.5 and 8.0 arcsec from the pulsar (Slane et al. 2004), indicating that the torus is inclined with respect to the plane of the sky. To derive an estimation of the termination shock radius requires the use of a geometrical model similar to the one presented by Ng & Romani (2004). As this is beyond the scope of the present modelling, the average value r_s = 5.25 arcsec is used.

Radio measurements show that the nebula is expanding at a velocity of V_pwn ∼ 910 ± 360 and ∼550 ± 220 km s⁻¹ along the major and minor axes of the nebula, respectively. The X-ray measurements were extracted from circular annuli centred on the pulsar, and the average expansion velocity V_pwn = 730 km s⁻¹ is thus used, implying σ = 5.55 × 10⁻⁴ and B_s = 8.4 μG for the KC84 model. Unfortunately normalized flux measurement is not available, and the fitting was therefore done based on the spatial variation of only Γ.

Fig. 11 shows the KC84 model can reproduce the Γ observations in the very inner regions of the PWN, but fails to reproduce the data for the largest part of the nebula. By contrast, the PKK14 model is able to predict the evolution of Γ. The value of |$\bar{B}=46\,\mu {\rm G}$| predicted by the PKK14 model is comparable to a previous estimate |$\bar{B}=40\,\mu {\rm G}$| derived by Tanaka & Takahara (2013), but larger than the value of |$\bar{B}=14\,\mu {\rm G}$| estimated by the MAGIC collaboration following their detection of 3C 58 at TeV gamma-ray energies (Aleksić et al. 2014). Note that the above-mentioned values are all smaller than the estimated equipartition value |$\bar{B}=80\,\mu {\rm G}$| (Green & Scheuer 1992).

Figure 11.

3C 58: fit to the data using the PKK14 and KC84 models. The X-ray data are taken from Slane et al. (2004) and correspond to the 2.2–8 keV energy range.

5 DISCUSSIONS AND CONCLUSIONS

We have adopted a MHD description of the PWN system with dynamics governed by particles with mean free paths much smaller than the system size. The overall success of MHD models in reproducing the dynamics, emission, and individual features of PWNe (e.g. Volpi et al. 2008; Camus et al. 2009; Olmi et al. 2014; PKK14) indicates that this assumption can be justified in retrospect. On top of the MHD flow, test particles representing the high-energy tail of the distribution functions are integrated, given the electric and magnetic fields of the fluid. It is worthwhile to note that in flat spectrum sources like PWN, the energetics is dominated by TeV particles (e.g. section of Beskin et al. 2015) with a Lorentz factor only an order of magnitude smaller than the lowest energy test particles we have considered. Because of the lack of self-consistency, we have focused only on the spatial particle transport and leave transport in momentum space (in situ particle acceleration) for future studies.

We find that turbulence in PWNe gives rise to high diffusive particle transport. This idea has been previously employed by Tang & Chevalier (2012), who found that the synchrotron spectral index of several young PWNe could be modelled well by relying on diffusive transport alone (i.e. neglecting advection). In the model of Tang & Chevalier (2012) the diffusion coefficient κ was left as a free parameter, whereas the present paper aims to constrain κ using three-dimensional MHD simulations of PWNe (PKK14). The MHD results are then used as input in a particle transport model, with the ultimate goal of modelling the X-ray emission from G21.5−0.9, the inner regions of the Vela PWN, and 3C 58.

Comparing the results from MHD simulations that contain an integrated test-particle component with the diffusion coefficient derived from a formalism based on correlations of the flow velocity field, we conclude that MHD turbulence gives rise to high diffusivities that are important for particle transport. For typical parameters of young PWNe, we obtain radial diffusion coefficients in excess of 10²⁶ cm² s^{− 1}. Motivated by these results, we suggest a simple scaling relation between the termination shock size L_s and the diffusion coefficient in the PWN: D_rr ∝ v_fL_s, where v_f ≈ 0.5c is a typical flow velocity at the driving scale corresponding to L_s (see equation 7). Both in our simulated PWN and in real systems, the driving is localized in the centre of the nebula, with the turbulence decaying as one moves away from the centre. This would lead to a lower radial transport in the outer regions of the nebula, which is however compensated for by the formation of a radial component in the magnetic field.

The transport is split into advection with the mean flow plus a diffusive contribution due to advection in the turbulent flow component. While both contributions are ‘advective’, only the mean flow is associated with adiabatic losses in the modelling. Not surprisingly, we do not obtain a significant energy dependence in the diffusion coefficient for particles with gyroradii smaller than L_s. Thus for average field strengths of B₀ < 300 μG and shock radii L_s > 0.1 Ly, we can safely neglect any energy dependence even for particles emitting MeV synchrotron radiation. The outward particle transport is clearly non-Bohm.

Diffusive transport is most efficient in the kink unstable polar jet where D_rr > 10²⁷ cm² s^{− 1}, marking good agreement between the test-particle approach and results based on correlations of the flow field alone. By contrast, the corresponding Bohm diffusion coefficient for X-ray emitting particles is orders of magnitude lower, on the order of 5 × 10²³ cm² s^{− 1} in a 100 μG field. In the outer equatorial and mid-latitude regions we find stronger discrepancies between the two approaches where the flow field underestimates the diffusive transport. We suggest that the development of a poloidal component to the magnetic field in these outer regions is responsible for the increased efficiency of diffusive transport.

To model the observed spatial evolution of the X-ray photon spectra of the three PWNe G21.5−0.9, the inner regions of Vela, and 3C 58, a steady-state, spherically symmetric transport model that includes convection, diffusion, adiabatic cooling, and synchrotron losses is used. This model has as input time averaged radial profiles for the magnetic field, velocity, and diffusion coefficient that are obtained from the MHD model. In order to better understand the results, the modelling was also done using the magnetic field and velocity profiles from the well-known MHD model of KC84. As diffusion is not a feature of the latter model, two scenarios were considered. In the first κ ∝ B⁻¹, and in the second κ ∝ B⁻². The model of KC84 is also useful to derive σ, the ratio of magnetic to particle energy in the nebula (in the MHD model of PKK14, σ = 1).

It was found that using the magnetic and velocity profiles of PKK14 allows one to find reasonable fits to the evolution of the photon index Γ, whereas using the profiles from the KC84 model leads to fits that are, in general, not as good. However, the model of KC84 does lead to somewhat better fits to the evolution of the X-ray synchrotron flux. However, taking into account the fits to both the spectral index and flux, the MHD model of PKK14 is the preferred model.

Using the PKK14 model, it is found that diffusion is the more important transport mechanism in G21.5−0.9 and 3C 58, whereas advection is the more important mechanism in the inner regions of the Vela PWN. However, for all three sources the Péclet number is close to unity, indicating that both advection and diffusion play an important role in particle transport.

One noteworthy point is that while the model of KC84 has trouble in fitting the spectral indices of G21.5−0.9 and 3C 58, the model does allow one to find a reasonable fit to the photon index and flux data of the inner regions of the Vela PWN. From these fits the value σ = 0.14 is derived, significantly higher than the value of σ = 3 × 10⁻³ derived for the Crab nebula (Kennel & Coroniti 1984b), with this latter value often taken as a canonical value for PWNe. However, the larger σ value for Vela is not entirely unexpected as Sefako & de Jager (2003) and Bogovalov et al. (2005) have estimated comparable values.

Recently, Cerutti, Philippov & Spitkovsky (2016) presented global particle in cell (PIC) simulations of pulsar magnetospheres with the aim to model the pulsar light curve. Although the scientific objective was very different, the rapid progress in global kinetic models prompts us to evaluate the feasibility of global PIC simulations of pulsar wind nebulae. The advantages of the PIC methodology are obvious as they allow to study particle acceleration and transport in collisionless plasma in a self-consistent way.

In general, explicit PIC algorithms need to resolve scales as small as the electron plasma skin depth |$c/\omega _{{\rm c}}=\sqrt{\Gamma _{\rm p}m_{{\rm e}}c^{2}/4\pi e^{2} n_{{\rm e}}}$| and preferentially the Larmor radius r_g of the particles (equation 11). The requirement on temporal scales is analogue, with the gyrofrequency ω_g = c/r_g of cold particles typically being the most severe constraint. The skin depth and Larmor radius can be expressed as |$(c/\omega _{{\rm p}})^{2}=r_{{\rm g}}^{2}\sigma$| (Sironi & Spitkovsky 2011), where σ denotes the usual wind magnetization. For a realistic simulation of the entire (downstream) Crab termination shock, we need to resolve the Larmor radius of TeV leptons that carry most of the energy. Thus for a 100 μG field strength, the termination shock spans approximately 10⁴r_g of this species, or |$10^{4} /\sqrt{\sigma }$| skin depths. Given that current simulations already cover hundreds of skin depths (e.g. Sironi & Spitkovsky 2014), global PIC simulations of relativistic particle transport downstream of the termination shock could soon become feasible. However, global models of particle acceleration at the cold wind shock seem far out of reach due to lacking dynamic range. Ideally, they would need to resolve also the magnetic stripes spaced one light cylinder apart (thus smaller than the termination shock by ∼10⁸).

For further studies with the MHD and test particle approach, we suggest performing simulations with a larger dynamic range between the driving scale, domain size, and the dissipation scale in order to obtain better statistics on the turbulent flow. It is also important to check the results against advanced cosmic ray turbulence theories developed for the Solar system (e.g. Breech et al. 2008; Zank et al. 2012) as the presence of subsonic motion in the PWN and large relativistic velocity fluctuations could render them quite different. Apart from spatial diffusion, future studies will aim to investigate diffusion in momentum space with special attention to Fermi-II type acceleration in the nebula proper.

This work was partially funded through the STFC under the standard grant ST/I001816/1, the ERC Synergy Grant ‘BlackHoleCam’ (Grant 610058), and the NSF grant number 1306672.

1

https://gitlab.com/mpi-amrvac/amrvac

2

In principle the integration time should be chosen as large as possible. However, we note that convergence is lost for larger integration times, most likely due to finite box effects.