A MODEL OF THE SPECTRAL EVOLUTION OF PULSAR WIND NEBULAE

In this paper, we assume that the PWN is a uniform sphere expanding at a constant velocity $v_{\rm {PWN}}$. This assumption of constant velocity is the easiest way to take into account the expansion of the PWN, although the real behavior of the expansion must be more complex. Gelfand et al. (2009) investigated the dynamical evolution of the PWN surrounded by SN ejecta. They showed that the PWN expands at an almost constant velocity at first and then it turns to shrink and bounces in a late phase of the evolution after about 10 kyr from birth. We consider that the age of the PWN is younger than 10 kyr in this paper and that constant velocity is a good assumption in this range. The radius of the PWN $R_{\rm {PWN}}(t)$ at a time t is given by

$$\begin{equation} R_{\rm PWN}(t) = v_{\rm PWN} t. \end{equation} \tag{ 1 }$$

For the components inside the PWN, we assume that the PWN is composed of the magnetic field and the relativistic electron–positron plasma, both of which are injected from the pulsar inside the PWN. The evolution of the spin-down power L(t) is given by

$$\begin{equation} L(t) = L_{\rm 0} \left(1+\frac{t}{\tau _{0}} \right)^{-\frac{n+1}{n-1}}, \end{equation} \tag{ 2 }$$

where L₀ is the initial spin-down power and τ₀ is the initial spin-down time. Both parameters are fixed, if the current pulsar period P, its time derivative $\dot{P}$, braking index n, and age of the pulsar $t_{\rm {age}}$ are known, assuming that the moment of inertia of the pulsar is 10⁴⁵ g cm².

We divide the energy injection from the pulsar into the magnetic field energy $\dot{E}_{\rm {B}}$ and the relativistic particle energy $\dot{E}_{{e}}$ using the time-independent parameter η (0 ⩽ η ⩽ 1). The fraction parameter η is the ratio of the magnetic field energy injection to the spin-down power, i.e.,

$$\begin{eqnarray} L(t) & = & \dot{E}_{{e}}(t) + \dot{E}_{\rm {B}}(t) \nonumber \\ & = & (1 - \eta) \cdot L(t) + \eta \cdot L(t). \end{eqnarray} \tag{ 3 }$$

The fraction parameter η in our model is similar to the magnetization parameter σ in the KC model, although they are not the same. The magnetization parameter σ is the ratio $\dot{E}_{\rm {B}} / \dot{E}_{{e}}$ at the pulsar wind region immediately upstream from the termination shock. On the other hand, the fraction parameter η pertains to the ratio $\dot{E}_{\rm {B}} / (\dot{E}_{\rm {B}} + \dot{E}_{{e}})$ into the PWN region.

For the particle injection, we also need to determine the injection spectrum of the relativistic particles. Following Venter & de Jager (2006), we assume that the injection spectrum of the relativistic particles Q_inj(γ, t) obeys a broken power law

$$\begin{equation} Q_{\rm inj}(\gamma, t) = \left\lbrace \begin{array}{@{}ll} Q_{\rm 0}(t) (\gamma /\gamma _{\rm {b}})^{-p_{\rm 1}} & \mbox{ for \quad$\gamma _{\rm min} \le \gamma \le \gamma _{\rm b}$,} \\ Q_{\rm 0}(t) (\gamma / \gamma _{\rm {b}})^{-p_{\rm 2}} & \mbox{ for \quad$\gamma _{\rm b} \le \gamma \le \gamma _{\rm max}$,} \end{array} \right. \end{equation} \tag{ 4 }$$

where γ is the Lorentz factor of the relativistic electrons and positrons. We introduce time-independent parameters γ_min, γ_b, γ_max, p₁, and p₂ which are the minimum, break, and maximum Lorentz factors and the power-law indices at the low and high energy ranges of the injection spectra, respectively. Below, we use the terms low/high energy particles to refer to the particles which have energy lower/higher than γ_b. Note that most PWNe have values p₁ < 2 and p₂>2, which means that the energy of the particle injection is dominated by the particles γ = γ_b.

We require that the normalization $Q_{\rm {0}}(t)$ satisfies the following equation:

$$\begin{eqnarray} &&(1-\eta)L(t) = \int _{\gamma _{\rm min}}^{\gamma _{\rm max}} Q_{\rm inj}(\gamma, t) \gamma m_{e} c^{2} d\gamma, \end{eqnarray} \tag{ 5 }$$

where m_e and c are the mass of an electron (or positron) and the speed of light, respectively. From Equations (1), (4), and (5), the normalization parameter Q₀(t) is given by

$$\begin{eqnarray} Q_{\rm 0}(t) &=& \frac{L_{\rm 0} \cdot (1-\eta)}{m_{e}c^{2}} \cdot \Bigg(1+\frac{t}{\tau _{\rm 0}} \Bigg)^{-\frac{n+1}{n-1}} \Bigg(\frac{\gamma _{\rm b}^{2}(p_{\rm 1}-p_{\rm 2})}{(2-p_{\rm 1})(2-p_{\rm 2})}\nonumber\\ &&+ \frac{\gamma _{\rm b}^{p_{\rm 2}} \cdot \gamma _{\rm max}^{2-p_{\rm 2}}}{2-p_{\rm 2}} - \frac{\gamma _{\rm b}^{p_{\rm 1}} \cdot \gamma _{\rm min}^{2-p_{\rm 1}}}{2-p_{\rm 1}} \Bigg)^{-1}. \end{eqnarray} \tag{ 6 }$$

Because the magnetic field lines are stretching and winding, it is difficult to model the evolution of the magnetic field B(t) in the context of the uniform PWN. We have to solve the relativistic magnetohydrodynamics (MHD) equations to determine the realistic magnetic field evolution (e.g., Del Zanna et al. 2004). In this paper, for simplicity, we assume that the magnetic field evolution is determined in the form of the magnetic field energy conservation, i.e.,

$$\begin{eqnarray} \frac{4\pi }{3} (R_{\rm {PWN}}(t))^{3} \cdot \frac{(B(t))^{2}}{8\pi } & = & \int _0^{t} \eta L(t^{\prime }) dt^{\prime } \nonumber \\ & = & \eta E_{\rm spin}(t), \end{eqnarray} \tag{ 7 }$$

where E_spin(t) is the integrated spin-down energy at a time t. From Equations (1), (2), and (7), the magnetic field at a time t is given by

$$\begin{equation} B(t) = \sqrt{\frac{3 (n-1) \eta L_0 \tau _0}{(R_{\rm PWN})^3} \left\lbrace 1 - \left(1+ \frac{t}{\tau _{\rm 0}} \right)^{-\frac{2}{n-1}} \right\rbrace }. \end{equation} \tag{ 8 }$$

The magnetic field approximately evolves as B(t) ∝ t^−1.5 for t>τ₀. Note that although this magnetic field evolution model may be ad hoc, its behavior is very similar to those adopted in other works as we compare our model with other representative treatments (Rees & Gunn 1974; Kennel & Coroniti 1984a; de Jager et al. 2009; Venter & de Jager 2006).

Rees & Gunn (1974) considered that the magnetic field evolution is determined by the number of turns of the central pulsar. They considered that the magnetic field in the PWN was built up by the winding of field lines because of the pulsar spin. The total number of the turns of the magnetic field line is given by

$$\begin{equation} N(t) = \int ^t_{t_{\min }} \frac{\Omega (t^{\prime })}{2 \pi } dt^{\prime }. \end{equation} \tag{ 9 }$$

This number is proportional to the magnetic flux, which means

$$\begin{equation} B(t) \propto N(t) (R_{\rm PWN}(t))^{-2}. \end{equation} \tag{ 10 }$$

For t>τ₀, we obtain B(t) ∝ t^n/(1−n) so that B(t) ∝ t^−1.5 for n = 3 and B(t) ∝ t^−5/3 for n = 2.5. Our choice of the magnetic field evolution is not so different.

Kennel & Coroniti (1984a) considered the steady state spatial structure of PWN. In the KC model, the magnetic field increases with the distance from the pulsar till the magnetic pressure dominates over the particle pressure. Although the assumption of the steady state PWN makes it difficult to compare with our model, it seems to be natural to regard that the magnetic energy is a constant fraction of the total energy (∼ the particle energy), which means

$$\begin{equation} (B(t))^2 (R_{\rm PWN}(t))^3 \propto \frac{E_{\rm spin}(t)}{R_{\rm PWN}(t)}, \end{equation} \tag{ 11 }$$

where we include the adiabatic cooling of the total energy, so that B(t) ∝ t⁻² for t>τ₀. Equation (11) is alternatively interpreted as the magnetic flux conservation. The magnetic field decreases more rapidly than our model, but it is still close to our model.

In de Jager et al. (2009), they mentioned that the spatially averaged magnetic field strength of the PWN decreases as B(t) ∝ t^−1.3 in their calculation of non-relativistic MHD equations. This is again close to our model.

Finally, Venter & de Jager (2006) assumed that the magnetic field of the PWN evolved as B(t) = B₀/(1 + (t/T)^k) and Zhang et al. (2008) used this model with k = 0.5, i.e., B(t) ∝ t^−0.5 for t>T = 500 yr. This is somewhat different from the other models presented above and from ours. Based on these comparisons, we regard that our model is close to the reality, although other choices are also possible.

We assume that the distribution of the particles in the PWN is isotropic, so that the particle distribution function can be easily volume integrated to be described by the energy distribution N(γ, t). The evolution of the particle distribution N(γ, t) is given by the continuity equation in the energy space,

$$\begin{equation} \frac{ \partial }{ \partial t} N(\gamma, t) + \frac{ \partial }{ \partial \gamma } \left(\dot{\gamma }(\gamma, t) N(\gamma, t) \right) = Q_{\mathrm{inj}}(\gamma, t). \end{equation} \tag{ 12 }$$

We consider the cooling effects of the relativistic particles $\dot{\gamma }(\gamma, t)$ including the synchrotron radiation, the inverse Compton scattering off the cosmic microwave background (CMB) radiation, and the adiabatic expansion, i.e.,

$$\begin{equation} \dot{\gamma }(\gamma, t) = \dot{\gamma }_{\mathrm{syn}}(\gamma,t) + \dot{\gamma }_{\mathrm{IC}}(\gamma) + \dot{\gamma }_{\mathrm{ad}}(\gamma,t). \end{equation} \tag{ 13 }$$

The inverse Compton cooling $\dot{\gamma }_{\mathrm{IC}}(\gamma)$ does not depend on time because we consider that the target photon field is only the CMB. As the cooling effect, we do not include the inverse Compton scattering off the synchrotron radiation field because it is never a more important cooling process than the synchrotron cooling and because it costs too much computer power. Note that we include it in the calculation of the radiation spectrum, as described in Section 2.4.

The synchrotron cooling $\dot{\gamma }_{\mathrm{syn}}(\gamma,t)$ is given by (e.g., Rybicki & Lightman 1979)

$$\begin{equation} \dot{\gamma }_{\mathrm{syn}}(\gamma,t) = - \frac{4}{3} \cdot \frac{\sigma _{\rm {T}}^{}}{m_{{e}} c} \cdot U_{\mathrm{B}}(t) \gamma ^2, \end{equation} \tag{ 14 }$$

where σ_T is the Thomson cross section and U_B(t) is the magnetic field energy density calculated from Equation (8).

The inverse Compton cooling $\dot{\gamma }_{\mathrm{IC}}(\gamma)$ is given by (e.g., Blumenthal & Gould 1970)

$$\begin{eqnarray} \dot{\gamma }_{\mathrm{IC}}(\gamma) &=& {-} \frac{3}{4} \cdot \frac{\sigma _{\mathrm{T}}^{} h}{m_{{e}} c} \cdot \frac{1}{\gamma ^2} \int _0^{\infty } \nu _{\mathrm{fin}} d\nu _{\mathrm{fin}}\nonumber\\ &&\times\int _0^{\infty } \frac{n_{\mathrm{CMB}}(\nu _{\mathrm{ini}})}{\nu _{\mathrm{ini}}} \cdot f(q, \Gamma _{\rm {\epsilon }})\nonumber\\ &&\cdot \theta (1-q) \cdot \theta (q-1/4\gamma ^2) d\nu _{\mathrm{ini}}, \end{eqnarray} \tag{ 15 }$$

where h is the Planck's constant, $\nu _{\rm {ini}}$ and $\nu _{\rm {fin}}$ are the frequency of the CMB photons and that of scattered photons, respectively, and $n_{\rm {CMB}}(\nu _{\rm {ini}})$ is the distribution of the CMB described in Equation (19) below. $\Gamma _{\rm {\epsilon }} = 4 \gamma h \nu _{\rm {ini}} / (m_{{e}} c^2)$, $q = h \nu _{\rm {fin}} /(\Gamma _{\rm {\epsilon }} (\gamma m_{{e}} c^2 - h \nu _{\rm {fin}}))$, $f(q, \Gamma _{\rm {\epsilon }}) = 2q \ln q + (1+2q)(1-q) + 0.5 (1-q) (\Gamma _{\epsilon }q)^2 / (1+\Gamma _{\epsilon }q)$, and θ is the step function.

Finally, the adiabatic cooling $\dot{\gamma }_{\mathrm{ad}}(\gamma, t)$ is given by

$$\begin{eqnarray} \dot{\gamma }_{\mathrm{ad}}(\gamma, t) & = & {-} \frac{v_{\mathrm{PWN}}}{R_{\rm PWN}(t)} \cdot \gamma \nonumber \\ & = & {-} \frac{\gamma }{t}, \end{eqnarray} \tag{ 16 }$$

where we use Equation (1). Note that the adiabatic cooling is independent of the expansion velocity v_PWN.

As was stated in Section 1, we consider that it is more reasonable to treat the adiabatic loss rather than the escape loss. Zhang et al. (2008) considered an escape of the particles instead of an adiabatic loss. They treated the escape of the particles based on the Bohm diffusion. Because higher energy particles have larger Lamor radii, the high energy particles suffer from the escape loss in their models. On the other hand, the adiabatic loss is the dominant cooling process for the low energy particles in our model, since the high energy particles suffer from stronger radiative loss.

When the evolutions of the magnetic field and the particle distribution are determined, we can calculate the spectral evolution of the PWN. To compare the calculated luminosity with the observed flux, we assume that the radiation is isotropic. The radiation processes that we consider in our model are the synchrotron radiation and the inverse Compton scattering off the CMB (IC/CMB) and the synchrotron radiation (SSC).

The synchrotron radiation luminosity $L_{\nu,\rm syn}(t)$ is given by (e.g., Rybicki & Lightman 1979)

$$\begin{equation} L_{\nu,\rm syn}(t) = \int _{\gamma _{\rm min}}^{\gamma _{\rm max}} N(\gamma, t) P(\nu, \gamma, B(t)) d\gamma, \end{equation} \tag{ 17 }$$

where $P(\nu, \gamma, B(t)) = \sqrt{3} e^3 B(t) / (m_{e} c^2) F(\nu / \nu _{\rm c}(t))$ is the emissivity of the synchrotron radiation per particle, ν_c(t) = 3eB(t)γ²/(4πm_ec), F(x) = x∫^∞_xK_5/3(y)dy with $x = \nu / \nu _{\rm {c}}$, and K_5/3 being the modified Bessel function of order 5/3.

The inverse Compton scattering luminosity $L_{\nu _{\rm fin},\rm IC}(t)$ is given by (e.g., Blumenthal & Gould 1970)

$$\begin{eqnarray} L_{\nu _{\rm fin},\rm IC}(t) &=& \frac{3}{4} \cdot \frac{\sigma _{\rm T}^{} h \nu _{\rm fin} }{m_{e} c} \int _{\gamma _{\rm min}}^{\gamma _{\rm max}} \frac{N(\gamma, t)}{\gamma } d\gamma \nonumber\\ &&\times\int _0^{\infty } \frac{n_{\rm ph}(\nu _{\rm ini})}{\nu _{\rm ini}} \cdot f(q, \Gamma _{\rm \epsilon }) d\nu _{\rm ini}, \end{eqnarray} \tag{ 18 }$$

where $n_{\rm {ph}}(\nu _{\rm ini})$ is the distribution of the target photon fields including the CMB and the synchrotron radiation. The distribution of the CMB n_CMB(ν) is given by the blackbody distribution:

$$\begin{equation} n_{\rm CMB}(\nu) = \frac{8\pi }{c^3} \cdot \frac{\nu ^2}{\exp {(h \nu / k_{\rm {B}} T_{\rm {CMB}})}-1}, \end{equation} \tag{ 19 }$$

where k_B is the Boltzmann's constant and T_CMB is the temperature of the CMB. The distribution of the synchrotron radiation n_syn(ν, t) is given by

$$\begin{equation} n_{\rm syn}(\nu, t) = \frac{L_{\nu,\rm syn}(t)}{4\pi R^2_{\rm PWN}(t) c} \cdot \frac{1}{h \nu } \cdot \overline{U}, \end{equation} \tag{ 20 }$$

where $\overline{U} \sim 2.24$ is the average of U(x) = (3/2)∫¹₀(y/x)ln((x + y)/(x − y))dy in a spherical volume and x is the ratio of the distance from the center of the PWN to the PWN radius R_PWN(t) (cf. Atoyan & Aharonian 1996).

In Equation (20), Zhang et al. (2008) did not use R_PWN(t), but a new parameter r_syn (<R_PWN(t)), considering that the spread of the synchrotron radiation field is smaller than the radius of the PWN. Because only the synchrotron photons whose frequencies are below optical wavelengths contribute to the SSC flux and because the observed size of the optical synchrotron nebula is comparable to the size of the PWN, we use R_PWN(t) in Equation (20).

Figure 1 shows the calculated current spectrum of the Crab Nebula with the current observational data. The adopted parameters are shown in Table 1. As seen in Figure 1, the SSC flux is stronger than the IC/CMB flux in γ-rays. Note that the IC/CMB flux is almost in the Thomson regime, while the SSC flux is largely affected by the Klein–Nishina effect.

Zoom In Zoom Out Reset image size

Table 1. Parameters Used to Reproduce the Current Observed Spectrum of the Crab Nebula

Adopted Parameter	Symbol	Value
Current period (s)	$\mbox{\it P}$	3.31 × 10−2
Current period derivative (s s−1)	$\dot{\mbox{\it P}}$	4.21 × 10−13
Braking index	$\mbox{\it n}$	2.51
Age (yr)	$\mbox{\it t}_{\rm age}$	950
Expansion velocity (km s−1)	vPWN	1800
Fitted parameter
Fraction parameter ($\dot{\mbox{\it E}}_{\rm B} / (\dot{\mbox{\it E}}_{\rm B} + \dot{\mbox{\it E}}_{e})$)	η	0.005
Low energy power-law index at injection	$\mbox{\it p}_1$	1.5
High energy power-law index at injection	$\mbox{\it p}_2$	2.5
Maximum energy at injection	γmax	7.0 × 109
Break energy at injection	γb	6.0 × 105
Minimum energy at injection	γmin	1.0 × 102

Download table as:  ASCIITypeset image

We fit the data with the parameters η = 0.005, γ_max = 7.0 × 10⁹, γ_b = 6.0 × 10⁵, γ_min = 1.0 × 10², p₁ = 1.5, and p₂ = 2.5. The fraction parameter η governs the absolute values of the fluxes and the flux ratio of the inverse Compton scattering to the synchrotron radiation, as discussed in more detail in Section 3.3. The KC model derived σ ≪ 1 from the viewpoint of the current dynamical structure of the Crab Nebula, while we determine η ≪ 1 from the viewpoint of the spectral evolution. The parameters γ_max, γ_b, p₁, and p₂ are fixed to reproduce the observed synchrotron spectral shape, such as the spectral breaks and the photon indices, while γ_min should be regarded as an upper limit to reproduce the radio flux at the lowest frequency. In Section 3.4, these fitted parameters characterizing particle injection are discussed in detail.

In our calculation, the current magnetic field strength of the Crab Nebula turns out to be B_now = 85 μG, which is smaller than ∼300 μG used by Atoyan & Aharonian (1996). This difference of the magnetic field strength can be explained as follows. Atoyan & Aharonian (1996) adopted B_KC ∼ 300 μG from the KC model and adjusted the particle number to reproduce the observations. They applied roughly half a spin-down power compared with the KC model to reproduce the spectrum and thus the other half is missing. On the other hand, all the injected spin-down power is divided between the magnetic field and the particle energies in our model. If we adopt B_now = B_KC ∼ 300 μG, the synchrotron flux and also the SSC flux increase by about an order of magnitude. Note that the relativistic MHD simulation by Volpi et al. (2008) also indicates a smaller value of the spatially averaged magnetic field strength ∼100 μG, which is close to our value B_now = 85 μG.

Our model can calculate the past and future spectra of the Crab Nebula. All the calculated results in this section shown in Figures 2–5 are with the use of the same parameters as in Figure 1.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Figure 2 shows the spectral evolution. As seen in Figure 2, the synchrotron flux decreases with time and the SSC flux also decreases with time in accordance with the synchrotron flux, while the IC/CMB flux decreases more slowly than the SSC flux. We discuss their time dependence in Section 3.3. An important result in Figure 2 is that the flux ratio of γ-rays to X-rays increases with time. This supports the view that old PWNe can be observed as γ-ray sources with no or a weak X-ray counterpart.

The radio/optical observations of the Crab Nebula have suggested that the radio/optical flux of the Crab Nebula is decreasing with time. The inferred rate of the radio flux decrease is −0.17% ± 0.02% yr⁻¹ and almost independent of the frequency for a range 86–8000 MHz (Vinyaikin 2007). Our model predicts the current rate ∼−0.16% yr⁻¹, which is almost consistent with the observation. The inferred rate of the optical continuum flux decrease is −0.55% yr⁻¹ calibrated at 5000 Å (Smith 2003). Our model predicts the current rate ∼−0.24% yr⁻¹, which is a factor of 2 smaller than the observation. However, the trend that the decreasing rate increases with frequency matches the observations. This is because the optical emitting particles suffer from stronger synchrotron cooling than the radio emitting particles.

To understand the detailed features shown in Figure 2, we examine features of the evolution of the particle energy distribution in Figure 3. For comparison, the injected spectrum of the particles till 10 kyr without the cooling effects is plotted by the dot-dashed line in Figure 3. The evolution of the particle distribution is characterized differently in four energy ranges. First, for γ>10⁸, the particle number increases with time. This increase of the high energy particles relates to an increase of the IC/CMB flux in Figure 2 after 1 kyr above 10 TeV. Second, for 10⁵ < γ < 10⁸, the behavior of the particle distribution is complex. The particle number and the power-law index in this energy range do not monotonically change. The radiation spectrum in the range from infrared to X-rays reflects this complex behavior. Third, for 10² < γ < 10⁵, the change of the particle distribution is small. The difference between all the lines for several epochs is less than a factor of 2. This leads to an important conclusion that the radio flux decrease is mainly due to the decrease of the magnetic field, so that the observations of the radio flux decrease support our model of the magnetic field evolution. Lastly, for γ < 10², there exist particles with $\gamma < \gamma _{\rm {min}}$. This is because the adiabatic cooling is still effective at low energy. We discuss more about the features of the particle γ>10² in later paragraphs, which relate to the spectral evolution above 10⁷Hz. As discussed below, these features indicate the importance of considering the evolution of the particle injection rate and the cooling time due to various processes.

Figure 4 plots the evolution of cooling time $\tau (\gamma, t) = \gamma / \vert \dot{\gamma }(\gamma, t) \vert$ and helps to understand features of the particles with γ>10² in Figure 3. As seen in Figure 4, at an age of 1 kyr, the synchrotron cooling is effective for the particles with γ>10⁶, while the adiabatic cooling is effective with γ < 10⁶. However, because of a rapid decrease of the magnetic field energy density, the adiabatic cooling dominates at the γ < 10⁸ at an age of 10 kyr. From Equations (8), (14), and (16), τ_syn ∝ γ⁻¹t² for t < τ₀ and τ_syn ∝ γ⁻¹t³ for t>τ₀, while τ_ad ∝ t at any time.

We examine the increase of the particle number for γ>10⁸ shown in Figure 3. The increase means that the particle injection dominates the cooling effect in this energy range. Equation (12) is approximately expressed as (∂/∂t ∼ t⁻¹ = τ⁻¹_dyn and $\partial \dot{\gamma } / \partial \gamma \sim \tau ^{-1}_{\rm cool} = (\tau ^{-1}_{\rm syn} + \tau ^{-1}_{\rm IC} + \tau ^{-1}_{\rm ad})$)

$$\begin{equation} N(\gamma, t) \sim \frac{\tau _{\rm dyn} \tau _{\rm cool}}{\tau _{\rm dyn} + \tau _{\rm cool}} \cdot Q_{\rm {inj}}(\gamma, t). \end{equation} \tag{ 21 }$$

Because synchrotron cooling dominates in this energy range τ_cool ∼ τ_syn < τ_dyn, Equation (21) becomes $N \sim \tau _{\rm syn} Q_{\rm {inj}}$. For t < τ₀, the injection term behaves as $Q_{\rm {inj}} \propto \gamma ^{- p_2}$, so that $N \propto \gamma ^{- p_2 - 1} t^2$. For t>τ₀, the injection term behaves as $Q_{\rm {inj}} \propto \gamma ^{- p_2} t^{-7/3}$ for n = 2.5, so that $N \propto \gamma ^{- p_2 - 1} t^{2/3}$. This is the reason why the particle number increases with time and the particle distribution is softer than the injection distribution for γ>10⁸.

We next examine the steadiness of the particle number for 10² < γ < 10⁵ shown in Figure 3. For t < τ₀, because τ_cool ∼ τ_ad ∼ τ_dyn in Equation (21), we have $N \propto \gamma ^{- p_1} t$. For t>τ₀, because the injection term more rapidly decreases with time than the adiabatic cooling term, Equation (12) is approximated by

$$\begin{equation} \frac{ \partial }{ \partial t} N(\gamma, t) \sim \frac{ \partial }{ \partial \gamma } \left(\frac{\gamma }{t} N(\gamma, t) \right). \end{equation} \tag{ 22 }$$

When we assume a form N ∝ γ^−αt^β, we obtain β = 1 − α. For α = p₁ = 1.5, we have β = −0.5. The particle number increases until t ∼ τ₀, and then it decreases slowly with time as N ∝ t^−0.5 for 10² < γ < 10⁵.

The feature of the particle distribution for 10⁵ < γ < 10⁸ shown in Figure 3 is complex because the dominant cooling process changes with time as seen in Figure 4. Until an age of a few thousand years, the feature is similar to that for γ>10⁸ (Equation (21) with τ_cool ∼ τ_syn), i.e., the particle number increases with time and the particle distribution is softer than the injection distribution. After that, the feature becomes similar to that for 10² < γ < 10⁵ for t>τ₀ (Equation (22)) and the particle number decreases with time. Because the time when τ_ad ∼ τ_syn is different for different energy, the particle distribution becomes harder than $\propto \gamma ^{-p_2 - 1}$ at an age of 10 kyr.

In Figure 5, we show the evolution of each energy component: the particle energy, the magnetic field energy, the radiated energies via the synchrotron radiation and the inverse Compton scattering, and the wasted energy via the adiabatic expansion, all of which are normalized by the integrated spin-down energy E_spin(t). As shown in Figure 5, although almost all the injected energy goes to the particle energy, the particle energy decreases with time by the cooling effects. The magnetic field energy is always a constant fraction of the injected energy ηE_spin(t) as assumed in Equation (7). Most of the particle energy is cooled by the synchrotron radiation in the early phase (<1 kyr) because the cooling break due to the synchrotron cooling is a little smaller than γ_b as seen in Figure 3. On the other hand, most of the injected energy goes to the adiabatic loss in the late phase (>1 kyr), which means that most of the injected energy goes to the kinetic energy of the SN ejecta. The inverse Compton cooling is not important in the case of the Crab Nebula for t < 10 kyr, which is also seen in Figure 4. Note that the integrated spin-down energy E_spin(t) also increases with time, for example, E_spin(100 yr) = 9.1 × 10⁴⁸ erg, E_spin(1 kyr) = 3.9 × 10⁴⁹ erg, and E_spin(10 kyr) = 5.5 × 10⁴⁹ erg.

We discuss the evolution of powers of the synchrotron radiation and the inverse Compton scattering and their dependence on the fraction parameter η, which will help to apply our model to other PWNe in future.

We assume t>τ₀ in this section aiming to discuss primarily old PWNe. For simplicity, we also assume that the particle energy is expressed as N(γ_b, t)γ²_bm_ec² = (1 − η)ξ(t)E_spin(t), where ξ(t) < 1 accounts for the cooling effects. As seen in Figure 5, $N(\gamma _{\rm b}, t) \gamma ^2_{\rm b} m_{e} c^2 \sim \rm const.$ and $\xi (t) \sim \rm const. \sim 0.1$ is a fairly good approximation. Moreover, we assume that powers of all radiation mechanisms are dominated by the emission from the particles γ = γ_b. This assumption is valid when the low energy power-law index <2 and the high energy power-law index >3 in the particle distribution.

The power ratio of the inverse Compton scattering to the synchrotron radiation is given by

$$\begin{equation} \frac{P_{\rm IC}(t)}{P_{\rm syn}(t)} \sim f_{\rm KN} \frac{U_{\rm ph}(t)}{U_{\rm B}(t)}, \end{equation} \tag{ 23 }$$

where f_KN < 1 represents the Klein–Nishina effect and U_ph(t) is the energy density of the target photon field. f_KN ∼ 1 for the IC/CMB and f_KN < 0.1 for the SSC.

The power of the synchrotron radiation is given by (e.g., Rybicki & Lightman 1979)

$$\begin{eqnarray} P_{\rm {syn}}(t) & \sim & \frac{4}{3} \sigma _{\rm {T}} c \gamma ^2_{\rm {b}} U_{\rm {B}}(t) \cdot \gamma _{\rm {b}} N(\gamma _{\rm {b}}, t) \nonumber \\ & = & \frac{\sigma _{\rm {T}} \gamma _{\rm b}}{\pi m_{{e}} c} \cdot \frac{\eta (1 - \eta)}{R^3_{\rm PWN}(t)} \xi (t) E_{\rm spin}^2(t)\nonumber\\ &\propto & \eta (1 - \eta) t^{-3}, \end{eqnarray} \tag{ 24 }$$

where γ_bN(γ_b, t) is the number of the particles around γ = γ_b and we use Equation (7). Note that although the power of the synchrotron radiation has been conventionally compared with the spin-down power L(t), the power of the synchrotron radiation relates to the integrated spin-down energy E_spin(t) rather than the instantaneous spin-down power L(t) in Equation (24).

Now Equation (23) becomes, in the case of the IC/CMB (U_ph(t) = U_CMB),

$$\begin{eqnarray} \frac{P_{\rm IC/CMB}(t)}{P_{\rm syn}(t)} & \sim & f_{\rm KN} \left(\frac{B(t)}{3\,\mu {\rm G}} \right)^{-2} U_{\rm CMB} \nonumber \\ & \propto & \frac{R^3_{\rm PWN}(t)}{\eta E_{\rm spin}(t)} \propto \eta ^{-1} t^3. \end{eqnarray} \tag{ 25 }$$

In the case of the SSC, we estimate the synchrotron radiation energy density U_syn(t) as

$$\begin{equation} U_{\rm syn}(t) \sim \frac{P_{\rm syn}(t)}{4 \pi R^2_{\rm PWN}(t) c} \propto \eta (1 - \eta) t^{-5}. \end{equation} \tag{ 26 }$$

Then, we obtain

$$\begin{equation} \frac{P_{\rm SSC}(t)}{P_{\rm syn}(t)} \sim \frac{f_{\rm KN} \xi (t)}{3} \cdot \frac{t}{\tau _{\rm syn}(\gamma _{\rm b}, t)} \cdot \frac{v_{\rm PWN}}{c} \cdot \frac{1 - \eta }{\eta } \propto (1- \eta) t^{-2}, \end{equation} \tag{ 27 }$$

where τ_syn(γ_b, t) = N(γ_b, t)γ²_bm_ec²/P_syn(t) is the synchrotron cooling time of the particles with γ = γ_b at time t.

All the powers P_syn(t) ∝ η(1 − η)t⁻³, P_IC/CMB(t) ∝ (1 − η), and P_SSC(t) ∝ η(1 − η)²t⁻⁵ depend on the fraction parameter η and time in different ways. These general characteristics of the radiation powers make it possible to estimate the fraction parameter η from the viewpoint of the spectral evolution. For example, the time dependence of P_SSC(t) and P_IC/CMB(t) can explain that the SSC flux is dominant in high energy γ-rays in Figure 1. Using the same parameters as in Figure 1, we find $P_{\rm IC/CMB}(t_{\rm {age}}) / P_{\rm syn}(t_{\rm {age}}) \sim 10^{-3}$ and $P_{\rm SSC}(t_{\rm {age}}) / P_{\rm syn}(t_{\rm {age}}) \sim 10^{-2}$.

In de Jager & Djannati-Ataï (2008), they suggested that some of the TeV γ-ray sources without an X-ray counterpart can be old PWNe because the X-ray emitting particles are cooled more rapidly than the TeV γ-ray emitting particles by the synchrotron cooling. In our model, old PWNe can also be the γ-ray sources with no or weak X-ray counterpart because the IC/CMB is almost time independent. However, the same result comes from a different reason. Because the number of the high energy particles is even increasing in our model as seen in Figure 3, this result comes from the rapid decrease of the magnetic field. PWNe which are old and have a small fraction parameter η can be recognized as the γ-ray source without an X-ray counterpart. Note that P_IC/CMB(t) slowly decreases with time, when we properly consider the cooling effects ξ(t), as seen in Figure 2 at 3 kyr and 10 kyr.

The fitted parameters p₁, p₂, γ_max, γ_b, and γ_min relate to the broken power-law injection in Equation (4). The above parameters include the information about the acceleration at the termination shock and the physics of the pulsar wind and the pulsar magnetosphere. We discuss them in the framework of our model.

For the low and high energy power-law indices at injection p₁ = 1.5 and p₂ = 2.5, we determine these values to reproduce the observed radio and X-ray photon indices, Γ_r and Γ_X. Radio photon index Γ_r is related to p₁, but X-ray photon index Γ_X is not simply related to p₂ because of the synchrotron cooling. Because it is difficult to make a large spectral break Γ_X–Γ_r = ΔΓ>0.5 from the single power-law injection (e.g., Reynolds 2009), we adopt the broken power-law injection.

For the maximum energy γ_max = 7.0 × 10⁹, we determine this value to reproduce the observed spectral break at ∼100 MeV. The maximum energy γ_max is conventionally given by comparing the acceleration time with the cooling time or the size of the acceleration site with the Lamor radius. Both conditions give γ_max ∼ 10¹⁰ using the magnetic field strength ∼100 μG and the size of the acceleration site ∼0.1 pc. The fitted value is near the limit of theoretical expectation.

For the break energy γ_b = 6.0 × 10⁵, we determine this value to reproduce the observed spectral break around optical wavelengths. Although the KC model related γ_b to the bulk Lorentz factor of the pulsar wind immediately upstream from the termination shock, our model does not allow this connection as discussed below.

Finally, for the minimum energy γ_min = 1.0 × 10², we determine this value to reproduce the flux of the observed minimum frequency at radio wavelengths. Because p₁>1, the particles around γ_min determine the injection of the particle number as $\dot{N}_e(t) \sim Q_{\rm inj}(\gamma _{\rm min}, t) \gamma _{\rm min}$. The particle number conservation outside the pulsar light cylinder that leads to the particle number flux into the PWN is much larger than the Goldreich–Julian number flux $\dot{n}_{\rm GJ}$, which is the particle number flux from the pulsar polar cap. In our model, the pair production multiplicity at an age of a thousand year $\kappa \sim (L(t) / \gamma _{\rm b} m_{e} c^2)(\gamma _{\rm b} / \gamma _{\rm min})^{p_1 - 1} \sim 10^6$ from Equations (4) and (5). Theoretically, the multiplicity is estimated as κ ∼ 10³–10⁵. Our value κ ∼ 10⁶ is somewhat large. The mean energy of the injected particles $L(t) / \dot{N}_e(t) \sim \gamma _{\rm b} m_{e} c^2(\gamma _{\rm b} / \gamma _{\rm min})^{- p_1 + 1}$ is significantly smaller than γ_bm_ec².

Atoyan & Aharonian (1996) did not use the broken power-law injection. They divided the particles inside the PWN into the low energy particles as the relic electrons and the high energy particles as the wind electrons. This assumption can reduce the multiplicity, but the origin of the relic electrons becomes another problem. It is difficult to discuss anything more about this problem from our spectral evolution model. One thing that we should note is that the relation between the radio flux decrease and the magnetic field evolution is always kept, i.e., the relation is independent of how and when the radio emitting particles are injected. This is because the low energy particles whose power-law index is p₁ do exist from the observation and their distribution is hardly changed by the cooling effect, as discussed in Section 3.2.

Our model of the magnetic field evolution is somewhat ad hoc. As discussed in Section 2.2, however, the time dependence of the magnetic field strength B ∝ t^−a for t>τ₀ seems to be in a range of 1.5 ⩽ a ⩽ 2.0 from other theoretical considerations and we adopt a = 1.5. Moreover, because our result for the radio flux decrease of the Crab Nebula is almost consistent with the observations, our model of the magnetic field evolution may be near the truth.

For the assumption of the uniform PWN, many non-uniform structures have been observed, such as the filamentary structures and the spatial variations of photon indices. However, for the calculation of the total spectrum of the PWN, the energetics of the PWN is important in the lowest order. We consider that the assumption of the uniform PWN is reasonable for the calculation of the total spectral evolution of the PWN.

For the injection spectrum of the particle distribution, the acceleration of the particles is an unsolved problem and we adopt the broken power-law injection. It should be noted that one of the important conclusions in our study, that old PWNe can be observed as γ-ray sources with no or weak X-ray counterpart, is not affected by the broken power-law assumption. This is because low energy particles do not contribute to X-ray and high energy γ-ray emissions.

The use of the time-independent parameters γ_min, γ_b, and γ_max can be improved as time-dependent parameters, because the physical condition of pulsar wind termination shock may change with the decrease of the spin-down power of the pulsar. Considering the time dependence of γ_max(t), Venter & de Jager (2006) use the condition r_L(t) < 0.5r_s with time-dependent magnetic field, where r_L is the electrons' Larmor radius and r_s is the radius of the termination shock. As discussed in Section 3.4, our model satisfies this condition. Both γ_min and γ_b are important parameters because these may include the information about the pulsar magnetosphere and the pulsar wind. However, the time dependences of them are uncertain. For simplicity, we used all of them as the time-independent parameters in this paper.

Constant velocity expansion is a good assumption for young PWNe, although the expansion of the PWN should be calculated by taking into account the environment of the PWN (e.g., Gelfand et al. 2009). In our model, the magnetic field evolution explicitly depends on the expansion velocity v_PWN (see Equation (8)). To understand a little more about the effects of the expansion evolution, we study how the Crab Nebula would be observed in the context of the constant velocity expansion. The Crab Nebula is one of the sources without an observable supernova remnant (SNR) shell. It may be because the surrounding interstellar medium is less dense than other PWNe with observable SNR shells. If the Crab Nebula were in different surroundings, the expansion velocity would also be affected. That is, if it were in less or more dense surroundings, the expansion velocity would be more rapid or slower.

An example of an expansion that is twice as rapid is shown in Figure 6. All the parameters except the expansion velocity are the same as in Figure 1. An example of a half-velocity expansion is shown in Figure 7 and all parameters except for the expansion velocity are the same as in Figure 1. Note that both of them are hypothetical PWNe, not the Crab Nebula itself. For the spectrum of the rapid expansion case, the absolute value of the flux becomes smaller and the flux ratio of the inverse Compton scattering to the synchrotron radiation is larger than the real Crab Nebula shown in Figure 1, and for the spectrum of the slow expansion case, vice versa. Comparing the particle distribution in Figure 6 with that in Figure 7, the low energy particles take the same distribution, but the high energy particle distribution in Figure 7 is steeper than that in Figure 6. These spectral behaviors against the expansion velocity can be understood from Equations (8), (14), and (16). In our model, the magnetic field becomes smaller when the expansion velocity becomes larger. This leads to the difference in the absolute flux and the synchrotron cooling which changes the high energy particle distribution. On the other hand, because the adiabatic cooling does not depend on the absolute value of the expansion velocity, the low energy particle distribution does not change.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Lastly, Atoyan & Aharonian (1996), included the infrared photons and the starlight for the target photon fields of the inverse Compton scattering. Although these soft photons can significantly contribute to the γ-ray flux of other PWNe (e.g., Porter et al. 2006), this is not the case for the Crab Nebula. Because the Crab Nebula is located far away from the galactic center (∼10 kpc) and galactic plane (∼200 pc), the energy density of these soft photon fields is less than solar neighborhood (∼8 kpc from the galactic center). Even if we assume that the energy density of the infrared photons and the starlight is the same as the solar neighborhood, the inverse Compton scattering off these photon fields contributes less than 30% of the current total γ-ray flux. Note that it also does not much affect the γ-ray spectrum, when the SSC flux decreases if the energy density of these soft photon fields is less than half that in the solar neighborhood since inverse Compton scattering off the CMB dominates there.

In this paper, we built a model of the spectral evolution of PWNe and applied the Crab Nebula as a calibrator of our model. We solved the equation for the particle distribution function considering adiabatic and radiative losses with a simple model of magnetic field evolution.

The flux decrease of the γ-rays is more moderate than radio to X-rays, because the magnetic field decreases rapidly, which implies that old PWNe can be observed as γ-ray sources with no or weak X-ray counterpart. de Jager & Djannati-Ataï (2008) obtained the same result but for a different reason: that the X-ray emitting particles are cooled more rapidly than TeV γ-ray emitting particles.

The current observed spectrum of the Crab Nebula is reconstructed when the fraction parameter has a small value η = 0.005. This is consistent with the prediction of the magnetization parameter σ ≪ 1 obtained by Kennel & Coroniti (1984a). They obtained σ ≪ 1 from the viewpoint of the current dynamical structure of the Crab Nebula, while we determine η ≪ 1 from the viewpoint of the spectral evolution.

The smaller value of the current magnetic field B_now = 85 μG is needed to reconstruct the observed spectrum of the Crab Nebula. This is consistent with that of Volpi et al. (2008) for the spatially averaged magnetic field strength ∼100 μG from their relativistic MHD simulation, but smaller than ∼300 μG in most other papers (e.g., Atoyan & Aharonian 1996).

Our model can predict the spectral evolution of the Crab Nebula, and the observed flux decrease of the Crab Nebula at radio wavelengths can be explained by our model. This conclusion does not depend on the assumption of the broken power-law injection and gives the validity of our magnetic field evolution model. The observed flux decrease of the Crab Nebula in optical wavelengths is somewhat larger than our model, but the trend that the decreasing rate increases with frequency matches observations.

The minimum energy γ_min is related to the pair production multiplicity in the pulsar magnetosphere, since low energy particles are assumed to be injected in the same way as high energy particles in our model. Our result for the minimum energy γ_min = 1.0 × 10² and the low energy power-law index p₁ = 1.5 means that the multiplicity κ ∼ 10⁶ is necessarily larger than other models which adopt a separate origin of low energy particles.