Modelling spectral evolution of pulsar wind nebulae inside supernova remnants

Abstract

1 INTRODUCTION

Pulsars (PSRs) are known to release most of their rotational energy in the form of a relativistic magnetized wind, whose particle component is mostly made of electron–positron pairs, with, possibly, a small number of ions, while the magnetic field is almost purely toroidal far away from the light cylinder (e.g. Goldreich & Julian 1969). The wind is initially confined by the slowly expanding supernova remnant (SNR), which leads to the formation of a termination shock (TS) at some distance from the PSR and of a bubble of relativistically hot fluid beyond it. This is the so-called pulsar wind nebula (PWN), which shines through synchrotron (SYN) and inverse-Compton (IC) emission from radio wavelengths up to gamma rays.

The typical energy that a PSR injects into the PWN during its entire lifetime is about 10⁻² of the typical energy of the supernova (SN) explosion (10⁵¹ erg). Therefore, the presence of an energetic PSR has little effect on the global evolution of the SNR. In contrast, the evolution of the PWN strongly depends on its interaction with the SNR (see Section 3).

Modelling of the evolution and dynamical properties of PWN-SNR systems has improved considerably since the pioneering 1D analytical works by Rees & Gunn (1974) and Kennel & Coroniti (1984a). In recent years, the impressive advances in X-ray observations, due to Chandra and XMM–Newton, have led to a renewed interest in these objects. A major effort has been devoted to explaining the observed X-ray properties of the Crab nebula and other young systems through multidimensional time-dependent approaches to their dynamics (Del Zanna, Amato & Bucciantini 2004; Komissarov & Lyubarsky 2004; Del Zanna et al. 2006). Several authors (Blondin, Chevalier & Frierson 2001; van der Swaluw et al. 2001; Bucciantini et al. 2004b; de Jager, Ferreira & Djannati-Ataï 2008) have presented numerical evolutionary models, extended to advanced ages (typically up to ∼10⁵ yr) and to the case of fast moving PSRs (van der Swaluw, Downes & Keegan 2004), but in fact fully multidimensional 3D and/or magnetohydrodynamics (MHD) models are still lacking in these cases. The reason for this is that detailed modelling of older systems is far more difficult and complex than for the young ones. To this one must add the fact that, with respect to a younger system, investigating older ones would require sampling a much larger parameter space. The latter are brighter, hence providing us with high-quality multifrequency data, and are often better constrained in terms of the most relevant physical quantities, such as age, distance and PSR energetics.

The overall evolution of old objects is much more complex than the smooth expansion appropriate to describe young systems. In general, even within the most simplistic approach, the long-term evolution of a PWN-SNR system depends on many parameters: the SN energy, the mass of the ejecta, their density structure, the density of the interstellar medium (ISM), the PSR luminosity and its spin-down history (Bucciantini et al. 2004b). A model for particle injection and spectral evolution must be added to all of this, if the aim is that of deriving information from direct comparison with observations. At present, multidimensional numerical simulations are not suitable to properly investigate the large parameter space that older systems may span.

The situation is actually delicate also for many young systems. Indeed, accurate multidimensional results, including emission maps to be directly compared with observations, have been so far presented only for the Crab Nebula. But this is a very special object, for which measurements are available throughout almost the entire electromagnetic spectrum, and almost all the aforementioned parameters are known, including the age and the PSR rotational history (initial spin frequency and braking index). For other young systems, such as 3C 58 or MSH 15−52, the quantity and quality of information available are much poorer than in the case of Crab and simplified evolutionary models are again the only possible tool for investigation, at least for a start.

Models for the evolution of the emission properties of PWNe have a long history, from the work by Pacini & Salvati (1973) to the most recent developments by de Jager et al. (2008) and Gelfand, Slane & Zhang (2009). These latter works have attempted at taking into account both the dynamics of PWN-SNR systems, as derived from accurate numerical simulations, and the spectral evolution of the particle distribution function. In particular, the work by Gelfand et al. (2009) is exemplary for showing the richness of behaviours that is found even within the one-zone treatment (where all quantities, from the dynamical ones, such as pressure and magnetic field, to the particle distribution function, are assumed to be uniform inside the nebula) and the degree of complexity that one can investigate using simple tools.

While constraining the unknown parameters is preliminary to detailed multidimensional modelling of the systems, which can only sample a much more limited range of parameters, the following question arises: how strong and reliable the results obtained with a one-zone approach are? The recent success of 2D models of PWNe at reproducing so many of the observed features has interesting consequences in this respect in that it lends support to the one-zone description over the 1D approximation. The latter inevitably leads to onion-like structures (Kennel & Coroniti 1984b), where particles injected at earlier times are located at larger radii, with a one-to-one correspondence between the properties of the particle distribution function and the distance from the central PSR. Comparison of multidimensional models with observations has shown that such ordered structure is not what is realized in PWNe, where an overall turbulent flow is injected at the TS and survives throughout most of the nebular extent (Camus et al. 2009). For the typical flow times, which are much smaller than the PWN age, this turbulent flow structure leads to efficient mixing, which is what a one-zone description is based on. This also finds some direct hints from observations where detailed data are available, as in the case of the Crab nebula: the remarkable homogeneity of the radio spectral index across the nebula and the extent of the X-ray emitting region which exceeds the predictions of any simple 1D flow model in the absence of particle diffusion far more efficient than the Bohm rate (Amato et al. 2000).

In the following we present a simplified evolutionary model, similar in spirit to the work by Gelfand et al. (2009). Our approach combines analytical results and numerical simulations for the description of the dynamics of the PWN-SNR system and the evolution of the particle distribution function within the PWN: the key features of the interaction between a PWN and the surrounding SNR are adequately reproduced during different evolutionary stages. The model is here applied to a few systems of various ages, both young and older ones with the aim of showing how this kind of modelling can help the extraction of information on the system properties from multiwavelength observations: we try to clarify how these properties depend on the different parameters, also discussing potential degeneracies and the observations that would best serve to disentangle them.

Special attention is devoted to the particularly interesting and difficult task of constraining the particle injection in PWNe. The relevance of this issue in terms of PSR physics and shock acceleration physics deserves a somewhat extended discussion, which will be the topic of Section 2.

The rest of the paper is organized as follows. Sections 3–5 describe, respectively, how the nebular evolution is implemented in our model, the scheme we adopt to follow the evolution of the particle distribution function and the emission model. In Section 6 we show the application of our model to several objects, and, finally, in Section 7 we discuss the implications of our findings.

2 PARTICLE INJECTION IN PWNe

40 yr of research on rotation powered PSRs have led to a conundrum: observations of young PWNe have clearly established large particle injection rates into the nebulae, well in excess of the electrodynamic minimum suggested by Goldreich & Julian (1969). The latter corresponds to

1

where is the magnetospheric potential, is the spin-down luminosity, with I being the neutron star's moment of inertia (I₄₅=I/10⁴⁵ g cm²), Ω = 2π/P is the stellar angular velocity and P is the rotation period (P₃₃=P/33 ms). The simplest estimates (Rees & Gunn 1974; Kennel & Coroniti 1984a) show that X-ray emitting PWNe (Kargaltsev & Pavlov 2008) with SYN cooling times of the X-ray emitting particles well less than the nebular ages have electron (and positron) injection rates much larger than shown in equation (1): values reported in the literature for the multiplicity are typically ∼10⁴, based on analysis of the X-ray emission, a result consistent with theoretically derived pair creation rates for the young, high voltage PSRs which have been subjected to such analysis (e.g. Kennel & Coroniti 1984a; Hibschman & Arons 2001; Gaensler et al. 2002). Indeed, the fact that the inferred injection rates are this large is one of the major pillars of support for the theoretical conclusion that PSRs have substantial outflows of e^±, this being the only known means through which such cool objects can have winds denser than the electrodynamic minimum.

Such high density flows with κ≫ 1 support the use of the force-free limit of MHD theory in modelling the torques on rotation-powered neutron stars (Spitkovsky 2006; Kalapotharakos & Contopoulos 2009). As already mentioned, MHD theory has been quite successful in modelling the multidimensional dynamics and appearance of young PWNe (Del Zanna et al. 2004, 2006; Komissarov & Lyubarsky 2004; Volpi et al. 2008). Such models confirm the early inference (Rees & Gunn 1974; Kennel & Coroniti 1984a) that the system behaves as if just upstream of the PSR wind TS, the plasma has low magnetization σ_w≡ (B²/8πγ_wn_±m_±c²) ≪ 1 (n_± is the number density of pairs; the particle density is 2n_±): typically, the average of σ_w in latitude with respect to the rotation axis is ∼0.02 for systems that show a jet–torus morphology (Del Zanna et al. 2004; Komissarov & Lyubarsky 2004). However, this analysis has been done in full only for the Crab nebula.

The nebular dynamics in the MHD model is insensitive to the specific value of the upstream 4-velocity u_w=cβ_wγ_w, as long as γ_w≫ 1. By applying their 1D MHD model to the optical and harder photon emission in the Crab, Kennel & Coroniti (1984b) inferred γ_w≈ 10^6.5, with σ_w≈ 0.005. More modern models (Del Zanna et al. 2006; Volpi et al. 2008) require latitude-averaged σ_w to be somewhat larger (so that magnetic hoop stress can create the jet component of the torus–jet structure), while account of the high-energy SYN emission requires particle spectra at the TS with parameters similar to the 1D Kennel and Coroniti model: again γ_w∼ 10⁶, although this is not explicitly stated since the distribution functions were not tied to the specifics of the MHD flow.

If the apparent low value of σ_w means that the wind just upstream of the TS is really weakly magnetized, σ_w and γ_w are closely tied to the pair multiplicity κ_±=κ/2, since when σ≪ 1 the wind carries the rotational energy lost from the PSR in the kinetic energy of the flow, ⁠. Since ⁠, and L=cΦ², then γ_w=eΦ/2 m_±c²κ_±. In the case of the Crab nebula and PSR, Φ = 4 × 10¹⁶ V, while extant theoretical models of pair creation, from polar caps, slot gaps or outer gaps (Hibschman & Arons 2001; Cheng 2004; Hirotani 2008; Harding, private communication), all yield κ_±∼ 10⁴, for young PSRs with Φ > 10¹⁵ V. Thus, theory also says γ_w∼ 10⁶– if one confines one's analysis to the high-frequency emission from PWNe, theory and observations appear to be in good accord.1

Using the X-ray emitting particles in PWNe, whose SYN cooling times are short, takes advantage of such PWNe being calorimeters for the contemporary particle injection rates, which gives the results some independence from the uncertainties of evolutionary models. However, radio and infrared (IR) emitting particles are much more numerous than the X-ray emitting ones, due to the rapid decline of SYN emissivity with declining particle energy (Kennel & Coroniti 1984b; Gallant et al. 2002). Models of the underlying PSR must be able to account for all the radiating relativistic particles found in the PWNe – the spectral continuity in several systems, demonstrated in the results described below, suggests that the lower energy particles are indeed injected by the PSR, rather than being accelerated out of the non-relativistic material often found embedded within PWNe. Thus, measuring the full relativistic particle content in PWNe is an important experimental input into modelling of PSRs, distinct from the modelling of the energy input into the nebulae. For the latter, the force-free model (Spitkovsky 2006; Kalapotharakos & Contopoulos 2009) provides the essential description, but unfortunately this is independent of the multiplicity, so long as κ≫ 1.

We study several young PWNe (Crab, 3C 58, B1509, Kes 75) for which the data can be used to reasonably constrain the pair injection rate – in principle, one would like to use nebulae with ages known, and reasonably complete (including near-IR and far-IR and millimetre) spectral energy distributions (SEDs) are required. We show that all these systems have spectral continuity from the radio through the IR to the X-ray bands, suggesting that a single source for the radiating particles is present in each system. We use these full SEDs to derive new estimates for the pair creation multiplicity. We discuss the constraints these results set on PSR pair production gap models – polar caps, slot gaps and outer gaps – concluding that no existing model adequately explains particle injection rates.

We also discuss alternate hypotheses that low-energy particles are picked up from thermal gas in the nebulae or are fossils left over from some unnamed acceleration process in the early history of the nebula (Atoyan & Aharonian 1996), or represent the effects of a second acceleration mechanism operating at low energy (Gallant et al. 2002), such as cyclotron acceleration stimulated by a heavy ion component of the wind (Hoshino et al. 1992; Amato & Arons 2006) or, more generally, any component with energy greater than m_±c²γ_w. We use the spectral continuity to argue that such a two-component model is unlikely, implying that some additional piece of physics needs to be added to the shock acceleration model to account for these broken power laws, in addition to the extra physics needed to account for the total injection rate.

3 NEBULAR EVOLUTION

There are three phases of interest in the evolution of a PWN. In this section, we briefly describe them and explain how we model each of them.

We consider objects of a relatively young age (less than 40 000 yr) with the PSR still inside the SNR, and assume spherical symmetry, neglecting the PSR proper motion [see van der Swaluw et al. (2004) for a discussion of this point]. We also assume, in deriving the PWN evolution, that radiation losses of the particles are negligible. Given the weak dependence of the PWN radius on the PSR luminosity, and the fact that radiation losses only affect the particle energy content and not the magnetic energy content, we do expect this to be a good approximation. This allows us to use analytical formulae and substantially reduce the computational requirements.

3.1 SNR evolution

In order to model the evolution of a PWN inside an SNR, we first need to solve for the evolution of the SNR itself. An excellent semi-analytic model for the evolution of an SNR without an embedded PWN was provided by Truelove & McKee (1999). Using pressure balance between the shocked ISM downstream of the outer forward shock and the shocked ejecta beyond the reverse shock, they solve for the evolution of both the outer forward shock and the inner reverse shock from the early ejecta-dominated phase (Hamilton & Sarazin 1984) to the later Sedov phase (Ostriker & McKee 1988). It is shown that the evolution can be cast in a dimensionless form by using the following characteristic variables:

2

3

where ρ₀ is the density of the ISM (assumed to be uniform), E_SN is the kinetic energy of the SN explosion, whose canonical value is 10⁵¹ erg, and M_ej is the mass of the ejecta. The only free parameter in the Truelove & McKee model is the mass distribution in the ejecta. Unfortunately, this is not directly constrained by observations, and choices in the literature are mostly based on theoretical assumptions; for example, Chevalier (1982) suggested a distribution with an inner plateau and an outer steep profile while Truelove & McKee (1999) focused on the case of self-similar ejecta with a density distribution ρ_ej∝r^−αt^α−3 and a velocity profile v_ej∝r. Many recent numerical works on young PWNe have indeed adopted the latter distribution, with great success in reproducing the observations (van der Swaluw et al. 2001; Bucciantini et al. 2003, 2004a; van der Swaluw 2003; Del Zanna et al. 2004, 2006).

Once the evolution of the SNR is known, it can be used to constrain the evolution of the PWN.

3.2 Free-expansion phase

The first phase of the PWN evolution is generally referred to as a free-expansion phase and has the PWN expanding inside the cold ejecta. This phase lasts for a few thousand years, until the PWN reaches the reverse shock. In this phase, the evolution of the PWN is independent of the evolution of the SNR shell, because no contact has yet been established between the two. In the case of self-similar ejecta and constant PSR luminosity L, a solution for the radius of the PWN as a function of time has been known for a long time (Chevalier & Fransson 1992; van der Swaluw et al. 2001):

4

In the more general case in which the PSR luminosity changes in time according to

5

Bucciantini et al. (2004a) have shown that it is possible to find an analytic solution and they provide a series expansion for it, showing that the first few leading terms of the series already provide an excellent approximation. In this work we adopt this approach, which allows us to properly include in our model the spin-down properties of the PSR.

The free-expansion phase lasts as long as the radius of the PWN is smaller than the radius of the reverse shock computed using the Truelove & McKee model. Once the PWN reaches the reverse shock, the reverberation phase begins.

3.3 Reverberation phase

Once the PWN has reached the SNR reverse shock, it is in contact with the shock-heated ejecta. The pressure of the hot ejecta is higher than the pressure of the relativistic material inside the PWN. As a consequence, the PWN expansion halts, and the system starts contracting (van der Swaluw et al. 2001).

An analytic model for this phase is not available; however, the evolution can be treated within the so-called thin-shell approximation (Giuliani 1982): the evolution of the PWN outer radius is described in terms of the evolution of a thin shell of material enclosing a mass M_sw equal to the mass of ejecta that has been swept up by the PWN up to the beginning of the reverberation phase. This shell is bounded on the inner side by a hot relativistic plasma with pressure P_in and on the outer side by the shock-heated ejecta with pressure P_out:

6

The value of the PWN pressure is computed numerically using energy conservation:

7

where t_r is the time at which the reverberation phase begins. It is less clear what P_out should be. A lower limit is given by the pressure that the shock-heated ejecta would have in the absence of a PWN. An upper limit corresponds to the pressure of the Sedov solution. Numerical simulations show that the interaction of the SNR with the PWN leads to additional heating of the ejecta, due to sound waves that are launched inside the SNR during the reverberation phase (see fig. 3 of Bucciantini et al. 2003). On the other hand, it takes longer for the SNR to relax to the Sedov solution, while during the first compression the value of P_out is found to be close to about 50 per cent of the Sedov value (Bucciantini et al. 2003). In our model, we use this fiducial fudge factor.

In 1D simulations (van der Swaluw et al. 2001), the reverberation phase is characterized by a series of compressions and expansions of the nebula, until the system relaxes to the Sedov–Taylor phase, which is finally established once the PWN reaches pressure equilibrium with the SNR. This behaviour, which resembles a damped oscillator, is an artefact of the 1D geometry. In a more realistic multidimensional regime (Blondin et al. 2001), the evolution of the PWN during the reverberation phase is subject to strong Rayleigh–Taylor (RT) instabilities and efficient mixing of the relativistic material with the SNR matter. From a dynamical point of view this mixing acts as a viscous term on the evolution of the nebula, and one might expect, instead of a series of oscillations, an almost complete relaxation to the Sedov–Taylor solution after the first compression.

Of course in the presence of mixing and the related clumpiness, the volume occupied by the relativistic plasma will not be directly related to the radial extent of the nebula. However, equation (7) will still hold if one interprets R(t) as an effective radius related to the total volume of the relativistic plasma, ⁠, rather than as the radial extent of the nebula (which in general might be larger because of clumpiness). In one-zone models, this effective volume is also all that matters for computing particle adiabatic losses.

In our model the reverberation phase ends after the first compression, once the PWN pressure reaches the value of the Sedov solution for the SNR.

3.4 Sedov–Taylor phase

Once the pressure inside the PWN reaches the value proper of the Sedov solution corresponding to the SNR forward shock, the Sedov–Taylor phase begins. This usually happens at an age of ∼10⁴ yr. By this time, because of spin-down, the PSR luminosity can in general be neglected in the evolution of the system. The speed of the forward shock, v_fs, and the post-shock pressure are both given by the Truelove & McKee (1999) model. The pressure inside the SNR is then assumed to be equal to the central pressure of the Sedov solution P_out∼ 0.5 ρ₀v²_fs.

The nebula evolves according to

8

where t_s is the time at which the Sedov–Taylor phase begins and, again, the radius is related to the total volume of the relativistic plasma rather than to the radial extent of the nebula.

One last comment is about the role of RT instabilities. While during the reverberation and Sedov–Taylor phase, reinterpretation of the radius appearing in equations (6)–(8) in terms of effective volume is necessary, during the initial free-expansion phase, RT and mixing are not as important: Bucciantini et al. (2004b) and Jun (1998) have shown that, in general, the radius derived from the 1D solution is at most 10–15 per cent smaller than the true radial extent of the nebula. We might then conclude that, at least for very young objects, the 1D radial model provides a reliable estimate both for the radial extent of the nebula and for the volume occupied by the relativistic fluid.

4 PARTICLE EVOLUTION

Once the evolution of the nebula is known, one can compute the evolution of the particle distribution function.

The energy of a particle is evolved according to

9

where B(t) and U(t) are the magnetic field and the background photon energy density in the nebula, respectively, σ_T is the Thompson cross-section and m is the particle mass. The magnetic field in the nebula is computed assuming that the ratio of magnetic to total energy in the nebula η_M is constant in time (0 < η_M < 1), while the total energy is computed considering the PSR injection and the adiabatic losses of the 1D model as described in Section 3. The background photon energy density includes different contributions: cosmic microwave background (CMB) and starlight, which are constant in time, and SYN emission, which is computed, instead, together with the nebular evolution. We have verified that IC losses are in general negligible with respect to SYN losses, and SYN–IC is important only in young compact objects, where the magnetic field is stronger and the SYN emissivity higher.

Given a particle injected at time t₀ with energy E₀, one can solve equation (9) to derive the energy E(t, t₀, E₀) of the particle at time t. An analytic solution of equation (9) is not, in general, available, neither is its inverse, and the energy evolution has to be computed numerically. Once the evolution of the energy is known, then it is possible to compute the evolution of the particle spectrum according to

10

where is the particle injection rate per unit energy interval.

Given the very short SYN lifetime of the particles at the high-energy end of the spectrum, we use a Lagrangian scheme in energy space. The energy space is originally divided into energy bins whose extrema are given by E_i(t). These evolve according to equation (9). When a bin moves to energies lower than 50 keV, it is removed from the sample. New bins are continuously added at the high-energy end of the distribution function, as the old ones evolve to lower energies. The time-step is then dictated by the requirement of sufficient energy resolution in the high-energy portion of the spectrum.

The total number of particles inside the ith bin, at each time, is then computed by integrating numerically the following equation:

11

where t₀ is the initial time at which injection in the bin begins. This approach guarantees conservation of particle number.

In principle, the functional form of can be arbitrarily chosen.

For the pairs, we assume a broken power law. Given that injection is due to the PSR wind, and that the only energy scale available in the wind is the PSR voltage, ε_v=eΦ, we have decided to use a scale-free model for a broken power law:

12

where γ₁ < 2 < γ₂, ε_c is the peak of the injected energy distribution and ε_m is the low-energy cut-off, usually found to be small compared to ε_c. The scale-free approximation implies that the ratios ν_e=ε_m/ε_v and μ_e=ε_c/ε_v, are constant in time. μ_e and ν_e are in principle free parameters of the model, and their value is obtained from a fit to the data. However, the condition ε_c≪ε_v translates into a requirement that μ_e≪ 1.

Once μ_e and ν_e are fixed, C₀(t) is given by energy conservation. We have

13

with η_e being the fraction of the PSR luminosity injected into pairs, and we have assumed ε_m≪ε_c and ε_c≪ε_v. The particle number flux in the wind is

14

and the average energy/particle in the spectrum of equation (12) is, for ε_m≪ε_c≪ε_v,

15

In expression (15), we have identified the average energy/pair in the spectrum with the upstream flow energy/pair, but we have not specifically used the shock jump conditions. This is because, as we will discuss in Section 7, the distribution of equation (12) implies some additional mechanism for energy redistribution, which may or may not be directly connected to the TS. From equation (15), it is clear how it is possible to estimate the PSR multiplicity through a fit of the model to the PWN emission.

If the wind is characterized by a single value of the Lorentz factor, then ⁠. However, the approximation of a single Lorentz factor might not be correct; in this case the ratio ε_m/ε_c cannot be a priori determined, but is found from requiring a fit to the entire spectrum. In particular, once an upper limit for ε_m is derived from fitting the low-frequency radio emission, this can be used to infer the particle number flux and to estimate the average Lorentz factor.

We consider the possibility that a minority by number of high-energy particles in the equatorial return current (either positrons, ions or electrons, depending on the return current sign) are part of the PSR wind and carry a fraction η_p of the energy flux. Focusing on the case when ions are the implied particles, one might look for signature of their presence. For simplicity, we consider them injected with a monoenergetic spectrum centred at an energy γ_wm_pc² (equivalent to assuming that they are moving at the same Lorentz factor as the pairs injected at ε_c). Ions are affected by different loss mechanisms compared to pairs: while radiation losses are negligible, p–p scattering and diffusion outside the nebula might be important. Following the approach of Amato, Guetta & Blasi (2003), we assume that both diffusion and p–p scattering remove particles from the distribution function. Due to these effects, the number of ions of energy E between time t and t+ d t changes by

16

where E(t) is given by equation (9), and where the characteristic p–p time-scale is

17

with σ₀ = 5 × 10⁻²⁶ cm² and n_t is the target proton density. The typical diffusion time-scale is

18

where r_L is the ions' Larmor radius.

5 EMISSION

There are three channels of non-thermal emission in PWNe. From radio frequencies to MeV photon energies, the dominant emission process is SYN radiation by the accelerated electrons and positrons gyrating in the nebular magnetic field. At higher energies, there are two possible contributions: IC emission by the high-energy pairs and gamma-ray emission from neutral pion decay, following p–p scattering between accelerated and target protons.

The SYN luminosity is computed using the ‘monochromatic’ approximation. The power per unit frequency emitted by a single particle with energy γm_ec² is

19

where σ_T is the Thompson cross-section, B is the nebular magnetic field and

20

is the characteristic emission frequency.

Equation (19) has then to be integrated over the pair distribution function (see the previous section) in order to derive the total nebular SYN luminosity at a given frequency.

IC is computed using the full Klein–Nishina cross-section (Jones 1968; Blumenthal & Gould 1970) and assuming scattering over three main target photon fields: (1) CMB, modelled as a perfect blackbody, with a typical temperature of 2.75 K; (2) SYN emission from the PWN, consistently computed as described above; (3) ‘Galactic background’, which in principle includes both Galactic starlight and any local contribution in the optical–IR band. This third term is modelled as a suppressed blackbody: both the temperature (in the range of a few hundreds/thousands Kelvin) and the suppression factor (of the order of 10⁻¹⁰) are, in principle, free parameters that can be adjusted. The role of this contribution will be discussed in more detail when dealing with specific objects.

Gamma-ray emission from neutral pion decay is computed following Amato et al. (2003). In the scaling regime, a monoenergetic proton distribution with energy E_p leads to the following pion injection spectrum:

21

where

22

and the constant K_π is found through the normalization condition:

23

where N(E_p) is the number of protons with energy E_p.

The number of photons emitted per unit time and energy interval is then

24

The total pion emissivity is then computed integrating over the total proton energy distribution. The density of target protons n_t that enter equation (17) is the main unknown: if the protons are well confined within the PWN, then a reasonable assumption would be to use the swept-up ejecta mass; if they can leak out, then one should consider the total ejecta mass. A realistic model for proton diffusion in PWN-SNR systems is far beyond the scope of this paper: we keep n_t as a free parameter and discuss the implications of its assumed value whenever relevant.

6 RESULTS

In this section we show the results obtained when applying our model to a few objects of different ages, both young and old. For each object, we will not attempt at obtaining the best possible fit of the data, which in our view would not be very significant given the intrinsic limitations of a one-zone model. In fact, our primary interest is to put some constraints on the unknown parameters and to clarify which observations would be more relevant to improve our understanding of these systems. At the same time, we will discuss the limitations of this approach and make clear which part of our results should be considered as really reliable.

6.1 The Crab Nebula

The Crab Nebula is the system for which we have the best constraints, regarding both the PSR injection properties and the luminosity at different wavelengths. Available data extend from low radio frequencies (Baldwin 1971; Baars 1972), through mm/IR (Mezger et al. 1986; Bandiera, Neri & Cesaroni 2002), optical/UV radiation (Hennessy et al. 1992; Veron-Cetty & Woltjer 1993), X-rays (Kuiper et al. 2001) and gamma rays from MeV to TeV energies (Aharonian et al. 2004; Albert et al. 2008; Abdo et al. 2010a). The source distance is estimated to be 2 kpc, and the nebular volume is 30 pc³ (Hester 2008), corresponding to a spherical radius of about 6 ly. This is the only system for which we know the exact age (950 yr). We also have a reliable estimate of the braking index (2.51 implying β = 2.33; Lyne, Pritchard & Graham-Smith 1993), which, when combined with the present luminosity L≃ 5 × 10³⁸ erg s⁻¹, allows us to derive an initial luminosity L(0) = 3.5 × 10³⁹ erg s⁻¹ and a characteristic spin-down time-scale τ≃ 730 yr.

Much more uncertain are the properties of the SNR which should be surrounding the nebula. Several attempts at detecting both the forward and the reverse shock have only provided upper limits. More information can be derived by studying the filamentary network surrounding the PWN and comparing the results with theoretical expectations based on hydrodynamical models of the interaction between the nebula and the surrounding ejecta (Hester et al. 1996; Jun 1998; Bucciantini et al. 2004b). The visible remnant contains at least 1–2 M_⊙ of He-rich line-emitting gas and the PWN expands at a speed of about 1.5–2 × 10³ km s^-1. Traditional models for the evolution of the Crab Nebula (Chevalier & Fransson 1992; van der Swaluw et al. 2001; Bucciantini et al. 2003; Del Zanna et al. 2004) have all assumed a flat density profile of the ejecta. Steeper ejecta (α≥ 1 in Table 1) give a mass in the filamentary network higher than what is observed. Assuming flat ejecta and a canonical SN energy of 10⁵¹ erg, in order to reproduce the proper size of the nebula a mass M_ej = 9.5 M_⊙ is needed, very close to the most recent estimate by MacAlpine & Satterfield (2008). The lack of any clearly detected SNR shell prevents any reliable estimate of the ISM density, and we simply assume a fiducial value of 0.1 cm⁻³.

The PSR and SNR parameters are listed in Table 1, together with the parameters that describe the particle energy distribution at injection. The low-energy spectral index γ₁ is fixed according to the radio data, while the high-energy index γ₂ is chosen in order to minimize the X-ray residuals. We want to stress here that the X-ray emission is concentrated in the central region of the PWN and it is strongly affected by the details of the flow dynamics just downstream of the TS (Del Zanna et al. 2004, 2006). It is unrealistic to expect that a simplified one-zone model can provide an accurate description of the high-energy spectrum. A one-zone model can at most provide an indication of the best power-law index that can fit the data, but more realistic multidimensional models are necessary to address the emission properties in the X-ray band. Two features are interesting to notice: (1) γ₂ is higher than the value of 2.23 typical of relativistic Fermi shock acceleration with isotropic scattering in the fluid frame and (2) there are indications that a single power law cannot reproduce the complete set of X-ray data points.

Several features are present at high energies, but it is not clear if they are intrinsic [a suggestion in this sense has been put forward by Volpi et al. (2008), where they present a multidimensional model of the emission that shows similar features], if they are simply due to calibration issues among different instruments or even for one instrument at different energies, or if they are due to temporal variability of the emission properties (Bucciantini 2008).

The present model can reproduce the SYN part of the spectrum within few per cent accuracy (greater discrepancies being present at high energies as discussed above). At present, our particle distribution function does not include any smooth high-energy cut-off, the reason being the very short lifetime of high-energy particles that prevents efficient numerical integration of their evolution. However, a good match with the data can be obtained by introducing a posteriori an exponential cut-off in the particle energy distribution. This reflects in a cut-off in the emission spectrum at a frequency of around ∼ 10²² Hz. Our model requires that 80 per cent of the PSR spin-down luminosity goes into accelerated pairs. This does not leave much energy to be carried by higher energy particles (i.e. ions; Amato et al. 2003). In principle, the presence of ions might have observable consequences in the TeV band; however, even assuming that the remaining 20 per cent of the spin-down energy is all carried by ions, their contribution to the gamma-ray emission turns out to be negligible. On the other hand, if the return current is made by a high-energy lepton beam (Arons, in preparation), there might be radiative consequences of this latter scenario that have not been investigated, however, again, they are not expected to give any appreciable contribution to gamma rays.

TeV emission is due to IC scattering by the pairs, for which the main target is the nebular SYN emission in this case. We estimate a magnetic field B≃ 200 μG, slightly lower than previous estimates (Aharonian et al. 2004; Hester 2008). Overall the broad-band spectrum is very well reproduced, with the largest discrepancies limited to the EGRET data points around 10 MeV. These points have large uncertainties and are likely affected by calibration issues (Abdo et al. 2009). Indeed, recent Fermi data (Abdo et al. 2010a) are consistent with our model curve (Fig. 1).

Figure 1

The Crab nebula integrated emission spectrum. Data points are from Baldwin (1971), Baars (1972), Mezger et al. (1986), Bandiera et al. (2002), Veron-Cetty & Woltjer (1993), Hennessy et al. (1992), Kuiper et al. (2001), Aharonian et al. (2004), Albert et al. (2008) and Abdo et al. (2010a) (triangles are EGRET points and squares are Fermi points). The solid line is the total luminosity, dashed line is the IC-CMB and dotted line is IC–SYN.

Open in new tabDownload slide

Baldwin (1971) has shown that the radio power-law spectrum of the Crab nebula extends down to the ionospheric cut-off at ∼30 MHz. One can use this piece of information to derive an upper limit on the minimum particle energy at injection. In order for the radio spectrum to extend to those energies as an uninterrupted power law, we need to assume ε_m/ε_c < 7 × 10⁻⁴. From equation (15), we then find for the wind Lorentz factor γ_w < 5 × 10⁴. This has to be compared with the typical Lorentz factor of the particles at ε_c, which is 7.4 × 10⁵, and with the minimum Lorentz factor of injected particles, which is <500. The estimated value of γ_w translates into a lower limit on the pair multiplicity of κ≳ 10⁶.

Our model allows us to compute a posteriori the energy radiation losses. We do this in order to verify to what degree the adiabatic approximation for the evolution of the PWN radius is correct. For the Crab Nebula, we find that about half of the total energy injected into the nebula in the form of pairs has been lost via SYN emission. However, the dependence of the radius on the PSR luminosity L(t) is in general weak [scales as L(t)^1/5 for constant luminosity], so we expect at most modification of the order of 20 per cent in the radial evolution and may be in the adiabatic losses. Such a value is well within the simplifications and the approximations of the model. In order to properly take into account radiation losses, one needs to abandon any analytic solution for the dynamics and solve the coupled system of equations for the dynamics and the emission simultaneously.

6.2 3C 58

3C 58 is another example of a young PWN, which shares many similarities with the Crab nebula. It has a typical non-thermal spectrum extending from radio to X-rays (Salter et al. 1989; Green & Scheuer 1992; Torii et al. 2000; Slane et al. 2004, 2008). It shows clear evidence of energy injection from the central pulsar PSR J0205+64 in the form of a jet–torus structure (Slane, Helfand & Murray 2002). Recent Spitzer measurements of the PWN IR luminosity have showed clear evidence for a possible injection break, as in the Crab (Slane et al. 2008). At a typical distance of 3.2 kpc (Chevalier 2005), its size is about 5 × 9 pc, equivalent to a volume of ∼140 pc³ (corresponding to a spherical radius of ∼3.3 pc). An association with SN 1181 has been proposed (Stephenson & Green 2002). However, recent measurements of the expansion of the filamentary structure surrounding the source (Bietenholz 2006; Rudie & Fesen 2007), interpreted with a simple model of the PWN-SNR system (van der Swaluw 2003; Chevalier 2005), have ruled out this possibility.

The PSR's present spin-down power is L≃ 2.7 × 10³⁷ erg s⁻¹, with a characteristic dipole age t_c = 5390 yr (Murray et al. 2002). Neither the true age nor the PSR braking index is known, and this prevents any reliable determination of the initial spin-down power. However, it is possible to obtain some constraint from the information available on the expansion of the PWN and surrounding SNR. In particular, Bocchino et al. (2001) measured a mass in the filaments of ∼ 0.1 M_⊙, and thermal emission corresponding to a typical velocity of the swept-up shell of ejecta V_SH∼ 340 km s⁻¹. Assuming that the filamentary structure corresponds to the swept-up shell of ejecta (in analogy with the case of the Crab Nebula), we can use our model for the system evolution and obtain an estimate of both the system age and the PSR braking index depending on the density profile of the ejecta.

In Table 2, we present a set of models derived for different values of the unknown braking index (3 and 2.5) and for different density profiles of the ejecta. In all the different models, the swept-up mass and the nebular volume are consistent with extant measures. The shock speed varies slightly, but is in all cases consistent with measurements, within the approximation of a spherical model. We again assume a canonical SN explosion energy of 10⁵¹ erg. As in the case of Crab the lack of any clearly detected SNR shell prevents any reliable estimate of the ISM density, and we simply assumed a fiducial value of 0.1 cm⁻³.

By looking at Table 2, one immediately realizes that the inferred age is not very sensitive to either the braking index or the assumed ejecta structure. The reason for this is that in all cases the age is much less than the characteristic age, and the PSR spin-down power in the past was not much larger than the present value. It is interesting to note that in all models the shock speed tends to be higher for steeper ejecta, favouring the case r⁻². However this would result in a very small mass of the ejecta, inconsistent both with estimates based on the observed filamentary knots (∼8 M_⊙; Rudie & Fesen 2007) and with stellar evolutionary models.

It is our opinion that cases with a density profile which is either flat or r⁻¹ are more likely, even if they are associated with a lower shock speed: this is because steeper ejecta lead to unreasonably small ejecta masses. In all cases, the age estimated for the system rules out any association with SN 1181. However the filamentary structure, which we have interpreted as mostly due to swept-up ejecta, is rather complex and with various velocity components (Rudie & Fesen 2007), so that it cannot be ruled out that part of this filamentary network might be associated with the 1181 event.

The question that naturally arises is whether spectral fitting allows us to discriminate between the models listed in Table 2. Unfortunately, the answer is negative: fits to the observations lead to similar results in terms of the required injection parameters for all of these models. The reason is that all of them correspond to similar underlying ages, similar integrated spin-down power (i.e. total injected energy) and the same size (same adiabatic losses). The parameters reported in Table 1 correspond to the case of ejecta with α = 1 and braking index 3. Spectral fitting is shown in Fig. 2. Very similar values of the parameters are found also for a braking index of 2.5 (β = 2.33).

Figure 2

3C 58. Data points are radio from Salter et al. (1989) and Green & Scheuer (1992), X-ray from Slane et al. (2004) and IR from Slane et al. (2008) (diamonds are IRAS and Spitzer points). The solid line is the total luminosity, dashed line is the IC-CMB and dash-dotted line is IC–SYN.

Open in new tabDownload slide

If we compare the particle spectrum at injection with that of Crab, we find that the spectral index is now flatter at low energies and steeper at high energies. Again, the required pair injection efficiency is very high ∼80 per cent and the ratio μ_e∼ 8.3 × 10⁻⁵ is only a factor of 2 lower than for Crab. The inferred magnetic field is 43 μG, higher than other recent estimates (Slane et al. 2008). The ratio of the magnetic energy to the total energy is ∼0.5: this is quite high and leads to a value for η_e+η_M > 1 which apparently violates the total energetics. The discrepancy is however ∼30 per cent, and we consider it to be within the approximations of our model (particularly critical in this respect are the simplifications related to the assumption of spherical geometry) and the uncertainties of the data (especially the estimated distance, and hence volume). Models with η_e+η_M = 1 always underpredict the radio emission, independently of the value chosen for η_M.

Within the present model, we consider our estimates of η_e and η_M rather solid: attempts at fitting the data assuming a lower magnetization and a higher particle content result in underproduction of the radio emission, while larger magnetization leads to the break frequency moving below 10¹¹ Hz, in contrast to observations. One final remark we want to make about η_M is that the idea of a relatively high magnetization in this PWN is also supported by its substantial elongation, created by hoop stresses of the toroidal magnetic field (Begelman & Li 1992).

The existence of radio measurements down to 100 MHz allows us to derive also in this case a lower limit on the pair multiplicity: we find that γ_w < 3 × 10⁴ to be compared with the typical Lorentz factor of the particles at ε_c, which is 1.2 × 10⁵, and with the minimum Lorentz factor of injected particles, which is <1600. This implies κ > 5 × 10⁵. All the latter values are very close to the corresponding estimates for Crab.

6.3 MSH 15−52

The system formed by the pulsar PSR B1509−58 and the SNR MSH 15−52 is striking for its relatively large size compared with its presumed age. Its distance is estimated to be 5.2 ± 1.4 kpc (Gaensler et al. 1999). Assuming the central value as a fiducial distance (which we will do in the following), the SNR radius turns out to be ∼ 83 × 10¹⁸ cm. PSR B1509−58 is a young PSR. Its current spin-down power is 1.8 × 10³⁷ erg s⁻¹, with a characteristic dipole age of 1550 yr (Kaspi et al. 1994; Livingstone et al. 2006). This is one of the very few PSRs for which the braking index is known: β = 2.087. Unfortunately, our ignorance of the true age of the system still prevents knowledge of the initial spin-down power, L₀: this can only be determined as a function of the assumed age, t. A maximum possible age of the system can be estimated from the braking index and the characteristic dipole age, and the result is t < 1690 yr. For such an age, the size of the SNR implies a high E_SN/M_ej ratio, typical of SN Ib/c (Gaensler et al. 1999; Mazzali, Iwamoto & Nomoto 2000; Iwamoto et al. 2000; Tsvetkov 2002; Mazzali et al. 2003; Tanaka et al. 2009). The size of the PWN is poorly constrained because of the high radio foreground from the SNR. Spectral modelling (Du Plessis et al. 1995; Chevalier 2005; Nakamori et al. 2008), however, suggests that the cooling frequency lies just below the Chandra band (see also the present discussion), which is also in agreement with the size of the nebula at TeV energies being comparable with the Chandra size. If cooling starts being important at frequencies no lower than the X-rays, then one can estimate the size of the PWN from Chandra images: the result is a PWN radius of ∼ 17 × 10¹⁸ cm, for the assumed distance.

Knowing the braking index, and assuming an age (t) and an ejecta density profile (α), we can use our model to determine which values of E_SN and M_ej return the assumed radius for both the PWN and SNR. The only other parameter that enters this calculation is the local density of the ISM. We assume ρ₀ = 0.001 cm^-3: a higher density would hint at even higher values of E_SN/M_ej and in any case does not allow to match both radii; a lower value marginally affects the results because the system turns out to be in the ejecta-dominated phase anyway.

If one accepts the hypothesis that B1509−58 was born in an SN Ib/c explosion, then one can further constrain the model requiring that the mass of the ejecta be in the range of 4–10 M_⊙ and the kinetic energy released be in the range of 5–20 × 10⁵¹ erg. In Table 3, we report the different age ranges. Flat ejecta tend to favour slightly younger systems. In spite of the relatively small differences between the ages reported in Table 3, spectral modelling can prove a useful tool to discriminate between different scenarios for this source. In the case of 3C 58 the inferred ages for all models were much smaller than the characteristic dipole age, so the injection properties of the PSR were almost the same, and it was not possible to discriminate using the SED. On the contrary, for B1509, the inferred ages are close to the characteristic dipole age, and different models correspond to different spin-down histories and different integrated injection energies. In such a situation the SED, and in particular the radio data, can be used to rule out some of the models allowed by the dynamics.

Observations cover the radio GHz band (Gaensler et al. 1999), although the quality of the data in this band is not very high, X-rays from Chandra to INTEGRAL (Gaensler et al. 2002; DeLaney et al. 2006; Forot et al. 2006), Fermi in the MeV–GeV range (Abdo et al. 2010b), Cangaroo-II and HESS in the TeV (Aharonian et al. 2005; Nakamori et al. 2008). Only models with L₀ > 10³⁹ erg s⁻¹ can reproduce the radio, while lower values of L₀ do not satisfy the energetic requirements. Table 3 shows the limits imposed on the model by SED fitting. Flat ejecta correspond to a very narrow range of possible ages and a value of μ_e∼ 1.2 × 10⁻⁶ smaller than what is found for both the Crab Nebula and 3C 58. Vice versa, ejecta with α = 1 correspond to spin-down properties and multiplicity that are closer to what is found in Crab and 3C 58. Compared with these two sources, we find a slightly harder high-energy spectrum at injection γ₂∼ 2.1–2.16. However, this is consistent with the hard spectrum observed in the MeV–GeV range by Fermi (Abdo et al. 2010b). Indeed, the photon index derived from X-ray observations (Gaensler et al. 2002; DeLaney et al. 2006) of the inner region is ∼ 1.4–1.6 and is much smaller than the previously used value of ∼2.7 which seems inconsistent with both observations and energetics (Nakamori et al. 2008). Moreover, the injection break is located between 5 and 50 GHz, at a significantly lower frequency than what was assumed in previous models (Nakamori et al. 2008).

Steeper ejecta give larger values of μ_e, as reported in Table 3. All possible models require a high acceleration efficiency, η_e∼ 0.7, and a strong magnetic field, η_M∼ 0.53, implying B > 20 μG. Again, as in the case of 3C 58, there is an energetic problem, but as in the previous case, this system is strongly magnetized, as also suggested by the presence of a strong X-ray jet [see Del Zanna et al. (2004) for a discussion of the dependence of the jet strength on the magnetization].

The lowest frequency radio point allows us to constrain the value of ν_e, as it was done in Crab and 3C 58. We find that even for α = 1, the average wind Lorentz factor is γ_w < 10⁴, quite small compared with other systems. This translates into a lower limit on multiplicity, which has to be greater than ∼ 2–3 × 10⁵.

Fig. 3 shows the spectrum derived from the model presented in Table 1. From spectral fitting it is clear that the Lorentz factor corresponding to ε_c has to be smaller than 5 × 10⁴ (μ_e has to be smaller than 4 × 10⁻⁶), otherwise it is not possible to reproduce the correct slope and intensity in the radio band (the injection break shifts below 1GHz). As a result η_M must be larger than 0.3, and the magnetic field cannot be smaller than 15 μG. Only by measuring the PWN emission in the range of 10¹⁰–10¹² Hz would it be possible to discriminate among various models and put better constraints on the multiplicity.

Figure 3

B1509−58; MSH 15−52. Data points are radio from Gaensler et al. (1999), X-ray from Gaensler et al. (2002), DeLaney et al. (2006) and Forot et al. (2006), and gamma ray from Nakamori et al. (2008), Aharonian et al. (2005), Kuiper et al. (1999) and Abdo et al. (2010b) (diamonds are EGRET points, squares are Fermi points and triangles are COMPTEL/HESS points). The solid line is the total luminosity, dashed line is the IC-CMB and dash–dotted line is the IC on the enhanced local IR background at 400 K.

Open in new tabDownload slide

As it has already been noted (Aharonian et al. 2005) it is not possible to reproduce the observed gamma-ray emission in the TeV range, either with IC on the CMB or with the average Galactic background. The contribution from IC–SYN is negligible. Various models have been presented to account for such discrepancy.

One suggestion is that of a lower value of the magnetic field in the nebula, which would lead to infer a higher particle content and consequently enhance the gamma over X-ray ratio (Aharonian et al. 2005). However, within this model this does not seem to be a viable solution: in order to increase the particle content up to the value required to fit the gamma rays as IC on the standard Galactic background, we would have to violate the energy fixed by the spin-down history by a factor of ∼5.

Another proposal is that of a possible contribution to the TeV flux of π^o decay, in the presence of relativistic hadrons (Nakamori et al. 2008). It was however immediately realized that this would need so much energy to be put in protons as to require millisecond magnetar conditions at birth for the PSR, which seems unlikely judging from the present PSR-SNR properties. In addition, the swept-up ejecta do not provide a sufficient target density of thermal protons.

The third possibility is that the local photon background, in particular the IR, could be much higher than the average galactic value. Indeed, the SNR itself could be the origin of this excess (Du Plessis et al. 1995; Nakamori et al. 2008). The fit presented in Fig. 3 assumes for the local IR background a blackbody with a temperature of ∼400 K suppressed by a factor of 3 × 10⁻⁷, which corresponds to an energy density about five to seven times higher than the average galactic background. The suggestion by Helfand (2007) that the correlation between the SNR detected by HESS and H ii regions might be due to enhancements in the local photon background is interesting in this regard.

6.4 Kes 75

The SNR Kes 75 (G29.7−0.3) is a shell-type remnant which hosts at its centre a young PWN powered by pulsar PSR J1846−0258. This is a 324-ms PSR, with a characteristic dipole age of 723 yr, a spin-down luminosity of 8.3 × 10³⁶ erg s⁻¹ and a measured value of the braking index β = 2.65, which provides an upper limit for the age of the system of ∼880 yr (Livingstone et al. 2006).

The estimated distance of the system suffers from large uncertainties. The first estimates pointed to a distance of 19–21 kpc (Becker & Helfand 1984). However, it was immediately realized that such a large value implied very peculiar conditions for the SN explosion (Helfand, Collins & Gotthelf 2003; Morton et al. 2007) and the ISM, together with an extremely high SYN efficiency of the PWN. Recent studies have revised the distance to a much smaller value, in the range of 5–7 kpc (Eisenhauer et al. 2005; Leahy & Tian 2008). We will assume a distance of 6 kpc in all our models.

Assuming a distance of 6 kpc, the radii of the PWN and of the SNR are, respectively, R = 1.8 and R_snr = 9.24 ly. The velocity of the forward SN shock is v_fs = 3700 km s⁻¹ (Helfand et al. 2006) as estimated via the Si X-ray line. This places an upper limit on the age of the system t < R_snr/v_fs∼ 800 yr. Various attempts have been made to measure the ejecta mass and the SN energy; however, the bulk of the SNR emission comes mainly from two regions of the SNR shell, and this tends to bias all estimates. In general, very high masses are inferred for both the swept-up ISM and the SN ejecta (Becker & Helfand 1984; Morton et al. 2007; Leahy & Tian 2008), hard to reconcile with theoretical expectations for an SN Ib/c and standard evolutionary models. It must be stressed, however, that mass estimates are very sensitive to the assumed volume of the emission region, which is not well known.

Proceeding as for previous cases, we will first attempt to constrain the possible SNR-PWN parameters through models for the dynamical evolution of the system. Given an assumed age t and an ejecta profile α, there is only one set of values of E_SN, M_ej and ρ₀ that gives the correct size of both the PWN and the SNR and the correct forward shock speed. Further constraints on the parameter space come from the properties of SN Ib/c which limit M_ej in the range 5 M_⊙≲M_ej≲ 16 M_⊙. Cases with α≥ 1 are acceptable only for a system older than 700 yr and are in general unlikely because they require a very low SN energy, E_SN≲ 0.3 × 10⁵¹ erg. Cases with flat ejecta are admissible for a larger range in age from 450 to 650 yr. The typical SN energy does not seem to depend on the age and is E_SN∼ 2 × 10⁵¹ erg; younger systems require lower values of M_ej, and larger values of ρ₀, as shown in Table 4. Cases with marginally steep ejecta, α = 0.5, suggest a larger age, a less energetic explosion, E_SN∼ 0.6 × 10⁵¹ erg, and lower ISM density.

The question arises again whether it is possible to discriminate among the different parameter sets for the dynamics based on the SED. Unfortunately, this is not the case: as it was already realized earlier (Chevalier 2004), Kes 75 is a particle-dominated system by large (more than 90 per cent of the PSR spin-down power seems to be converted into accelerated particles; see below), so the differences in particle content due to the different energetics of the various models can be easily compensated by small changes in the magnetization to give the same SYN emission.

However, fitting the model to the multiband emission spectrum is non-trivial. The radio emission spectrum implies a low-energy particle spectral index at injection γ₁ = 1.7 (Becker & Helfand 1984; Salter et al. 1989; Bock & Gaensler 2005). The average photon index in the Chandra band is 1.9 while deep X-ray images have shown that the photon index in the vicinity of the PSR is ∼1.8 (Blanton & Helfand 1996; Collins, Gotthelf & Helfand 2002; Helfand et al. 2003; Ng et al. 2008).

Evidence for a possible spectral break below 100 GHz was presented by Bock & Gaensler (2005), but these data seem to be inconsistent with previous measurements (Salter et al. 1989) and a single power law cannot be ruled out.

Indeed in our model an injection break below 100 GHz would require a very hard high-energy injection, inconsistent with spectral information at X-ray frequencies. IR data do not provide good constraints (Morton et al. 2007). INTEGRAL data above 15 keV show a particularly hard spectrum, not fully consistent with an expected cooling break in the Chandra band (Terrier et al. 2004), but the PSR might contribute to the flux above 40 keV and be responsible for the excess emission at high energies.

If the X-ray emission in the Chandra band corresponds to freshly injected high-energy particles, one must infer γ₂ > 2.6. However, we found that it is not possible to fit the overall SED when adopting such a soft high-energy injection. The reason is simple: given that the average spectral index is close to the one measured at the injection, the SYN break frequency must be above the Chandra band, but this implies a magnetic field B≲ 10 μG. Even assuming η_e = 1, with this small field all the models under-produce the radio emission by at least a factor of 5. We find that it is possible to fit the SYN spectrum, including the average slope in the X-rays, only if 2.2 ≤γ₂≤ 2.4.

As to the other parameters, we note the following: the value of η_M depends on age, with younger systems (t∼ 450 yr) requiring a field B∼ 30 μG and older ones (t∼ 650 yr) requiring B∼ 20 μG, but it does not depend on γ₂; the value of μ_e depends on γ₂, ranging from 2.5 × 10⁻⁴ for γ₂ = 2.4 to 2.5 × 10⁻⁵ for γ₂ = 2.2, but it is independent of age. All the models require a very high injection efficiency η_e∼ 1, to reproduce the INTEGRAL points. In Fig. 4, we show the spectrum derived according to the model in Table 1.

Figure 4

Kes 75; SNR G29.7-0.3. Data points are radio and IR from Becker & Helfand (1984), Salter et al. (1989) and Bock & Gaensler (2005), X-ray from Collins et al. (2002) and Terrier et al. (2004), and gamma ray from Terrier et al. (2008) and Djannati-Ataĭ et al. (2008). The solid line is the total luminosity, dashed line is the IC-CMB, dash–dotted line is the IC on the enhanced local IR background at 1000 K and dotted line is the IC–SYN.

Open in new tabDownload slide

The Lorentz factor corresponding to ε_c ranges between 2 × 10⁵ and 1.5 × 10⁶, very close to the values derived in Crab. The one corresponding to the model of Fig. 4 is 7 × 10⁵. The minimum Lorentz factor is lower than ∼ 5 × 10⁴, and the average Lorentz factor is ∼ 7 × 10⁴, with an associated multiplicity that has to be greater than 10⁵.

Kes 75 has been detected in gamma rays by HESS (Djannati-Ataĭ et al. 2008; Terrier et al. 2008). The particles responsible for TeV gamma-ray emission, in a purely leptonic model, are the ones emitting SYN radiation in the Chandra band. However the photon index in the TeV range is found to be ∼2.3, far steeper than the average X-ray photon index of 1.9. IC scattering on local SYN radiation gives a negligible contribution of TeV gamma rays. Despite this being a young system, PSR J1846−0258 has a low spin-down luminosity, which results in a particle content about 2 orders of magnitude lower than in Crab and a correspondingly lower energy density of the SYN photon field. IC scattering on CMB photons is about 50 times stronger. However, even this contribution is about a factor of 10 below the observed flux. In addition, scattering on the CMB is expected to occur in the Thompson regime and cannot account for the steepening of the spectral index in gamma rays. A suggestion that has been made is that this steeper spectrum hints at a much warmer radiation background, for which IC should take place in the Klein–Nishina regime. This requirement places the average blackbody temperature of the seed background photons in the range of 1000–2000 K.

With a blackbody spectrum at 1000 K, suppressed by a factor of 5 × 10⁻⁹, it is possible to reproduce the correct gamma-ray luminosity. This contribution corresponds to a local enhancement of the IR background which is only a factor of a few above the Galactic average. Interestingly, as in the case of B1509, Kes 75 is also surrounded by a bright SNR shell, which might be responsible for the enhancement.

The question arises if the TeV gamma rays can be explained by a hadronic component, and in particular by π^o decay. Given the very low magnetization, inferred from the SYN spectrum, it can be shown that only weak constraints can be put on η_p; in particular, one can obtain a good fit of the radio and Chandra X-ray data, using η_e∼ 0.5 (which however underpredicts INTEGRAL data) and a slightly higher nebular field, without violating the energetics. However even by assuming that half of the spin-down energy goes into protons, in order to have a significant contribution to the emission in the 10-TeV range, at least 10 M_⊙ of target thermal protons are needed. This is far in excess of the swept-up ejecta mass of ∼ 0.1 M_⊙, even if it is consistent with the total ejecta mass (in this case, one needs to assume that protons escape from the PWN and interact with the SNR). The main problem in this case is to properly reproduce the observed TeV spectrum. If we assume that its origin is purely hadronic, we predict a peak at 10²⁸ Hz instead of 10²⁶ Hz as is observed. This shift is due to our choice of injection energy for the protons, which is tied to γ_w. Fermi observations should be able to distinguish between a leptonic and a hadronic model, by constraining the emission below 1 TeV. To properly investigate whether π^o decay is a viable possibility of explaining the gamma-ray data, one would need a model for the diffusion of protons outside the nebula, which at the moment we do not have. The simplest possible estimate of the diffusion time of protons outside the nebula, and in the ejecta, gives values that are about an order of magnitude larger than the age of the nebula.

6.5 Old objects

In this section, we discuss two relatively old objects. The code has been developed for the investigation of systems also beyond the free-expansion phase. However the late dynamics, especially if the PSR kick velocity is important, can be quite complex. In addition, old systems are usually very poorly constrained from the observational point of view. The age is generally known only as an order of magnitude estimate, and quite often either the central PSR is not observed (as in the case of G327.1 and IC 443) or the SNR is not observed.

The following discussion highlights the difficulties that one has to face when trying to model sources for which both the quantity and the quality of the data are really poor.

6.5.1 W44

The SNR W44 is known to contain an old PWN, associated with pulsar PSR B1843+01. This is a 267-ms PSR, with a characteristic dipole age of 20 380 yr and a spin-down luminosity of 4.3 × 10³⁵ erg s⁻¹ (Wolszczan Cordes & Dewey 1991). The value of the braking index is not known, while the distance of the PSR is estimated to be 3.1 kpc (Wolszczan Cordes & Dewey 1991).

The SNR W44 (3C 392) has an elongated shape with axes of 25 × 11 arcmin². Its distance is estimated from H i absorption to be 2.6 kpc, corresponding to a typical SNR radius of 11–13 pc (Cox et al. 1999). The small discrepancy between the estimated distance to the PSR and to the SNR is within the uncertainties, and the association is considered a secure one. The nebula is thought to be a transition object from the spherical shape of young systems to later bow-shock-like morphologies (van der Swaluw et al. 2004). The central pressure is estimated from X-ray observations to be ∼ 1.4 × 10⁹ dyn cm⁻² (Cox et al. 1999). Modelling the SNR (Cox et al. 1999; Shelton et al. 1999) leads to the following estimates for the relevant dynamical parameters: E_SN = 10⁵¹ erg, M_ej = 5 M_⊙, ρ₀∼ 6 cm⁻³; for the expansion velocity of the SNR, H i emission gives v_fs∼ 150 km s⁻¹ (Koo & Heiles 1995). However, all existing models assume the spin-down time of the PSR as the age of the nebula; in addition, the quoted value of v_fs is inferred from an H i emitting ring structure which does not trace exactly the SNR, and estimates can be as high as 330 km s⁻¹ (Koo & Heiles 1995).

The PWN is observed in radio (Frail et al. 1996; Giacani et al. 1997) with a typical luminosity of 200 mJy. It has a distinctive cometary shape with the PSR located at the tip of a protruding finger of emission. It is assumed that the PSR has been displaced by its proper motion with respect to the core of the radio emission. X-rays have been detected with both Chandra and XMM–Newton from the head of the cometary nebula, and the flux is measured to be ∼ 2.7 × 10⁻¹³ erg cm^-2 s⁻¹ in the 2–10 keV band (Petre, Kuntz & Shelton 2002; Harrus et al. 2006). The radio emission is about twice as extended as the X-rays, which has been interpreted as due to effective cooling in the nebula, for which, however, there is no indication from the spatial behaviour of the spectral index: this is found to be ∼ 2.2 ± 0.3 with no appreciable variations with the distance from the PSR (Harrus et al. 2006). No IR or gamma-ray emission is detected.

The PWN is quite weak compared to the PSR spin-down energy and also appears to be quite small in size. At the distance of the SNR, a typical value for the volume is ∼0.5 pc³, corresponding to a radius of ∼1.5 ly. It is interesting to note that spectral models of the PWN have often assumed a much smaller age of <5000 yr (Petre et al. 2002), which clearly contradicts SNR models. However such models are usually developed assuming that the spectrum is similar to that in the Crab Nebula, which we think is unjustified.

The major problem in modelling this system, as was already recognized by Petre et al. (2002), is the relatively low radio luminosity compared to the integrated PSR spin-down power. This suggests a much younger system: in particular, we find that a braking index close to 2 and an age of about 10 000–15 000 yr are required to reproduce the observed radio emission. The average magnetic field must be well below equipartition (η_M has to be <0.002, corresponding to a nebular field below 10 μG) in order to avoid overproducing the radio emission.

This small value gains some support from the lack of spectral steepening with distance, which suggests that cooling is not important. Indeed, the difference in size between the radio and X-ray nebula could be explained by the presence of a much stronger magnetic field in the head of the nebula than in the body, resulting in the suppression of X-ray emission in the latter. In Table 5, we list the values of the parameters used to produce the curve in Fig. 5. The value we find for the central pressure is ∼ 1.6 × 10⁹, in agreement with X-ray observations, but the expansion speed is found to be higher, v_fs∼ 300 km s⁻¹, as a consequence of the younger assumed age.

Figure 5

W44. Data points are radio from Frail et al. (1996) and Giacani et al. (1997), X-ray from Petre et al. (2002) and Harrus et al. (2006), upper limits on gamma ray from Aharonian et al. (2002) and Abdo et al. (2009). The solid line is the total luminosity, dashed line is the IC-CMB, dash–dotted line is the IC on the local IR background (same as for Crab nebula) and dotted line is IC–SYN.

Open in new tabDownload slide

It is evident that the quantity and quality of existing data (only four measurements are found in the literature) do not allow us to constrain the PWN model as much as it is possible for younger systems. However, by repeating the same analysis we did in the previous section, we find that the parameter μ_e in the model is 8 × 10⁻³, which corresponds to a Lorentz factor at ε_c equal to 2 × 10⁵, while the average wind Lorentz factor is found to be γ_w≲ 10⁴, corresponding to a multiplicity κ≳ 10⁵.

One thing to note is that the spectrum barely extends to X-ray frequencies. As pointed out above, this is due to the low value of the magnetization we adopt. If the magnetization in the head is higher, as could be expected due to compression in the bow shock, this might enhance the X-ray emission.

Another striking feature is the bump at around 10¹⁵ Hz: this is a typical feature of post-reverberation systems. During the reverberation phase, there is a typical energy at which SYN losses are balanced by adiabatic compression gains. Particles tend to accumulate at this energy.

6.5.2 K2/3 Kookaburra

The ‘Kookaburra’ is a complex of compact and extended radio/X-ray and gamma-ray sources. Diffuse X-ray emission and point-like sources have been detected by ASCA (Roberts et al. 1999; Roberts, Romani & Johnston 2001), XMM–Newton and Chandra (Ng Roberts & Romani 2005). We are here interested in the north-east wing of the complex, where at radio and X-ray frequencies a nebula is found, hosting a young and energetic pulsar, PSR J1420−6048. The PSR location is also coincident with TeV gamma-ray emission detected by HESS (Aharonian et al. 2006).

PSR J1420−6048 is a 68.2-ms PSR, with a characteristic dipole age of 13 050 yr and a spin-down luminosity of 10³⁷ erg s⁻¹ (D'Amico et al. 2001). The braking index is not known. The distance of the PSR is estimated to be 5.5 ± 0.8 kpc (Aharonian et al. 2006).

Identifying the PWN associated with PSR J1420−6048 is, however, rather problematic. The wing-like structure hosting the PSR is usually referred to as K2. In coincidence with the PSR an enhancement of radio emission, usually referred to as K3, is also observed. K2 has an extent of ∼ 15 × 10 arcmin², a total flux at 20 cm of ∼1 Jy and a spectral index of 0.2 ± 0.2, while K3 has an extent of about 3 arcmin, a total excess flux at 20 cm of ∼20 mJy and a spectral index of 0.4 ± 0.5 (Roberts et al. 1999). No IR detection has been reported.

X-rays have been detected by ASCA in the 2–10 keV band, with an extent of about 7 arcmin, an integrated flux of 4.8 × 10⁻¹² erg cm^-2 s⁻¹ and a spectral index of 1.4 (Roberts et al. 2001). A much more compact (0.5 arcmin) but extended source has been detected by Chandra around the PSR, with a total flux, extrapolated to the 2–10 keV band, of 1.3 ± 0.14 × 10⁻¹² erg cm^-2 s⁻¹ and with a spectral index of 2.3 ± 0.9 (Ng Roberts & Romani 2005).

Gamma rays were first detected by EGRET (Thompson et al. 1996) and more recently by HESS (Aharonian et al. 2006). The total HESS luminosity in the 0.4–20 TeV band is 5.1 × 10³⁴ erg s⁻¹, for the assumed distance, and the spectral index is found to be ∼2.2. The 2σ angular extension of the gamma-ray emission (∼7 arcmin) implies a nebular size of ∼11 pc.

One immediately realizes that these data are somewhat confusing. We think that the standard interpretation, according to which the small K3 region is the PWN while the K2 region is the SNR (Roberts et al. 2001), seems unlikely. Indeed, within the framework of our model, this assumption leads to unphysical values for the SN energy and ejecta mass and to overpredict the radio emission.

The larger extent of the ASCA source, compared to the Chandra one, goes against the expectation that the nebular size should shrink at higher frequencies. The spectral index of the TeV emission (interpreted as IC-CMB) is consistent with the steeper Chandra spectrum, but not with the flat ASCA one.

Even the radio data are not conclusive: if one assumes that the PWN contributes also to the K2 emission, rather than just to K3, the radio emission changes by about 2 orders of magnitude. Indeed, as for W44, we will present here just a simple model fit to the data.

The major difficulty in accounting for what is observed relates to the low efficiency in both X-rays and radio, together with the ratio L_γ/L_X > 1, which is suggestive of a large nebular content of low-energy particles. A decent fit to the data can only be achieved if one assumes a relatively young system, with a steep high-energy spectrum (similar to 3C 58) and a weak magnetization (comparable with that in Kes 75). Fig. 6 shows the result of the model corresponding to the parameters in Table 5. The nebular radius is found to be ∼5 pc in agreement with the HESS size. The magnetic field is inferred to be ∼ 5 μG. To get the steeper Chandra spectrum, one needs to assume a cooling break around 10¹⁶ Hz. This is also necessary in order to reproduce the correct TeV spectrum. As already suggested by Aharonian et al. (2006), in the TeV range the IC-CMB dominates over the IC on the standard galactic background. Indeed, our model underpredicts the TeV emission by a factor of 4–5. The local photon background must be higher if the observed TeV flux has to be explained as IC emission: the TeV flux shown in Fig. 6 was obtained by assuming an IR background in the form of a blackbody at 200 K, with a photon density about a factor of 4–5 higher than the standard galactic background.

Figure 6

K3/2 Kookaburra. Data points are radio from Roberts et al. (1999), X-ray from Roberts et al. (2001), and gamma ray from Ng Roberts & Romani (2005) and Aharonian et al. (2006). The solid line is the total luminosity, dashed line is the IC-CMB, dash–dotted line is the IC on the local IR background (assumed as a blackbody at 200 K) and dotted line is IC–SYN.

Open in new tabDownload slide

The parameter μ_e in the model is found to be 5 × 10⁻⁵, which corresponds to a Lorentz factor at ε_c equal to 2 × 10⁵, while the average wind Lorentz factor is γ_w≲ 10⁴, corresponding to a multiplicity κ≳ 10⁵.

7 DISCUSSION AND IMPLICATIONS

Let us briefly review here our findings and what general conclusions can be drawn from the results of our attempt at modelling different objects, both young and old.

7.1 Summary of the results

First, it is remarkable that our one-zone model, despite its simplifications, has worked surprisingly well for the observed SEDs in all the cases we considered. This is already an interesting result, since the assumption of efficient mixing that is at its base is likely a poor approximation at high energies (in the X-ray band), where the flow pattern and the magnetic field structure in the inner nebula are expected to play a dominant role.

The agreement between the model fit and the data is extremely good in the case of Crab, where the unknown parameters are reduced to a minimum (essentially the SNR properties). Our model seems to work rather well for young systems in general, but the quality of the data and the limited spectral coverage lead to larger uncertainties. In the case of strongly magnetized sources, such as 3C 58 and the nebula associated with PSR B1509, we find that in order to reproduce the data our model requires an energy input in the PWN which exceeds the PSR release by 20–30 per cent, judging from the current spin-down power. This conclusion does not seem to be affected by uncertainties in the distance to the sources: moving them closer to (further from) us would decrease (increase) the estimated brightness, but at the same time it would make them smaller and younger (larger and older), leaving the discrepancy between the radiated and accumulated energy almost unaltered. We think it likely that this discrepancy is due to the simplifying assumptions of the one-zone model, which become progressively more important at large magnetizations. In particular, as also shown by recent multidimensional simulations (Volpi et al. 2008), the radio and X-ray emitting particles might sample different magnetic fields. At the same time, it is worth keeping in mind that the estimate of the PSR spin-down power is based on the assumed canonical moment of inertia for neutron stars: I = 10⁴⁵ g cm². A 20 per cent discrepancy in two systems is compatible with present uncertainties.

Our model seems to constrain with good accuracy the value of ε_c. It is not possible to vary the value of ε_c by more than a factor of a few and still reproduce the radio and X-ray data, without violating the overall energy by more than an order of magnitude. At the same time, our model allows us to derive, for each object, an upper limit on the wind Lorentz factor and a lower limit on the pair multiplicity. The limits we put on γ_w and κ are strictly valid only if the wind has a unique Lorentz factor and the energy scales ε_m and ε_c respect the assumed scaling with time, proportional to the PSR voltage. It is not easy to predict how a possible latitude dependence of the wind Lorentz factor, such as that included in axisymmetric models of PWNe (e.g. Del Zanna et al. 2004), or a different dependence on time of the energy scales would affect our conclusions.

7.2 Implications for the wind composition

A controversial issue in the study of PWNe, and more broadly PSR winds, is the possible presence of ions or high-energy particles. Given that a direct signature of their presence in the wind is not available, all efforts have focused on the search for possible observables in the PWNe that could be indicative of their presence. Let us briefly recall why ions where invoked in the first place, what is the present status of those arguments and what the present work tells us about their contribution to the PWN emission.

The idea that a particle population with a larger Larmor radius (ions or leptons in the return current) could be carrying most of the wind energy was made attractive by the following three main reasons: these particles could help solve the problem of particle acceleration at the PSR wind TS; they could explain the wisps variability in the Crab nebula; finally, in the case of ions, they could solve the discrepancy between the predictions of 1D MHD models and the gamma-ray flux observed from Crab.

As far as particle acceleration is concerned, at a relativistic pair shock, efficient diffusive shock acceleration is shown to be possible only if the shock magnetization is extremely low (Spitkovsky 2008). The 2D MHD models of PWNe (Del Zanna et al. 2006; Camus et al. 2009) use a simple model of field reversal across the equator where a current sheet occurs, such that the appropriate conditions for acceleration would be realized only in a small sector of the PSR wind TS, with no more than a few per cent of the total wind energy flux flowing through it. On the other hand, complete dissipation of the ‘striped’ component of a PSR's wind (Coroniti 1990) would create a weakly magnetized current sheet occupying the full rotational latitude sector |λ| < i, i = angle between the rotation axis and the direction of the PSR's magnetic moment. Then the fraction of the energy flow carried in the current sheet would be sin i, which can be quite large. The truth would be between these extremes. One, however, needs to remember that the presence of a jet puts a lower limit on the average magnetization, and full dissipation of a large striped wind zone might fail to produce the observed structures.

On the other hand, the presence of an energetically significant ion component in the PSR wind is found to lead to efficient acceleration of the pairs (Hoshino et al. 1992; Amato & Arons 2006). This effect is the consequence of the large Larmor radius of the ions, whose gyration introduces long wavelength turbulence through cyclotron instability. The same result could come from the presence of high-energy leptons, accelerated as a consequence of runaway dynamics in the equatorial current layer (Arons, in preparation).

As to the wisps' variability, Spitkovsky & Arons (2004) showed that this could be explained as the result of compression waves associated with the high-energy current carriers' (ions, in their case) gyration in the shock region; in order to reproduce the observed brightness contrast, one would again need the current carriers (the ‘beam’) to be energetically significant: γ_beamm_beam/γ_wm_± > 2 κ_±. An alternate explanation of the wisp variability, which does without the kinetic effects of the high-energy particles' beam, has been recently provided by Camus et al. (2009), who have shown that it can also be recovered within the framework of pure MHD, due to global instabilities of the TS.

One final consideration concerns the possible detection of signatures of high-energy protons from gamma rays. Modelling of the TeV emission from the Crab Nebula, based on the 1D Kennel & Coroniti flow (Atoyan & Aharonian 1996), was shown to underpredict the observed flux by a factor of ∼5, leading them to suggest a possible contribution from the decay of neutral pions produced in nuclear collisions of relativistic protons. More recently, the presence of relativistic protons as the source of high-energy gamma-ray emission was also suggested in the Vela PWN (Horns et al. 2006). In the case of Vela, however, a later determination of the density of the thermal protons (that would serve as a target for nuclear collisions) resulted in too low a value and led to substantial questioning of the initial claim (LaMassa, Slane & de Jager 2008). At the same time, Volpi et al. (2008) showed that within a multidimensional model for the flow structure in Crab, the estimated IC-TeV emission is easily overproduced, for the same parameters that allow one to better fit the lower energy SYN emission. This suggests that the results are indeed strongly model-dependent and that the discrepancies are likely due to problems with the originally adopted MHD model.

This is also the conclusion we reach in this work. We can reproduce the TeV emission from the Crab and find results in agreement with Fermi observations without any need for a proton contribution. In the case of both B1509 and 3C 58, our pair+magnetic field energetic is already exceeding the estimated PSR energy input by 20 per cent, and fitting the gamma rays through π⁰ decay requires that protons alone carry far more energy than the nebula is currently estimated to store. Finally, in the case of Kes 75, in order to ascribe the excess TeV emission to protons, one needs to assume that protons can effectively diffuse out of the nebula and experience the whole ejecta mass as a target.

We also find that, as one can easily see from Table 1, given the already high value of the quantity η_e+η_M, there is little energy (∼20 per cent at most) left to be carried by any other higher energy component. We can thus conclude that modelling of the SED and of the PWN–SNR evolution suggests that, if ions are present, they contribute a negligible fraction of the energy in the nebula. There is still the unsolved problem of particle acceleration, which neither this work nor the MHD simulations can address. A definitive answer would probably require a merging of models for shock acceleration and MHD models for the global evolution of the PWN.

7.3 Interpretation of spectral breaks

Our study has reinforced a long-standing puzzle: the electron (and positron) accelerator at work in PWNe knows how to create a spectrum convex in energy space, well represented by a broken power law. As was advanced by Kennel & Coroniti (1984b), the spectral steepening observed between optical/soft X-ray energies and the harder (ɛ > 10 keV) spectrum is well understood as the effect of SYN energy losses, with the pre-cooled spectrum being a power law from Table 1. The steepening between mid-IR and optical, when observed, requires energy space structure in the accelerator (an ‘intrinsic break’). The high-energy spectrum can be attributed to diffusive acceleration at the TS, at least qualitatively – Fermi acceleration in relativistically very low σ shocks in the test particle limit yields a spectrum N (E) ∝E^−2.2 (Keshet & Waxman 2005). The radio spectrum requires a much harder distribution of the lower energy particles, ⁠.

The standard model (Kennel & Coroniti 1984a,b), developed to account for the radiation from Crab at near-IR and shorter wavelengths, assigns the conversion of the PSR wind energy in the particle spectra that emit the observed SYN radiation to diffusive Fermi acceleration at the wind TS. Diffusive shock acceleration, with isotropic particle scattering in the upstream and downstream flow frames, always shows particle spectra with a Maxwellian at low energy, plus a power-law supra-thermal tail at high energy, with the temperature of the Maxwellian set by the shock jump conditions (except when the acceleration of the tail is very efficient). These properties are well exhibited by particle-in-cell (PIC) simulations of relativistic shocks in unmagnetizede^± plasma or, if magnetized, in upstream quasi-parallel flow geometry (Spitkovsky 2008; Spitkovsky & Sironi 2009). Our models identify the transition between the soft, high-energy spectrum and the hard, low-energy spectrum as being at the energy ε_c, typically 10⁵ - 10⁶m_±c² (see Table 1). In the Kennel & Coroniti (1984b) model, this energy was identified as the lower cut-off of the shock-accelerated power law and as the ‘temperature’≈γ_wm_±c² of the flow downstream of the TS, thus giving rise to the belief that PSR winds have an upstream flow Lorentz factor γ_w∼ 10⁶. The question of how a transverse relativistic shock with significant magnetization (σ_w > 10⁻³) can efficiently accelerate a non-thermal particle spectrum remains unanswered – existing PIC simulations suggest that the required large amplitude turbulence does not form.

The Kennel and Coroniti model and its descendants deliberately neglected the radio emission from the Crab Nebula (and, by extension, other PWNe). The large (by number) population of particles with E < ε_c suggests γ_w to be much smaller, if the radio and mid-IR to far-IR emitting particles come from the PSR. The fact that the PSRs are the most likely source of the low-energy particles in each of the nebulae gains support from the existence of the radio ‘wisp’ features in the Crab (Bietenholz et al. 2004) closely associated with the similar time variable structures seen in optical and X-ray imaging, as well as similar structure seen in 3C 58 (Bietenholz 2006) – these structures are clearly coincident with the TS. Then, in the regions of the flow populated by the large particle flux feeding the low-energy population, γ_w∼ 10⁴ or smaller.

The low-energy particle spectrum is definitely not a single temperature relativistic Maxwellian, with T∼ε_c. Since the TS is not spherical, becoming more oblique in higher rotational latitude regions (Del Zanna et al. 2004; Komissarov & Lyubarsky 2004), the post-shock temperature declines with increasing latitude, suggesting the possibility that the low-energy particles enter in the polar regions of the outflow. However, replicating the observed spectrum as the envelope of a sequence of Maxwellians requires a factor of ∼100 variation in temperature, which requires most of the mass flux to be nearest to the rotation poles. There is no sign of pole to equator asymmetry in the radio emission near the PSRs, except for the radio wisp features in the Crab, which are components of the immediate post-shock flow, as judged from the optical and harder photon emission. In addition, PSR pair production models show no sign of such strong polar enhancements of the mass flux. However we cannot rule out that mixing and effects related to integration along the line of sight might prevent the detection of the implied variations, even if these will have to be highly efficient.

In any case, let us assume in the following that also the low-energy particle spectrum results from some acceleration process, rather than from the convolution of different thermal distributions, and let us speculate on the nature of such a process.

Turbulence and associated Fermi II acceleration in and around the TS is one appealing possibility. The ‘wisp’ motion, observed from the radio through X-rays [see Hester (2008) for a comprehensive review], has recently been interpreted as the result of the strongly variable TS structure found in high-resolution MHD simulations of the Crab nebula (Camus et al. 2009). The shock instability implied is a TS variant of the standing accretion shock instability, with outer scale variable velocity δv∼v∼ 0.6c and length-scale δr∼r∼ 0.5–1 light-year. The magnetized motions observed in the simulations (which do well in replicating the time variable spatial structure observed in the nebular ‘wisps’) can act as an accelerator through the Fermi II process (Kardashev 1962; Stawarz & Petrosian 2008), creating electron and positron spectra N(E) =N₀(E/E₀)^−s.

In quasi-linear models of Fermi II acceleration in isotropic small amplitude Alfvén turbulence (e.g. Stawarz & Petrosian 2008), ⁠, where ε=T_accel/T_escape. These models include an analogue of scattering from large amplitude magnetic inhomogeneities (Fermi's original model), in the case where the wave energy spectrum scales in proportion to k⁻², with k being the wavenumber. This case is germane to acceleration due to pair interaction with large-scale moving magnetic fluctuations: the acceleration time-scale for interaction with large-scale moving ‘eddies’ is

25

where λ₀ is the outer scale of the turbulence, comparable to the shock radius in our case (Camus et al. 2009), v_eddy the turbulent velocity and B and δB the magnetic field and its fluctuations, respectively. The high downstream flow velocity suggests that particles escape the turbulence zone through advective loss rather than diffusive escape (at energies much less than a PeV, microscopic diffusion across B, even at the Bohm rate, is negligible compared to the advective losses). Then at flow speed ∼c/3, T_escape∼ 3 R_shock/c, leading to ε∼ (1/3) (B/δB)² (λ₀/R_shock). Taking the parameters to be unity yields ε∼ 1/3, in which case the accelerated spectrum is very hard: N(E) ∝E^−1.3 which is similar to the spectrum inferred for the radio emitting pairs in PWNe. Because of the strong radial mixing observed in the MHD models of PWNe, the spectrum created in this turbulent acceleration zone will fill the whole body of the SYN emitting nebulae, which provides a natural explanation of the observed lack of gradients in the radio spectral index of the Crab (Bietenholz & Kronberg 1992). Thus, we revive, in a modern form, the suggestion of Barnes & Scargle (1973) that the wisp motions are responsible for a part of the particle acceleration required to account for PWNe SYN emission.

Turning this idea into a full physical model may require extending the MHD models to 3D. All models to date have been 2D, axisymmetric, with an exclusively toroidal magnetic field. Particles then cannot scatter radially, in order to sample adjacent moving eddies, but are confined to toroidal flux tubes, which move in and out in radius around average positions in the outflow (and inflow, at higher latitudes) in the turbulent region. It is possible that magnetic pumping in these flux tubes might substitute for particle scattering in quasi-isotropic turbulence as an acceleration mechanism (e.g. Melrose 1969; Kuijpers et al. 1997). However, polarization studies have already demonstrated that there are substantial poloidal magnetic fields in PWNe inner regions, at least in the Crab. The TS instability identified by Camus et al. (2009) is not likely to be restricted to 2D – perhaps the assumption of quasi-isotropic turbulence, with forcing scale ∼R_shock, is the most natural starting point for further pursuit of this idea.

It is clear that a key point for the model to succeed is that the proposed turbulently accelerated spectrum at low energy merges smoothly into the shock-accelerated spectrum at high energy. Since the macroscopic turbulence is closely associated with the shock, a smooth merger is at least thinkable. In this context, the meaning of ε_c changes completely with respect to the Kennel and Coroniti interpretation: this should not be identified with the post-shock temperature (which according to the present model is about an order of magnitude lower), but rather with the energy above which the low-energy accelerator (Fermi II in the turbulent flow, in the suggestion made here) fades out and relativistic diffusive shock acceleration (DSA) takes over.

7.4 How do we move forward?

Multi-D high-resolution MHD and PIC simulations of the TS region separating the freely expanding wind from the subsonically expanding nebulae will shed light on the possibly turbulent flow and on whether that turbulence can act as the low-energy accelerator discussed previously. A useful first step would be to couple the diffusion–advection equation in energy space with the MHD calculations. Observationally, IR measures of the young PWNe, using the Herschel space telescope (Pilbratt 2009) as well as near-IR instruments on the ground and in space, will greatly improve our understanding of the physics behind the broken power laws in energy space inferred here for the injected particle distributions. These instruments will have more than sufficient angular resolution to resolve the PWNe, thus emphasizing the need for multidimensional models to quantitatively interpret the observations.

For all the objects studied, the plasma injection rates, here inferred from the evolutionary models, not simply taken from the average over the nebular lifetimes, exceed 10⁵ times the Goldreich–Julian rate (Table 6). Such high rates can be understood only as the result of pair creation within the PSRs' magnetospheres; given the spectral continuity in the SEDs, models that rely on the radio and far-IR emitting particles, which dominate the injection rates, being accelerated from the thermal plasma in and around the nebulae, are less likely than injection from the PSRs. Our results support and extend a long-held suspicion (Gallant et al. 2002), previously based only on analysis of the Crab nebula, that the particle loss rates from young, high voltage PSRs are substantially in excess of the inferences of particle outflow from all known magnetospheric pair creation models. The answer may lie in intrinsic time dependence of pair creation within PSR polar caps (e.g. Timokhin 2009) or in the so far unexplored outer magnetosphere accelerators associated with the boundary layer return currents separating the closed and open field lines (Arons 2009). Fermi gamma-ray observations are and will be useful in constraining magnetospheric pair creation models.

NB was in part supported by NASA through Hubble Fellowship grant HST-HF-01193.01-A, awarded by the Space Telescope Science Institute, which is operated by the Association of Universities for Research in Astronomy, Inc., for NASA, under contract NAS 5-26555; and from NORDITA. JA acknowledges the support of the US National Science Foundation (AST-0507813), the NASA Astrophysics Theory Program (NNG06108G) and the taxpayers of California.

REFERENCES