A SELF-CONSISTENT AND SPATIALLY DEPENDENT MODEL OF THE MULTIBAND EMISSION OF PULSAR WIND NEBULAE

Pulsar wind nebulae (PWNe) are generally believed to be the very high energy γ-ray sources, which are created by the interaction of the pulsar wind with the nonrelativistic ejecta of the ambient supernova remnant (SNR; Goldreich & Julian 1969; Kennel & Coroniti 1984a). A termination shock is formed in the interaction and then the particles are accelerated to extreme relativistic energy at the shock (Rees & Gunn 1974; Kennel & Coroniti 1984b). From the termination shock, the high energy particles are injected into the nebula, and emit photons ranging from radio to TeV γ-rays during the evolution of the PWN. In recent years, the High Energy Stereoscopic System (H.E.S.S.) survey of the galactic plane has detected about 30 PWNe as the γ-ray sources (see a recent status report of Carrigan et al. 2013), which probably provided more information to study the radiation processes of PWNe.

To explain the multiband spectral energy distribution (SED) of a PWN, various models have been proposed. Usually, the radio and X-rays are attributed to the process of relativistic electron's synchrotron radiation, and the very high energy γ-rays are produced by the inverse Compton scattering with various soft photons (Atoyan & Aharonian 1996; Gaensler & Slane 2006). In addition, the relativistic particles coming directly from the pulsar magnetosphere are usually assumed to be uniform (e.g., Zhang et al. 2008; Fang & Zhang 2010; Tanaka & Takahara 2010; Martín et al. 2012; Zhu et al. 2015). In their models, the distribution of particles is spatially independent and the diffusion of particles is often approximated as an escape timescale or neglected. In this case, the multiband observational results of a PWN can be well reproduced by these models. However, the observations implied that the photon index in the X-ray band is a function of position (e.g., Slane et al. 2000; Bocchino et al. 2001, 2005; Mangano et al. 2005; Schöck et al. 2010; Holler et al. 2012). Moreover, Tang & Chevalier (2012) used a model with diffusion to study the spectral index profiles of Crab, 3C 58, and G21.5-0.9, and found that the advection-diffusion model can fit the spectral index profiles better than the pure advection model, which means that the diffusion probably plays an important role in changing the particle spectral index. Vorster & Moraal (2013) calculated the particle (e.g., electron) spectrum as a function of distance by solving the Fokker–Planck transport equation, and studied the effect of diffusion on the particle spectrum in PWN. In their calculations, the spatial distributions of magnetic field and convection velocity in PWNe are approximated as $B\propto 1/{r}^{(1-\alpha )}$ and $V\propto 1/{r}^{\alpha }$, respectively, which are qualitatively similar to the results calculated by Kennel & Coroniti (1984a), but will break down in the case of $\alpha \geqslant 1.0$. They expected that particle diffusion is important for reducing the amount of synchrotron losses, and then explained the variation in the photon index of the X-ray synchrotron radiation as a function of position. Recently, Porth et al. (2016) investigated the particle diffusion in PWNe using three-dimensional magnetohydrodynamic (MHD) simulations. In their model, an eddy particle diffusion is adopted, and the distributions of the magnetic field and the convection velocity in the PWN are described as the Gaussian functions, respectively, which are derived from fitting the averaged data predicted by the MHD model, and then these Gaussian functions can easily be incorporated into the kinetic equation. In the steady-state system, for simplicity, this model under the assumption of spherical symmetry (i.e., diffusion coefficient only has radial component) can well reproduce the spatial variation in the photon index of X-ray synchrotron emissions of G21.5-0.9, Vela, and 3C 58. However, in this model, an adiabatic heating effect will appear in the inner regions of the PWN due to the effect of the rapid decrease in velocity in these regions (Porth et al. 2016).

Motivated by the above discussions, we describe a self-consistent and spatially dependent model to study the multiband emission of PWNe. In this model, the magnetic field in the PWN is assumed to be toroidal and the spatial variations of magnetic field and convection velocity are described by the semi-analytic model of Kennel & Coroniti (1984a); the dynamical evolution of the PWN is described by using the model of Bucciantini et al. (2011); the temporal evolutions of electron and photon spectra are self-consistently solved in this system. In our treatment, the temporal spectral evolution of electrons is described as a Fokker–Planck equation in spherical coordinates, where the electrons are injected with a time-dependent broken power law (e.g., Atoyan & Aharonian 1996; Venter & de Jager 2006), and the processes of convection, diffusion, and adiabatic and synchrotron radiation losses (Vorster & Moraal 2013), as well as inverse Compton scattering loss and photon–photon pair production are included; meanwhile, the temporal spectral evolution of photons is described as a spatially dependent conservation equation. This model is applied to MSH 15–52 to explain the multiband observed data and the spatial variations of photon index and surface brightness in radio, X-ray, and γ-ray bands.

The paper is organized as follows. In Section 2, we describe the self-consistent spatial model and numerical method. In Section 3, we apply the model to the MSH 15–52. Finally, in Section 4, we give the conclusions and discussions.

A realistic PWN is a dynamical system that has significantly and complicatedly expanded (e.g., Chevalier & Fransson 1992; Blondin et al. 2001; Van der Swaluw et al. 2001; Gelfand et al. 2009; Bucciantini et al. 2011) and has a quite complicated magnetic field structure shown in MHD simulations (e.g., Del Zanna et al. 2004; Porth et al. 2014). In the self-consistent and spatially dependent model, for simplicity, under the assumption that the system is spherically symmetric and the magnetic field in the PWN is toroidal, the dynamical evolution of the PWN is simulated using the model of Bucciantini et al. (2011) and the spatial variations of magnetic field and convection velocity are calculated with the semi-analytic model of Kennel & Coroniti (1984a). The self-consistent and spherical symmetric model can only apply to the young PWN before the reverse shock crushes the nebula, because the structure of the nebula could be very complicated after the crushing phase and a bow shock nebula could possibly form due to asymmetric compression in this phase (e.g., Chevalier 1998; van der Swaluw et al. 1998, 2004). In such a system, a central pulsar has a spin-down luminosity L(t), which is given by (Pacini & Salvati 1973)

$$\begin{eqnarray}&&L(t)={L}_{0}{\left(1+\displaystyle \frac{t}{{\tau }_{0}}\right)}^{-\tfrac{n+1}{n-1}},\end{eqnarray} \tag{ 1 }$$

where L₀ is the initial spin-down power, n is the breaking index, and ${\tau }_{0}$ is the initial spin-down timescale of the pulsar, which is determined by the current pulsar period P, its time derivative $\dot{P}$, age of the pulsar ${T}_{\mathrm{age}}$, and pulsar moment of inertia I (e.g., Gaensler & Slane 2006; Slane 2008). It is believed that a major fraction of the spin-down luminosity of a central pulsar will be converted into the mechanical luminosity of a pulsar wind. Relativistic electrons of such a wind will propagate through a surrounding medium and be accelerated at the termination shock of the pulsar wind, and then injected into the PWN (Kennel & Coroniti 1984a; Aharonian et al. 1997). Meanwhile, non-thermal photons are produced through various mechanisms of the relativistic electrons. Below, we describe the evolutions of relativistic electrons and non-thermal photons in the system and our numerical method.

2.1. Basic Equations

On the one hand, the electron's evolution is generally described as the Fokker–Planck equation (e.g., Parker 1965), but in the spherically symmetric system, the equation can be simplified (e.g., Vorster & Moraal 2013). On the other hand, non-thermal photons' evolution can be described by the photon conservation equation (Mastichiadis & Kirk 1995). Assuming the number density of the electrons with Lorentz factor γ is ${n}_{{\rm{e}}}\equiv {n}_{{\rm{e}}}(r,\gamma ,t)$ (i.e., the number of particles in the energy interval ($\gamma ,\gamma +d\gamma$) at time t and radial distance r per unit volume), and the photon's number density is ${n}_{\gamma }\equiv {n}_{\gamma }(r,x,t)$ with the dimensionless photon energy $x=h\nu /({m}_{{\rm{e}}}{c}^{2})$, we can express them as

$$\begin{eqnarray}\displaystyle \frac{\partial {n}_{{\rm{e}}}}{\partial t} & = & \kappa \displaystyle \frac{{\partial }^{2}{n}_{{\rm{e}}}}{\partial {r}^{2}}+\left[\displaystyle \frac{1}{{r}^{2}}\displaystyle \frac{\partial }{\partial r}({r}^{2}\kappa )-V\right]\displaystyle \frac{\partial {n}_{{\rm{e}}}}{\partial r}-\displaystyle \frac{1}{{r}^{2}}\displaystyle \frac{\partial }{\partial r}[{r}^{2}V]{n}_{{\rm{e}}}\\ & & +\displaystyle \frac{\partial }{\partial \gamma }[\dot{\gamma }{n}_{{\rm{e}}}]+{Q}_{{\rm{e}}}^{\gamma \gamma }+{Q}_{\mathrm{inj}},\end{eqnarray} \tag{ 2 }$$

and

$$\begin{eqnarray}&&\displaystyle \frac{\partial {n}_{\gamma }}{\partial t}+\displaystyle \frac{{n}_{\gamma }}{{\tau }_{\mathrm{esc}}}={Q}_{\gamma }^{\mathrm{syn}}+{Q}_{\gamma }^{\mathrm{ics}}+{L}_{\gamma }^{\mathrm{ssa}}+{L}_{\gamma }^{\gamma \gamma },\end{eqnarray} \tag{ 3 }$$

respectively. The physical quantities in Equations (2) and (3) are described as follows.

In Equation (2), the momentum diffusion in spherical coordinates is neglected, $\dot{\gamma }$ is the sum of energy losses and is given by

$$\begin{eqnarray}&&\dot{\gamma }={\dot{\gamma }}_{\mathrm{ad}}+{\dot{\gamma }}_{\mathrm{syn}}+{\dot{\gamma }}_{\mathrm{ics}},\end{eqnarray} \tag{ 4 }$$

where ${\dot{\gamma }}_{\mathrm{ad}}$ is the adiabatic loss, which is

$$\begin{eqnarray}&&{\dot{\gamma }}_{\mathrm{ad}}(r,\gamma ,t)=\displaystyle \frac{1}{3}{\rm{\nabla }}\cdot {\boldsymbol{V}}\gamma =\displaystyle \frac{1}{3}\left[\displaystyle \frac{2V}{r}+\displaystyle \frac{\partial V}{\partial r}\right]\gamma ,\end{eqnarray} \tag{ 5 }$$

where V is the convection velocity; ${\dot{\gamma }}_{\mathrm{syn}}$ is the synchrotron loss and is given by (Rybicki & Lightman 1979)

$$\begin{eqnarray}&&{\dot{\gamma }}_{\mathrm{syn}}(r,\gamma ,t)=\displaystyle \frac{4}{3}\displaystyle \frac{{\sigma }_{{\rm{T}}}}{{m}_{{\rm{e}}}c}{\gamma }^{2}{U}_{{\rm{B}}}(r,t),\end{eqnarray} \tag{ 6 }$$

where ${\sigma }_{{\rm{T}}}$ is the Thomson cross-section and ${U}_{{\rm{B}}}(r,t)\,={B}^{2}(r,t)/8\pi$ is the magnetic field energy density; ${\dot{\gamma }}_{\mathrm{ics}}$ is the IC loss and is given by (Blumenthal & Gould 1970)

$$\begin{eqnarray}{\dot{\gamma }}_{\mathrm{ics}}(r,\gamma ,t) & = & \displaystyle \frac{3}{4}{\sigma }_{{\rm{T}}}c\displaystyle \frac{1}{{\gamma }^{2}}{\displaystyle \int }_{0}^{\infty }{x}_{f}{{dx}}_{f}{\displaystyle \int }_{0}^{\infty }\displaystyle \frac{{n}_{\gamma }(r,{x}_{i},t)}{{x}_{i}}\\ & & \times f(q,{\rm{\Gamma }})H(1-q)\left(q-\displaystyle \frac{1}{4{\gamma }^{2}}\right){{dx}}_{i},\end{eqnarray} \tag{ 7 }$$

where x_i and x_f are the initial and final energies of the scattered photons, respectively, H is the Heaviside step function and ${n}_{\gamma }(r,{x}_{i},t)$ is the distribution of target photon fields, ${\rm{\Gamma }}=4\gamma {x}_{i}$, $q={x}_{f}/({\rm{\Gamma }}(\gamma -{x}_{f}))$, and $f(q,{\rm{\Gamma }})=2q\mathrm{ln}q+(1+2q)(1-q)\,+$ $0.5(1-q){({\rm{\Gamma }}q)}^{2}/(1+{\rm{\Gamma }}q)$.

In the above expressions of Equations (4)–(7), the electron's energy losses depend on the PWN evolution. As mentioned above, we calculate the dynamical evolution of the PWN using the model of Bucciantini et al. (2011). In their model, the dynamical evolution of a PWN is divided into three phases: initial free-expansion, reverberation, and Sedov–Taylor phases. The object is assumed to be a spherical symmetry and the proper motion of the pulsar is neglected and the radiation losses of the particles are negligible. For the young objects, such as Crab, 3C 58, and MSH 15–52, dynamic modeling indicated that they are expanding in the free-expansion phase (Bucciantini et al. 2011). In this case, combined mass and momentum conservation with the energy conservation, an analytic solution with a series expansion is derived to calculate the radius of the PWN in the free-expansion phase

$$\begin{eqnarray}{R}_{{\rm{w}}}(t) & = & {\left(\displaystyle \frac{{E}_{\mathrm{sn}}^{3}{L}_{0}^{2}}{{M}_{\mathrm{ej}}^{5}}\right)}^{1/10}{\left(\displaystyle \frac{t}{1+t/{\tau }_{0}}\right)}^{6/5}\\ & & \times \displaystyle \frac{1}{1-s}\displaystyle \sum _{i=0}^{\infty }{c}_{i}{s}^{i},\end{eqnarray} \tag{ 8 }$$

where ${E}_{\mathrm{sn}}$ is the energy of the supernova explosion, ${M}_{\mathrm{ej}}$ is the mass of supernova ejecta, $s=(t/{\tau }_{0})/(1+t/{\tau }_{0})$, and the expressions of the coefficients c_i are shown in Bucciantini et al. (2004). Moreover, the radius of the termination shock can be estimated as

$$\begin{eqnarray}&&{r}_{0}=\sqrt{\displaystyle \frac{L(t)}{4\pi {{cP}}_{\mathrm{ts}}(t)}},\end{eqnarray} \tag{ 9 }$$

where ${P}_{\mathrm{ts}}(t)$ is the pressure of the gas at the termination shock (see Equations (7) and (18) of Bucciantini et al. 2004).

Assuming that V₀ is the convection velocity at r₀, the convection velocity at any radius r, which depends on the magnetization parameter σ and decreases with the increase of radial radius, can be written as

$$\begin{eqnarray}&&V(r)={V}_{0}v(r),\end{eqnarray} \tag{ 10 }$$

where v(r) is the radial profile of the velocity and it is a dimensionless quantity (see Equations (4.11) and (5.12c) of Kennel & Coroniti 1984a). Following Kennel & Coroniti (1984a), the convection velocity at the termination shock is given by

$$\begin{eqnarray}{V}_{0}^{2} & = & \displaystyle \frac{8{\sigma }^{2}+10\sigma +1}{16\sigma +1}+\displaystyle \frac{1}{16\sigma +1}\\ & & \times {[64{\sigma }^{2}{(\sigma +1)}^{2}+20\sigma (\sigma +1)+1]}^{\tfrac{1}{2}}.\end{eqnarray} \tag{ 11 }$$

Based on the MHD restrictions shown in Kennel & Coroniti (1984a), i.e., ${VBr}={V}_{0}{B}_{0}{r}_{0}=\mathrm{constant}$, where B₀ is the magnetic field strength at r₀ and is determined by the spin-down power L(t), so the temporal evolution of the radial magnetic field can be expressed as

$$\begin{eqnarray}&&B(r,t)={B}_{0}(t)\displaystyle \frac{{r}_{0}}{{rv}(r)},\end{eqnarray} \tag{ 12 }$$

where ${B}_{0}(t)$ is calculated by ${W}_{{\rm{B}}}(t)={\int }_{{r}_{0}}^{{R}_{{\rm{w}}}(t)}{B}^{2}(r,t){r}^{2}{dr}/2$ where ${W}_{{\rm{B}}}(t)$ is the magnetic energy. In an expanding system, the magnetic energy is balanced by the injection rate of electromagnetic energy and the adiabatic loss. Thus, the magnetic energy can be obtained by the equation (Torres et al. 2013)

$$\begin{eqnarray}&&\displaystyle \frac{{{dW}}_{{\rm{B}}}(t)}{{dt}}=\eta L(t)-\displaystyle \frac{{W}_{{\rm{B}}}(t)}{{R}_{{\rm{w}}}(t)}\displaystyle \frac{{{dR}}_{{\rm{w}}}(t)}{{dt}},\end{eqnarray} \tag{ 13 }$$

where $\eta \equiv {\dot{E}}_{{\rm{B}}}(t)/L(t)$ is the magnetic energy fraction and ${\dot{E}}_{{\rm{B}}}(t)$ is the rate of electromagnetic energy (e.g., Gelfand et al. 2009). Thus, the magnetization parameter can be expressed as $\sigma =\eta /(1-\eta )$.

Once the radial magnetic field is determined, we can estimate the diffusion coefficient κ, which can be written as $\kappa ={\kappa }_{0}{\kappa }_{{\rm{r}}}{\kappa }_{\gamma }$, where ${\kappa }_{0}$ is the diffusion coefficient at r₀ (Lerche & Schlickeiser 1981; Caballero-Lopez et al. 2004). Following Caballero-Lopez et al. (2004), the radial dependence of diffusion ${\kappa }_{{\rm{r}}}$ can be modeled as ${\kappa }_{{\rm{r}}}\propto 1/B(r,t)$ and the energy dependence of diffusion ${\kappa }_{\gamma }$ is given by ${\kappa }_{\gamma }\propto \gamma$. Therefore, the diffusion coefficient is given by

$$\begin{eqnarray}&&\kappa (r,\gamma ,t)={\kappa }_{0}(t)\left(\displaystyle \frac{\gamma {m}_{{\rm{e}}}{c}^{2}}{1\,\mathrm{TeV}}\right)\left[\displaystyle \frac{{rv}(r)}{{r}_{0}}\right].\end{eqnarray} \tag{ 14 }$$

In Equation (2), there are two injection terms of the electrons, one represents the electrons derived from the central pulsar, which are accelerated at the termination shock of the pulsar wind and then injected into the PWN (Kennel & Coroniti 1984a; Aharonian et al. 1997), the other represents the injection of electrons produced in the photon–photon pair production process. For the first injection term, the multiband observations of photon emissions from the PWN indicated that the electron spectrum can be divided into two components (e.g., Weiler & Panagia 1978; Gaensler & Slane 2006). Therefore, the particle evolution models often use a broken power law to describe the spectrum of the injected relativistic electrons (e.g., Venter & de Jager 2006; Zhang et al. 2008; Tanaka & Takahara 2010; Martín et al. 2012; Zhu et al. 2015), and then the injection rate is assumed to be

$$\begin{eqnarray}{Q}_{\mathrm{inj}}(\gamma ,t)={Q}_{0}(t)\,\left\{\begin{array}{ll}{(\gamma /{\gamma }_{{\rm{b}}})}^{-{\alpha }_{1}} & \mathrm{for}\,\,\gamma \lt {\gamma }_{{\rm{b}}}\\ {(\gamma /{\gamma }_{{\rm{b}}})}^{-{\alpha }_{2}} & \mathrm{for}\,\,\gamma \geqslant {\gamma }_{{\rm{b}}}\end{array}\right.,\end{eqnarray} \tag{ 15 }$$

where ${\gamma }_{{\rm{b}}}$ is the break energy and the parameters ${\alpha }_{1}$ and ${\alpha }_{2}$ are the spectral indices. ${Q}_{0}(t)$ is the normalization coefficient and is estimated by

$$\begin{eqnarray}&&(1-\eta )L(t)={\int }_{0}^{\infty }\gamma {m}_{{\rm{e}}}{c}^{2}{Q}_{\mathrm{inj}}(\gamma ,t)d\gamma .\end{eqnarray} \tag{ 16 }$$

The second injection term is given by (Coppi & Blandford 1990)

$$\begin{eqnarray}{Q}_{{\rm{e}}}^{\gamma \gamma }(r,\gamma ,t) & = & 4{\sigma }_{{\rm{T}}}{{cn}}_{\gamma }(r,2\gamma ,t){\displaystyle \int }_{0}^{\infty }{n}_{\gamma }(r,x^{\prime} ,t)\\ & & \times {R}_{\gamma \gamma }(2\gamma x^{\prime} ){dx}^{\prime} ,\end{eqnarray} \tag{ 17 }$$

where

$$\begin{eqnarray}&&{R}_{\gamma \gamma }(\omega )=0.652\displaystyle \frac{{\omega }^{2}-1}{{\omega }^{3}}\mathrm{ln}\omega H(\omega -1),\end{eqnarray} \tag{ 18 }$$

with ${\boldsymbol{\omega }}=2{\boldsymbol{\gamma }}{\boldsymbol{x}}^{\prime}$ and H is the usual Heaviside function.

We now turn to Equation (3), ${\tau }_{\mathrm{esc}}={R}_{{\rm{w}}}(t)/c$ is the photon escape timescale. In the right side of Equation (3), ${Q}_{\gamma }^{\mathrm{syn}}$ is the emissivity of the synchrotron radiation and is given by (Blumenthal & Gould 1970)

$$\begin{eqnarray}&&{Q}_{\gamma }^{\mathrm{syn}}(r,x,t)={\int }_{0}^{\infty }{n}_{{\rm{e}}}(r,\gamma ,t){P}_{\mathrm{syn}}(r,x,\gamma ,t)d\gamma ,\end{eqnarray} \tag{ 19 }$$

where ${P}_{\mathrm{syn}}(r,x,\gamma ,t)$ is the single particle synchrotron emissivity

$$\begin{eqnarray}&&{P}_{\mathrm{syn}}(r,x,\gamma )=\displaystyle \frac{\sqrt{3}{e}^{3}B(r,t)}{{{hxm}}_{{\rm{e}}}{c}^{2}}F\left(\displaystyle \frac{x}{{x}_{{\rm{c}}}}\right),\end{eqnarray} \tag{ 20 }$$

where ${x}_{{\rm{c}}}=3{heB}(r,t){\gamma }^{2}/(4\pi {m}_{{\rm{e}}}^{2}{c}^{3})$, $F(y)=y{\int }_{y}^{\infty }{K}_{5/3}({\rm{z}})d{\rm{z}}$ with $y=x/{x}_{{\rm{c}}}$, and ${K}_{5/3}$ is a modified Bessel function of the order of 5/3. ${Q}_{\gamma }^{\mathrm{ics}}$ is the emissivity of IC scattering in the Klein–Nishina regime and is given by (Blumenthal & Gould 1970)

$$\begin{eqnarray}&&{Q}_{\gamma }^{\mathrm{ics}}(r,x,t)={\int }_{0}^{\infty }{n}_{{\rm{e}}}(r,\gamma ,t){P}_{\mathrm{IC}}(r,\gamma ,x,t)d\gamma ,\end{eqnarray} \tag{ 21 }$$

where ${P}_{\mathrm{IC}}(r,\gamma ,x,t)$ is the scattered photon emissivity per electron

$$\begin{eqnarray}{P}_{\mathrm{IC}}(r,\gamma ,x,t) & = & \displaystyle \frac{3}{4}\displaystyle \frac{{\sigma }_{T}c}{{\gamma }^{2}}{\displaystyle \int }_{0}^{\infty }\displaystyle \frac{{n}_{\gamma }(r,{x}_{i},t)}{{x}_{i}}\\ & & \times f(q,{\rm{\Gamma }}){{dx}}_{i},\end{eqnarray} \tag{ 22 }$$

where x_i is the energy of target photons, and x is the final energy of photons, the soft photons in IC scattering include the cosmic microwave background (CMB), the galactic far-infrared background (FIR), the near-infrared and optical photon field due to the stars (NIR), and synchrotron photons. The exact value of each photon field will be defined in Section 3. ${L}_{\gamma }^{\mathrm{ssa}}$ is the loss term for the photons due to synchrotron self-absorption and can be written as (Rybicki & Lightman 1979; Mastichiadis & Kirk 1995)

$$\begin{eqnarray}{L}_{\gamma }^{\mathrm{ssa}}(r,x,t) & = & {\sigma }_{{\rm{T}}}c\displaystyle \frac{\pi }{6{\alpha }_{f}}{x}^{-1/2}{b}^{-1/2}{n}_{\gamma }(r,x,t)\\ & & \times {\left[\displaystyle \frac{\partial }{\partial \gamma }\left(\displaystyle \frac{{n}_{{\rm{e}}}}{{\gamma }^{2}}\right)\right]}_{\gamma ={(x/b)}^{1/2}},\end{eqnarray} \tag{ 23 }$$

where ${\alpha }_{f}$ is the fine-structure constant, and $b=B(r,t)/{B}_{\mathrm{cr}}$ with ${B}_{\mathrm{cr}}={m}_{{\rm{e}}}^{2}{c}^{3}/(e{\hslash })=4.414\times {10}^{13}{\rm{G}}$ is the critical magnetic field. Finally, ${L}_{\gamma }^{\gamma \gamma }$ is another loss term of the photons due to photon–photon pair production, and is expressed as (Coppi & Blandford 1990)

$$\begin{eqnarray}{L}_{\gamma }^{\gamma \gamma }(r,x,t) & = & -\,{\sigma }_{{\rm{T}}}{{cn}}_{\gamma }(r,x,t){\displaystyle \int }_{0}^{\infty }{n}_{\gamma }(r,x^{\prime} ,t)\\ & & \times {R}_{\gamma \gamma }({xx}^{\prime} ){dx}^{\prime} .\end{eqnarray} \tag{ 24 }$$

2.2. The Numerical Method

To solve the particle spectrum evolution in Equation (2), an Alternating Direction Implicit method (Douglas 1962) is used and the boundary conditions of both energy and radius must be determined. In the energy direction, following Vorster & Moraal (2013), the free-escape boundary conditions are imposed, i.e., ${n}_{{\rm{e}}}(r,{\gamma }_{\min },t)=0$ and ${n}_{{\rm{e}}}(r,{\gamma }_{\max },t)=0$. The Larmor radius ${r}_{{\rm{L}}}$ of the particles must be smaller than the termination shock radius, which leads to the maximum energy of electron satisfying

$$\begin{eqnarray}&&{\gamma }_{\max }(t)=\displaystyle \frac{\zeta \varepsilon e}{{m}_{e}{c}^{2}}\sqrt{\displaystyle \frac{\sigma L(t)}{(1+\sigma )c}},\end{eqnarray} \tag{ 25 }$$

where ζ is the compression ratio at the shock, and for strong shocks ($\sigma \ll 1$), $\zeta =3.0$ (Venter & de Jager 2006). The parameter ε is the shock radius fraction and is defined as $\varepsilon ={r}_{{\rm{L}}}/{r}_{0}$. In the radial direction, the particles in the outer boundary should escape from the system. Thus, the outer boundary condition is given by ${n}_{{\rm{e}}}({R}_{{\rm{w}}},\gamma ,t)=0$. In the inner boundary, which is located at the termination shock, the number of the particles should be equal to the injected particles. Thus, the inner boundary condition can be written as (e.g., Ng & Gleeson 1975)

$$\begin{eqnarray}&&{V}_{0}{n}_{{\rm{e}}}-\kappa ({r}_{0},\gamma ,t)\displaystyle \frac{\partial {n}_{{\rm{e}}}}{\partial r}=\displaystyle \frac{{Q}_{\mathrm{inj}}(\gamma ,t)}{4\pi {r}_{0}^{2}}.\end{eqnarray} \tag{ 26 }$$

Equation (3) is numerically solved using the Crank–Nicolson method (Crank & Nicotson 1996) and the free-escape boundary conditions are also imposed, i.e., ${n}_{\gamma }(r,{x}_{\min },t)=0$ and ${n}_{\gamma }(r,{x}_{\max },t)=0$. The spatial average radiation luminosity at different positions can be calculated by

$$\begin{eqnarray}&&L(r,{E}_{\gamma },t)={x}^{2}{m}_{{\rm{e}}}{c}^{2}{\int }_{{r}_{0}}^{r}4\pi {r}_{1}^{2}{n}_{\gamma }({r}_{1},x,t)\displaystyle \frac{{{dr}}_{1}}{\tau {^{\prime} }_{\mathrm{esc}}},\end{eqnarray} \tag{ 27 }$$

where $\tau {^{\prime} }_{\mathrm{esc}}=(r-{r}_{0})/c$ and ${E}_{\gamma }$ is the photon energy. The surface brightness at different positions are calculated according to Holler et al. (2012).

In this section, we apply the model to MSH 15–52 (G320.4-1.2 or HESS 1514-591). This PWN is a composite SNR that is associated with the energetic pulsar PSR B1509-58 (or PSR J1513-5908). The pulsar has a 150 ms rotation period (Seward & Harden 1982; Trussoni et al. 1996). As one of the the youngest and most energetic known in the galaxy, it has a period derivative of $1.5\times {10}^{-12}\,{\rm{s}}\,{{\rm{s}}}^{-1}$, a characteristic age of ∼1600 years, and a measured braking index $n=2.839\pm 0.003$ (Livingstone et al. 2005). Based on the HI absorption measurement, the value of the distance to system is estimated to be 5.2 ± 1.4 kpc (Gaensler et al. 1999), and the age of MSH 15–52 is likely lower than ∼1700 years from the observation of HI emission (Dubner et al. 2002).

In the radio band, Gaensler et al. (2002) estimated that the flux density at both 0.8 and 1.4 GHz is 2 ± 1 Jy based on the data from Whiteoak & Green (1996) and Gaensler et al. (1999). In the X-ray band, the observations of the Chandra satellite have revealed that the PWN has a power-law photon index of 2.05 ± 0.04 in the energy band of $0.5\mbox{--}10\,\mathrm{keV}$ (Gaensler et al. 2002). At hard X-ray energies, the PWN has been observed with BeppoSAX and INTEGRAL-IBIS telescopes. The power-law fit of the emission from the PWN gave a photon index of 2.1 ± 0.2 and a flux of $1.2\times {10}^{-10}\,\mathrm{erg}\,{{\rm{s}}}^{-1}\,{\mathrm{cm}}^{-2}$ in the 20–200 keV energy range (Mineo et al. 2001; Forot et al. 2006). In the GeV γ-ray band, the PWN was detected by the Fermi-LAT, and the LAT spectrum can be well described by a power law with a photon index of $1.57\pm {0.17}_{\mathrm{stat}}\pm {0.13}_{\mathrm{syst}}$ in the 1–100 GeV energy band (Abdo et al. 2010). In the very high energy γ-ray band, the PWN was measured by H.E.S.S., and the energy spectrum can be well fitted by a power law with a spectral index of $2.27\pm {0.03}_{\mathrm{stat}}\pm {0.02}_{\mathrm{syst}}$ in the energy range of 280 GeV to 40 TeV (Aharonian et al. 2005).

The observations revealed that MSH 15–52 has two outflow jets in the southeast and northwest directions and has a large thermal shell structure (RCW 89) in the north of the PWN, especially bright and extended in X-rays (Gaensler et al. 2002; Yatsu et al. 2005; Forot et al. 2006; Yatsu et al. 2009; An et al. 2014). At very high energies, the source also shows a similar morphology as that shown in X-rays (Aharonian et al. 2005; Nakamori et al. 2008). These observations show that the geometric structure of MSH 15–52 is very complicated and deviates from the spherical symmetry. As the first-order approximation, we apply the self-consistent and spherical symmetry model to MSH 15–52, which can help us to study the effects of radial variation of azimuthally integrated spectra in the PWN.

In our calculation, we used the parameters show in Table 1. The initial sin-down timescale ${\tau }_{0}=290\,\mathrm{years}$ is determined by the observations (e.g., Gaensler & Slane 2006). With an age of ${T}_{\mathrm{age}}=1450\,\mathrm{years}$ and an initial spin-down power of ${L}_{0}=3.9\times {10}^{39}\,\mathrm{erg}\,{{\rm{s}}}^{-1}$, the current spin-down luminosity of the pulsar is $9.2\times {10}^{37}\,\mathrm{erg}\,{{\rm{s}}}^{-1}$, which is consistent with the value derived from Livingstone et al. (2005) due to the uncertainty of the moment of inertia. The SN explosion energy of ${E}_{\mathrm{sn}}=1.5\times {10}^{52}\,\mathrm{erg}$ is higher than that shown in Torres et al. (2014) and Fang & Zhang (2010) because the radius of the PWN in our model is larger. However, our modeling result ${R}_{{\rm{w}}}({T}_{\mathrm{age}})=312\,\mathrm{arcsec}$ is consistent with the observed results of $300\,\mathrm{arcsec}$ (e.g., Gaensler et al. 2002; Schöck et al. 2010; An et al. 2014). Moreover, the SN explosion energy is consistent with the value of $4\mbox{--}20\times {10}^{51}\,\mathrm{erg}$ adopted by Bucciantini et al. (2011). The soft photon fields in the calculation are the same as that adopted in Torres et al. (2014), i.e., the $2.7\,{\rm{K}}$ CMB has an energy density of $0.25\,\mathrm{eV}\,{\mathrm{cm}}^{-3}$, the $20\,{\rm{K}}$ FIR has a energy density of $4\,\mathrm{eV}\,{\mathrm{cm}}^{-3}$, and the $400\,{\rm{K}}$ NIR has a energy density of $20\,\mathrm{eV}\,{\mathrm{cm}}^{-3}$.

Table 1.  Values of Parameters for MSH 15–52

Input Parameters	Symbol	Value
Ejected mass (${M}_{\odot }$)	${M}_{\mathrm{ej}}$	4.0
SN explosion energy (${10}^{51}\,\mathrm{erg}$)	${E}_{\mathrm{SN}}$	15.0
Period ($\mathrm{ms}$)	P	150
Period derivative (${\rm{s}}\,{{\rm{s}}}^{-1}$)	$\dot{P}$	$1.5\times {10}^{-12}$
Initial spin-down power ($\mathrm{erg}\,{{\rm{s}}}^{-1}$)	L0	$3.9\times {10}^{39}$
Initial spin-down timescale (years)	${\tau }_{0}$	290
Braking index	n	2.839
Age (years)	${T}_{\mathrm{age}}$	1450
Distance (kpc)	d	5.2
Shock radius fraction	ε	0.4
Fitted Parameters
Magnetic fraction	η	0.04
Low energy power-law index	${\alpha }_{1}$	1.65
High energy power-law index	${\alpha }_{2}$	2.45
Break energy	${\gamma }_{{\rm{b}}}$	$4.0\times {10}^{5}$
Initial diffusion coefficient(${\mathrm{cm}}^{2}\,{{\rm{s}}}^{-1}$)	${\kappa }_{0}$	$1.7\times {10}^{20}$

Download table as:  ASCIITypeset image

We now study the variations of both electron and photon spectra with radial distance. We calculate the present-day electron spectra and the corresponding non-thermal photon spectra at different positions, and the results are shown in Figure 1. The results indicate that both electron and photon spectra of the PWN are all a function of distance and decrease with increasing the distance. The spectrum shift to lower energies is attributed to adiabatic losses and the spatial variation of the spectral shape is dominated by the magnetic field and diffusion coefficient (e.g., Vorster & Moraal 2013; Lu et al. 2016).

Zoom In Zoom Out Reset image size

As mentioned in Section 2, non-thermal photons are produced by various radiation mechanisms. So it is worth investigating the different photon components with the distance. As an example, Figure 2 shows the spectra of non-thermal photons at different positions of $r=28.1\,\mathrm{arcsec}$ and $r=279.9\,\mathrm{arcsec}$, respectively. It is obvious that the radio and X-rays emission from the PWN are produced by the synchrotron radiation and γ-rays are mainly produced via the inverse Compton scattering with the FIR and NIR photons. Moreover, FIR and NIR are the major contributors to the high energy spectra of MSH 15–52 at all of the positions, while the contribution of self-synchrotron Compton (SSC) becomes less and less important as the distance increases.

Zoom In Zoom Out Reset image size

In order to study the effect of photon–photon pair production on the SED and the spatial photon index distribution of MSH 15–52, we calculate the present-day timescales of the photon escape and the photon–photon pair production changing with the photon energy, respectively, and the results are shown in Figure 3, where the photon escape timescale is calculated by ${\tau }_{\mathrm{esc}}(t)={R}_{{\rm{w}}}(t)/c$ and the photon–photon pair production timescale is given by

$$\begin{eqnarray}&&{\tau }_{\gamma \gamma }(r,x,t)={\left[{\sigma }_{{\rm{T}}}c{\int }_{0}^{\infty }{n}_{\gamma }(r,x^{\prime} ,t){R}_{\gamma \gamma }({xx}^{\prime} ){dx}^{\prime} \right]}^{-1}.\end{eqnarray} \tag{ 28 }$$

Compared to the photon escape, the process of photon–photon pair production is negligible because of the timescale of ${\tau }_{\mathrm{esc}}\ll {\tau }_{\gamma \gamma }$ in all of the positions. Therefore, the process of photon–photon pair production has little effect on the SED and the spatial photon index distribution of MSH 15–52. In fact, the photon–photon pair production process can affect the photon spectrum only when some specific conditions are satisfied (more details are shown in Petropoulou & Mastichiadis 2011).

Zoom In Zoom Out Reset image size

We also calculate the distributions of the surface brightness at different positions of the PWN in MSH 15–52 and the results are shown in Figure 4, where the different color lines correspond to the surface brightness at different positions. It is clear that the surface brightness of the PWN is also a function of distance. Especially in the X-rays, the photon index strongly depends on distance because the electron spectral index in the high energy varies with the distance (see Figure 1).

Zoom In Zoom Out Reset image size

We use our model to fit the observed SED of MSH 15–52 and the result is shown in Figure 5. The ${\chi }^{2}$ value of the model results with the observed data is 2.86. Note that an upper limit arrow presented around 0.1 MeV, which is obtained from Forot et al. (2006) based on the INTEGRAL-IBIS telescopes, implies a sharp declining in the energy range. Forot et al. (2006) pointed out that the possible spectral break near 160 keV can be explained by a simple jet scenario with efficient synchrotron cooling, and it is beyond the discussion in our paper. From the model, a present-day magnetic field ${B}_{0}({T}_{\mathrm{age}})=26.6\,\mu {\rm{G}}$ at the termination shock is derived. The value of the magnetic fraction $\eta =0.04$, i.e., the magnetization parameter $\sigma =0.042$, implies that the PWN in MSH 15–52 has a strong shock. This also means that the convection velocity has a value of ${V}_{0}=0.41c$ at the shock. For a strong shock with $\sigma =0.042$, a present-day spatial average magnetic field of $\bar{B}({T}_{\mathrm{age}})\,=14.9\,\mu {\rm{G}}$ is derived. This result is slightly smaller than the value of $17\,\mu {\rm{G}}$ obtained by Aharonian et al. (2005). With an initial diffusion coefficient of ${\kappa }_{0}=1.7\times {10}^{20}\,{\mathrm{cm}}^{2}\,{{\rm{s}}}^{-1}$, the present-day diffusion coefficient at the shock is estimated to be ${\kappa }_{0}({T}_{\mathrm{age}})=6.6\times {10}^{24}\,{\mathrm{cm}}^{2}\,{{\rm{s}}}^{-1}$ and a spatial average diffusion coefficient $\bar{\kappa }({T}_{\mathrm{age}})=1.4\times {10}^{25}\,{\mathrm{cm}}^{2}\,{{\rm{s}}}^{-1}$ of the PWN is obtained, which is smaller than the value that was estimated by An et al. (2014) using the advection-diffusion model shown in Tang & Chevalier (2012), i.e., the diffusion coefficient is estimated to be $5.6\mbox{--}8.2\times {10}^{25}\,{\mathrm{cm}}^{2}\,{{\rm{s}}}^{-1}$ for the magnetic field strength $B=8\mbox{--}17\,\mu {\rm{G}}$. It should be pointed out that the fitting around ${10}^{5}\,\mathrm{MeV}$ is not very good, which may be caused by using the same soft photon fields as those of Torres et al. (2014) in this paper. We believe that choosing different densities and temperatures of the soft photon fields can get a good fitting.

Zoom In Zoom Out Reset image size

Using the parameters obtained by fitting the observed SED, we calculate the variations of the spectral index and surface brightness with the distance for the PWN. Schöck et al. (2010) integrated the spectrum over the full azimuthal angle and analyzed the changes of the spectral index and surface brightness in the energy range of $0.5\mbox{--}9.0\,\mathrm{keV}$ by using the XMM-Newton data under the assumption of spherical symmetry. Their results clearly showed that the photon index increasing with the increase of distance while the surface brightness decreases as the distance increases in the X-ray band. An et al. (2014) studied the radial variation of the photon index and surface brightness in the energy range of $3\mbox{--}10\,\mathrm{keV}$ based on the NuSTAR data. Different from Schöck et al. (2010), they analyzed the northern nebula and the jet separately because the spectral index has a significant variation in the azimuthal direction. Their results also showed that the photon index increases but the surface brightness decreases with the distance in X-rays in the northern nebula. The comparisons of our results with the observed data in the energy range of $0.5\mbox{--}9.0\,\mathrm{keV}$ given by Schöck et al. (2010) and in the energy range of $3\mbox{--}10\,\mathrm{keV}$ obtained from An et al. (2014) are shown in Figure 6. It is clear that the radial variations of the photon index and surface brightness in X-rays detected by the XMM-Newton and the radial variation of the photon index in X-rays in the northern nebula observed by NuSTAR can be reproduced in the frame of our model. However, the radial variation of the surface brightness in X-rays in the northern nebula detected by the NuSTAR cannot be reproduced well in our model. A possible reason may be that the radial surface brightness data obtained from An et al. (2014) only include the annular regions in the northern nebula and the southern part is neglected, while we just consider the case that the spectrum is integrated over the full azimuthal angle.

Zoom In Zoom Out Reset image size

Finally, we predict the changes of the spectral index and surface brightness with the distance in the radio, X-ray, and γ-ray bands. Our calculating results are shown in Figure 7, where the above parameters are used. It can be seen that the spectral indexes depend on the distance in almost all of the energy ranges except in the radio band. Especially, in the energy ranges of $0.5\mbox{--}10\,\mathrm{keV}$ and $20\mbox{--}200\,\mathrm{keV}$, the spectral indexes have an obvious increase with the increase of distance. While in the energy of $1\mbox{--}100\,\mathrm{GeV}$, the spectral index only slightly varies as the distance increases. In the energy range of $0.2\mbox{--}40\,\mathrm{TeV}$, the variation of spectral index also is not very obvious, except in the region close to the outer boundary of the PWN. Although the spectral indexes in these energy ranges all increase with the increase of distance, the changing characteristics are not identical (see Figure 7). Note that the variation of the spectral index in the TeV γ-rays is not consistent with that predicted by Schöck et al. (2010), because their radially symmetric model does not include the effect of particle diffusion and the dynamical evolution of the PWN. The results also show that the surface brightness of MSH 15–52 decreases with the increase of distance in all of the energy ranges, but the spatial changing characteristics are also not identical. Especially, the surface brightness has more obvious change in X-ray band (see Figure 7).

Zoom In Zoom Out Reset image size

In this paper, a self-consistent and spatially dependent model is presented to study the multiband emissions of PWNe. In the model, the temporal evolution of electrons is described as a Fokker–Planck equation and the processes of convection, diffusion, adiabatic loss, radiation loss, and photon–photon pair production are included. The temporal evolution of non-thermal photons is given by a spatially kinetic equation with the processes of synchrotron radiation, inverse Compton scattering off soft photons, synchrotron self-absorption, and the pair production including. Numerically solving these two equations, we studied the spatial variations of electron spectrum and multiband radiative properties of the PWN in MSH 15–52 with appropriate parameters.

In the modeling of the SED of a PWN, it is usually assumed that the particles in the system are uniform and isotropic. However, the observations implied that the photon index in the X-ray band is a function of position (e.g., Slane et al. 2000; Bocchino et al. 2001; Mangano et al. 2005; Schöck et al. 2010; Holler et al. 2012). Furthermore, Vorster & Moraal (2013) and Tang & Chevalier (2012) pointed out that the diffusion of particles plays an important role in changing the particle spectral index. In our calculations, we found that the SEDs of both electrons and non-thermal photons are all a function of distance (see Figure 1). We note that the photon indexes in the X-ray band have a more obvious increase with the increase of distance (see Figure 7). The main reason is that the spatial variation of electron SED in the higher energy part is dominated by the competition between synchrotron loss and diffusion (Lu et al. 2016), so the photon index is also a function of distance. The spatial variations in the X-ray band are especially more obvious.

With the appropriate parameters, we calculated the photon SED for the PWN in MSH 15–52. The emission from radio to X-rays is produced by the synchrotron radiation while the TeV γ-rays are dominated by the inverse Compton scattering with the FIR and NIR photons at all of the positions, the contribution of the SSC process is negligible (see Figure 2). The process of photon–photon pair production has little effect on the SED of MSH 15–52 due to the timescale ${\tau }_{\mathrm{esc}}\ll {\tau }_{\gamma \gamma }$ in all of the positions (see Figure 3). The surface brightness is also a function of distance and the photon indexes in the X-ray and very high energy γ-ray bands strongly depend on the position (see Figure 4). The observed SED of MSH 15–52 can be well reproduced in our model (see Figure 5). With the fitting results, we obtained that the magnetic fraction is $\eta =0.04$, which is slightly smaller than the value 0.05 adopted by Torres et al. (2014) and much smaller than the value 0.4 shown in Schöck et al. (2010). Meanwhile, the observed variations of both spectral index and surface brightness with positions based on the azimuthally integrated can also be well reproduced in our model (see Figure 6). Finally, we predicted the spatial changing characteristic of photon index and surface brightness from radio to TeV γ-rays, which should be tested in future observations.

We would like to thank an anonymous referee for very constructive comments. This work is partially supported by the National Natural Science Foundation of China (NSFC 11433004 and 11173020), Top Talents Program of Yunnan Province, the Natural Science Foundation of Yunnan Province (2013FB063), the Young Teachers Program of Yuxi Normal University, and the Program for Innovative Research Team (in Science and Technology) in University of Yunnan Province (IRTSTYN).