On the Gamma-Ray Nebula of Vela Pulsar. I. Very Slow Diffusion of Energetic Electrons within the TeV Nebula

TL;DR

This study models the slow diffusion of TeV electrons in the Vela pulsar's nebula, using analytical solutions to match observed spectral and spatial data, revealing diffusion rates much lower than in the interstellar medium.

Contribution

The paper provides the first analytical modeling of TeV electron diffusion in Vela X PWN, constraining the diffusion coefficient and suggesting slow diffusion is common in PWNe.

Findings

Diffusion coefficient for 10 TeV electrons is ~1×10^{26} cm^2/s.

Diffusion rate is over three orders of magnitude lower than in the ISM.

Results align with recent HAWC observations of Geminga nebula.

Abstract

High-energy particle transport in pulsar wind nebulae (PWNe) plays an essential role in explaining the characteristics revealed in multiwavelength observations. In this paper, the TeV-gamma-ray-emitting electrons in the Vela X PWN are approximated to be injected impulsively when the cocoon is formed due to the interaction between the SNR reverse shock and the PWN. By solving the diffusion-loss equation analytically, we reproduce the broadband spectral energy distribution and surface brightness profile simultaneously. The diffusion coefficient of TeV electrons and positrons, which is well constrained by the spectral and spatial properties of the TeV nebula, is thus determined to be $1 \times 10^{26}$\,cm$^{2}$\,s$^{-1}$ for 10\,TeV electrons and positrons. This coefficient is more than three orders of magnitude lower than that in the interstellar medium, in agreement with a constraint recently obtained from HAWC observations of a TeV nebula associated with the Geminga pulsar. These results suggest that slow diffusion of high-energy particles might be common in PWNe.

Tables

Table 1. Table 1: Fitting Parameters

Parameter	Quantity
$T_{\text{age}}$ (yr)	12000
$Q_{\text{0}}$ (cm-3 s-1)	$6.7\times {10}^{36}$
$d_{\text{Vela}}$ (pc)	287
${\tau }_{\text{s}}$ (yr)	7500
$B_{\text{eff}}$ ($\mu$G)	6
$D_0$ (cm2 s-1)	${10}^{26}$
$\alpha$	1.7

Equations

$$\frac{\partial }{\partial t}f(\gamma ,r,t)=\frac{D(\gamma )}{r^2}\frac{\partial }{\partial r}{r}^2\frac{\partial }{\partial r}f(\gamma ,r,t)+\frac{\partial }{\partial \gamma }\left(Pf\right)+Q_{\text{inj}}(\gamma ,t),$$

$$Q_{\text{inj}}(\gamma ,t)=Q_0\delta (t-{\tau }_{\text{s}}){\left(\frac{\gamma }{10^7}\right)}^{-\alpha },$$

$$\frac{\partial \mathcal{F}}{\partial z}=D_1(z)\frac{{\partial }^2\mathcal{F}}{\partial r^2}$$

$$\displaystyle f(\gamma,r,t)=\left\{\begin{array}[]{ll}\frac{\gamma^{2}_{t}Q_{\textup{inj}}(\gamma_{t},t)}{\gamma^{2}\pi^{3/2}r^{3}_{\textup{dif}}}\exp\left({-\frac{r^{2}}{r^{2}_{\textup{dif}}}}\right)&\ \ \gamma\leq\gamma_{{}_{\textup{max}}},\\ 0&\ \ \gamma>\gamma_{{}_{\rm max}}.\\ \end{array}\right.$$

$$r_{\text{dif}}(\gamma ,t)=2\sqrt{D(\gamma )t\frac{1-(1-\gamma /{\gamma }_{{}_{\text{max}}}{)}^{1-\delta }}{(1-\delta )\gamma /{\gamma }_{{}_{\text{max}}}}},$$

$$F_{\text{LS}}={\int }_{-\infty }^{\infty }fdl={\int }_{-\infty }^{\infty }\frac{{\gamma }_t^2{Q}_{\text{inj}}({\gamma }_t)}{{\gamma }^2{\pi }^{3/2}{r}_{\text{dif}}^3}\exp \left(-\frac{l^2+R^2}{r_{\text{dif}}^2}\right)dl=\frac{{\gamma }_t^2{Q}_{\text{inj}}({\gamma }_t)}{{\gamma }^2\pi r_{\text{dif}}^2}\exp \left(-\frac{R^2}{r_{\text{dif}}^2}\right).$$

$${\int }_0^{T_{\text{age}}}{\int }_{R_{\text{in}}}^{R_{\text{out}}}{F}_{\text{LS}}2\pi RdRdt=\frac{{\gamma }_t^2{Q}_0{\left(\gamma /10^7\right)}^{-\alpha }}{{\gamma }^2}\left[\exp \left(-\frac{D_{\text{in}}^2}{r_{\text{dif}}^2(\gamma ,T_{\text{age}}-{\tau }_s)}\right)-\exp \left(-\frac{D_{\text{out}}^2}{r_{\text{dif}}^2(\gamma ,T_{\text{age}}-{\tau }_s)}\right)\right],$$

Full text

ISM: supernova remnants — ISM: individual objects (Vela X) — diffusion — (ISM:) cosmic rays

††journal: ApJ

1 INTRODUCTION

The unexpected excess of high-energy positrons ($\geq 10$ GeV) observed by PAMELA (Adriani et al., 2013), AMS-02 (Accardo et al., 2014), VERITAS (Staszak & VERITAS Collaboration, 2015), Fermi-LAT (Ackermann et al., 2012; Abdollahi et al., 2017), and DAMPE (DAMPE Collaboration et al., 2017) has been explored extensively by many researchers, yet its origin is still a matter of debate. While dark matter particles are believed to be a possible positron source, pulsars and pulsar wind nebulae (PWNe) are more natural astrophysical contributors owing to their high pair production capability (Hooper et al., 2009). Pulsars inside PWNe produce ultra-relativistic wind of electron-positron pairs (hereafter, “electrons” refers to both electrons and positrons). A shock forms as the wind interacts with slow-moving supernova ejecta, and the electrons are thermalized and/or reaccelerated by the shock (Gaensler & Slane, 2006). The resulting high-energy electrons may escape from PWNe, and propagate in the nearby interstellar medium (ISM). Considering radiative energy loss experienced by high-energy electrons while propagating, TeV positron sources must be located within hundreds of parsecs from Earth if they are to contribute to the positron excess.

However, this scenario is questioned by the recent HAWC’s detection of extended TeV $\gamma$-rays from the Geminga PWN and PSR B0656+14. Although this observation confirms the existence of TeV electrons in PWNe, the diffusion coefficients, constrained by the $\gamma$-ray surface brightness profile, are too small to account for the positron flux observed at Earth (Abeysekara et al., 2017). Two-zone diffusion models seem to partially remedy this issue (Hooper et al., 2017; Fang et al., 2018). Nevertheless, the local PWNe origin of the positron excess is oppugned by non-detection of the Geminga PWN in Fermi observations (Xi et al., 2018), and the nature of positron excess is still under debate (Shi & Liu, 2019).

The Vela pulsar, lying at a distance of 287 pc to Earth (Dodson et al., 2003a), is one of the positron contributor candidates in the TeV energy range. Thanks to its proximity, the Vela X PWN is spatially resolved in multiwavelength observations from the radio to $\gamma$-rays. The extended radio nebula (ERN) of Vela X, which is spatially correlated with the GeV nebula seen by Fermi, is much larger than the X-ray “cocoon” (Markwardt & Ögelman, 1995; Frail et al., 1997). A counterpart of the cocoon was later revealed in TeV $\gamma$-rays (Aharonian et al., 2006). Although there are extensive observational and theoretical studies, the formation of the complex frequency-dependent morphology of the PWN remains to be addressed (Tibaldo et al., 2018). Two electron components have been assumed to explain the frequency dependent morphology (de Jager et al., 2008). However, challenges arise when TeV emission beyond the cocoon (extended to $\sim 1.2^{\circ}$) is detected (Abramowski et al., 2012; Reynolds et al., 2017). We hereafter refer to the TeV-$\gamma$-ray-emitting plasma as the TeV nebula. Hinton et al. (2011) argued that diffusive escape of high-energy electrons must be taken into consideration to explain the soft GeV spectrum of the ERN while the TeV nebula formed over the past few hundreds of years, so that high-energy electrons are still trapped in a cocoon structure. The escaped electrons, if diffused to Earth, should contribute to the local cosmic ray positron spectrum, suggesting that the Vela X PWN may be one of the positron contributors (Hinton et al., 2011). However, AMS-02 measurement shows that the positron fraction saturates at $\sim 300$ GeV and likely drops above $\sim 500$ GeV (Recchia et al., 2018; Aguilar et al., 2019), disfavoring the pulsars and PWNe as the origin of the positron excess.

Tibaldo et al. (2018) found that the TeV nebula of the Vela X PWN has a very hard spectrum above 100 GeV. This spectrum suggests that the emission is dominated by electrons at tens of TeV. Such a distribution is a natural consequence of a hard electron spectrum with a power-law index of less than 2 being subjected to synchrotron and IC radiative energy losses. The spatial offset of the cocoon from pulsar suggests that it was likely caused by interaction with the reverse shock of the supernova remnant (SNR) thousands of years ago (Slane et al., 2018). We therefore propose that the TeV nebula formed via impulsive injection of a hard power-law distribution of energetic electrons before the PWN was displaced by interaction with the reverse shock several thousands year ago. The consequent parallel diffusion and radiative energy loss of these high-energy particles in a large-scale magnetic field can naturally account for the spatial and spectral properties of the TeV nebula (Dodson et al., 2003b).

Based on the escape model proposed by Hinton et al. (2011), Huang et al. (2018) argued that diffusion coefficient in the Vela X PWN should be less than $10^{28}$ cm2 s*-1* for 10 TeV electrons to avoid excessive cosmic ray electron flux above 10 TeV. Here we suggest that the extended TeV emission beyond the X-ray cocoon arises from diffusion of high-energy electrons and use the spatial and spectral properties of the TeV nebula to constrain the diffusion coefficient in this nebula directly. The electrons are injected as a power-law function of energy, and the diffusion-loss equation can be solved analytically. Since the effective magnetic field strength for radiative energy loss can be constrained by the high-energy cutoff of the TeV spectrum and the $\gamma$-ray emission is dominated by the electrons at tens of TeV, the diffusion coefficient is essentially the only key parameter one can adjust to fit the TeV $\gamma$-ray radial brightness profile and the spectra of the inner region ($0^{\circ}$–$0.8^{\circ}$, slightly larger than the X-ray cocoon) and an outer ring ($0.8^{\circ}$–$1.2^{\circ}$, diffused emission beyond the cocoon; Abramowski et al. 2012). Our model is described in §2. The application of this model to the TeV nebula of the Vela pulsar is presented in §3. The discussion and conclusions are provided in §4. In a follow-up paper (Bao & Chen 2019, hereafter Paper II), we explain the spectrum of the ERN and the associated GeV $\gamma$-ray emission.

2 MODEL DESCRIPTION

We assume the TeV surface brightness profile stems from the diffusion and radiation losses of high-energy electrons. For simplicity, a simple diffusion-loss model is adopted to constrain the diffusion coefficient within the TeV nebula of Vela pulsar. The transport of electrons can be described by

[TABLE]

where $f(\gamma,r,t)$ is the electron distribution function, $r$ is the radial distance to the center of the TeV nebula, $D(\gamma)=D_{0}\left(\gamma/\gamma_{{}_{\textup{10\,TeV}}}\right)^{\delta}$ is an energy-dependent diffusion coefficient (with $\gamma_{{}_{\textup{10\,TeV}}}$ the Lorentz factor of 10 TeV electrons), $P$ is the radiation energy loss rate, and $Q_{\textup{inj}}$ is the electrons’ injection rate associated with the pulsar’s spin-down. Since the slow-diffusion regions ($\sim$ 30–50 pc in radii) in the two-zone diffusion models (Fang et al., 2018; Huang et al., 2018) are apparently larger than the Vela SNR ($\sim 20$ pc in radius), we assume a spatially independent diffusion coefficient. The $\gamma$-ray spectral indices are almost a constant throughout the cocoon (Abramowski et al., 2012), therefore the diffusion coefficient is expected to have a weak energy dependence; we thus fix $\delta=1/3$ (Kolmogorov diffusion, as is assumed in Abeysekara et al. 2017). Although a broken power law is usually used to explain the broadband spectral energy distribution (SED) of PWNe, a single power law is sufficient to explain the broadband spectrum of the Vela X PWN (see Paper II), therefore the injection rate is assumed to be

[TABLE]

where $Q_{0}$ is the injection constant for the electrons with $\gamma=10^{7}$. Here we assume that most of plasma injected before the reverse shock-PWN interaction is compressed into the ERN, and approximate the injection to be impulsive for the following reasons. The total energy in the cocoon is estimated to be $E_{\textup{cocoon}}\sim 1.5\times 10^{46}$ erg (Abdo et al., 2010), much lower than the energy in the ERN ($5\times 10^{48}$ erg; Abdo et al. 2010). Hydrodynamic simulations have revealed that the cocoon was formed $\sim 4$ kyr ago (Slane et al., 2018), when the SNR reverse shock disrupted the PWN, creating a tail (the cocoon) to the south, which contains a small fraction of fresh plasma released by the pulsar. Meanwhile, the relatively faint TeV emission and bright X-ray emission in the very vicinity of the Vela pulsar imply a strong magnetic field formed after the reverse shock interacted with the PWN (see, e.g., Hinton et al., 2011) that can cool the electrons (injected after the X-ray cocoon was created) to low energy in a short time ($\sim 10^{2}(B/400\ \mu\textup{G})^{-2}(E/700\ \rm GeV)^{-1}$ yr). Therefore, injection after the cocoon was formed does not contribute to TeV electrons. The cooling timescale of 70 TeV electrons is $\sim 4.4\times 10^{3}(B/6\mu\textup{G})^{-2}$ yr (Hinton et al., 2011), hence the very high-energy electrons injected $\sim 4$ kyr ago have been cooled down, giving rise to the observed TeV spectral cutoff.

Given the very hard GeV spectrum of the TeV nebula (Abramowski et al., 2012; Tibaldo et al., 2018), $\alpha$ needs to be less than 2, giving rise to an energy distribution piling up where the radiative cooling timescale is equal to the age of the injected electrons (see Figure 1). To get the structure of the nebula, we assume that the nebula is spherical symmetric111 Although the X-ray cocoon and its TeV counterpart is hardly spherically symmetric in the central region ($\sim 0.8^{\circ}$ in length), spherical symmetry could be a good approximation for the spatially integrated spectra and the larger extension of TeV emission out to $2^{\circ}$ in radius. On the other hand, the asymmetry in the central region would only affect the first two data points in the surface brightness profile. Also, the radial TeV brightness profile is given from circular annuli (Abramowski et al., 2012), so spherical symmetry could be a convenient technical treatment to fit the data. Moreover, the diffusion coefficient is derived from fitting the relative surface brightness, therefore the absolute magnitude of the injected energy is of less importance. The asymmetry does not affect our constraint on the diffusion coefficient obtained. The SED, however, depends on the total number of the $\gamma$-ray-emitting electrons injected in the TeV nebula after interaction with the reverse shock, which can be adjusted by changing the injection constant $Q_{0}$. Here we only consider the synchrotron and inverse Compton (IC) losses, and the loss rate can be described by $d\gamma/dt=-p_{2}\gamma^{2}$, where $p_{2}=3.23\times 10^{-20}(B/5\,\mu G)^{2}\ \textup{s}^{-1}$. Following the procedure described in Atoyan et al. (1995), we define $\mathscr{F}=(d\gamma/dt)rf(\gamma,r,t)$, hence Equation 1 can be written as

[TABLE]

(see Atoyan et al. 1995 for the definition of $z$ and $D_{1}(z)$); consequently, the electron distribution function $f(\gamma,r,t)$ can be obtained:

[TABLE]

where $\gamma_{\textup{max}}=1/(p_{2}t)$, $\gamma_{t}=\gamma/(1-p_{2}t\gamma)$ is the initial energy of the electrons that are cooled down to $\gamma$ after time $t$. Atoyan et al. (1995) showed, $r_{\textup{dif}}\approx 2\sqrt{D(\gamma)t}$ for $\gamma<0.5\gamma_{{}_{\rm max}}$. Then we integrate the electron distribution function along the line of sight $l$ (noting that $l^{2}+R^{2}=r^{2}$, with $R$ being the projection distance):

[TABLE]

According to Abramowski et al. (2012), the radial brightness profile is extracted from annuli with the same width of 12’. To obtain the $\gamma$-ray fluxes observed now, $F_{\textup{LS}}$ is further integrated over each annulus and pulsar lifetime:

[TABLE]

where $R_{\textup{in}}$ and $R_{\textup{out}}$ represent the inner and the outer radius of an annulus, respectively, and $T_{\textup{age}}$ represents the age of the Vela PWN.

3 Application to Vela X

The Vela pulsar and its associated SNR have been extensively studied thanks to the pulsar’s short distance of 287 pc obtained using VLBI parallax (Dodson et al., 2003a). The age of the pulsar derived from its proper motion is $9000$–$27000$ yr (Aschenbach et al., 1995). It is the first H.E.S.S. TeV source showing a prominent high-energy spectral cutoff (Abramowski et al., 2012). Aharonian et al. (2006) showed that the magnetic field strength in the X-ray cocoon is $\sim$ 4 $\mu$G based on its X-ray flux. The formation of the X-ray cocoon is attributed to the reverse shock-PWN interaction thousands of years ago (Slane et al., 2018), and the TeV counterpart can be explained by the electron diffusion from the cocoon. We focus on the $\gamma$-ray data and consider the cosmic microwave background radiation (CMB) and the far infrared radiation (FIR) as the background photons for calculation of the $\gamma$-rays via the IC process. We assume blackbody spectrum with temperature of 25 K and energy density of 0.2 eV cm*-3* to approximate the FIR radiation field presented in Porter et al. (2008). We fit the $\gamma$-ray spectra of both the inner and outer regions as well as the surface brightness profile given in Abramowski et al. (2012). The model parameters are listed in Table 1, where $d_{\textup{Vela}}$ is the distance to the Vela pulsar, $B_{\textup{eff}}$ is the effective magnetic field strength ($B_{\textup{eff}}^{2}=B^{2}+8\pi u_{{}_{\textup{CMB}}}+8\pi u_{{}_{\textup{FIR}}}$, with $u_{{}_{\textup{CMB}}}$ and $u_{{}_{\textup{FIR}}}$ being the energy density of CMB and FIR photons, respectively). Because the energy losses are dominated by synchrotron radiation and IC off low-energy photons, the energy losses are assumed to be quadratic in energy. However, the Klein-Nishina effect is incorporated into our calculation of SEDs. Because the electron spectrum is very hard and unbroken below 70 TeV (see Figure 1), IC off starlight is severely suppressed by the Klein-Nishina effect. Meanwhile, the energy density of starlight is as low as 0.3 eV cm*-3*, with a temperature of $\sim 3000$ K (Porter et al., 2008). The IC off starlight contributes $\sim 1\%$ of the $\gamma$-ray emission and is thus neglected.

Figure 1 shows the spatially integrated electron distribution of the TeV nebula for parameters given in Table 1. The spectrum is cut off at $\sim 70$ TeV due to the radiative cooling. Figure 2 shows the spectral fit to the inner and outer rings (left panel) and fit to the normalized brightness profile (right panel). The SED can be fit well with $\alpha=1.7$, which is in agreement with the radio spectrum of the ERN (Paper II). Since the magnetic field is weak in the X-ray cocoon, the synchrotron lifetime of 10 TeV electrons is $\approx 3\times 10^{4}(B/6\mu\textup{G})^{-2}$ yr (Hinton et al., 2011), larger than the age of the Vela SNR, therefore the corresponding Fermi $\gamma$-ray spectrum (Tibaldo et al., 2018) is a single powerlaw, consistent with theoretical expectations for $\gamma$-ray-emitting electrons with a power-law distribution. The spectrum of the outer ring is slightly harder than the inner ring due to the increase of the diffusion coefficient with energy. The spectrum of the outer ring is slightly harder than the inner ring due to the increase of the diffusion coefficient with energy. In Figure 3, we show the spatial distribution of the $\gamma$-ray indices for different $\delta$. The $\gamma$-ray indices are calculated for nine annuli centered in the central position of the cocoon, with the same width of $0.1^{\circ}$. It can be seen that, for $\delta=1$, the $\gamma$-ray indices have apparent spatial variation, whereas for $\delta=1/3$, the spatial distributions of $\gamma$-ray indices are consistent with the indices extracted from sectors along the cocoon. The dependence of the TeV $\gamma$-ray brightness profile and spectra on the diffusion coefficient is plotted in Figure 4. It shows that the diffusion coefficient is well constrained ($D_{0}\approx 1\times 10^{26}$ cm2 s*-1*) by the observed spectra and brightness profile.

4 discussion and conclusions

The discovery of a high-energy spectral cutoff in combination with a recent detection of a very hard GeV spectrum indicate that the $\gamma$-ray emission of the TeV nebula of the Vela X PWN is dominated by electrons with energies of tens of TeV via the IC process. The weak magnetic field strength derived from X-ray observations and the high-energy cutoff are consistent with the radiative cooling of TeV electrons injected upon the formation of the cocoon, and the small total energy in the TeV nebula is consistent with the impulsive injection. In this paper, with a simple diffusion and radiative loss model, we demonstrate that in the scenario in which the TeV cocoon formed via interaction of the PWN with the reverse shock several thousand years ago, the $\gamma$-ray spectra and radial brightness profile of the TeV nebula can be reproduced. Since the magnetic field strength is constrained by the injection time and the electron index $\alpha$ can be constrained by the radio SED of the ERN, the diffusion coefficient $D_{0}$ is the only key parameter that can be well constrained by the brightness profile and radial spectral change of the TeV nebula. The diffusion coefficient of TeV electrons and positrons is thus determined to be $1\times 10^{26}$ cm2 s*-1* for 10 TeV electrons, which is more than three orders of magnitude lower than the typical value in the ISM. Our results, in combination with an early result from HAWC observations of the Geminga PWN suggest that slow diffusion might be common in TeV PWNe.

At 10 TeV, the diffusion coefficient is almost the same as the Bohm diffusion ($8.3\times 10^{25}$ cm2 s*-1*) in a magnetic field of $4\,\mu$G, which implies sub-Bohm diffusion at even higher energies. Since $r_{\textup{dif}}\sim 2\sqrt{DT_{\textup{age}}}\sim 6$ pc for 100 TeV electrons, such slow diffusion could hardly transport the positrons below 100 TeV out of the Vela SNR. One possible explanation for such a slow diffusion in the cocoon is efficient trapping of electrons in a magnetic mirror. Such trapping needs to be extended to the whole TeV nebula. Multiwavelength observations can be used to probe the magnetic field structure (Abramowski et al., 2012). More detailed exploration of such a scenario will be presented in a future paper.

The energy of high-energy electrons inside the cocoon is estimated to be $E_{\textup{cocoon}}\sim 1.5\times 10^{46}$ erg (Abdo et al., 2010), which is extremely low, but much higher than the energy of $\sim 10^{44}$ erg carried by the magnetic field. The cocoon is expected to form due to the reverse shock-PWN interaction (Slane et al., 2018), hence the bulk of the pulsar spin-down energy (and bulk of TeV electrons) is deposited into the ERN. The SED of the ERN can be explained by intense radiation loss in the reverberation phase for the interaction of the nebula with reverse shock; the interpretation of the broadband spectrum of the ERN is delegated to an ensuing paper (Paper II). Moreover, if the magnetic field strength in the ERN is high, the soft GeV spectrum may be attributed to hadronic processes.

Slow diffusion around SNRs has been suggested for years (see e.g., Fujita et al., 2009; Gabici et al., 2010; Li & Chen, 2010; Yan et al., 2012). Theoretical works also show that cosmic rays may diffuse slowly due to self-generation waves and wave-wave turbulent cascading from a scale, which is in concordance with the size of the SNR (see e.g., Blasi et al., 2012), and may also work inside PWNe with weak magnetic fields.

We thank the anonymous referee for the helpful comments. This work is supported by the National Key R&D Program of China under grants 2018YFA0404203, 2015CB857100, and 2017YFA0402600, NSFC under grants 11773014, 11633007, 11851305, U1738122, and 11761131007, and the International Partnership Program of Chinese Academy of Sciences under grant 114332KYSB20170008.