Particle acceleration model for the broadband baseline spectrum of the Crab nebula

TL;DR

This paper presents a one-zone model for the Crab nebula's broadband spectrum, successfully explaining observed synchrotron and inverse Compton emissions with a single electron distribution, offering a new interpretation of spectral deviations.

Contribution

It introduces a simple, analytical model using a log-parabola electron injection spectrum that fits the entire broadband spectrum of the Crab nebula, contrasting with more complex multi-power-law models.

Findings

Model accurately reproduces the Crab nebula's observed spectrum.

Single log-parabola electron distribution explains spectral features.

Provides a framework applicable to other astrophysical outflows.

Abstract

We develop a simple one-zone model of the steady-state Crab nebula spectrum encompassing both the radio/soft $X$-ray and the GeV/multi-TeV observations. By solving the transport equation for GeV-TeV electrons injected at the wind termination shock as a log-parabola momentum distribution and evolved via energy losses, we determine analytically the resulting differential energy spectrum of photons. We find an impressive agreement with the observed spectrum of synchrotron emission, and the synchrotron self-Compton component reproduces the previously unexplained broad $200$-GeV peak that matches the Fermi/LAT data beyond $1$ GeV with the MAGIC data. We determine the parameters of the single log-parabola electron injection distribution, in contrast with multiple broken power-law electron spectra proposed in the literature. The resulting photon differential spectrum provides a natural interpretation of the deviation from power-law customarily fit with empirical multiple broken power-laws. Our model can be applied to the radio-to-multi-TeV spectrum of a variety of astrophysical outflows, including pulsar wind nebulae and supernova remnants, as well as to interplanetary shocks.

Figures

Equations

$$Q(>\gamma )=N_0(\gamma /{\gamma }_0{)}^{-[s-1+r\log (\gamma /{\gamma }_0)]}$$

$$\frac{\partial \bar{N}(\gamma ,t)}{\partial t}-\frac{\partial }{\partial \gamma }(b{\gamma }^2\bar{N}(\gamma ,t))=Q(\gamma ,t),$$

$$Q(\gamma )=\frac{N_0}{{\gamma }_0}{\left(\frac{\gamma }{{\gamma }_0}\right)}^{-s-r\log \left(\frac{\gamma }{{\gamma }_0}\right)}\left|s-1+2r\log \left(\frac{\gamma }{{\gamma }_0}\right)\right|,$$

$$Q(\gamma ,t)=Q(\gamma )(1+a)\frac{t^a}{{\tau }_0^{1+a}}\Theta \left({\tau }_0-t\right)\Theta \left(t\right),$$

$$N(\gamma )$$

$$=\frac{1+a}{b{\gamma }^2{\tau }_0^{1+a}}{\int }_{\gamma }^{{\gamma }_{\max }}d{\gamma }^{\prime }Q({\gamma }^{\prime }){\left[{\tau }_0+\frac{\gamma -{\gamma }^{\prime }}{b\gamma {\gamma }^{\prime }}\right]}^a$$

$$=(1+a)N_0\frac{{\tau }_{\mathrm{syn}}}{{\tau }_0}\frac{1}{\gamma }{\int }_1^{x_{\max }}dx{\left[1+\frac{(1-x){\gamma }_0{\tau }_{\mathrm{syn}}}{\gamma x{\tau }_0}\right]}^a$$

$$\times {\left(\frac{\gamma x}{{\gamma }_0}\right)}^{-s-r\log \left(\frac{\gamma x}{{\gamma }_0}\right)}\left|s-1+2r\log \left(\frac{\gamma x}{{\gamma }_0}\right)\right|,$$

$${\gamma }_{\max }=\{\begin{array}{ll}\frac{\gamma }{1-b\gamma {\tau }_0}& \text{mbox}if\gamma <\frac{1}{b{\tau }_0}\\ \infty & \text{mbox}otherwise\end{array}$$

$$x_{\max }=\{\begin{array}{ll}\frac{1}{1-b\gamma {\tau }_0}& \text{mbox}if\gamma <\frac{1}{b{\tau }_0}\\ \infty & \text{mbox}otherwise\end{array}.$$

$$\nu F_{\nu }^{\mathrm{syn}}(ϵ)=\frac{{\sigma }_{T}c{B}^2}{12\pi 4\pi d^2}{\left(\frac{ϵ}{{ϵ}_0}\right)}^{3/2}N\left(\sqrt{\frac{ϵ}{{ϵ}_0}}\right),$$

$$\nu F_{\nu }^{\mathrm{syn}}(ϵ)=$$

$$\times {\int }_1^{x_{\max }}dx{\left[1+\sqrt{\frac{{ϵ}_s}{ϵ}}\frac{{\tau }_{\mathrm{syn}}}{{\tau }_0}\frac{1-x}{x}\right]}^a$$

$$\times {\left(x\sqrt{\frac{ϵ}{{ϵ}_s}}\right)}^{-s-r2\log \left(x\sqrt{\frac{ϵ}{{ϵ}_s}}\right)}\left|s-1+2r\log \left(x\sqrt{\frac{ϵ}{{ϵ}_s}}\right)\right|,$$

$$x_{\max }=\{\begin{array}{ll}{\left(1-\sqrt{\frac{ϵ}{{ϵ}_s}}\frac{{\tau }_0}{{\tau }_{\mathrm{syn}}}\right)}^{-1}& \text{mbox}ifϵ<{ϵ}_s{\left(\frac{{\tau }_{\mathrm{syn}}}{{\tau }_0}\right)}^2\\ \infty & \text{mbox}otherwise\end{array}.$$

$$n_s(ϵ)=\frac{4\pi }{chϵ}{J}_{\nu }(ϵ);J_{\nu }(ϵ)=\frac{3R}{4}{j}_{\nu }(ϵ)$$

$$\Rightarrow n_s(ϵ)=\frac{3\pi R}{chϵ}{j}_{\nu }(ϵ),$$

$$\nu F_{\nu }^{\mathrm{syn}}(ϵ)=\frac{4\pi }{3}\frac{R^3}{d^2}\nu j_{\nu }(ϵ)$$

$$n_s(ϵ)$$

$$=\frac{3{\sigma }_T{B}^2}{64{\pi }^2{R}^2}\frac{(1+a)N_0{\tau }_{\mathrm{syn}}}{{\tau }_0{ϵ}_s^2}{\gamma }_0^2{\left(\frac{ϵ}{{ϵ}_s}\right)}^{-1}$$

$$\times {\int }_1^{x_{\max }}dx{\left[1+\sqrt{\frac{{ϵ}_s}{ϵ}}\frac{{\tau }_{\mathrm{syn}}}{{\tau }_0}\frac{1-x}{x}\right]}^a{\left(x\sqrt{\frac{ϵ}{{ϵ}_s}}\right)}^{-s-r\log \left(x\sqrt{\frac{ϵ}{{ϵ}_s}}\right)}$$

$$\times \left|s-1+2r\log \left(x\sqrt{\frac{ϵ}{{ϵ}_s}}\right)\right|.$$

$$\frac{d{N}_{\gamma ,e}}{dtd{ϵ}_1}=\frac{3{\sigma }_{T}c}{4{\gamma }^2}\frac{n_{\mathrm{tot}}(ϵ)}{ϵ}dϵG(q,\bar{\Gamma }),$$

$$G(q,\bar{\Gamma })=2q\ln q+(1+2q)(1-q)+\frac{{\bar{\Gamma }}^2{q}^2(1-q)}{2(1+\bar{\Gamma }q)}$$

$$\bar{\Gamma }=\frac{4ϵ\gamma }{m_e{c}^2},q=\frac{{ϵ}_1}{\bar{\Gamma }(m_e\gamma c^2-{ϵ}_1)}.$$

$${\Pi }_{\mathrm{IC}}(\gamma ,{ϵ}_1)$$

$$\nu F_{\mathrm{IC},\nu }^{\mathrm{tot}}({ϵ}_1)$$

$$={ϵ}_1^2\frac{3c{\sigma }_T}{16\pi d^2}\int dϵ\frac{n_{\mathrm{tot}}(ϵ)}{ϵ}\int d\gamma \frac{N(\gamma )}{{\gamma }^2}G(q,\bar{\Gamma })$$

$$\nu F_{\mathrm{IC},\nu }^{\mathrm{ssc}}({ϵ}_1)={ϵ}_1^2\frac{27{\sigma }_T}{64\pi R^2}\int dϵ\frac{\nu F_{\nu }^{\mathrm{syn}}(ϵ)}{{ϵ}^3}\int d\gamma \frac{N(\gamma )}{{\gamma }^2}G(q,\bar{\Gamma }).$$

$$\frac{d\mathcal{N}}{d\epsilon dAdt}(\epsilon )\propto (\epsilon /{\epsilon }_0^M{)}^{-\alpha +\beta \cdot \log (\epsilon /{\epsilon }_0^M)}$$

Full text

††pubyear: 2017††pagerange: Particle acceleration model for the broadband baseline spectrum of the Crab nebula–Particle acceleration model for the broadband baseline spectrum of the Crab nebula

1 Introduction

Current theoretical models of the overall spectral energy distribution (SED) of the Crab Nebula fail in reproducing both the peak in the very-high energy (hereafter VHE) range and the radio/$X$-ray data with a single and physically justified population of accelerated particles. The photon differential spectrum from ground-based observatories, i.e., HEGRA Aharonian et al. (2004), MAGIC (Aleksić et al., 2015), HESS II (Holler et al., 2015) and VERITAS (Meagher, K. for the VERITAS Collaboration, 2015) is consistent with a log-parabola distribution in the range $\sim 50$ GeV up to $100$ TeV. However, the joint photon differential spectrum matching the VHE observations of MAGIC ($0.05-30$ TeV) with the Fermi/LAT spectrum in the range $1-200$ GeV exhibits a broad and flat peak that is fit (Aleksić et al., 2015) by a modified log-parabola with an additional ad hoc free-parameter. Such a VHE spectrum cannot originate from an electron-positron spectrum evolved from the standard power-law distribution of shock-accelerated energetic particles; however, no attempt has been reported so far to reconstruct the electron-positron source spectrum.

The VHE emission of the Crab nebula is usually modelled with Inverse Compton (hereafter IC) radiation from a non-thermal population of electron-positron pairs via scattering of three possible photon targets: synchrotron radiation emitted by the same electron-positron population within the nebula (Synchrotron self-Compton) in the spectral range from the optical to $X$ rays, thermal far-infrared radiation (typically associated with dust) or cosmic microwave background (CMB) radiation. A variety of radiation processes proposed to explain the broadband spectrum of the Crab pulsar wind nebula (PWN) seems to fail in matching the VHE emission with the emission at smaller wavelengths from radio to $X$-rays, as discussed in Aleksić et al. (2015): the MHD spherical flow model of the downstream pulsar wind structure (Kennel & Coroniti, 1984; Meyer et al., 2010) fits well only the synchrotron part, not the VHE part of the spectrum; a simplified model with uniform magnetic field throughout the nebula and two electron populations (relic of the pulsar wind, that cooled to produce synchrotron radio, and wind) fails to reproduce the $200$ GeV peak; a time-dependent model solving the diffusion-loss equation that includes synchrotron, IC and bremsstrahlung losses fails to reproduce the size of the nebula inferred at smaller wavelengths. Thus, the determination of the energetic particles distribution leading to the observed $\gamma$-ray spectrum and compatible with the low-energy spectrum (radio and soft $X$-rays) is an open problem.

Within the PWN standard scenario, the magnetized wind from the pulsar advects electron-positron pairs presumably (Arons & Tavani, 1994) accelerated to a power-law distribution by the wind termination shock (hereafter TS) via a diffusive shock-acceleration á la Fermi. An alternative scenario based on shock-driven reconnection in the striped wind was proposed by Pétri & Lyubarsky (2007). In this paper we adopt the former scenario and show that a physically motivated energetic particles distribution can explain the broadband observations.

Energetic particle spectra measured in-situ at interplanetary shocks exhibit in several cases deviation from a single power-law. A multispacecraft analysis of the energetic ions spectra (Mewaldt et al., 2012) reported that all the $16$ Ground Level Events during solar cycle $23$ can be fit by broken power laws or by the Band function (Band et al., 1993) in the range $0.1$ to $500-700$ MeV with break energies ranging from $\sim 2$ to $\sim 46$ MeV. This result suggests that no-power-law spectra over $2$ or $3$ energy decades are perhaps a commonplace for shocks.

In this paper we investigate the microscopic origin of the broadband photon spectrum of the Crab nebula. We show that a single population of electrons evolved via energy losses from a log-parabola injection spectrum leads to: 1) quasi-log-parabola photon differential spectrum produced via synchrotron self-Compton scattering and IC scattering over CMB in agreement with the overall $\gamma-$ray observations between $1$ GeV and $100$ TeV; 2) low-energy synchrotron spectrum in the range $(10^{-5}-10^{4})$ eV in agreement with radio/$X$-rays observations (also reported in Aleksić et al. (2015)). The injection of the energetic electrons is found to decrease at a very slow rate over the entire lifetime of the nebula ($\sim 10^{3}$ years). We provide a probabilistic microscopic interpretation of the underlying process of particle acceleration at the Crab nebula TS.

This paper is organised as follows. In Sect. 2 we outline the basic assumptions of the microscopic model for the energetic particles. In Sect. 3 the model for the synchrotron radiation and the IC scattering emission is described in detail. In Sect. 4 we compare the predicted photon spectrum with the broadband observations, contrast the log-parabola electron injection distribution with asymptotic power-laws. In Sect. 5 we derive the relation between the parameters of the electrons distribution and the observed photon spectrum and estimate the diffusion-advection scale upstream of the TS. Sect 6 summarises our findings.

2 Outline of the model

For a population of energetic particles at an infinitely planar collisionless shock wave, the process of acceleration is self-sustained if individual particles, advected downstream after shock-crossing, can contrast the flow advection and return to the shock to undergo further acceleration. The reservoir of far downstream particles can be phenomenologically treated by introducing the probability of return, $P(v)$, i.e., the probability that an individual particle with speed $v$ in the frame of the downstream fluid can return to the shock111By definition, $P(v)$ is the fraction of the particle flux that crosses the shock from upstream to downstream and return to the shock over the total flux of particles crossing from upstream to downstream.. In contrast, for the upstream fluid, all particles penetrating into the upstream are overtaken by the shock, unless turbulence is unable to confine them and they escape to a far-away observer; the escape can be mathematically modelled by requiring a free escape boundary whose location upstream is generally treated as a free parameter.

Numerical simulations have extensively used $P(v)$ by implementing the scattering off the magnetic fluctuations with a Monte-Carlo method (Ellison et al., 1996) or by synthesizing the magnetic turbulence with a prescribed power spectrum (Giacalone, 2005; Fraschetti & Giacalone, 2015). For flows crossed by a non-relativistic shock wave (e.g., interplanetary or supernova remnant shock), if $U_{2}$ is the speed of the downstream plasma in the shock rest frame, such probability reads $P(v)=[(v-U_{2})/(v+U_{2})]^{2}$, increasing as a function of $v$ and independent of energy for $v\rightarrow c$ (Jones & Ellison, 1991, e.g.). The definition of $P(v)$ requires the isotropy of the pitch-angle distribution in the local plasma frame, or in other terms that $v$ is so much greater than the flow speed $U$ that second order terms in $U/v$ can be disregarded. We note that, if the condition $v\gg U$ is not satisfied, i.e., at supra-thermal energy, particles are not fast enough to isotropize their pitch-angle distribution in the downstream plasma frame close to the shock. 222An example is the emergence of spikes in the intensity profiles of supra-thermal particles measured across interplanetary shocks (Lario et al., 2003), explained as result of ion reflection by the magnetic barrier at quasi-perpendicular shocks (Gieseler et al., 1999) or pile-up of particles streaming along turbulent field lines crossing the shock surface multiple times (Fraschetti & Giacalone, 2015, and references therein). As a result, the supra-thermal particle spectrum downstream close to the shock is harder than the one-dimensional prediction of the diffusive shock acceleration model. The same expression of $P(v)$ holds for ultra-relativistic shocks (Peacock, 1981) and was extended to any relativistic shock speed (Summerlin & Baring, 2012), provided that the particle distribution is isotropic in the downstream plasma frame.

The momentum gain relates differently to the flow speed between ultra-relativistic flows and non-relativistic flows. For an isotropically distributed incoming population of relativistic electrons crossing a quasi-perpendicular ultra-relativistic shock, such as a pulsar wind TS at large-scale, the momentum gain upon the first cycle of crossing and re-crossing the shock is $p_{f}/p_{i}\simeq\Gamma^{2}/2$, where the initial ($p_{i}$) and final ($p_{f}$) momentum are calculated in the upstream fluid frame and $\Gamma$ is the Lorentz factor of the upstream fluid in the shock frame. For the subsequent shock-crossings, $p_{f}/p_{i}$ is limited by geometrical constraints to smaller values (Gallant & Achterberg, 1999). For a non-relativistic fluid, the momentum gain upon first cycle of shock crossing is smaller than for ultra-relativistic flows. Upon any cycle of shock crossing, including the first one, $\Delta p=p_{f}-p_{i}$ can be written, by averaging over an isotropic distribution, as $\langle\Delta p/p\rangle=4/3\,\Delta U/v$, where $\Delta U$ is the change of flow speed across the shock, and $v$ is the particle speed ($v$ and $p$ are computed in the local downstream plasma frame). This expression for $\Delta p/p$ is valid only if higher-order terms in $U/v$ are neglected.

Several effects can hamper the energy gain of charged particles interacting with shock waves. In the case of a large-scale quasi-perpendicular magnetic obliquity, such as the Crab TS, the acceleration proceeds as long as the transport in the direction perpendicular to the average field due either to field line meandering (Jokipii, 1966) or to the motion across the local perturbed field lines (Fraschetti & Jokipii, 2011) allows particles to re-cross the shock. Ultra-relativistic shocks, such as presumably a pulsar wind TS, have predominantly perpendicular obliquity leading to the high compression of the magnetic component perpendicular to the average direction of the shock motion; moreover, the particle speed and the shock speed are very close, in contrast with non-relativistic shocks, so that the isotropy in the downstream medium can be attained only far from the shock. Therefore, particles are more likely to be lost in the downstream of a shock at ultra-relativistic than at non-relativistic flow speeds.

In this paper we propose that the Crab wind TS accelerates to the multi-TeV range a single population of electrons that produces the observed baseline differential spectrum between the radio and the VHE band. We consider here the realistic case that the shock cannot be assumed infinitely planar at the typical spatial scale of particle motion (i.e., gyroradius $r_{g}$) of the highest-energy particles ($\sim$ TeV) as a result of inherent geometrical effects limiting the particle acceleration, for instance the corrugation at various scales of the shock surface or its large-scale curvature; the gyroscale is usually neglected as typically much smaller than the other scales. As the particle speed grows, such geometrical effects modify the probability that a particle remains in proximity of the shock introducing a mild energy-dependence, that accounts for the highest-energy particles with $r_{g}$ relatively large to leak out of the system, thereby leading to a softening of the photon spectrum.

Throughout the paper we shall assume for simplicity that the energy gain per cycle is a constant fraction of the particle energy, i.e. the energy gain is linear in energy. Then, we consider the probability ${\cal P}(\gamma)$, where $\gamma$ is the particle Lorentz factor in the downstream plasma frame where the observed radiation is emitted, as the probability of remaining in the acceleration region, both upstream or downstream, that incorporates the probability of return to the shock from downstream, the geometrical factors limiting the acceleration listed above and the particle escape upstream. For such VHE electrons we use the microscopic ansatz that ${\cal P}(\gamma)$ depends on $\gamma$ as ${\cal P}(\gamma)=g/\gamma^{q}$, where $g$ is a normalisation factor and $0<q<1$ a parameter. Such an expression does not assume isotropy in the local plasma frame. The anisotropy expected at ultra-relativistic flows in the downstream local frame behind and very close to the shock, as compared to the gyroscale of $\simeq 1$ TeV electrons downstream, does not contradict the ansatz ${\cal P}(\gamma)=g/\gamma^{q}$. Test-particle simulations show that isotropy in the local downstream plasma frame is recovered (Summerlin & Baring, 2012).

It has been shown that such ansatz leads to a log-parabola spectrum of energetic electrons (Massaro et al., 2004) that produce a quasi-log-parabola synchrotron spectra, modulo a logarithmic factor, in good agreement with the steady $X$-ray spectra observed from blazars. Moreover, a significant fraction of photon spectra of extended sources listed in the 3FGL catalog are better represented by a log-parabola functional form than with a single power-law (Acero et al., 2015).

We define $Q(>\gamma)$ as the number density of energetic electrons accelerated to a total energy higher than $m_{e}\gamma c^{2}$, where $m_{e}$ is the electron mass, and $R=\gamma_{i}/\gamma_{i-1}$ as the ensemble-averaged increment in energy at the $i^{th}$ shock-crossing, given by the ratio of the $i^{th}$, $\gamma_{i}$, to the previous, $\gamma_{i-1}$, particle Lorentz factor. The form ${\cal P}(\gamma)=g/\gamma^{q}$ yields333If the probability of return ${P}(v)$ is independent of the particle energy, the power-law differential spectrum $\sim N_{0}\gamma^{-p}$ is re-obtained, where the index $p=-{\rm ln}P(v)/{\rm ln}R$ can be re-written in terms of the shock density compression ratio in the test-particle regime. If the backreaction of energetic particles grows as to modify the upstream turbulence, the spectrum will be concave, instead of a simple power-law, as particles with different energy experience different shock compression. (Massaro et al., 2004)

[TABLE]

where $N_{0}$ is a normalization constant, $\gamma_{0}$ is the electron injection factor and the spectral parameters $s$ and $r$ are related to the microscopic statistical parameters $g,q,R$ as follows: $s=-{\rm log}(g/\gamma_{0}^{q})/{\rm log}(R)-(q-2)/2$, $r=q/2{\rm log}(R)$. It is noteworthy clarifying that $\gamma_{0}$ is the injection threshold into the process of acceleration limited by geometrical factors in the VHE range, and not the injection energy from the supra-thermal pool of particles that undergo diffusive shock acceleration.

3 Theoretical broadband flux from log-parabola electron distribution

We assume that the non-thermal population of energetic electrons/positrons pairs is produced with a log-parabola injection spectrum in Eq. 1, accounting for a slowly decreasing rate, at the pulsar wind TS up to energies $\sim 100$ TeV, evolves under the influence of energy losses, and feeds the observed baseline photon emission from the radio up to the multi-TeV range. We calculate analytically in this section the total synchrotron and IC fluxes including energy losses during the entire nebula life-time. All quantities are calculated in the reference frame of the fluid downstream of the TS, i.e., the nebula, where synchrotron photons are emitted by the multi-TeV electrons and, in part, upscattered to VHE range.

3.1 Time dependence of the electron spectrum

In this sub-section we provide an analytic expression for the spectrum of electrons that results from injection at the TS with the spectrum given by Eq.1. The GeV-TeV electrons are assumed to be injected in a spatially uniform fashion at the source, i.e., Crab PWN TS, so that the differential electron number spectrum $\bar{N}(\gamma,t)$ is independent of location. Therefore, $\bar{N}(\gamma,t)$ satisfies the simplified transport equation

[TABLE]

with an energy loss rate $-b\,\gamma^{2}$ and a differential injection rate of the electrons $Q(\gamma,t)$ at time $t$.

From Eq. 1, the differential electron production rate $Q(\gamma)=dQ(>\gamma)/d\gamma$ is given by

[TABLE]

that is a log-parabola multiplied by a pre-factor logarithmically dependent on $\gamma/\gamma_{0}$. The injection term in Eq. 2 is taken to have the form

[TABLE]

where $\tau_{0}\simeq 10^{3}$ years denotes the current age of the remnant, and $a<0$ is a free parameter modulating the electron injection rate; $\Theta(x)$ indicates the Heaviside function.

Defining the synchrotron-loss time scale as $\tau_{\mathrm{syn}}=1/b\gamma_{0}$, we find for the differential electron number spectrum at the present time, $\bar{N}(\gamma,\tau_{0})\equiv N(\gamma)$

[TABLE]

where $\delta(\cdot)$ in the first line of Eq. 5 is the Dirac-delta function and we have used $x=\gamma^{\prime}/\gamma$. The upper limit of integration is

[TABLE]

or

[TABLE]

3.2 Synchrotron radiation

The synchrotron power emitted by a single electron, averaged over an isotropic electron distribution, is given in the local plasma frame by $P_{\mathrm{syn}}=(\sigma_{T}c/6\pi)\gamma^{2}B^{2}$, where $\sigma_{T}=({8\pi/3})({e^{4}/m_{e}^{2}c^{4}})$ is the Thomson cross-section, and $B$ is the external magnetic field assumed to be uniformly embedded in the nebula. The total synchrotron flux at Earth from a source at distance $d$, namely $\nu F_{\nu}^{\mathrm{syn}}$, is found by folding $P_{\mathrm{syn}}$ with the differential electron spectrum in Eq. 5: $\nu F_{\nu}^{\mathrm{syn}}={1\over 4\pi d^{2}}\int d\gamma P_{\mathrm{syn}}N(\gamma)$. We use the monochromatic approximation, i.e., the electron power is concentrated around the characteristic synchrotron frequency $\nu_{s}=0.29(3e\gamma^{2}B)/(4\pi m_{e}c)=\gamma^{2}\nu_{s}^{0}$. Thus, the total flux can be recast as $\nu F_{\nu}^{\mathrm{syn}}\simeq{1\over 4\pi d^{2}}\int d\gamma P_{\mathrm{syn}}N(\gamma)|_{\nu=\nu_{s}}$ which leads to

[TABLE]

where we have used ${\epsilon_{0}}=h\nu_{s}^{0}$ with the Planck constant, $h$. Inserting the expression for $N(\gamma)$ in Eq. 5 we find for the total synchrotron flux

[TABLE]

where we have used $\epsilon_{s}=\gamma_{0}^{2}\,\epsilon_{0}$. To retain consistency in notation we write

[TABLE]

Equation 9 contains a log-parabola distribution in $x\,\sqrt{{\epsilon}/{\epsilon_{s}}}$, first factor in the third line of Eq. 9, inside the $x$-integrand that takes into account the energy losses due to synchrotron radiation during the entire nebula lifetime. The factor in $|\cdot|$ in the third line of Eq. 9 is usually neglected in the fit of blazars $X$-ray spectra, as only logarithmically increasing; such factor controls the width of the spectral peaks at $\sim 1$ eV and $\sim 200$ GeV. The flux given in Eq. 9 will be compared in the next section with radio and $X$-ray data for the Crab nebula baseline spectrum.

3.3 Synchrotron self-Compton and inverse Compton scattering from CMB radiation

In this sub-section we analytically calculate the predicted IC flux due to scattering off the synchrotron photons determined in Sect. 3.2 and off CMB radiation within the nebula. The Lorentz-invariant differential density of synchrotron photons inside a spherical source of radius $R$, $n(\epsilon)=dn/(dV\,d\epsilon)$, is for optically thin systems related to the mean intensity, $J_{\nu}$, and that to the emission coefficient, $j_{\nu}$, as

[TABLE]

where $\epsilon$ is the target photon energy in the local plasma frame. The total synchrotron flux at Earth depends on the emission coefficient as

[TABLE]

Thus, the differential photon density of synchrotron photons is

[TABLE]

where in the last step we have used Eq. 9. The total differential photon density includes $n_{s}(\epsilon)$ and the differential CMB photon density $n_{\mathrm{CMB}}$: $n_{\rm tot}(\epsilon)=n_{s}(\epsilon)+n_{\mathrm{CMB}}$.

For an ambient photon population isotropic in the local plasma frame with a differential density $n(\epsilon)$, the spectrum of the scattered photons with initial energy between $\epsilon$ and $\epsilon+d\epsilon$ and final energy $\epsilon_{1}$ per single electron is (Blumenthal & Gould, 1970)

[TABLE]

with

[TABLE]

and

[TABLE]

where Eqs (14, 15) are derived under the only restriction $\gamma\gg 1$. The calculation of the IC photon number emitted per unit time by a single electron in the local plasma frame over the target photon distribution $n(\epsilon)$, namely $\Pi_{\mathrm{IC}}(\gamma,\epsilon_{1})$, leads to

[TABLE]

We note that in the $\epsilon$-integral of Eq.17 we calculate the contribution of target synchrotron photons at all frequencies with no monochromatic approximation. Since the synchrotron photon energies extend over a wide range ($10^{-5}$ eV up to $100$ keV), Eq. 17 entails scattering at all regimes between Thomson ($\bar{\Gamma}\ll 1$) and Klein-Nishima ($\bar{\Gamma}\gg 1$). The total IC photon number emitted per unit time by the electron population is found by folding Eq. 17 with differential distribution $N(\gamma)$ in Eq.5 $\Pi_{\mathrm{IC}}(\epsilon_{1})=\int{d\gamma}\Pi_{\mathrm{IC}}(\gamma,\epsilon_{1})N(\gamma)$.

The observed IC energy flux is related to $\Pi_{\mathrm{IC}}(\epsilon_{1})$ by

[TABLE]

where the term due to synchrotron self-Compton can be re-written as

[TABLE]

We emphasize that the energy fluxes in the radio/$X$-ray band (Eq. 9) and in the VHE band (Eq. 18) use the same differential energetic electron distribution (Eq. 5). Those energy fluxes will be compared in the next section with the baseline broadband spectrum of the Crab nebula.

4 Comparison with broadband observations

Figure 1 compares the theoretical differential spectrum of synchrotron photons given by Eq. 9 with the observed radio and soft $X$-ray spectra ($10^{-5}$ eV up to $\sim 10$ MeV), and likewise the theoretical IC radiation as in Eq. 18 with the joint Fermi/LAT and MAGIC spectra. The best-fit involves $4$ free parameters for the overall spectrum: $s=1.75$, $r=0.08$, $N_{0}=2.6\times 10^{49}$, and $a=-0.05$ (see also Fig. 5) for an energy scale $\gamma_{0}=2\times 10^{4}$. We have used a magnetic-field strength $B$ constant and uniform throughout the nebula of the order of the typical value, i.e., $B=140\,\mu$G (cfr. the best-fit value for the constant $B$-field model in Aleksić et al. (2015), $B=143\,\mu$G). We notice that for this magnetic-field strength the synchrotron losses dominate for $\gamma^{\prime}\gtrsim 10^{2}\gamma_{0}=2\times 10^{6}$, leading to $\tau_{\mathrm{syn}}\lesssim\tau_{0}$; at smaller electron energy the loss time is larger than the age of the Nebula itself. As a consequence the synchrotron flux peaks at a photon energy $\epsilon=\gamma^{\prime 2}\epsilon_{0}=2.7$ eV and diminishes at higher energy (see Fig. 1). Finally, we have used a size of the nebula $R=1$ pc and a distance from the observer $d=1.9$ kpc.

The VHE region of the spectrum is zoomed-in in Fig.2 and compared with reconstructed spectra obtained by a number of VHE ground-based observatories. The overlaid log-parabola fit of the MAGIC data only (red in Fig. 2) reported in (Aleksić et al., 2015) has the form

[TABLE]

with index $\alpha=2.47\pm 0.01$, curvature $\beta=-0.24\pm 0.01$ and scale $\varepsilon_{0}^{M}=1$ TeV (note that $\varepsilon_{0}^{M}\simeq m_{e}c^{2}\gamma^{\prime}$).

Figure 2 shows that MAGIC’s empirical log-parabola in Eq. 20 and our solution in Eq. 18 equally well reproduce the MAGIC data down to the IC peak at $\simeq 200$ GeV. However, a pure log-parabola curve fails to explain the lower energy (Fermi/LAT) spectrum down to $\sim 1$ GeV, that Aleksić et al. (2015) proposes to fit with an empirically modified log-parabola relying on an additional fitting parameter. In contrast, we show that the theoretical physically motivated curve derived here is in satisfactory agreement with the Fermi/LAT rising part of the spectrum, the broad $200$ GeV peak and the VHE spectrum. We note that the value $\tau_{0}=1,000$ years used here is a crude approximation and might be smaller due to the time required to accelerate electrons to the TeV range ($~{}50-100$ years). The range $10-10^{3}$ MeV of the spectrum is overestimated by an extrapolation of our model, and we stop plotting the synchrotron spectrum beyond 10 MeV. We note that the standard synchrotron emissivity looses validity in this energy band as one approaches the radiation-reaction limit at $\approx 150\ \mathrm{MeV}$ (Guilbert et al., 1983). A possible explanation is also that the acceleration processes operating in the Crab nebula can not reach electron energies higher than a PeV444PeV electrons radiate at $10\,\mathrm{MeV}$ in a $140$-$\mathrm{\mu G}$ field.. An additional feature not reproduced here is the re-brightening in the range $0.01-0.1$ eV (see Fig. 1) possibly due to dust contribution (Meyer et al., 2010).

Figure 3 shows the overall theoretical spectrum with three different values of $r$ ($0.2,0.08,0$). In particular, the $r=0$ curve, corresponding to a single power-law electron-injection spectrum (see Eq. 1), cannot reproduce the spectral curvature in the synchrotron range. To be noted from the figure is the mismatch of the $r=0$ curve both below and above the cooling break at $2.7$ eV, implying that at least two breaks in the injection spectrum, if not more, are required to reproduce the SED of the Crab nebula. As each break requires 2 free parameters, a combination of multiple broken power-laws energy spectra with ad-hoc energy breaks is an unlikely explanation for the synchrotron $\nu F_{\nu}$ spectrum, regardless of the consistency with the VHE spectrum. The main result shown by Fig. 3 is that a single electron distribution accounts naturally for the curvature of the spectrum across several orders of observed photon energy with no need of energy breaks, provided one allows for log-parabola injection spectra.

Figure 4 shows the steepening of the synchrotron photon spectrum as a result of the steepening of electron spectrum. The slope of the $\nu F_{\nu}$ spectrum in the radio range is expected to be the closest to that resulting from a single-power-law electron distribution as the logarithmic factor in Eq. 9 is very small ($\epsilon\simeq\epsilon_{s}\simeq 0.27$ meV and $x\simeq 1$). Figure 1 indicates a $\nu F_{\nu}$ spectral index of $\sim 0.64$ for the radio band, corresponding to an electron slope $\sim 1.72$, very close to the best-fit value of $s=1.75$. In Fig. 5 we explore the observational consequences of a possible inherent time-dependence of the particle injection by varying the parameter $a$ (see Eq. 4) over the range $-0.5,-0.05,0$. We can conclude that observations are consistent with a constant ($a=0$) or a slowly declining ($a=-0.05$) rate of particle injection. Figure 6 compares the best-fit electron distribution $N(\gamma)$ from Eq. 5 with two asymptotic power laws. The overlaid differential electron production rate $Q(\gamma)$ from Eq. 3 shows the effect of the energy losses.

The parameter $\gamma_{0}$ is poorly constrained, because apart from a logarithmic factor any variation in $\gamma_{0}$ in Eq. 3 can be compensated by adjustments of other parameters ($s$, $r$, or $N_{0}$). Technically, $\gamma_{0}$ may be interpreted the Lorentz factor of electron that return to the upstream region for the first time, and hence for highly relativistic shocks $\gamma_{0}$ can be very large. We note that the radio spectrum is best reproduced if we allow for curvature in the synchrotron spectrum down to $100\,\mathrm{MHz}$, corresponding to an electron Lorentz factor of about $10^{3}$. Our figures demonstrate that $\gamma_{0}=2\cdot 10^{4}$ is also an acceptable choice.

5 Discussion

The model of the broadband Crab nebula spectrum presented in the previous section relies on the ansatz of an energy dependence of the probability ${\cal P}$ that the multi-TeV energetic electrons undergoing acceleration at the wind TS can be confined in the proximity of the shock less efficiently as the particle energy grows. As mentioned above, such an assumption is also supported by the modelling of blazars $X$-ray spectra during the past decade, by the 3FGL catalog of extended astrophysical sources, and by in-situ measurements of interplanetary shocks.

By using a crude assumption, we can derive from the best-fit parameters of the observed photon spectrum an estimate of the microscopic parameter $q$. The electron-positron pair wind, possibly loaded by ions (Hoshino et al., 1992), impinges on the TS with a speed customarily regarded as ultra-relativistic, although the order of magnitude of $\Gamma$ is substantially model-dependent and controversial. The value $\Gamma\simeq 10^{6}$ early proposed by Kennel & Coroniti (1984), with a single electron population injected at $1$ TeV, could not account for radio-to-IR observations. We adopt here $\Gamma\simeq 10^{2}$, consistent with a cold, magnetically dominated upstream wind, outward to the TS, that appears point-like in optical band (Kirk et al., 2009, for a discussion). For such a shock, the relation $r=q/2{\rm log}(R)$ yields $q\simeq 2r{\rm log}(\Gamma^{2}/2)$ for the first shock crossing ($R\simeq\Gamma^{2}/2$) and $q\simeq 2r{\rm log}(2)$ thereafter ($R\simeq 2$). The assumed value $\Gamma\simeq 10^{2}$ leads to $q=0.59$. Such a small value of $q$ has nevertheless a dramatic impact on the dependence of ${\cal P}$ on the particle energy, indicating an extremely efficient depletion of the shock region. Figure 7 depicts ${\cal P}(\gamma;\Gamma)=g/\gamma^{q(\Gamma)}$ where best-fit values are used for $g=\gamma_{0}^{q}\times 10^{-(s+(q-2)/2){\rm log}R}$. For completeness we show in Fig. 7 also $g/\gamma^{q}$ for shock-crossings after the first one although the small values of ${\cal P}$ for first shock-crossing make subsequent crossings very unlikely to take place. For a non-relativistic strong shock (assuming a density compression equal to $4$) the relation $r=q/2{\rm log}(R)$ yields $q=2rU_{1}/v$, where we have used $R=1+\Delta p/p$ and $U_{1}$ is the upstream fluid speed. Since the energy gain per first crossing for non-relativistic shocks is smaller than for ultra-relativistic shocks, the value of $q$ is likely to be typically smaller ($U_{1}/v\ll 1$), i.e. $q<0.1$.

In contrast with the infinitely planar shock case, the probability ${\cal P}$ drops dramatically with the flow $\Gamma$. Monte Carlo test-particle simulations for infinitely planar shocks in Lemoine & Pelletier (2003) found a negligible dependence of the probability of return, $P(v)$, on $\Gamma$ and a mild dependence on $\gamma$. By introducing the ansatz ${\cal P}=g/\gamma^{q}$, we argue that at VHE the non-planar large-scale structure of the shock likely affects the capability of the shock of containing the most energetic particles.

In terms of macroscopic interpretation, we consider the case that electrons undergo scattering in the fluid upstream of the shock (Achterberg et al., 2001; Pelletier et al., 2009) so that the relevant physical length-scale is the scale of confinement of particles upstream or the diffusion-advection scale, defined as $D/c$, where $D$ is the transport coefficient in the average direction of the flow impinging on the shock, that is a radial outflow at the TS to large-scale; $D/c$ has to be larger than the scattering mean free path in order to guarantee particle confinement, hence acceleration. In the downstream fluid the turbulence has to be strong enough to enable the return of the electrons back to the shock. At the wind TS, the large-scale magnetic field is quasi-perpendicular to the flow direction555At non-relativistic quasi-perpendicular shocks, the acceleration is very fast and efficient (Jokipii, 1987) as confirmed by test-particle simulations at fluid-scale (Fraschetti & Giacalone, 2015, and reference therein) and by hybrid simulations at ion kinetic scale (Giacalone, 2005) for interplanetary shocks travelling into pre-existing turbulence. Such result is also confirmed via hybrid simulations for high-Mach number supernova remnant shocks travelling, e.g., in partially ionized fluid (Ohira, 2016).; thus, $D$ is well approximated by $D_{\perp}$ at large-scale. In a quasi-perpendicular configuration we do not expect self-generated waves ahead of the ultra-relativistic shock to enhance the scattering of particles and thus keeping them close to the shock, although streaming of particles could be relevant in regions of the TS with a quasi-parallel magnetic field (Vainio et al., 2004); thus, $D_{\perp}$ is mainly due to pre-existing wind turbulence. We can compare $D_{\perp}/c$ with the typical TS scales. We adopt for the sake of simplicity the Bohm scaling, i.e., namely scattering mean free path equals $r_{g}$ determined by the unperturbed magnetic field. The Bohm scaling generally applies to weak field regime where the power of the magnetic fluctuations along all three space directions is nearly isotropic and corresponds to large-scale diffusion coefficients comparable in the direction parallel and perpendicular to the field. In the Bohm regime $D/c=r_{g}/3$; thus, the scattering mean free path ($\sim r_{g}$) becomes comparable with $D/c$, the scattering becomes inefficient, and particles are likely to leave the system. The Bohm assumption sets an upper limit on $D_{\perp}$: if the turbulence is weaker (Fraschetti & Giacalone, 2012, e.g.), $D_{\perp}$ will be smaller and $D_{\parallel}$ will be larger respectively leading to smaller values of the diffusion-advection scale and larger values of the the scattering mean free path, thereby enhancing the particle leak-out.

We have shown that the synchrotron spectrum is produced by electrons in the GeV - TeV range orbiting around a large-scale magnetic field with strength $B_{0}=140\,\mu$G. The value of $B_{0}$ probes the region behind and close to the shock where cooling length-scale of TeV electrons is small; the field $B_{0}$ might be the result of local amplification due to the inhomogeneities of the pulsar wind advected through the shock (see Mizuno et al. (2014) for relativistic and Fraschetti (2013) for non-relativistic flows) and the average nebular field might be smaller than $B_{0}$. An electron emerging upstream will have an energy in the upstream frame $E^{\prime}=\Gamma E/\sqrt{2}$, where $\Gamma/\sqrt{2}$ is the Lorentz factor of the downstream fluid relative to the upstream frame (Blandford & McKee, 1976); if the magnetic field $B_{0}$ is only the result of shock compression from the jump conditions at a quasi-perpendicular shock front, the upstream field in the upstream frame is $B^{\prime}\sim B_{0}/{\sqrt{r^{2}-1}\Gamma}$, where the “$\sim$” accounts for the local change of magnetic obliquity along the shock surface and $r=3$ is the density compression for high Mach number ultra-relativistic shocks. Thus, by taking $B_{0}=140\,\mu$G, the upstream electron gyroradius in the upstream frame is given by $r^{\prime}_{g}=E^{\prime}/eB^{\prime}\simeq 2\sqrt{2}\Gamma^{2}E/eB_{0}\sqrt{2}\simeq 1.5\times 10^{-1}$ pc for $E\simeq 1$ TeV, where $e$ is the electron charge and $\Gamma\simeq 10^{2}$ (shrinking to $r^{\prime}_{g}\simeq 1.5\times 10^{-3}$ pc at $E=10$ GeV). The resulting diffusion-advection length in the upstream wind frame at $1$ TeV is $D/c=r_{g}^{\prime}/3=5\times 10^{-2}$ pc. Particles in the cold wind upstream are confined within a region much smaller than the radius of the nebula TS ($\sim 0.1$ pc), typically identified with an equatorial $X$-ray ring comprised of a collection of knots (Hester et al., 2002).

6 Conclusions

We have developed a simple analytic one-zone model for the broadband baseline photon spectrum of the Crab nebula. We find an impressive agreement with observation over the entire range $(10^{-5}-10^{14})$ eV, except the MeV-range and an IR re-brightening at $\sim 0.01-0.1$ eV. The energetic electrons, injected at $\sim 10$ GeV, are assumed to be accelerated to a log-parabola distribution at the nebula wind termination shock. The electron spectrum is then obtained as solution to the time-dependent continuity equation accounting for energy losses and a possibly time-dependent injection rate. This model can account for the broad and flat IC peak matching the Fermi/LAT emission beyond $1$ GeV with the VHE emission up to tens of TeV and the synchrotron emission in the radio up to soft $X$-ray range. We show that the resulting synchrotron emission accounts naturally for the curvature of the broadband observed photon spectrum. Our broadband model is consistent with a uniform magnetic field throughout the nebula $B_{0}=140\,\mu$G, in agreement with previous MHD models (Meyer et al., 2010), although this field strength can also result from local downstream amplification. Simple assumptions on the large-scale perpendicular topology of the magnetic field at the nebula termination shock allow to infer the energy dependence of the probability that the energetic electrons remain in the proximity of the shock. The model hereby proposed is applicable to a broad variety of sources (pulsar wind nebulae, supernova remnants) over an unprecedentedly broad energy range.

Acknowledgements

The work of FF was supported, in part, by NASA under Grants NNX13AG10G, NNX15AJ71G. We thank the referee for useful comments. We thank D. Mazin, N. Otte and M. Holler for sharing low-energy, VERITAS and HESS II data respectively. FF acknowledges J. R. Jokipii, J. Giacalone and J. Kóta for continuing discussions and E. J. Summerlin for email exchange. FF thanks for the hospitality the Harvard/Smithsonian Center for Astrophysics where part of this work was performed.