Broadband Photon Spectrum and its Radial Profile of Pulsar Wind Nebulae

TL;DR

This study tests the standard 1-D steady model for Pulsar Wind Nebulae against observations, finding it struggles to simultaneously reproduce spectra and spatial profiles, and suggests including additional physical processes for improvement.

Contribution

The paper critically evaluates the 1-D steady model for PWNe, demonstrating its limitations and proposing the need for incorporating spatial diffusion and other effects.

Findings

The 1-D model cannot simultaneously fit the spectrum and surface brightness.

Parameters constrained by the spectrum lead to smaller nebula sizes than observed.

Ignoring radio and optical data allows some fits but results in unphysical advection times.

Abstract

The observed radial profiles of the X-ray emission from Pulsar Wind Nebulae (PWNe) have been claimed to conflict with the standard one-dimensional (1-D) steady model. However, the 1-D model has not been tested to reproduce both the volume-integrated spectrum and the radial profile of the surface brightness, simultaneously. We revisit the 1-D steady model and apply it to PWNe 3C 58 and G21.5-0.9. We find that the parameters of the pulsar wind, the radius of the termination shock $r_{\rm s}$ and magnetization $\sigma$, greatly affect both the photon spectrum and radial profile of the emission. We have shown that the parameters constrained by the entire spectrum lead to a smaller X-ray nebula than observed one. We have also tested the case that reproduces only the observations in X and gamma-rays, ignoring the radio and optical components. In this case, there are parameter sets that reproduce both the spectrum and emission profile, but the advection time to the edge of the nebula becomes much smaller than the age. Our detailed discussion clarifies that the standard 1-D steady model has severe difficulty to reproduce both the volume-integrated spectrum and the surface brightness simultaneously. This implies that the model should be improved by taking into account extra physical processes such as spatial diffusion of particles. Additionally, we calculate the surface brightness profile of the radio, optical and TeV gamma-rays. The future observations in these wavelengths are also important to probe the spatial distributions of the relativistic plasma and the magnetic field of PWNe.

Tables

Table 1. Table 1 Obtained parameters for the test calculation with σ = 10 − 4 𝜎 superscript 10 4 \sigma=10^{-4} shown in Fig. 3 . See Fig. 3 for the other parameters. In these parameter sets, the magnetic field always reaches its maximum at r = r N 𝑟 subscript 𝑟 N r=r_{\rm N}

Radius of termination shock (pc)	$0.05$	$0.1$	$0.2$	$0.4$
Total number of pairs (${10}^{50}$)	$9.5$	$3.4$	$0.97$	$0.25$
Advection time (yr)	$6800$	$2400$	$690$	$180$
Total pressure at $r=r_N$ (${10}^{-10}$dyn ${\text{cm}}^{-2}$)	$36$	$16$	$5.4$	$1.5$
Maximum magnetic field ($\mu$G)	$250$	$100$	$32$	$8.7$
Averaged magnetic field ($\mu$G)	$240$	$88$	$26$	$6.8$
$r_{\mathrm{eq}}/r_N$	$1.4$	$2.9$	$5.8$	$12$

Table 2. Table 2 Obtained parameters for the test calculation with r s = 0.1 subscript 𝑟 s 0.1 r_{\rm s}=0.1 pc shown in Fig. 4 . See Fig. 3 for the other parameters.

Magnetization Parameter	${10}^{-2}$	${10}^{-3}$	${10}^{-4}$	${10}^{-5}$	${10}^{-6}$
Total number of pairs (${10}^{50}$)	$0.37$	$1.6$	$3.4$	$4.0$	$4.2$
Advection time (yr)	$280$	$1200$	$2400$	$2900$	$2900$
Total pressure at $r=r_N$ (${10}^{-10}$erg cm-3)	$2.3$	$8.9$	$16$	$21$	$23$
Maximum magnetic field ($\mu$G)	$130$	$130$	$100$	$43$	$14$
Magnetic field at $r=r_N$ ($\mu$G)	$69$	$120$	$100$	$43$	$14$
Averaged magnetic field ($\mu$G)	$86$	$120$	$87$	$34$	$11$
$r_{\mathrm{eq}}/r_N$	$0.29$	$0.91$	$2.89$	$9.13$	$28.9$

Table 3. Table 3 Parameters in our calculations.

		3C 58		G21.5-0.9
Given parameters	Symbol	broadband	alternative	broadband	alternative
Spin-down luminosity (erg s-1)	$L_{\mathrm{sd}}$	$3.0\times {10}^{37}$		$3.5\times {10}^{37}$
Distance (kpc)	$D$	${2.0}^a$		${4.8}^b$
Radius of the nebula (pc)	$r_N$	$2.0$		$0.9$
Fitting parameters
Break energy (eV)	$E_b$	$4.1\times {10}^{10}$	$1.1\times {10}^{11}$	$2.6\times {10}^{10}$	$3.1\times {10}^{12}$
Low energy power-law index	$p_1$	$1.26$		$1.1$
High energy power-law index	$p_2$	$3.0$	$3.2$	$2.3$	$3.0$
Radius of the termination shock (pc)	$r_s$	$0.13$	$0.26$	$0.05$	$0.1$
Magnetization parameter	$\sigma$	$1.0\times {10}^{-4}$	$1.0\times {10}^{-3}$	$2.0\times {10}^{-4}$	$3.0\times {10}^{-2}$
Obtained parameters
Initial bulk Lorentz factor	${\gamma }_u$	$7.3\times {10}^3$	$5.5\times {10}^5$	$2.1\times {10}^4$	$1.8\times {10}^7$
Pre-shock density (cm-3)	$n_u$	$1.1\times {10}^{-11}$	$5.1\times {10}^{-16}$	$1.1\times {10}^{-11}$	$3.5\times {10}^{-18}$
Pre-shock magnetic field ($\mu$G)	$B_u$	$0.79$	$1.2$	$3.1$	$19$
Maximum Energy (eV)	$E_{\max }$	$9.5\times {10}^{13}$	$3.0\times {10}^{14}$	$1.4\times {10}^{14}$	$1.7\times {10}^{15}$
Advection time (yr)	$t_{\mathrm{adv}}$	$1500$	$330$	$800$	$38$
Averaged magnetic field ($\mu$G)	$B_{\mathrm{av}}$	$31$	$21$	$120$	$61$
Total pressure at $r=r_N$ (erg cm-3)	$p_{\mathrm{tot}}$	$3.7\times {10}^{-10}$	$8.3\times {10}^{-11}$	$2.1\times {10}^{-9}$	$1.2\times {10}^{-10}$
$r_{\mathrm{eq}}/r_N$		$3.8$	$2.4$	$2.3$	$0.37$

Table 4. Table 4 Parameters in the calculations for the Crab Nebula.

Given parameters a	Symbol	Crab
Spin-down luminosity (erg s-1)	$L_{\mathrm{sd}}$	$5\times {10}^{38}$
Distance (kpc)	$D$	2.0
Radius of the nebula (pc)	$r_N$	1.8
Fitting parameters a
Break energy (eV)	$E_b$	$2.5\times {10}^{11}$
Low energy power-law index	$p_1$	$1.6$
High energy power-law index	$p_2$	$2.4$
Radius of the termination shock (pc)	$r_s$	$0.1$
Magnetization parameter	$\sigma$	$5.0\times {10}^{-3}$
Obtained parameters
Initial bulk Lorentz factor	${\gamma }_u$	$3.2\times {10}^3$
Pre-shock density (cm-3)	$n_u$	$1.7\times {10}^{-9}$
Pre-shock magnetic field ($\mu$G)	$B_u$	$30$
Maximum Energy (eV)	$E_{\max }$	$2.7\times {10}^{15}$
Advection time (yr)	$t_{\mathrm{adv}}$	$380$
Averaged magnetic field ($\mu$G)	$B_{\mathrm{av}}$	$234$
Total pressure at $r=r_N$ (erg cm-3)	$p_{\mathrm{tot}}$	$2.0\times {10}^{-9}$
$r_{\mathrm{eq}}/r_N$		$0.45$

Equations

$$L_{\mathrm{sd}}=4\pi r_{\mathrm{s}}^2{n}_{\mathrm{u}}{u}_{\mathrm{u}}{\gamma }_{\mathrm{u}}{m}_{\mathrm{e}}{c}^3\left(1+\sigma \right),$$

$$\sigma \equiv \frac{B_{\mathrm{u}}^2/4\pi }{n_{\mathrm{u}}{u}_{\mathrm{u}}{\gamma }_{\mathrm{u}}{m}_{\mathrm{e}}{c}^2}.$$

$$n_{\mathrm{d}}=\frac{n_{\mathrm{u}}{u}_{\mathrm{u}}}{u_{\mathrm{d}}},$$

$$u_{\mathrm{d}}^2=\frac{8{\sigma }^2+10\sigma +1+\sqrt{64{\sigma }^2{\left(\sigma +1\right)}^2+20\sigma \left(\sigma +1\right)+1}}{16\left(\sigma +1\right)},$$

$$P_{\mathrm{d}}=\frac{n_{\mathrm{u}}m{c}^2{u}_{\mathrm{u}}^2}{4{\gamma }_{\mathrm{d}}{u}_{\mathrm{d}}}\left[1+\sigma \left(1-\frac{{\gamma }_{\mathrm{d}}}{u_{\mathrm{d}}}\right)\right],$$

$$B_{\mathrm{d}}=B_{\mathrm{u}}\frac{{\gamma }_{\mathrm{d}}}{u_{\mathrm{d}}},$$

$$\sqrt{1+u^2(r)}\left(\delta +\frac{\left(u_{\mathrm{d}}^2/\sigma \right)-\frac{1}{2}}{u_{\mathrm{d}}^2+\frac{1}{4}}{\left(\frac{u(r)r^2}{u_{\mathrm{d}}{r}_{\mathrm{s}}^2}\right)}^{-\frac{1}{3}}+\frac{u_{\mathrm{d}}}{u(r)}\right)$$

$$={\gamma }_{\mathrm{d}}\left(\delta +\frac{\left(u_{\mathrm{d}}^2/\sigma \right)-\frac{1}{2}}{u_{\mathrm{d}}^2+\frac{1}{4}}+1\right),$$

$$n_{\mathrm{tot}}(r)$$

$$B(r)$$

$$P(r)$$

$$r\simeq r_{\mathrm{eq}}\equiv r_{\mathrm{s}}/\sqrt{3\sigma },$$

$$n_{\mathrm{tot}}(r)=\int n(E,r)dE.$$

$$\displaystyle n(E,r_{\rm s})=\left\{\begin{array}[]{ll}\displaystyle\frac{n_{0}}{E_{\rm b}}\left(\frac{E}{E_{\rm b}}\right)^{-p_{1}}&\left(E_{\rm min}<E<E_{\rm b}\right)\\ \displaystyle\frac{n_{0}}{E_{\rm b}}\left(\frac{E}{E_{\rm b}}\right)^{-p_{2}}&\left(E_{\rm b}<E<E_{\rm max}\right)\end{array}\right.,$$

$$C_1\equiv [\frac{1}{p_1-1}\left\{{\left(\frac{E_{\min }}{E_{\mathrm{b}}}\right)}^{1-p_1}-1\right\}$$

$${+\frac{1}{p_2-1}\left\{1-{\left(\frac{E_{\max }}{E_{\mathrm{b}}}\right)}^{1-p_2}\right\}]}^{-1}.$$

$$E_{\max }=e{B}_{\mathrm{u}}{r}_{\mathrm{s}}=\sqrt{\frac{e^2}{c}\frac{L_{\mathrm{sd}}\sigma }{1+\sigma }}.$$

$$P_{\mathrm{d}}=\frac{1}{3}\int En(E,r_{\mathrm{s}})dE,$$

$$n_{\mathrm{d}}=\frac{3L_{\mathrm{sd}}}{16\pi r_{\mathrm{s}}^{2}c{u}_{\mathrm{d}}{\gamma }_{\mathrm{d}}{E}_{\mathrm{b}}\left(1+\sigma \right)}\frac{C_2}{C_1}\left[1+\sigma \left(1-\frac{{\gamma }_{\mathrm{d}}}{u_{\mathrm{d}}}\right)\right],$$

$$C_2=[\frac{1}{2-p_1}\left\{1-{\left(\frac{E_{\min }}{E_{\mathrm{b}}}\right)}^{2-p_1}\right\}$$

$${+\frac{1}{p_2-2}\left\{1-{\left(\frac{E_{\max }}{E_{\mathrm{b}}}\right)}^{2-p_2}\right\}]}^{-1}.$$

$${\gamma }_{\mathrm{u}}=\frac{4}{3}\frac{E_{\mathrm{b}}}{m{c}^2}{\gamma }_{\mathrm{d}}\frac{C_1}{C_2}{\left[1+\sigma \left(1-\frac{{\gamma }_{\mathrm{d}}}{u_{\mathrm{d}}}\right)\right]}^{-1},$$

$$B_{\mathrm{u}}={\left[\frac{L_{\mathrm{sd}}}{c{r}_{\mathrm{s}}^2}\frac{\sigma }{1+\sigma }\right]}^{1/2},$$

$$n_{\mathrm{u}}=\frac{9L_{\mathrm{sd}}}{64\pi r_{\mathrm{s}}^{2}c{\gamma }_{\mathrm{d}}^2{m}_{\mathrm{e}}{c}^2\left(1+\sigma \right)}{\left(\frac{C_2}{C_1}\frac{m_{\mathrm{e}}{c}^2}{E_{\mathrm{b}}}\right)}^2$$

$$\times {\left[1+\sigma \left(1-\frac{{\gamma }_{\mathrm{d}}}{u_{\mathrm{d}}}\right)\right]}^2.$$

$$t_{\mathrm{adv}}={\int }_{r_{\mathrm{s}}}^{r_{\mathrm{N}}}\frac{dr}{cu(r)}.$$

$$u(r)\frac{\partial n(E,r)}{\partial r}$$

$$F_{\nu }=\frac{1}{D^2}{\int }_{r_{\mathrm{s}}}^{r_{\mathrm{N}}}{j}_{\nu }\left(r\right)r^{2}dr,$$

$$B_{\nu }\left(s\right)=2{\int }_{\max \left(r_{\mathrm{s}},s\right)}^{r_{\mathrm{N}}}\frac{j_{\nu }\left(r\right)rdr}{\sqrt{r^2-s^2}},$$

$${\nu }_{\mathrm{b}}$$

$$E_{\mathrm{cut}}(r)=\frac{E_{\max }}{1+\frac{1}{5}\frac{E_{\max }}{E_{\mathrm{bof}}}\left((\frac{r}{r_{\mathrm{s}}}{)}^5-1\right)},$$

$$E_{\mathrm{bof}}$$

$${\nu }_{\mathrm{c}}$$

Full text

Our detailed discussion clarifies that the standard 1-D steady model has severe difficulty to reproduce both the volume-integrated spectrum and the surface brightness simultaneously. This implies that the model should be improved by taking into account extra physical processes such as spatial diffusion of particles.

Additionally, we calculate the surface brightness profile of the radio, optical and TeV gamma-rays. The future observations in these wavelengths are also important to probe the spatial distributions of the relativistic plasma and the magnetic field of PWNe.

1. INTRODUCTION

Pulsar Wind Nebulae (PWNe) are extended sources around a rotation powered pulsar. They show a broadband spectrum from radio to $\gamma$-rays so that they would contain very high energy non-thermal particles (Gaensler & Slane, 2006; Kirk et al., 2009; Bühler & Blandford, 2014). The central pulsar releases its rotational energy as the pulsar wind, which is a highly relativistic magnetized electron-positron outflow, and plays a role of the energy source of the PWN. The strong termination shock, which is formed by the interaction between the pulsar wind and external supernova remnant (or interstellar medium), has been supposed to produce the non-thermal electrons and positrons, and they emit the synchrotron radiation (Rees & Gunn, 1974) and the inverse Compton emission (de Jager & Harding, 1992).

PWNe are characterized by a center-filled morphology caused by confinement in a supernova remnant (SNR), which is associated with the progenitor of the central pulsar. Most PWNe are detected as spatially extended sources in radio and X-rays. While the radio spectral index is almost spatially homogeneous (e.g., Bietenholz et al., 1997; Bietenholz & Bartel, 2008), the X-ray spectral index increases with the distance from the pulsar (e.g., Bocchino & Bykov, 2001; Slane et al., 2004; Schöck et al., 2010). Recently, PWNe are found to be very bright sources in also TeV $\gamma$-rays (Kargaltsev et al., 2013). One-zone models of PWNe, which are invented by Pacini & Salvati (1973) and developed by several authors (e.g., Bednarek & Bartosik, 2003; Chevalier, 2005; Tanaka & Takahara, 2010; Bucciantini et al., 2011; Vorster et al., 2013), can reproduce broadband spectra of entire nebulae well. The one-zone models do not describe the spatial distribution of the emission (e.g., Amato et al., 2000), thus it is indispensable to invest a model which includes the spatial structure of nebulae.

The model of Kennel & Coroniti (1984a, b, hereafter KC84s) has been established as a standard model of PWNe. They constructed a steady-state and 1-D magnetohydrodynamic model of the Crab Nebula. Also, assuming the particle acceleration at the termination shock, they calculated the evolution of non-thermal pairs along the flow and the synchrotron emission from advected particles. Adopting the model of KC84s, Atoyan & Aharonian (1996) succeeded in reproducing the entire photon spectrum including the inverse Compton component. However, as Reynolds (2003) suggested, it is unclear whether KC84s can apply to general PWNe other than the Crab Nebula. Furthermore, Slane et al. (2004) showed that the model of KC84s disagrees with the observed radial profile of the X-ray spectral index in 3C 58. They suggested that the radial profile of the X-ray spectral index in the model should change more rapidly, and the X-ray nebula size becomes more compact than the observation. Note that the entire photon spectrum was not taken into account as a model constraint in Slane et al. (2004).

This ignited the revision of KC84s. Tang & Chevalier (2012) introduced the effect of the spatial diffusion of the particles, and reproduced the X-ray radial profile. Porth et al. (2016) supported this idea via 3-D magnetohydrodynamic and test-particle simulations. In each of the studies, simultaneous verification of the entire spectrum and the spatial profile of the emission is not discussed, thus there is no consensus on the spatial structure in PWN models so far. In order to advance the study, it is essential to clarify controversial points in the simple steady 1-D model before introducing nontrivial effects such as the particle diffusion.

We choose two objects, 3C 58 and G21.5-0.9, to examine the 1-D steady model. Both of the PWNe show the feature that the extent of the X-ray emission is the same as the radio one, in contrast to the Crab Nebula, in which the observed size shrinks with increasing frequency. Furthermore, the 1-D steady model has been never applied to those two PWNe. Our purpose is to make the validity of the 1-D steady model clear for general PWNe, so that those two PWNe are suitable targets for testing.

In this paper, we revisit the 1-D steady model and calculate the photon spectrum and its radial profile numerically. In Section 2, we review the 1-D steady model of PWNe based on KC84s. The parameter dependence of this model is investigated in Section 3. The application to the two observed sources (3C 58 and G21.5-0.9) is presented in Section 4.

We discuss the 1-D steady modeling in Section 5 and summarize our results in Section 6.

2. Model

In this paper, we adopt the one-dimensional magnetohydrodynamical (MHD) model constructed by KC84s. We numerically solve the energy distribution of the electron–positron pair plasma along the outflow obtained by MHD equations. From the resultant spectral distribution of the pairs, we calculate the photon emission from those pairs, resultant surface brightness, and volume-integrated spectrum. In this section, we review the the 1-D steady model and present the method to calculate the radial evolution of the energy distribution of electrons and positrons.

Here, we consider that the PWN is a steady outflow and has a radius $r_{\rm N}$. The relativistic magnetized wind emitted from the central pulsar forms a strong termination shock at a radius $r_{\rm s}$. Assuming that the pre-shock plasma is cold, the wind property is regulated by three quantities, number density in the comoving frame $n$, bulk Lorentz factor $\gamma$, and magnetic field in the lab frame $B$ at the shock. Almost all of the pulsar spin-down luminosity $L_{\rm sd}$ is converted to the wind luminosity as

[TABLE]

where $u\equiv\sqrt{\gamma^{2}-1}$, the subscript u denotes values at just upstream of the shock, and $\sigma$ is the ratio of the magnetic energy flux to the particle energy flux at the upstream of the shock,

[TABLE]

The magnetic field of the wind is dominated by the toroidal component (e.g., Goldreich & Julian, 1969), and the upstream plasma is highly relativistic ($u_{\rm u}/\gamma_{\rm u}\simeq 1$), which means the downstream temperature is relativistic (adiabatic index $4/3$). The Rankine-Hugoniot jump conditions provide the values in the downstream (KC84s) as

[TABLE]

[TABLE]

[TABLE]

[TABLE]

where the subscript d denotes the values at just downstream of the shock, and $P$ is the thermal pressure. For $\gamma_{\rm u}\gg 1$ and $\sigma\ll 1$, we obtain $u_{\rm d}/\gamma_{\rm d}\simeq 1/3$, which coincides with the well-known result in the relativistic hydrodynamics.

As boundary conditions, we adopt Equations (3)–(6) at the radius $r=r_{\rm s}$, and solve the steady state and spherical symmetric MHD equations. Under the toroidal field approximation and the adiabatic assumption, the MHD equations are integrable. After some algebra with introducing $\delta\equiv u_{\rm d}/(\sigma u_{\rm u})$, we obtain (KC84s)

[TABLE]

from which we obtain the radial profile of the four velocity $u(r)$ (or equivalently $\gamma(r)$). In the strong shock approximation, $u_{\rm d}$ is a function of only $\sigma$ as shown in Equation (4). If $\delta\ll 1$ is established, the above flow Equation (7) is depicted by only one parameter $\sigma$ independently of $n_{\rm u}$ and $u_{\rm u}$. We calculate the radial profile of $u(r)$ numerically, while KC84s neglected $\delta$ and adopted $\gamma_{\rm d}\simeq 1$ in the downstream. Then, the MHD conservation laws provide the other quantities as follows:

[TABLE]

where $n_{\rm tot}(r)$ is the comoving number density in the wind.

In Figure 1, the radial profiles of $u(r)$ and $B(r)$ in our test calculations are shown for $\sigma=10^{-6}$, $10^{-5}$, $10^{-4}$, $10^{-3}$ and $10^{-2}$. The results are almost the same as the behavior shown in KC84s because $\delta\ll 1$ for all the cases. For $\sigma\ll 1$, at a small radius, the pressure ratio $\beta_{\rm pl}\equiv B^{2}/8\pi P$ ($\ll 1$ at $r=r_{\rm s}$) gradually increases with radius as $\beta_{\rm pl}\propto r^{2}$. In that regime, the outflow behaves as $u\propto r^{-2}$ (equivalently $n_{\rm tot}\propto r^{0}$), $B\propto r$ and $P\propto r^{0}$. At the radius

[TABLE]

$\beta_{\rm pl}$ becomes unity, namely the magnetic pressure starts to dominate. Outside this radius, the radial four speed is approximately constant with $B\propto r^{-1}$ and $P\propto r^{-8/3}$. Thus, the magnetic field has a maximum value at $r\simeq r_{\rm eq}$ as shown in Figure 1.

The energy distribution of electron–positron pairs $n(E,r)$ is calculated consistently with the MHD model,

[TABLE]

At the inner boundary $r=r_{\rm s}$, $n(E,r_{\rm s})$ is assumed to have a broken power-law energy distribution at injection as following:

[TABLE]

where the parameters are the break energy $E_{\rm b}$, minimum energy $E_{\rm min}$, maximum energy $E_{\rm max}$, and two power-law indices $p_{1}$ and $p_{2}$ for low and high energy regions, respectively. The normalization $n_{0}=C_{1}n_{\rm d}$ is determined by Equation (12) as

[TABLE]

While the origin of the radio spectral component may be different from that for the X-ray and optical components as discussed in KC84s, in this paper, we use a broken power-law distribution that has been adopted by the one-zone studies (e.g., Tanaka & Takahara, 2010).

The observed radio spectral index almost uniquely gives the index $p_{1}$, which is generally lower than 2 (e.g., Salter et al., 1989). In this case, the particles with energies $\sim E_{\rm min}$ dominate the particle number. For simplicity, we fix the minimum energy as $E_{\rm min}=10m_{\rm e}c^{2}$ and leave the pair density problem (c.f. Tanaka & Takahara, 2013a). The particles above $E_{\rm b}$ may be produced via the shock acceleration (Spitkovsky, 2008). The maximum energy $E_{\rm max}$ is determined by the same method in KC84s; the energy at which a gyro radius is equal to the shock radius provides

[TABLE]

The pressure obtained by $n(E,r_{\rm s})$,

[TABLE]

should satisfy Equation (5). We thus obtain

[TABLE]

where

[TABLE]

Notify that $\gamma_{\rm d}$ or $u_{\rm d}$ is already given as a function of only $\sigma$ (see Equation (4)), so that the quantities in the upstream are written with the six parameters $L_{\rm sd}$, $\sigma$, $r_{\rm s}$, $E_{\rm b}$, $p_{1}$, and $p_{2}$ as

[TABLE]

[TABLE]

and

[TABLE]

From the above relations, the functional shape described in Equation (15) is also written with the six parameters.

In this 1-D model, we have a unique parameter $r_{\rm s}$, which is not in the one-zone models, and $r_{\rm s}$ significantly affects the results as will be shown in Section 3. The flow solution given by Equation (7) provides the advection time as

[TABLE]

While the age of the PWN is an important parameter in the one-zone time-dependent models such as Tanaka & Takahara (2010), the parameter $r_{\rm s}$ in our steady model adjusts the advection time, which may be close to the age of the PWN.

The radial evolution of $n(E,r)$ and photon emission are calculated with the numerical code used in Sasaki et al. (2015), which are based on the time-dependent code in Asano & Mészáros (2011) (see also Asano & Mészáros, 2012). Taking into account the Klein–Nishina effect on the inverse Compton (IC) cooling, the code can follow the temporal evolution of the energy distribution along the stream. Transforming the elapsed time into radius as $dr=cu(r)dt$, our calculation is practically equivalent to solve the steady transport equation (e.g. Parker (1965); Ginzburg & Syrovatskii (1964))

[TABLE]

where $\dot{E}_{\rm syn}$ and $\dot{E}_{\rm IC}$ are the energy loss rates due to synchrotron radiation and IC scattering, respectively. The three terms of the right hand side of Equation (25) represent the effects of the radiative cooling, adiabatic cooling and volume expansion, respectively.

Here, we have used the solution for $u(r)$ with adiabatic approximation, namely the radiative cooling is assumed not to affect the dynamics of the flow. This approximation is valid when the cooling time for the particles with $E\sim E_{\rm b}$ is longer than the advection time (KC84s). Most of the results shown in this paper safely satisfy this condition.

The spectral emissivity $j_{\nu}(r)$ per unit volume is obtained consistently with the energy distribution $n(E,r)$, the magnetic field profile $B(r)$, and the interstellar radiation field (ISRF) with the Klein–Nishina effect. The model of the ISRF is taken from GALPROP v54.1 (Vladimirov et al., 2011, and references therein), in which the results of Porter & Strong (2005) are adopted, as shown in Figure 2. The ISRF is assumed as uniform and isotropic in the PWN, and we neglect the synchrotron self-Compton, which significantly contribute for only limited cases like the Crab Nebula (Torres et al., 2013). We neglect the contribution of bremsstrahlung as well, because the density of pairs is low enough (Atoyan & Aharonian, 1996).

Neglecting the emission from the pre-shock wind, the photon spectrum of the entire nebula $F_{\nu}$ is given by

[TABLE]

where $D$ is the distance to the PWN from us. The surface brightness $B_{\nu}$ is also given by

[TABLE]

where $s$ is the distance perpendicular to the line of sight from the central pulsar.

Our numerical code is checked by reproducing the result of Atoyan & Aharonian (1996) for the Crab Nebula. If we adopt the same assumption as Atoyan & Aharonian (1996), the resultant spectrum agree with the observed one. Note that our more conservative model leads to a slightly dimmer flux for the inverse Compton component as shown in Appendix.

3. Parameter Dependence in the 1-D Steady Model

In this section, we discuss how the entire spectrum $L_{\nu}\equiv 4\pi D^{2}F_{\nu}$ and the surface brightness $B_{\nu}$ depend on the parameters in this 1-D steady model. There are some previous studies that discussed the parameter dependence of the 1-D model. KC84s have already discussed how the parameters $u_{\rm u}$ and $\sigma$ change the total synchrotron luminosity $\int L_{\nu}d\nu$ (not the spectral distribution). While Schöck et al. (2010) have studied the X-ray spatial profile for different $r_{\rm s}$ assuming that the flow velocity decreases as a power-law of $r$ independently of $\sigma$, we investigate the parameter dependence consistently with the MHD flow solution. We focus on the dependence on the parameters $r_{\rm s}$ and $\sigma$, which largely affect the spatial structure of the emission. In this section, the nebula size $r_{\rm N}=2.0\ {\rm pc}$ is fixed. The external photon field is taken from the model for G21.5-0.9 in Figure 2. For four parameters out of the six parameters in our model, we adopt a parameter set as $L_{\rm sd}=10^{38}\ {\rm erg}~{}{\rm s}^{-1}$, $E_{\rm b}=10^{5}m_{\rm e}c^{2}$, $p_{1}=1.1$, and $p_{2}=2.5$, and change $r_{\rm s}$ or $\sigma$ below.

3.1. Characteristic frequencies and energies of the model

First, we introduce some typical particle energies and corresponding photon frequencies, and their dependence on the model parameters are discussed. In this subsection, we limit the discussion within the case of $r_{\rm N}<r_{\rm eq}$ (i.e., $r_{\rm N}/r_{\rm s}<(3\sigma)^{-1/2}$). In this case, the magnetic field has a maximum value $3B_{\rm u}r_{\rm N}/r_{\rm s}$ at the edge of the nebula. The cooling effect for pairs with energies $E=E_{\rm b}$ is also found to be negligible. We find the first typical frequency, the intrinsic break frequency,

[TABLE]

In the case of $B\propto r$ or equivalently $u\propto r^{-2}$, the energy of pairs injected with $E=E_{\rm max}$ decreases with $r$ via the synchrotron cooling as

[TABLE]

where the burn-off energy

[TABLE]

(KC84s). At $r=r_{\rm N}$, the maximum energy becomes $E_{\rm cut}(r_{\rm N})\simeq 5E_{\rm bof}(r_{\rm s}/r_{\rm N})^{5}$. Then, the cooling frequency at the outer boundary is given by

[TABLE]

Above $\nu_{\rm c}$, the entire spectrum should show the softening behavior.

The maximum particle energy decreases following Equation (29), while the magnetic field increases as $B(r)=3B_{\rm u}(r/r_{\rm s})$. The typical synchrotron frequency $\propto B(r)E_{\rm cut}(r)^{2}$ peaks at

[TABLE]

where we have assumed $E_{\rm max}\ll 5E_{\rm bof}$. With $r_{\rm pk}$, the cut-off frequency in the synchrotron spectrum is obtained as

[TABLE]

above which the flux decreases exponentially.

3.2. Shock radius dependence

One-zone models do not include the shock radius as a parameter. The value $r_{\rm s}$ is a characteristic parameter in the 1-D model. For the Crab Nebula (Schweizer et al., 2013), Vela (Helfand et al., 2001), and MSH 15-52 (Yatsu et al., 2009), a possible shock structure (inner ring) is detected with X-ray observations. In most of PWNe, however, the shock radii are not observationally constrained well. In Figure 3, we show the $r_{\rm s}$-dependences of the entire spectrum $L_{\nu}$ and the X-ray surface brightness with $\sigma=10^{-4}$. The synchrotron component dominates below $10^{20}$ Hz ($\sim 400$ keV) and the IC component dominates in $\gamma$-rays. Note that in the case for $r_{\rm s}=0.05$ pc, $E_{\rm cut}(r)<E_{\rm b}$ beyond $r\gtrsim 1.5$ pc , i.e., the adiabatic approximation is invalid. In Table 1, for various $r_{\rm s}$, we summarize the advection time, the volume-averaged magnetic field $B_{\rm av}$ as,

[TABLE]

and the maximum magnetic field $B_{\rm max}=B(r_{\rm N})$ because we have $r_{\rm eq}>r_{\rm N}$ from Equation (11) with $\sigma=10^{-4}$.

First, as a benchmark case, we take up the case of $r_{\rm s}=0.1$ pc (the red solid line in Figure 3), which leads to $r_{\rm eq}=5.7$ pc. There are two breaks in the synchrotron spectrum: the intrinsic break at $\nu_{\rm b}\sim 2.2\times 10^{12}$ Hz corresponding to $E_{\rm b}$ and the cooling break at $\nu_{\rm c}\sim 3.6\times 10^{14}\ {\rm Hz}$. Adopting the average magnetic field and the advection time, the cooling break energy of pairs is

[TABLE]

The corresponding cooling break frequency is written as

[TABLE]

The above value obtained with a one-zone-like treatment roughly agree with our results.

The analytical descriptions of the spectral indices $\alpha$ ($F_{\nu}\propto\nu^{-\alpha}$) are $(p_{1}-1)/2$ below $\nu_{\rm b}$, and $(p_{2}-1)/2$ for $\nu_{\rm b}<\nu<\nu_{\rm c}$. The calculated spectrum agrees with those values. In one-zone models, the index above $\nu_{\rm c}$ steepened by $1/2$. However, in the 1-D models following $B\propto r$ and $u\propto r^{-2}$, the spectral change $\Delta=(p+7)/18$ is slightly different from 1/2 (Reynolds, 2009), which also agrees with our result.

Let us remark the flux at $\nu=\nu_{\rm cut}$. When the synchrotron cooling is efficient for particles of $E=E_{\rm max}$, almost all of those energies are released by photon emission until $r=r_{\rm pk}$. Since the energy density of pairs, which have an energy $E=E_{\rm max}$ at $r=r_{\rm s}$, is estimated as $L_{\rm sd}(E_{\rm max}/E_{\rm b})^{2-p_{2}}$, the synchrotron luminosity at $\nu\sim\nu_{\rm cut}$ is calculated as $\nu_{\rm cut}L_{\nu_{\rm cut}}\sim 1.7\times 10^{-2}L_{\rm sd}$ for $\sigma=10^{-4}$ and $p_{2}=2.5$, which seems consistent with the flux in Figure 3. The above estimate does not depend on $r_{\rm s}$, which agrees with the results for $r_{\rm s}\leq 0.2$ pc, where $\nu_{\rm c}<\nu_{\rm cut}$.

A particle of energy $E$ emits synchrotron photons of frequency $\nu\propto E^{2}B$, and power $p_{\rm syn}\propto E^{2}B^{2}$, then the spectral emissivity $j_{\nu}$ is proportional to $n(E)E^{2}B^{2}(dE/\nu)$. When we can assume $\sigma\ll 1$ and $n_{\rm tot}\propto r^{0}$, $j_{\nu}\propto n_{0}E_{\rm b}^{p-1}\nu^{-(p-1)/2}B(r)^{(p+1)/2}$ ($p$ is index of the pair energy distribution) at the energy range where the cooling effect is negligible. The constant density implies that the magnetic field behaves as $B\sim B_{\rm d}r/r_{\rm s}\propto L_{\rm sd}^{1/2}r_{\rm s}^{-2}\sigma^{1/2}r$. Since we can treat $C_{2}$ as a constant for $p_{1}<2<p_{2}$ and $E_{\rm min}\ll E_{\rm b}\ll E_{\rm max}$, we obtain $n_{0}\propto L_{\rm sd}r_{\rm s}^{-2}E_{\rm b}^{-1}$. Finally, we obtain the entire specific luminosity ($L_{\rm\nu}\sim 4\pi\int_{r_{\rm s}}^{r_{\rm N}}drr^{2}j_{\nu}$) as

[TABLE]

where $p$ is $p_{1}$ for $\nu<\nu_{\rm b}$, and $p_{2}$ for $\nu_{\rm b}<\nu<\nu_{\rm c}$. At $\nu=\nu_{\rm b}\propto L_{\rm sd}^{1/2}\sigma^{1/2}E_{\rm b}^{2}r_{\rm s}^{-2}r_{\rm N}$ (Equation (28)),

[TABLE]

where the absolute value, which is difficult to evaluate analytically, is provided by our numerical results. Above $\nu=\nu_{\rm c}\propto L_{\rm sd}^{-3/2}\sigma^{-3/2}r_{\rm s}^{10}r_{\rm N}^{-9}$ (Equation (31)), using the above formula, the spectrum behaves as $\nu L_{\nu}\sim\nu_{\rm c}L_{\nu_{\rm c}}(\nu/\nu_{\rm c})^{-5(p_{2}-2)/9}$ so that

[TABLE]

for $\nu_{\rm c}<\nu<\nu_{\rm cut}$. Those formulae well agree with our numerical results.

The cooling frequency strongly depends on $r_{\rm s}$. The frequencies $\nu_{\rm b}$ and $\nu_{\rm c}$ have similar values for $r_{\rm s}=0.05$ pc, while $\nu_{\rm c}$ and $\nu_{\rm cut}$ merge for $r_{\rm s}=0.2$ pc. For $r_{\rm s}=0.4$ pc, $\nu_{\rm c}$ becomes higher than $\nu_{\rm cut}$. Namely, both of the radiative and adiabatic cooling effects are negligible in this case. This leads to the lower flux at $\nu=\nu_{\rm cut}$ for $r_{\rm s}=0.4$ pc than the fluxes for smaller $r_{\rm s}$.

The value of $E_{\rm cut}(r_{\rm N})$ increases with $r_{\rm s}$, because the synchrotron cooling becomes less effective. Reflecting this, the IC spectra show a soft-to-hard evolution with $r_{\rm s}$. The high-energy cut-off of the IC component is determined by the maximum energy of pairs. Electron–positron pairs expend only a small fraction of their energy in the IC emission. Since the $L_{\rm sd}$ and $E_{\rm b}$ are common for the examples in Figure 3, the IC flux in the low energy range ($10^{20}-10^{22}{\rm Hz}$) is basically proportional to the total number of corresponding low-energy particles in the nebula. As shown in the left panel in Figure 1, a flow with a small ratio of $r_{\rm N}/r_{\rm s}$ reaches the edge of the nebula before significant deceleration. Consequently, the advection time becomes shorter for a smaller $r_{\rm N}/r_{\rm s}$ as shown in Table 1. If we can neglect the cooling effect, the total particle number $\propto t_{\rm adv}L_{\rm sd}E_{\rm b}^{p_{1}-2}$ decreases with $r_{\rm s}$. This effect is seen as the flux growth with increasing $t_{\rm adv}$ below $\sim 10^{22}$ Hz. For $r_{\rm s}=0.05$ pc, the synchrotron cooling is crucial ($E_{\rm cut}(0.7r_{\rm N})\lesssim E_{\rm b}$), which practically reduces the particle number above $E_{\rm b}$ in the nebula. Therefore, the IC flux in this case does not follow the aforementioned trend of the flux growth.

The surface brightness profile in X-rays (see the right panel of Figure 3) is regulated by the synchrotron cooling. For a smaller $r_{\rm s}$, the stronger magnetic field results in a compact X-ray profile. The X-ray extent is proportional to $r_{\rm s}^{10/9}$ for $r_{\rm s}\leq 0.2$ pc as explained as follows. When the cooling effect is significant, $E_{\rm cut}\propto E_{\rm bof}(r_{\rm s}/r)^{5}\propto\sigma^{-1}r_{\rm s}^{6}r^{-5}$, while the magnetic field behaves as $B\propto\sigma^{1/2}r_{\rm s}^{-2}r$. For a given frequency $\nu\propto BE^{2}$, the maximum radius to emit photons of $\nu$ is proportional to $\sigma^{-1/6}r_{\rm s}^{10/9}$. This supports the X-ray size growth with $r_{\rm s}^{10/9}$ for the case where the cooling effect is significant. For $r_{\rm s}=0.4$ pc, the synchrotron cooling is negligible effect for the X-ray profile. In order to reconcile the fact that the X-ray extent is comparable to the radio nebula, a large $r_{\rm s}$ is preferable, though the synchrotron component becomes very hard and dim in this case.

As shown in Table 1, the total pressure $P_{\rm tot}\equiv 4Pu^{2}+P+B^{2}/8\pi$ at the outer boundary, which may balance with the pressure outside the nebula, decreases with $r_{\rm s}$ by roughly an order of magnitude. Since the uncertainty of the current observation of the outside pressure is larger than this variance, it may be difficult to constrain the value of $r_{\rm s}$ directly.

3.3. $\sigma$ dependence

Figure 4 shows the $\sigma$ dependences of the entire spectrum and the X-ray surface brightness. While the shock radius is fixed to $r_{\rm s}=0.1$ pc and the other parameters are the same as those in the previous subsection, the value of $\sigma$ changes from $10^{-6}$ to $10^{-2}$. As is shown in the figure, a complicated behavior appears in the spectral shape with increasing $\sigma$. The change of $\sigma$ modifies the profiles of the emission through two processes: the strength of magnetic field (see Equation (22)) and the deceleration profile as shown in Fig. 1. The magnetic field strength affects the typical synchrotron frequency and the cooling efficiency. The flow velocity profile adjusts the radius $r_{\rm eq}$, where the magnetic field becomes maximum, and the advection time, which controls the total energy in the nebula and the cooling efficiency (the ratio of the cooling time to the advection time).

For $\sigma\leq 10^{-4}$, the radius $r_{\rm eq}$ is outside $r_{\rm N}$ (see Table 2, note that $r_{\rm N}/r_{\rm s}$ corresponds to 20 in Figure 1), so that the behaviors of the characteristic frequencies are well explained by Equations (28), (31), and (33) as $\nu_{\rm b}\propto\sigma^{1/2}$, $\nu_{\rm c}\propto\sigma^{-3/2}$, and $\nu_{\rm cut}\propto\sigma^{11/10}$. In the case of $\sigma=10^{-6}$, the frequency $\nu_{\rm c}$ is much higher than $\nu_{\rm cut}$, so that the power-law portion for $\nu_{\rm c}<\nu<\nu_{\rm cut}$ is absent. Below $\nu_{\rm cut}$, the formulae of Equations (37) and (38) well represent the spectral behavior for $\sigma\leq 10^{-4}$. In the frequency range of $\nu_{\rm c}<\nu<\nu_{\rm cut}$, the spectrum practically follows Equation (39) for the cases of $\sigma=10^{-5}$ and $10^{-4}$. As $\sigma$ increases, the resultant stronger magnetic field leads to a large efficiency of the energy release during the advection time. The peak of the synchrotron flux at $\nu=\nu_{\rm c}$ grows with $\sigma^{p_{2}-2}$ (see Equations (31) and (37)) for $\sigma\leq 10^{-4}$, accompanying the shift of $\nu_{\rm c}$ to a lower frequency.

The spectral behavior deviates from the above trends for $\sigma\geq 10^{-3}$. This is because the radial evolution of the magnetic field can be no longer approximated by $B\propto r$. As shown in Table 2, the radius $r_{\rm eq}$ is inside the nebula radius $r_{\rm N}$ in this parameter region unlike the discussion in Section 3.1. The magnetic pressure prevents the flow from deceleration at $r>r_{\rm eq}$, and the adiabatic cooling starts to play its role. The decline of the advection time (see Table 2) leads to the reduction of the total energy in the nebula. As a result, above $\sigma=10^{-4}$, the synchrotron peak flux turns into decreasing and the cooling frequency turns into increasing. Therefore, we cannot decrease $\nu_{\rm c}$ extremely low in the 1-D steady model differently from the one-zone models.

In the case of $r_{\rm eq}<r_{\rm N}$, the contribution to the entire spectrum is dominated by the emission from non-thermal pairs at $r<r_{\rm eq}$. At $r>r_{\rm eq}$, electrons/positrons lose their energies via adiabatic cooling as $E\propto(r/r_{\rm eq})^{-2/3}$ and the magnetic field decays as $\propto r^{-1}$. Thus, the energy loss rate rapidly decreases as $r^{-10/3}$ so that the emission beyond $r=r_{\rm eq}$ is almost negligible. From Equation (29), assuming $E_{\rm max}\gg E_{\rm bof}$, we obtain $E_{\rm cut}(r_{\rm eq})\simeq 5E_{\rm bof}(r_{\rm s}/r_{\rm eq})^{5}\simeq 5E_{\rm bof}(3\sigma)^{5/2}$. In this case, the cooling frequency for enough large $\sigma$ that $r_{\rm eq}$ is smaller than $r_{\rm N}$ may be estimated with $E_{\rm cut}$ at $r=r_{\rm eq}$ as

[TABLE]

Although the estimate of the value of $\nu^{\rm(eq)}_{\rm c}$ almost agrees with our numerical results, the dependence on $\sigma$ is a little bit strong. This is because the case for $\sigma=10^{-2},10^{-3}$ are in a marginal situation ($r_{\rm eq}\sim r_{\rm N}$). Nevertheless, it is certain that the cooling frequency for $\sigma>10^{-4}$ ($r_{\rm eq}<r_{\rm N}$) becomes large with increasing $\sigma$.

In the frequency range of $\nu^{\rm(eq)}_{\rm c}<\nu<\nu_{\rm cut}$ for a larger $\sigma$, the spectrum is harder than the analytical estimate $\alpha=(5p_{2}-1)/9$ for $r_{\rm eq}<r_{\rm N}$. On the other hand, since $r_{\rm pk}\sim r_{\rm s}$, the frequency $\nu_{\rm cut}$ well agrees with the analytical estimate of Equation (33) even for a larger $\sigma$. For $\sigma\geq 10^{-4}$, as we have mentioned in the previous subsection, the luminosity around $\nu_{\rm cut}$ is almost independent of $\sigma$, while the peak luminosity at $\nu=\nu_{\rm c}$ decreases with increasing $\sigma$ following the decline of the total energy in the nebula. Those complicated effects lead to the spectral hardening between $\nu_{\rm c}$ and $\nu_{\rm cut}$.

The peak flux of the IC component declines monotonically with increasing $\sigma$. Below the spectral break frequency $10^{23}$ Hz, which corresponds to the photon energy emitted by particles of $E_{\rm b}$ interacting with dust photons, all the model curves for $\sigma\leq 10^{-4}$ almost overlap each other. In this range of $\sigma$, the flow profiles are almost the same, so that the IC emission processes, which does not directly depend on the magnetic field, are common. On the other hand, above $10^{23}$ Hz, the softening of the photon spectrum with increasing $\sigma$ is seen. The softening is caused by the decline of the cut-off energy $E_{\rm cut}(r_{\rm N})$ due to synchrotron cooling, though $E_{\rm max}$ is higher for a larger $\sigma$. Above $\sigma=10^{-4}$, the short advection time leads to the reduction of the IC flux. The drop of the cooling efficiency due to the short advection time causes the spectral hardening of the IC component. While the flux decrease of the synchrotron component due to the reduction of the advection time is mitigated by the magnetic field growth, the IC component more rapidly falls with $\sigma$ than the synchrotron one.

In the right panel of Figure 4, the X-ray surface brightness profile is shown. The size of the X-ray nebula contracts with increasing $\sigma$. The dependence of $\sigma^{-1/6}$ obtained in the previous subsection is a reasonable approximation.

4. Application to Observed Sources

Here, we apply our model to the observed sources, 3C 58 and G21.5-0.9, for which significant data sets are available to constrain the model parameters. Moreover, in both the two PWNe, the extent of the X-ray emission is close to the radio. We argue both the entire spectra and spatial profiles for those objects.

The images of 3C 58 were obtained in the radio wavelengths (e.g., Reynolds & Aller, 1988) and X-ray band (e.g Slane et al., 2004). The radial profile of photon index in the X-ray band was also obtained. The extents of the radio and X-ray images are similar as $\sim 5^{\prime}\times 9^{\prime}$. From the distance $D\simeq 2$ kpc (Kothes, 2013), we adopt $r_{\rm N}=2$ pc (Tanaka & Takahara, 2013b). The spin period and its time derivative for the central pulsar of 3C 58 (PSR J0205+6449) are $65.7\ {\rm ms}$ (Murray et al., 2002; Camilo et al., 2002) and $1.94\times 10^{-13}\ {\rm s\ s^{-1}}$ (Livingstone et al., 2009), respectively. The spin-down luminosity is estimated as $2.7\times 10^{37}\ {\rm erg\ s^{-1}}$, assuming $10^{45}\ {\rm g\ cm^{2}}$ for the moment inertia of the pulsar.

G21.5-0.9 shows spherical structures in the radio (Bietenholz & Bartel, 2008) and X-ray (Matheson & Safi-Harb, 2005, 2010; Camilo et al., 2006) images. The radio and X-ray sizes of the PWN ($\sim 40^{\prime\prime}$ in radius) are almost the same again. Adopting the distance $D=4.8\ {\rm kpc}$ (Tian & Leahy, 2008), we have $r_{\rm N}=0.9$ pc. PSR J1833-1034, the central object of G21.5-0.9, has a spin period $61.9$ ms (Gupta et al., 2005) and its derivative $2.02\times 10^{-13}\ {\rm s\ s^{-1}}$ (Roy et al., 2012), from which we obtain the spin-down luminosity $L_{\rm sd}=3.5\times 10^{37}\ {\rm erg\ s^{-1}}$.

The parameters to fit the spectra for the two PWNe are summarized in Table 3. See Figure 2 for the ISRFs taken from GALPROP v54.1 for the two cases. First, we discuss the parameter sets denoted with “broadband” in the Table 3 (hereafter we call them broadband model). The resultant radial profiles of $u(r)$ and $B(r)$ are shown in Figure 5. In both the two PWNe, $r_{\rm eq}$ is outside the nebula radius $r_{\rm N}$ in our parameter sets. The particle spectra in Figure 6 show the evolution of $E_{\rm cut}$ as discussed in Section 3.1. The volume averaged spectra (dashed lines) $\bar{n}(E)$ are well expressed by broken power-laws. The differences of the particle spectral index above the cooling break energy are 0.67 for both the cases. Given the particle energy $E$, Equation (29) implies the maximum radius those particles survive $r_{E}\propto E^{-1/5}$. The differences in the index seem consistent with a naive estimate $\bar{n}(E)\propto E^{-p_{2}}r_{E}^{3}$.

Figure 7 shows the volume-integrated photon spectra for the two PWNe. Our models roughly reproduce the entire structures of the spectra. In 3C 58, the data points obtained with Fermi (Abdo et al., 2013) may contain large systematic errors due to the emission from the central pulsar so that we treat those data as upper limits. The model spectrum in the X-ray range is apparently softer than the X-ray data in 3C 58 due to the cooling effect. As discussed in Section 3.3, $\nu_{\rm c}$ can be higher than the X-ray energy range by adopting a larger $\sigma$ or conversely lower $\sigma$. In such cases, the X-ray model spectrum may be as hard as observed one. However, when we adopt a lower $\sigma$ to make $\nu_{\rm c}$ above the X-ray frequency, the synchrotron component does not extend to the X-ray energy as Equation (33) indicates. We also do not find a consistent high-$\nu_{\rm c}$ model with a very large $\sigma$ or a slightly large $r_{\rm s}$, for which the radio and X-ray fluxes are hard to be reproduced simultaneously. One may suppose that a smaller $p_{2}$ can agree with the observed X-ray spectral index, even if $\nu_{\rm c}$ is below the X-ray frequency. The extrapolation from the X-ray data requires $\nu_{\rm c}<10^{14}$ Hz to make $\nu_{\rm b}$ above the radio data points. Such a low $\nu_{\rm c}$ is hard to be realized in this model (see the $\nu_{\rm c}$-turnover in Figure 4); the radio or X-ray flux becomes inconsistent for such extreme parameter sets. Therefore, our model spectra cannot be reconciled with the X-ray spectral index.

A similar problem to the case in 3C 58 arises in the X-ray spectrum of G21.5-0.9. When we fit a model with $p_{2}=2$, which yields a flat spectrum ($\nu L_{\nu}\propto\nu^{0}$) above $\nu_{\rm c}$, a much lower $\sigma$ is required to adjust the X-ray flux. For such a low $\sigma$, $\nu_{\rm cut}$ becomes lower than the X-ray band.

In those two objects, we are forced to have $\nu_{\rm c}$ below the X-ray band. As a result, the X-ray spectra show softer shapes than the observed ones. The X-ray extents are more compact than the radio images (Figure 8). The radial profiles of photon indices in 0.5-10.0 keV range (Figure 9) also deviate from the observed data. However, the discrepancy in 3C 58 is not so prominent compared to the model curve by Slane et al. (2004) based on Reynolds (2003). Note that the radial profiles of photon indices in optical (3944-4952$\AA$) and radio (4.75 GHz) band for two objects do not depend on angular distance from the pulsar. Since the $\nu_{\rm c}$ is higher than the frequencies of these bands in the parameters of the broadband models, the pairs can emit radio and optical photons without the cooling effect all over the nebula, thus the emission has same spectral index in each radial position.

The advection time of G21.5-0.9 well agrees with the age $870$ yr (Bietenholz & Bartel, 2008). However, Wang et al. (2006) argued that this object associates with BC48 guest star and its age is thus about $2000$ yr. In this case, $t_{\rm adv}$ in the broadband model becomes less than half of the age. In 3C 58, if this object associates with SN 1181 (Stephenson, 1971), $t_{\rm adv}$ is about 2 times larger than the age. The characteristic ages of these objects are $5370$ yr for 3C 58 and $4850$ yr for G21.5-0.9, respectively. Notify that the characteristic age tends to be longer than the actual pulsar age, especially for young pulsars. Meanwhile, the previous one-zone time dependent models have obtained the ages. For 3C 58, Bucciantini et al. (2011), Torres et al. (2013) and Tanaka & Takahara (2013b) obtained $t_{\rm age}\sim 2000$ yr, which is comparable with $t_{\rm adv}$ in our broadband model. However, Bucciantini et al. (2011) and Torres et al. (2014) adopted the different value of the distance to the object, and these three studies did not include the data of Aleksić et al. (2014). A direct comparison of the age with our advection times does not seem meaningful so much. On the other hand, for G21.5-0.9, $t_{\rm age}$ was estimated to be $870$ yr (Vorster et al., 2013; Torres et al., 2014) or 1000yr (Tanaka & Takahara, 2011) in one-zone models. Those are close to our estimate. These studies adopted the same condition for the distance and the observed flux, so that the coincidence of the age and $t_{\rm adv}$ encourages our 1-D model.

Next, we discuss the case, where the radio/IR/optical emission can be treated separately as an additional component. In the broken power-law spectrum for the pair injection, the low-energy portion dominates the number. In the case of the Crab Nebula, the required particle number to reconcile the radio flux is much larger than the theoretically expected value (Tanaka & Takahara, 2010, 2011). Atoyan & Aharonian (1996) and Olmi et al. (2015) treated the low-energy component as a different component from the wind particles in their calculations. Thus, as an alternative model, we assume that the low-energy particle component responsible for the radio/IR emission has a different origin from the high-energy component. In the alternative models, we incorporate only the high energy particles above $E_{\rm b}$ neglecting the emission below the optical band..

Blue dashed lines in Figure 7-9 are for the alternative models, whose model parameters are summarized in Table 3. We adopt a slightly large $r_{\rm s}$ and large $\sigma$, which lead to X-ray extents consistent with observation as shown in Figure 8. The difference of the X-ray profiles in the two PWNe is attributed to the effect of the adiabatic cooling in G21.5-0.9 as shown in Figure 6. Since we have adopted a larger $\sigma$ enough to establish $r_{\rm eq}<r_{\rm N}$ for G21.5-0.9, a signature of adiabatic cooling appears. Although those models seem to reproduce the observed X-ray surface brightness and X-$\gamma$ fluxes, the resultant advection times become very short (see Table 3).

5. Discussion

As shown in Section 4, we have fitted the entire spectra of 3C 58 and G21.5-0.9. The obtained $\sigma$ by fitting the entire spectrum of nebula in the broadband models is about 10 times smaller than the conceivable value in the Crab Nebula and obtained $r_{\rm s}$ is similar to the value of the Crab Nebula. Our 1-D model has difficulty to reproduce both of the hard spectra and the size of PWNe in X-rays. As discussed in Section 4, our model is hard to avoid the spectral softening due to the cooling effect in X-ray range. As a result, the X-ray nebula size becomes more compact than the observed extents. The 1-D model should be improved by introduction of possible physical processes, such as the spatial diffusion of high energy particles, reacceleration by turbulences, and amplification/dissipation of the magnetic field.

The 1-D model tends to lead a lower $\sigma$ than the values derived from one-zone models. One-zone time-dependent models resulted in $\sigma\sim 0.03$–$0.5$ (Bucciantini et al., 2011; Tanaka & Takahara, 2013b; Torres et al., 2013) and $0.01$–$0.2$ (Tanaka & Takahara, 2011; Vorster et al., 2013; Torres et al., 2014) for 3C 58 and G21.5-0.9, respectively, while the value in our paper is $\sigma\sim 10^{-4}$. Since the magnetic field increases with radius in the 1-D model, the average magnetic field $B_{\rm av}$ is consistent with the previous one-zone models (see also Equations (35) and (36)). On the other hand, $\nu_{\rm cut}$ is determined by the magnetic field near the shock rather than $B_{\rm av}$, in contrast to the cooling break $\nu_{\rm c}$. Therefore, the one-zone model overestimates the maximum synchrotron frequency about $B_{\rm av}/B_{\rm d}$ times higher than the 1-D model. Note that hard X-ray observations around $\nu_{\rm cut}$ are also interesting to investigate how the maximum energy of non-thermal pairs is determined. Although $E_{\rm max}$ is constrained by the size of $r_{\rm s}$ in our case as shown in Equation (17), $\nu_{\rm cut}\sim 100$ MeV in the Crab Nebula implies that $E_{\rm max}$ is determined by the balance of the acceleration and cooling times (e.g., de Jager et al., 1996).

The fitted shock radius $r_{\rm s}$ of 3C 58 is twice as large as that of G21.5-0.9. The total pressures at $r=r_{\rm N}$ in our models are $p_{\rm tot,3C~{}58}\sim 3.7\times 10^{-10}$ erg cm*-3* and $p_{\rm tot,G21.5}\sim 2.1\times 10^{-9}$ erg cm*-3*. From the pressure balance, the large $p_{\rm tot}$ implies the large plasma pressure of the surrounding remnants of their supernovae. Slane et al. (2004) obtained $kT_{e}\sim 0.23$ keV and $n_{\rm SNR}\sim 0.38$ cc*-1* for 3C 58 and then its pressure $\sim 1.4\times 10^{-10}$ erg cm*-3*. Matheson & Safi-Harb (2010) obtained $kT_{e}\sim 0.3$ keV and $n_{\rm SNR}\sim 0.63$ cm*-3* for G21.5-0.9 and then its pressure $\sim 3.0\times 10^{-10}$ erg cm*-3*. While the pressure values may be not so robust, this indicates that the surrounding SNR of G21.5-0.9 has higher pressure than 3C 58. In addition, the fact that the bright shell-like SNR is clearly seen in G21.5-0.9 also supports that the pressure for G21.5-0.9 would be higher than that for 3C 58.

In the broadband models, although we reproduce the flux levels of the entire spectrum, the X-ray spectral indices disagree with the observations. We argue the alternative models, in which the emission in radio and optical is assumed to be different from the direct emission from the pulsar wind (Atoyan & Aharonian, 1996; Olmi et al., 2014, 2015). As shown in Table 3, the obtained $\sigma$ in the alternative models tends to be larger than the value in the broadband models. This tendency is similar to some one-zone models (Torres et al., 2013; Vorster et al., 2013). The time-dependent model of Torres et al. (2013) introduced an order of magnitude larger energy density of the ISRF (i.e., the lager magnetic field strength) than ours in order to reproduce the X-ray spectral index of 3C 58. In the model of Vorster et al. (2013), the X-ray spectrum of G21.5-0.9 was also reproduced by a strong magnetic field ($230\mu$G), and hard spectral index ($p_{2}=2.0$). Note that their predicted GeV flux seems above the Fermi upper-limit (Ackermann et al., 2011). The cooling break was set well below keV range in those models, while it is hard to set $\nu_{\rm c}$ low enough in our model, because the high $\sigma$ leads to a shorter advection time than the cooling time scale of low-energy particles as discussed in Section 3.3. The temporal evolution of the magnetic field in one-zone models causes the gradual hardening of the particle spectrum (see also Tanaka & Takahara, 2010), which is favorable to fit the X-ray data differently from our steady model.

Since a larger $\sigma$ is required in the alternative models, the resultant short advection time prevents high-energy particles from cooling before reaching the edge of the nebula. However, such a short advection time may contradict the age of the PWNe. To validate $t_{\rm age}\gg t_{\rm adv}$, the efficient particle escape at the nebula surface should be required. Although the large amount of the escaped high-energy particles should emit photons outside the PWNe, such a signature outside PWNe has not been claimed. For example, the model of Holler et al. (2012), in which the radial velocity profile is artificially tuned, also implies $\sim 100$ yr for the advection time in G0.9+0.1, though the age is more than kyr. We should carefully note the advection time in modeling the outflow property (see Equation (24)). Even for the models including the effect of the spatial diffusion (Tang & Chevalier, 2012; Porth et al., 2016), a short diffusion time scale may be required to reproduce the X-ray surface brightness. This leads to the same problem in the short advection time case. If the outer supernova ejecta efficiently confines the PWN, the fast outflow implied in the high $\sigma$ model should be decelerated near the edge of the PWN, and should induce turbulence inside the PWN. As a result, the wind material may be efficiently mixed inside the PWN (e.g., Porth et al., 2014). In this case, the one-zone approximation may be rather adequate.

As shown in Figure 8, we calculated the radial radio profile in the broadband model. They appear to be almost uniform profile in radius, and seem to agree with the observational facts (e.g., Bietenholz & Bartel, 2008; Bietenholz et al., 2013). Additionally, the radiation in $\nu<\nu_{\rm c}$ ($\sim 10^{15}$ Hz) shows the similar profile as radio one. This is because the radial evolution of the density is common for the non-cooled particles. If a clear distinction of surface brightness between the radio and optical is detected, this would be the strong basis of the hypothesis that the radio emission of PWNe has a different origin from the optical and X-rays. We also calculate the radial profile of surface brightness in 0.8–2.0 TeV and 8–10 TeV. The extents of 3C 58 and G21.5-0.9 in 0.8–2 TeV are larger than X-rays in the broadband models, because the photon of $\sim 1$ TeV is emitted by non-thermal pairs with lower energies than the energy of particles emitting the synchrotron radiation at $\nu\sim\nu_{\rm cut}$ (c.f., de Jager & Djannati-Ataï, 2009). Since the extent of 3C 58 is about $200^{\prime\prime}$, CTA will be able to present the spatially resolved $\gamma$-ray map of 3C 58 in $\sim 1$ TeV.

6. Conclusion

We have revisited the 1-D steady model, and applied to the pulsar wind nebulae, in order to find a parameter set consistent with both the entire photon spectrum and surface brightness profile. It is still controversial whether the simple 1-D model reproduces observed properties of the PWNe other than the Crab Nebula or not. As we have shown in Section 3, both the entire photon spectrum and surface brightness profile largely depend on the parameters, the uncertain shock radius $r_{\rm s}$ and the magnetization parameter $\sigma$. The flux of inverse Compton component becomes dim with increasing $\sigma$. In contrast, the synchrotron component is not a monotonic function of $\sigma$. For the dependence on $r_{\rm s}$, while the synchrotron component becomes dim with increasing with $r_{\rm s}$, the IC component shows complicated behaviors. The X-ray size of a PWN becomes large with increasing $r_{\rm s}$ and decreasing $\sigma$.

We have fitted the entire spectrum of two observed sources 3C 58 and G21.5-0.9. Calculating the radial profile of the surface brightness for those models, we show that the resultant X-ray extents are significantly smaller than the observed sizes. Furthermore, we have performed another parameter set called “alternative” model, where we treat the radio and optical emissions as extra components. The alternative models successfully reproduce the observed X-ray surface brightness and the X-ray and $\gamma$-ray fluxes. However, those models imply too short advection time. In summary, the 1-D model constructed by KC84s has severe difficulty to reproduce both the spectrum and spatial emission profile of PWNe consistently. The model should be improved by taking some possible physical processes into consideration, such as spatial diffusion of non-thermal particles, reacceleration by turbulences.

Acknowledgments

First, we appreciate the anonymous referee for valuable comments. We are grateful to Kento Sasaki for providing the calculation code. This work is supported by Grants-in-Aid for Scientific Research Nos. 15K05069 (TT and KA), 25400227, 16K05291 (KA), and 2510447 (SJT) from the Ministry of Education, Culture, Sports, Science and Technology (MEXT) of Japan.

In order to check our model, we calculate the entire spectrum and the X-ray profile of the Crab Nebula and compare the result with Atoyan & Aharonian (1996) (hereafter AA96). We treat the SSC approximately by assuming that the synchrotron photon exist in the nebula homogeneously. The other differences of our model from AA96 are as follows. In our model, the low-energy particles responsible for the radio emission are supplied from the pulsar wind, while AA96 considered them as another component. In AA96, the maximum energy of pairs was treated as a parameter of the model (but see Equation (17) for our model). Finally, AA96 introduced the correction factor $\kappa$, which is a parameter to adjust the ratio of the synchrotron flux to the SSC flux. The parameter $\kappa$ (cf. AA96 adopted $\kappa\sim 0.5$) may represent the effects of deviation from the spherical symmetry or inhomogeneity inside the PWN.

In Figure 10, the entire photon spectrum and the X-ray radial profile for the Crab Nebula are shown. All the parameters to calculate the spectrum for the Crab Nebula are same as AA96 without “Obtained parameter”, and are summarized in Table 4. The flux of inverse Compton becomes a little smaller than the value observed and calculated by AA96. The difference in the SSC flux is due to the additional parameters $\kappa$ and $E_{\rm max}$ in AA96. Even if we take the smaller $\sigma$ to enhance the SSC flux, the maximum energy of synchrotron emission becomes much lower than the cut-off energy of observed spectrum (see Equation (33)). Within our conservative model assumption for the maximum energy of particles, we are not able to fit the spectrum around the cut-off of the spectrum.

The right panel of Figure 10 shows the radial profile of surface brightness in 3.0-5.0 keV range. The extent of the X-ray nebula calculated by our model is comparable with observations. Comparing with the cases of 3C 58 and G21.5-0.9, we do not have so strong motivation to improve the 1-D steady model in the case of the Crab Nebula.