Modeling the broadest spectral band of the Crab nebula and constraining the ions acceleration efficiency

It is a concensus that pulsar wind nebulae (PWNe) are efficient accelerators for the extremely relativistic leptons (electrons/positrons). A pulsar located in a PWN releases its rotational energy by driving an ultra-relativistic magnetized wind composed of electron-positron pairs. The pulsar wind interacts with the ambient medium and slows down at the so-called termination shock where the acceleration of the pulsar wind leptons occurs (Gaensler & Slane, 2006; Kargaltsev et al., 2015). It is generally believed that the accelerated leptons generate the electromagnetic emission from radio to very-high-energy $\gamma$-rays via synchrotron and inverse Compton (IC, probably including synchrotron self-Compton process) mechanism, namely the leptonic model (e.g., Kennel & Coroniti, 1984; Venter & de Jager, 2007; Zhang et al., 2008; Tanaka & Takahara, 2010; Bucciantini et al., 2011; Martín et al., 2012; Torres et al., 2014).

Although the leptonic model can well explain the observational features, an open question still remains whether, or how large a fraction of, protons (or more generally, ions) may exist in pulsar wind and are further accelerated in PWN. Theoretically, it has been suggested that protons can also be extracted from the surface of a neutron star and injected into the pulsar wind (Cheng et al., 1990; Arons & Tavani, 1994). Protons carrying a substantial amount of the wind energy were also predicted by the resonant cyclotron absorption model (Hoshino et al., 1992; Gallant & Arons, 1994). Moreover, ion acceleration in pulsars/PWNe was explored in some recent works (e.g. Fang et al., 2013; Chen & Beloborodov, 2014; Philippov & Spitkovsky, 2014; Kotera et al., 2015; Lemoine et al., 2015; Guépin et al., 2020), suggesting that pulsars/PWNs may be the sources of Cosmic Rays. Especially, Guépin et al. (2020) concluded that the energetic protons accelerated in the pulsar magnetosphere via the magnetic reconnection can take away the spin-down energy with a fraction of 0.2 to 4.0 per cent.

These theoretical works on proton acceleration in pulsars or PWNe indicate that the $\gamma$-ray emission of a PWN may contain a contribution from hadrons. Due to the difficulties of differentiating the leptonic and hadronic origin of the $\gamma$-ray emission, the most direct method to test this is the neutrino experiments. With the current stacked IceCube data toward 35 $\gamma$-ray-bright PWNe, it can only be inferred that PWNe are not the hadronic-dominated sources (Aartsen et al., 2020). For some individual sources, however, the hadronic component was indeed invoked to explain the $\gamma$-ray emission, e.g. Crab (Atoyan & Aharonian, 1996), Vela X (Horns et al., 2006; Zhang & Yang, 2009), PWN G54.1+0.3 (Li et al., 2010), and DA 495 (Coerver et al., 2019), although the purely leptonic model can also work, especially for the Crab nebula.

The Crab nebula, powered by pulsar PSR J0534+2200, is the best studied prototype PWN. Based on the Kennel & Coroniti (1984) 1D flow model, Atoyan & Aharonian (1996) calculated the Crab’s spectrum extending to the $\gamma$-ray band and showed that the pure IC emission underpredicts the observed flux, indicating possible additional contribution from protons. Nontheless, in a multidimensional simulation for the flow structure in the Crab, Volpi et al. (2008) found that the leptonic model can explain the EGRET data. With the updated GeV data obtained by Fermi (Buehler et al., 2012), the broadband spectrum of the Crab was widely modeled using the leptonic model (e.g. Zhang et al., 2008; Tanaka & Takahara, 2010; Bucciantini et al., 2011; Martín et al., 2012) without any need of a proton contribution.

Recently, $>100$ TeV photons from the Crab were detected by the Tibet AS+MD (Amenomori, 2019) and HAWC (Abeysekara et al., 2019) experiments. These new data have helped presenting the most complete spectrum across the electro-magnetic emission window, and may provide information on the proton acceleration. In this paper, we model the broadest spectral data of the Crab nebula and explore the possibility that the very-high-energy $\gamma$-ray emission contains the contribution from the energetic protons. The spectral evolution model is presented in section 2. In section 3, the hadronic component of $\gamma$-rays is constrained via the Markov Chain Monte Carlo (MCMC) method. Some discussions and conclusions are given in section 4.

The global multiband radiative properties of a PWN are generally described by the “1-Zone model" (see Bucciantini, 2014, for a review). The evolution of the lepton spectrum in the emission zone is given by the continuity equation in the energy space

where $\gamma$ is the Lorentz factor of the leptons, $N(\gamma,t)$ is the lepton distribution function, $\dot{\gamma}(\gamma,t)$ is the summation of the energy losses including the synchrotron radiation, the IC process, and the adiabatic loss (see Li et al., 2010), and $\tau(\gamma,t)$ represents the escape time which can be estimated via Bohm diffusion (e.g. Zhang et al., 2008). $Q(\gamma,t)$ represents the injection rate of leptons from shock acceleration process into the PWN emission zone per unit energy at a certain time and is generally assumed to be a broken power-law function

where $\gamma_{\rm bre}$ denotes the break energy, and the parameters $\alpha_{1}$ and $\alpha_{2}$ are spectral indices. The normalization is determined from the total energy in leptons ($W_{\rm e}$) extracted from the spin-down energy of a pulsar ($L_{\rm sd}$) via $W_{\rm e}=\eta_{\rm e}L_{\rm sd}=(1-\eta_{\rm B})L_{\rm sd}$, where $\eta_{\rm e}$ and $\eta_{\rm B}$ are the leptonic and magnetic energy fraction, respectively. The maximum energy of the injected leptons is obtained by introducing a parameter $\varepsilon$, indicating that the Larmor radius of the lepton must be less than the termination shock radius (see Eq. (8) in Martín et al., 2012).

Solving equation 1, we can obtain the lepton distribution at a certain age. Then we can calculate the multiwavelength nonthermal photon spectrum via the synchrotron and IC processes. We apply this evolutionary model to Crab nebula and fit its observational multiband data including the radio (Baldwin, 1971; Macías-Pérez et al., 2010), infra red (IR, Ney & Stein, 1968; Grasdalen, 1979; Green et al., 2004; Temim et al., 2006), optical (Veron-Cetty & Woltjer, 1993), X-ray (Hennessy et al., 1992), soft $\gamma$-ray, (Kuiper et al., 2001), GeV $\gamma$-ray (Fermi: Buehler et al. 2012), and TeV $\gamma$-ray (HEGRA: Aharonian et al. 2004; H.E.S.S.: Aharonian et al. 2006; MAGIC: Aleksić et al. 2015; ARGO-YBJ: Bartoli et al. 2015; Tibet AS+MD: Amenomori 2019) data. For the IC process, the seed photon fields include (i) the 2.7-K cosmic microwave background (CMB) radiation; (ii) the excess FIR radiation with a temperature of 70 K and an energy density of 0.5 ${\rm eV\ cm^{-3}}$ (Marsden et al., 1984); (iii) the NIR radiation with a temperature of 5000 K and an energy density of 1.0 ${\rm eV\ cm^{-3}}$ (Aharonian et al., 1997); and (iv) the synchrotron radiation photon.

In the calculation, we adopt the distance $d=2.0$ kpc, the radius $R_{\rm PWN}=2.1$ pc, the initial spin-down luminosity $L_{0}=3.1\times 10^{39}\ \ {\rm erg}\ {\rm s}^{-1}$, and the breaking index $n=2.509$ (Trimble, 1968; Lyne et al., 1988; Taylor et al., 1993). The eye-fit parameters are $\alpha_{1}=1.55$, $\alpha_{2}=2.5$, $\gamma_{\rm bre}=8.5\times 10^{5}$, $\eta_{\rm B}=0.012$, and $\varepsilon=0.3$, which are similar to previous results (e.g. Tanaka & Takahara, 2010; Martín et al., 2012). The corresponding spectral energy distribution (SED) is plotted in Figure 1, which includes the $\gamma$-ray absorption based on the cross section of $\gamma$-$\gamma$ interaction (Gould & Schréder, 1967) and the interstellar radiation field in Galaxy (Shibata et al., 2011). As can be seen, the IC process can well explain the $\gamma$-ray data up to $\sim$500 TeV and confirm the leptonic origin for the broadband emission. According to the fitted results, at the Crab’s age of 950 yr, the average magnetic field strength in the emission zone and the maximum energy of injected leptons are about $100\ \mu{\rm G}$ and 4 PeV, respectively. If protons can also be injected into the acceleration process, their maximum energy may be expected to reach similar order. Moreover, protons with this energy can be well confined in the nebula at the magnetic field of $100\ \mu{\rm G}$ according to the Larmor formula, and may contribute to the $\ga 100$ TeV photons via p-p interaction. In the next section, we study the probability that protons contribute some fraction of $\gamma$-ray emission from Crab.

Although the broadband emissions from the Crab are explained by the leptonic process, the $\gamma$-ray emission may have hadronic component which can be constrained by the $>$100 TeV photons. In order to quantitatively investigate this possibility, one should add the hadronic process to the above evolutionary model, and then constrain the model parameters by using a statistical method. However, for avoiding such time-costing calculation, we can do the estimation in a simplified way. We directly assume that the lepton distribution at the current age in the emission zone has a smooth broken power-law form with a high-energy cutoff:

where $\gamma_{\rm cut}$ is the cut-off energy, $s$ the smooth parameter, and $\beta$ the cutoff shape parameter. The normalization $A_{\rm e}$ can be calculated from the total energy in leptons $W_{\rm e}$. To check its accurateness, this distribution is firstly employed to fit the electrons in the evolutionary model with the above eye-fit parameters but for three different ages $t=500$, 950, and 5000 yr, which are shown in the black dotted lines in left panel of Figure 2. As it can be seen, the electrons in the evolutionary model can almost be described by the exponential cutoff smooth broken power law distribution.

We now apply this simplified leptonic model to refit the broadband data of the Crab and use the MCMC approach to constrain the model parameters (Model A). There are eight parameters in total: $\alpha_{1}$, $\alpha_{2}$, $\gamma_{\rm bre}$, $\gamma_{\rm cut}$, $B_{\rm PWN}$, $W_{\rm e}$, $s$, and $\beta$. The best-fit parameters are listed in Table 1 and the corresponding SED is plotted in Figure 3, in which the result of the evolutionary model is also displayed in magenta. The comparison of the electron distribution in the evolution model and that in Model A is also shown in Figure 2. We can see that this leptonic model is good enough to explain the broadband data, especial for the $\gtrsim$ 10 TeV data which will play an important role in constraining the hadronic component. Thus, we will use it instead of the evolutionary model to constrain the contribution of the protons. In addition, we can directly see that the maximum energy of the leptons at current age needs to be $\sim 3$ PeV to explain the data. This hints that PWNe are one of potential sources of Galactic CRs.

Next, we add hadronic photons resulted from the neutral pion-decay process to the simplified leptonic model. For the proton spectrum, we adopt the exponential cutoff power-law distribution,

where $\alpha_{\rm p}$ and $E_{\rm c,p}$ are the the proton index and the cutoff energy, respectively. The normalization $A_{\rm p}$ can be calculated from the total energy in protons $W_{\rm p}$. For p-p process, we use the analytic photon emissivity developed by Kelner et al. (2006), including the enhancement factor of 1.84 due to contribution from heavy nuclei (Mori, 2009). The density of target gas interacted by relativistic protons is not well constrained. The mean gas density in the nebula can be estimated from the total nebular mass (gas plus dust) of $7.2\pm 0.5M_{\odot}$ (Owen & Barlow, 2015) and the nebular radius $R_{\rm PWN}=2.1$ pc, giving a mean density $n\sim 5\ {\rm cm}^{-3}$. Considering the fact that relativistic particles were partially captured in the dense filaments in the Crab nebula, the effective gas density for interactions of relativistic protons may be much higher than the mean density (Atoyan & Aharonian, 1996). Here, we adopt the density of the target gas $n_{\rm t}=10\ {\rm cm}^{-3}$ as a fiducial value. In addition, the nonthermal bremsstrahlung process can be neglected at this density level (e.g., Atoyan & Aharonian, 1996).

Due to lack of constraint on the hadronic component, we fix the proton index and the cutoff energy. We set $E_{\rm c,p}=3$ PeV and consider two cases for the proton index: $\alpha_{\rm p}=2.0$ (Model B, generally assumed particle index) and $\alpha_{\rm p}=1.55$ (Model C, the same as the lepton index below the break). For this lepto-hadronic model, there are nine parameters in total: $\alpha_{1}$, $\alpha_{2}$, $\gamma_{\rm bre}$, $\gamma_{\rm cut}$, $B_{\rm PWN}$, $W_{\rm e}$, $s$, $\beta$, and $W_{\rm p}$. The best-fit parameters are listed in Table 1 and the corresponding broadband SED is plotted in Figure 4. The fitted values of all parameters except $W_{\rm p}$ in Models B and C are almost the same as in Model A. This means the contribution from hadronic process is negligible, which also can be seen from the result that the hadronic component (the blue line in Figure 4) is far below the observed $\gamma$-ray data. In order to void exceeding the $>100$ TeV data detected by Tibet AS+MD, it gives a constraint on the energy in protons as $W_{\rm p}\lesssim 1.9\times 10^{47}$ erg and $W_{\rm p}\lesssim 1.0\times 10^{47}$ erg at $2\sigma$ level for $\alpha_{\rm p}=2.0$ (Model B) and $\alpha_{\rm p}=1.55$ (Model C), respectively.

In the calculation, we see that the residual is mainly contributed by the radio to soft $\gamma$-ray data. To bring the $\gamma$-ray data into play, we will only use the GeV-TeV data (i.e., above 1 GeV) to constrain the model parameters. In doing this, we first fit the data from radio to soft $\gamma$-ray band by using a smooth broken power law function, which will describe the seed photons in the SSC process (see the gray line in left panel of Figure 5). Two cases $\alpha_{\rm p}=2.0$ (Model D) and $\alpha_{\rm p}=1.55$ (Model E) are also considered.

Just using $\gamma$-ray data, some of parameters do not converge very well if we do the same MCMC-fitting routine, and hence we make the following attempts. First, the cutoff shape parameter $\beta$ is the most serious one because it is fully determined by the $\sim$100 MeV Fermi data. Considering the results in Models B and C, we round it as $\beta=3$. In fact, due to the Klein-Nishina effect, this parameter will not obviously affect the spectral shape around 100 TeV. Secondly, different from the narrow peak of the 1D probability in Models B and C, the smooth parameter $s$ is in a large range 0.3-0.6 and is obviously smaller that in Model B and C. This means that the break smoothness of the leptonic spectrum is essentially determined by the IR-to-X-ray data, and the IC/$\gamma$-ray peak is broader than the synchrotron peak. It may be the reason that the IC peak is broader than suggested by several previous models which is pointed out by MAGIC group (Aleksić et al., 2015; MAGIC Collaboration et al., 2020). The broader $\gamma$-ray peak not only impacts the convergence of $s$ but also $\alpha_{1}$, $\alpha_{2}$, and $\gamma_{\rm bre}$. In order to constrain the parameters, we fix $\alpha_{1}=1.55$ for the index below the break and take $s=0.5$¹¹1We have also tried the case of $s=0.6$ which gives similar results. as a example for the smooth parameter.

Then there are five parameters needed to be determined in Models D and E: $\alpha_{2}$, $\gamma_{\rm bre}$, $\gamma_{\rm cut}$, $W_{\rm e}$, and $W_{\rm p}$. The best-fit parameters are listed in Table 1 and the corresponding SED is plotted in Figure 5. Different from Models B and C, the hadronic component dominates the $\gamma$-ray spectrum above 100 TeV (see the blue lines in left panels of Figure 5). This implies that the total energy in protons is constrained by the Tibet AS+MD data as $W_{\rm p}=2.6\times 10^{48}$ erg and $W_{\rm p}=1.4\times 10^{48}$ erg for $\alpha_{\rm p}=2.0$ (Model D) and $\alpha_{\rm p}=1.55$ (Model E), respectively.

As shown above, the broadband spectra of the Crab nebula can be fitted with the leptonic model, in which the $\gamma$-ray emission is dominated by the IC process. In this picture, we also constrain the contribution to $\gamma$-rays from the hadronic process. Fitting the radio to $\gamma$-ray data via the MCMC method (namely, Models B and C), we found the energy in protons confined in the Crab should be less than $1.9\times 10^{47}$ erg at $2\sigma$ level. It is corresponding to the energy conversion efficiency of $\eta_{\rm p}\sim 0.5\ (n_{\rm t}/10\ {\rm cm}^{-3})^{-1}$ per cent, considering the total spin-down energy of $\sim 4\times 10^{49}$ erg released by Crab during its whole life and the long lifetime of the energetic protons in the Crab’s environment. This fraction is far lower than the required value in the hadronic dominated models for some PWNe in which the spin-down energy are dominantly carried by ions, e.g., Vela X (Horns et al., 2006) and G54.1+0.3 (Li et al., 2010). Of course, the energy carried by protons is dependent to the target gas density $n_{\rm t}$, but it can be inferred that the hadronic emissions play a minor part in explaining the $\gamma$-ray emission from the Crab in Models B and C.

On the other hand, in Models D and E in which only $\gamma$-ray data were used alone as an extreme case, we find that the $\sim 100$ TeV $\gamma$-rays are dominated by the hadronic process. This requires that the energy stored in protons needs to be 1.4 – $2.6\times 10^{48}$ erg, indicating that protons steal a substantial spin-down energy with a fraction of $\eta_{\rm p}\sim 7\ (n_{\rm t}/10\ {\rm cm}^{-3})^{-1}$ per cent. This value is close to the upper limits of $<11\ (n_{\rm t}/10\ {\rm cm}^{-3})^{-1}$ per cent derived by Di Palma et al. (2017) considering the fact that IceCube has not detected any neutrinos from the Crab²²2It was scaled to the target density of $10\ {\rm cm}^{-3}$ based on the material in PWN with the mass of $10\ M_{\odot}$ adopted in the original paper.. We also calculate the neutrino spectrum and find that the neutrino fluxes around 30 TeV are easy to exceed the IceCube’s differential sensitivity as the increase of $\eta_{\rm p}\cdot n_{\rm t}$. Our predicted neutrino fluxes at 30 TeV for $\alpha_{\rm p}=2.0$ are displayed in Figure 6, which also shows the dependence of the neutrino flux on the product of $\eta_{\rm p}\cdot n_{\rm t}$. It should be noted that our results are obtained on the precondition of $E_{\rm c,p}=3$ PeV which approximately equals leptons’ maximum energy. Unlike the leptons, the energy losses are negligible for protons at current age. So their maximum energy is likely much larger than 3 PeV. Then, in order to avoid over-estimating the flux at $\sim 100$ TeV, one should decrease the energy in protons. Thus the conversion fraction $\eta_{\rm p}\sim 7$ per cent should be treated as an upper limits.

Recently, Xin et al. (2019) performed an analysis of the Fermi data toward the TeV source VER J2227+608 associated with a SNR-PSR complex. Considering the fact that the GeV-TeV $\gamma$-ray spectra have a hard spectral index ($1.90\pm 0.04$) and no high-energy cutoff with energy upto several dozens of TeV, they suggested that the $\gamma$-ray emission originate from the energetic protons accelerated by a PWN not a SNR. Based on the hadronic model, PSR J2229+6114 that is in charge of this PWN needs to convert its spin-down energy to protons with a amount of $6.0\times 10^{47}(n_{\rm t}/10\ {\rm cm}^{-3})^{-1}$ erg during its life $\sim 10^{4}$ yr. Considering the low spin-down luminosity of $2.2\times 10^{37}\ {\rm erg\ s^{-1}}$ of PSR J2229+6114 (Halpern et al., 2001), the corresponding energy conversion efficiency of protons is $\eta_{\rm p}\approx 10$. This value is consistent with our results given by Models D and E, or by Models B and C if the target density in Crab is down to $\sim 0.5\ {\rm cm}^{-3}$.

In addition, for Models D and E the population of leptons determined by the $\gamma$-ray data also can produce radio to soft-$\gamma$-ray emission via synchrotron process. In order to match the radio data, we obtain $B_{\rm PWN}\approx 63\ \mu{\rm G}$ (see red dashed lines in the left panel of Figure 5). As can be seen, this population underestimates the fluxes in the IR-to-soft-gray bands. It implies that another population of leptons are needed to fully explain the data, which is consistent with the two-component leptonic models (Zhu et al., 2015; Luo et al., 2020).

In conclusion, the broadband spectra of Crab can be explained by the leptonic model, giving a constraint on the energy conversion fraction of protons $\eta_{\rm p}<0.5$ per cent. With this small fraction, however, it can not fully rule out that PWNe are PeVatrons of protons based on the current data.

X.Z. is indebted to Siming Liu for the helpful discussions. We thank the support of National Key R&D Program of China under nos. 2018YFA0404204 and 2017YFA0402600 and NSFC grants under nos. U1931204, 11803011, 11851305, 11773014, 11633007, 11873065, 11673041, and 11533007.

No new data were generated or analysed in support of this research.