Time-Dependent Electron Acceleration in Pulsar Wind Termination Shocks: Application to the 2007 September Crab Nebula Gamma-Ray Flare

TL;DR

This paper applies a time-dependent electrostatic acceleration model to explain the 2007 Crab Nebula gamma-ray flare, successfully reproducing observed spectra and light curves, and demonstrating that shock-induced reconnection energizes electrons beyond classical limits.

Contribution

The study introduces a novel time-dependent self-similar analytical model for electrostatic acceleration in pulsar wind shocks, explaining both wave and sub-flare features of the Crab Nebula gamma-ray flares.

Findings

Model reproduces the gamma-ray spectrum and light curve of the 2007 flare.

Electrostatic acceleration can energize electrons beyond the synchrotron limit.

Shock-induced magnetic reconnection drives electron acceleration on sub-Larmor timescales.

Abstract

In 2007 September, the Crab Nebula exhibited a bright gamma-ray flare in the GeV energy range that was detected by AGILE. The observed emission at >160 MeV indicates that the radiating electrons had energies above the classical synchrotron radiation-reaction limit, thus presenting a serious challenge to classical models for electron acceleration in astrophysical environments. In this paper, we apply our recently developed time-dependent self-similar analytical model describing electrostatic acceleration in the explosive reconnection region around the pulsar wind termination shock to the 2007 September flare. This event was unique in that it displayed both long-duration "wave" and short-duration "sub-flare" features. The unusual temporal variation makes this flare an especially interesting test for our model. We demonstrate that our model can reproduce the time-dependent gamma-ray spectrum for this event, as well as the associated gamma-ray light curve, obtained by integrating the spectrum for photon energies >100 MeV. This establishes that our time-dependent electrostatic acceleration model can explain both wave and sub-flare transients, which lends further support to the theoretical framework we have developed. We also further examine the validity of the self-similar electric and magnetic field evolution implied by our model. We conclude that strong electrostatic acceleration driven by shock-induced magnetic reconnection is able to power the Crab Nebula gamma-ray flares by energizing the electrons on sub-Larmor timescales.

Tables

Table 1. Table 1: Model Free Parameters

event	$J_0$	$\frac{E_{\ast }}{B_{\ast }}$	$\hat{S}$	$\hat{C}$	$\mu$	$\sigma$	$\alpha$	$\theta$	$t_{\mathrm{adv}}\left(s\right)$	$t_{\ast }\left(\mathrm{MJD}\right)$
$W1$	$1.86\times {10}^{39}$	$0.06$	$3.16\times {10}^{-20}$	$0.2$	${10}^5$	$1.00\times {10}^{10}$	$5.81$	$4.80$	$1.01\times {10}^6$	$54364.0$
$F2$	$3.24\times {10}^{39}$	$0.04$	$4.27\times {10}^{-20}$	$0.2$	${10}^5$	$8.91\times {10}^9$	$6.85$	$6.90$	$1.08\times {10}^5$	$54379.8$

Table 2. Table 2: Derived Parameters

event	${𝒩}_0$	${\overline{\gamma }}_0$	$R_b\left(\mathrm{cm}\right)$	$B_{\mathrm{pk}}\left(\mu G\right)$	$\frac{E_{\mathrm{pk}}}{B_{\mathrm{pk}}}$	$A_{\ast }\left(s^{-1}\right)$	$w_{\ast }$	$w_{\mathrm{pk}}$	$t_{\mathrm{pk}}\left(s\right)$
$W1$	$3.97\times {10}^{38}$	$2.36\times {10}^9$	$1.01\times {10}^{16}$	$472$	$1.10$	$27.3$	$83.3$	$4.56$	$5.52\times {10}^5$
$F2$	$6.68\times {10}^{38}$	$2.74\times {10}^9$	$1.08\times {10}^{15}$	$713$	$1.23$	$16.3$	$125$	$4.07$	$2.38\times {10}^5$

Table 3. Table 3: Asymptotic Energy Values

event	${ℰ}_{\mathrm{inj}}$ (erg)	${ℰ}_{\mathrm{elec}}$ (erg)	${ℰ}_{\mathrm{synch}}$ (erg)	${ℰ}_{\mathrm{esc}}$ (erg)
$W1$	$8.89\times {10}^{41}$	$5.65\times {10}^{41}$	$2.40\times {10}^{41}$	$1.10\times {10}^{42}$
$F2$	$1.17\times {10}^{42}$	$5.45\times {10}^{41}$	$2.12\times {10}^{41}$	$1.52\times {10}^{42}$

Equations

$$\frac{\partial f}{\partial t}=\frac{-1}{p^2}\frac{\partial }{\partial p}\left\{p^2\left[A(t)m_{e}c-S(t)\frac{p^2}{m_{e}c}\right]f\right\}-\frac{f}{t_{\mathrm{esc}}(p,t)},$$

$$N_{\mathrm{tot}}(t)={\int }_0^{\infty }4\pi p^{2}f(p,t)dp.$$

$$t_{\mathrm{esc}}=t_{\mathrm{SRE}}=\frac{p}{C(t)m_{e}c},$$

$$t_{\mathrm{esc}}=t_{\mathrm{ad}}=\frac{R_{\mathrm{b}}}{v_{ds}},$$

$$A(t)=\frac{qE(t)}{m_{e}c}=A_{\ast }h(t),S(t)=\frac{{\sigma }_{{}_{\mathrm{T}}}{B}^2(t)}{6\pi m_{e}c}=S_{\ast }h(t),$$

$$h(t)=\{\begin{array}{ll}{e}^{\alpha t/t_{\mathrm{pk}}},& t\le t_{\mathrm{pk}},\\ e^{\alpha }{e}^{-\theta (\frac{t}{t_{\mathrm{pk}}}-1)},& t\ge t_{\mathrm{pk}},\end{array}$$

$$A_{\ast }=\frac{q{E}_{\ast }}{m_{e}c},S_{\ast }=\frac{{\sigma }_{{}_{\mathrm{T}}}{B}_{\ast }^2}{6\pi m_{e}c}.$$

$$E(t)=E_{\ast }h(t),B(t)=B_{\ast }\sqrt{h(t)},$$

$$E(t)\propto B^2(t).$$

$$v_{\mathrm{rec}}\propto t^{1/2}.$$

$$E\propto v_{\mathrm{rec}}B.$$

$$E\propto t,$$

$$x\equiv \frac{p}{m_{e}c},$$

$$\frac{\partial N}{\partial y}=\frac{{\partial }^2}{\partial x^2}\left(\frac{1}{2}\frac{d{\sigma }^2}{dy}N\right)-\frac{\partial }{\partial x}\left(\frac{dx}{dy}N\right)-\frac{1}{A_{\ast }h(y)t_{\mathrm{esc}}(x,y)}N,$$

$$y(t)\equiv A_{\ast }{\int }_{t_{\ast }}^{t}h(t^{\prime })d{t}^{\prime },$$

$$N(x,t)\equiv 4\pi (m_{e}c{)}^3{x}^{2}f(x,t).$$

$$\frac{1}{2}\frac{d{\sigma }^2}{dy}=0,\frac{dx}{dy}=1-\hat{S}{x}^2,$$

$$N(x,y)\Big{|}_{y=0}=\dfrac{J_{0}}{\sigma\sqrt{2\pi}}e^{\frac{-(x-\mu)^{2}}{2\sigma^{2}}}\ ,\qquad\frac{1}{\sqrt{\hat{S}}}\geq x\geq 0\ ,$$

$${\mathcal{N}}_0\equiv J_0{\int }_0^{1/\sqrt{\hat{S}}}\frac{1}{\sigma \sqrt{2\pi }}{e}^{\frac{-(x-\mu {)}^2}{2{\sigma }^2}}dx=\frac{J_0}{2}\left[\mathrm{Erf}\left(\frac{\mu }{\sqrt{2}\sigma }\right)-\mathrm{Erf}\left(\frac{\mu -1/\sqrt{\hat{S}}}{\sqrt{2}\sigma }\right)\right].$$

$${\mathcal{E}}_0=J_0{\int }_0^{1/\sqrt{\hat{S}}}\frac{m_e{c}^2\sqrt{x^2+1}}{\sigma \sqrt{2\pi }}{e}^{\frac{-(x-\mu {)}^2}{2{\sigma }^2}}dx.$$

$${\bar{\gamma }}_0=\sqrt{{\bar{x}}_0^2+1}=\frac{{\mathcal{E}}_0}{{\mathcal{N}}_0{m}_e{c}^2},$$

$$N_{\mathrm{rise}}(x,y)=\frac{J_0}{\sigma \sqrt{2\pi }}{\left[\frac{x_0(x,y)}{x}\right]}^{\hat{C}}{\left[\frac{1-\hat{S}{x}_0^2(x,y)}{1-\hat{S}{x}^2}\right]}^{1-\frac{\hat{C}}{2}}\exp \left\{-\frac{[\mu +x_0(x,y){]}^2}{2{\sigma }^2}\right\},x_{\min }(y)<x<\frac{1}{\sqrt{\hat{S}}},$$

$$x_{\min }(y)=\frac{1}{\sqrt{\hat{S}}}\tanh \left(y\sqrt{\hat{S}}\right),$$

$$x_0(x,y)=\frac{1}{\sqrt{\hat{S}}}\tanh \left[{\tanh }^{-1}(x\sqrt{\hat{S}})-y\sqrt{\hat{S}}\right].$$

$$N_{\mathrm{decay}}(x,y_{\mathrm{pk}})=N_{\mathrm{rise}}(x,y_{\mathrm{pk}}),$$

$$\begin{gathered} N_{\mathrm{decay}}(x,y)=\frac{J_0{e}^{-(t-t_{\mathrm{pk}})/t_{\mathrm{ad}}}}{\sigma \sqrt{2\pi }}\left[\frac{1-\hat{S}{x}_0^2(x,y)}{1-\hat{S}{x}^2}\right]\exp \left\{-\frac{{\left[\mu +x_0(x,y)\right]}^2}{2{\sigma }^2}\right\} \\ \times {\left[\frac{1-\hat{S}{x}_0^2(x,y-y_{\mathrm{pk}})}{1-\hat{S}{x}_0^2(x,y)}\right]}^{\hat{C}/2}{\left[\frac{x_0(x,y-y_{\mathrm{pk}})}{x_0(x,y)}\right]}^{-\hat{C}},x_{\min }(y)<x<\frac{1}{\sqrt{\hat{S}}},y\ge y_{\mathrm{pk}}, \end{gathered}$$

$$N(x,y)=\{\begin{array}{ll}{N}_{\mathrm{rise}}(x,y),& y_{\mathrm{pk}}\ge y\ge 0,\\ N_{\mathrm{decay}}(x,y),& y\ge y_{\mathrm{pk}}.\end{array}$$

$$Q_{\nu }(\nu ,\gamma )=\frac{\sqrt{3}{q}^{3}B}{m_e{c}^2}R\left(\frac{\nu }{{\gamma }^2{\nu }_s}\right)\propto \mathrm{erg}{\mathrm{s}}^{-1}\mathrm{H}{\mathrm{z}}^{-1},$$

$${\nu }_s\equiv \frac{3qB}{4\pi m_{e}c},$$

$$R(x)=\frac{x^{2}}{2}K_{4/3}\Big{(}\frac{x}{2}\Big{)}K_{1/3}\Big{(}\frac{x}{2}\Big{)}-\frac{3x^{3}}{20}\Big{[}K^{2}_{4/3}\Big{(}\frac{x}{2}\Big{)}-K^{2}_{1/3}\Big{(}\frac{x}{2}\Big{)}\Big{]}\ .$$

Full text

1 INTRODUCTION

The Crab Nebula has been observed extensively since its progenitor star went supernova in 1054 CE, and in recent decades its high-energy emission has displayed a rich variety of behaviors which have been documented by a series of space-based $\gamma$-ray and X-ray observatories. It is a relatively nearby and young pulsar wind nebula that, until recently, was thought to produce stable synchrotron emission from electrons and positrons (hereafter, electrons) in the remnant’s magnetosphere, and thus has been used as a calibration source (for a review see Bühler & Blandford, 2014). The nebula is powered by the rapidly rotating pulsar which is gradually slowing down. A portion of the resulting spin-down power, $\approx 5\times 10^{38}\,\rm erg\,s^{-1}$ is deposited into the outward-flowing electron wind. A standing termination shock forms at the radius at which ram pressure balances the gas pressure of the nebula, corresponding to a distance $r_{t}\sim 10^{17}\,$cm from the pulsar (Rees & Gunn, 1974; Montani & Bernardini, 2014).

The relativistic electrons in the upstream portion of the wind interact with the ambient magnetic fields and produce synchrotron radiation as they advect outwards. Particles accelerated at the termination shock can explain the broad energy distribution of the observed quiescent emission (Gaensler & Slane, 2006). The accelerated particles diffuse and advect downstream from the shock, producing the synchrotron emission observed in the outer region of the nebula. Thus, it is not surprising that the apparent angular size of the Crab Nebula increases with decreasing energy (Abdo et al., 2011). The oldest electrons continue to advect far downstream of the shock; the emission characterizes the radio synchrotron nebula out to a radius of $\sim 10^{18}\,$cm.

The 2007 September $\gamma$-ray flare event was observed by AGILE and analyzed by Striani et al. (2013). They define a classification scheme for episodes of enhanced $\gamma$-ray emission in which “sub-flares” are very short-duration events lasting a day or less, and “waves” are outbursts that have longer durations and peak brightnesses about half that of the sub-flares. The light curve plotted in Figure 1 from Striani et al. (2013) displays a clearly-defined and isolated wave event, designated W1, which takes place roughly between 2007 September 25-30 (MJD 54368-54373) and is followed by a second wave, W2. However, in the case of W2, after the wave reaches its peak, it is interrupted by a series of bright and short-lived sub-flares which are designated as F1, F2 and F3. The authors present a time-averaged spectrum for the W1 event over the dates in which it is defined, and a time-averaged spectrum for the W2 event up to its peak only. They also provide a spectrum for the peak of the F2 sub-flare event, integrated over 12-hour time bins. Based on the time constraints associated with the various spectral features that they define, we find it convenient to focus on analysis of the F2 sub-flare and the W1 wave event, since these have the highest quality spectral data.

Kroon et al. (2018, hereafter Paper 1) developed a time-dependent self-similar model describing the acceleration and energy losses experienced by a population of magnetically-confined relativistic electrons in the vicinity of the pulsar-wind termination shock. The observed synchrotron radiation produced by these electrons is determined by the local magnetic field strength, combined with the energy distribution of the radiating electrons. They assumed a Gaussian initial energy distribution for the incident electrons, which is a reasonable approximation for the conditions in the wind plasma just upstream from the shock (Cerutti et al., 2012, Paper 1). By combining this initial condition with a rigorous transport equation, they derived an exact analytical expression for the time-evolving energy distribution of the electron population, subject to the acceleration and losses encountered as they pass through the region of explosive magnetic reconnection surrounding the termination shock. The associated synchrotron spectrum can be computed exactly, and the model can be tested via comparisons with the $\gamma$-ray data. The synchrotron spectra computed using the model are self-consistent since synchrotron energy losses are included in the transport equation. In Paper 1 we used the time-dependent model to analyze the 2011 April Crab Nebula $\gamma$-ray event observed using the LAT detector onboard Fermi. Furthermore, an earlier, time-independent version of the model was used to analyze all of the available $\gamma$-ray spectra produced during flares observed between 2007-2013 (Kroon et al., 2016).

The paper is organized as follows. In Section 2 we present an abbreviated review of the time-dependent model developed in Paper 1. In Section 3 we apply the model to the analysis of the 2007 September event and discuss the resulting physical parameters. In Section 4 we provide our conclusions and discuss their astrophysical significance.

2 Particle Transport Formalism

In Paper 1 we modeled the 2011 April $\gamma$-ray flare using a time-dependent particle transport formalism, based on an analytical solution to a first-order transport equation, which includes terms describing electrostatic acceleration, synchrotron losses, and particle escape. The theoretical spectra computed at various time intervals compare quite favorably with the Fermi spectral data for the 2011 April event. Furthermore, the theoretical light curve for the event also agrees with the observations (Bühler et al., 2012). Here, we will apply the time-dependent model developed in Paper 1 to the analysis and interpretation of the 2007 September flare event observed by AGILE. The primary goal of this paper is to determine whether the same model can successfully account for the long-duration “wave” events observed during that transient (Striani et al., 2013).

Our model is similar to the “blob” picture proposed by Zrake (2016). In this scenario, the GeV synchrotron emission observed during the flares is produced by magnetically-confined populations of relativistic electrons that are processed via particle acceleration, losses, and escape occurring in the active reconnection region surrounding the pulsar-wind termination shock. Individual sub-flares and waves are attributed to separate blobs of plasma that independently encounter the shock. Since little is known about the precise nature of the upstream pulsar wind, we will assume that the energy distribution of the electrons in the blob as it impinges on the shock can be approximated using a broad Gaussian function whose mean ($\mu$) and standard deviation ($\sigma$) are treated as free parameters (Cerutti et al., 2012, Paper 1). Our time-dependent model allows us to follow the detailed evolution of the electron energy distribution as the particles experience strong electrostatic acceleration and synchrotron losses in the vicinity of the shock. In the remainder of Section 2 we provide a brief review of our model, and refer interested readers to Paper 1 for details.

2.1 Particle Transport Equation

We can model the evolution of the energy distribution of the blob electrons during a wave or sub-flare event using a spatially-averaged, time-dependent transport equation that includes terms describing electrostatic acceleration, synchrotron losses, and particle escape,

[TABLE]

where $p$ is the electron momentum, $t$ is time, $m_{e}$ is the electron mass, $c$ is the speed of light, and $f$ denotes the electron momentum distribution function. The mean time for particles to escape from the blob is represented by $t_{\rm esc}$, which can in principle depend on both $p$ and $t$. This is further discussed in Section 2.2. The (time-dependent) total number of electrons contained in the blob is computed from the distribution function $f$ using

[TABLE]

Note that particle escape tends to decrease the number of particles in the blob over time since we assume that there is no continual particle injection once the population of electrons is initialized.

2.2 Particle Escape Mechanisms

In the one-zone model considered here, the electron population represents a spatial average over the active region surrounding the shock, where explosive reconnection and particle acceleration occurs. Following Steinacker & Schlickeiser (1989), we characterize each region according to the dominant mechanism governing the escape of particles from the active region. On the upstream side, the process known as “shock-regulated escape” (SRE) is thought to dominate particle escape (see Paper 1). In this scenario, higher-energy particles (with large Larmor radii) have a higher probability of being recycled back to the upstream side of the shock through the shock. Conversely, lower-energy particles will almost with small Larmor radii are swept downstream with the bulk flow. The mean escape timescale for the SRE process is given by (Kroon et al., 2016)

[TABLE]

where $p$ is the particle momentum and $C(t)$ is a time-dependent parameter that controls the strength of the SRE process based on physical conditions such as the shock obliquity.

Once the plasma has moved through the termination shock, the dominant mechanism governing the escape of particles from the active reconnection zone surrounding the shock switches to an advection-dominated mode. In this region, even the particles with large Larmor radii are no longer recycled through the shock for additional acceleration, and instead they are swept away by the tangled magnetic field into the outer region of the nebula. Hence, on the downstream side of the shock, the escape timescale should be interpreted as the mean time for electrons to be advected away from the active zone around the termination shock. In this case, the escape timescale is independent of both energy and time, and we can write

[TABLE]

where $R_{\rm b}$ is the blob radius and $v_{ds}=c/3$ denotes the downstream flow velocity (Achterberg et al., 2001). Our model accounts for the variation in the dominant particle escape mechanism by incorporating two different formulations for the escape timescale, with Equation (3) utilized in the upstream region and Equation (4) utilized in the downstream region. This allows us to more accurately model the evolution of the electron distribution during the $\gamma$-ray flares within the confines of a spatially-averaged model. The upstream and downstream regions are temporally associated with the rising and decaying phases of the sub-flare (or wave), respectively.

2.3 Profile Function and Self-Similarity

In our self-similar time-dependent model (Paper 1), the electrostatic acceleration and synchrotron loss processes are parameterized by the functions $A(t)$ and $S(t)$ in Equation (1), which are related to the time-dependent electric and magnetic fields, $E(t)$ and $B(t)$, respectively, via

[TABLE]

where $q$ is the magnitude of the electron charge, the subscript “$*$” denotes the initial value of a quantity measured at the beginning of the $\gamma$-ray flare (at time $t=t_{*}$), and $h(t)$ represents the self-similar “profile function,” which parameterizes the variation of the electric and magnetic fields due to reconnection at the termination shock. We assume that the profile function, $h(t)$, varies according to (Bühler et al., 2012)

[TABLE]

where $\alpha$ and $\theta$ are the rising and decaying time constants, respectively, and $t$ and $t_{\rm pk}$ are measured with respect to the origin of the event at time $t_{*}$. This functional form for $h(t)$ was successfully used to fit the spectral data for the 2011 April flare from the Crab Nebula, and it will also allow us to reproduce the 2007 September event of interest here. It follows that the minimum value of the profile function, at the beginning of the event, is $h(t_{*})=1$, and the maximum value, which occurs at the peak of the sub-flare or wave, is $h(t_{\rm pk})=e^{\alpha}$. Based on Equations (5), we can express the initial values of the functions $A(t)$ and $S(t)$ in terms of the initial electric and magnetic fields, $E_{*}$ and $B_{*}$, respectively, obtaining

[TABLE]

Combining Equations (5) and (7), one finds that in the self-similar model, the time variations of the electric and magnetic fields are given by

[TABLE]

which implies that

[TABLE]

Some useful insight on the validity of this self-similar field variation is provided by the work of Lyutikov et al. (2017, 2018) and (Cerutti et al., 2014). These authors perform particle-in-cell (PIC) simulations to study the effect of explosive reconnection in plasmas with $\sigma\gg 1$. We note that Figure 7 from Lyutikov et al. (2017) indicates that the reconnection velocity $v_{\rm rec}$ increases approximately in proportion to the square-root of time during the early phase of the reconnection event, and therefore we can write

[TABLE]

Furthermore, it follows from Equation (4.2) in Lyutikov et al. (2017) that

[TABLE]

Figure 7 from Lyutikov et al. (2017) also indicates that

[TABLE]

early in the flare’s evolution. Combining Equations (10), (11), and (12), we find that $B\propto t^{1/2}$, and consequently $E\propto B^{2}$, in agreement with Equation (9). This supports the plausibility of the self-similar relationship between $E$ and $B$ implied by our model.

2.4 Fokker-Planck Equation

Moving forward with analysis of the particle transport equation, it is convenient to transform from the momentum $p$ to the dimensionless parameter

[TABLE]

so that $x=\sqrt{\gamma^{2}-1}$, where $\gamma$ is the electron Lorentz factor. In the case of the ultra-relativistic electrons that power the Crab Nebula $\gamma$-ray flares, $\gamma\gg 1$, and therefore $x\to\gamma$. We can now express Equation (1) in Fokker-Planck form by writing

[TABLE]

where the dimensionless time

[TABLE]

and the electron number distribution

[TABLE]

The “broadening” and “drift” coefficients in Equation (14) are given by

[TABLE]

respectively. Since our model does not include momentum diffusion, the broadening coefficient vanishes in our application. Solving Equation (14) requires the specification of an initial condition and also a functional form for the escape timescale $t_{\rm esc}(x,y)$. As explained in Section 2.2, we use Equation (3) to compute $t_{\rm esc}$ in the upstream region and Equation (4) in the downstream region. As discussed in Sections 2.5 and 2.6, we also employ two different initial conditions when solving Equation (14) for the electron number distribution $N$ during the rising and decaying phases of each $\gamma$-ray sub-flare (or wave).

2.5 Rising Phase

During the rising phase of the sub-flare (or wave), one can obtain an exact solution for the electron distribution in the plasma blob by solving Equation (14) subject to a Gaussian initial condition and setting the escape time $t_{\rm esc}=t_{\rm SRE}$ using Equation (3). The Gaussian initial condition represents the (incident) electron distribution in the upstream pulsar wind; that is

[TABLE]

where $\mu$ and $\sigma$ represent the mean and standard deviation, respectively, and $J_{0}$ is a normalization factor. The total number of electrons initially contained in the plasma blob is given by

[TABLE]

The corresponding total initial energy of the electrons, $\mathscr{E}_{0}$, is computed by numerically evaluating the integral

[TABLE]

By combining these relations, we find for the mean initial Lorentz factor for the electrons in the plasma blob is given by

[TABLE]

where $\bar{x}_{0}$ denotes the corresponding initial value for the mean dimensionless momentum.

Solving Equation (14) subject to the initial condition given by Equation (18) yields the exact solution for the electron distribution during the rising phase of the transient. The result obtained is

[TABLE]

where the minimum energy at dimensionless time $y$ is given by

[TABLE]

and the injection energy for particles with current energy $x$ at time $y$ is computed using

[TABLE]

The parameter $\hat{C}$ appearing in Equation (22) is related to the shock-regulated escape process (see Equation (32) from Paper 1). The exact solution for the electron distribution during the rising phase of the event given by Equation (22) can be used to compute a series of time-dependent theoretical $\gamma$-ray synchrotron spectra for comparison with a sequence of observational spectral data for a given sub-flare or wave, up to the peak time, $t=t_{\rm pk}$, or equivalently, $y=y_{\rm pk}$. The synchrotron spectra thus obtained are self-consistent since synchrotron energy losses are included in the electron transport equation. In order to complete the theoretical picture, we must also employ the solution for the electron distribution during the decaying phase of the event, as discussed below.

2.6 Decaying Phase

The determination of the electron distribution during the decaying phase of the event requires us to solve Equation (14), combined with an appropriate initial condition, while also setting the escape time $t_{\rm esc}=t_{\rm ad}$ using Equation (4). The initial condition is given by the rising phase solution (Equation (22)) evaluated at the peak of the sub-flare or wave; that is,

[TABLE]

where $y=y_{\rm pk}$ at the peak of the event. Equation (25) can be combined with Equation (14) to obtain the exact solution

[TABLE]

where $x_{0}(x,y)$ and $x_{0}(x,y-y_{\rm pk})$ are computed using Equation (24). Taken together, Equations (22) and (26) provide the exact solution for the electron number distribution during the rising and decaying phases of a single sub-flare or wave. The global solution is therefore given by

[TABLE]

2.7 Synchrotron Emission

The synchrotron spectra and light curves computed using the electron distribution represented by Equation (27) is self-consistent since synchrotron losses are included in the transport equation that we have solved. We assume here that the electron distribution is isotropic in the frame of the bulk flow, so that the synchrotron emission can be computed by convolving the electron number distribution function (Equation (27)) with the single-particle synchrotron emission function (e.g., Becker, 1992; Kroon et al., 2016)

[TABLE]

where

[TABLE]

and (Crusius & Schlickeiser, 1986)

[TABLE]

The functions $K_{4/3}(x)$ and $K_{1/3}(x)$ represent modified Bessel functions of the second kind. The corresponding total synchrotron power per unit frequency emitted by the isotropic electron distribution in the frame of the plasma is obtained by integrating over the electron number distribution, yielding

[TABLE]

with the corresponding observed flux density

[TABLE]

where $D$ is the distance to the source and $P_{\nu}(\nu,t)$ is evaluated using Equation (31). The theoretical light curve can be computed by integrating the flux with respect to the photon frequency $\nu$.

3 Application to 2007 September Event

In this Section, we apply our model to the 2007 September $\gamma$-ray flare observed by AGILE and refer to Striani et al. (2013) for the various data sets. Following Striani et al. (2013), we will use the terms “wave” and “sub-flare” to describe the two different types of temporal events observed during the 2007 September $\gamma$-ray flare. Waves are defined as events lasting about one week, whereas sub-flares last about a day or less. Hence these two types of events are qualitatively different, and therefore it is important to test our theoretical model by determining if it is capable of accounting for the spectral data and the light curves for both types of transients. Striani et al. (2013) provide spectra and light curves for the W1 wave and the sub-flares F1, F2, and F3. In this paper, we will focus on the W1 and F2 events since these have the highest quality of the available spectral data.

By utilizing our time-dependent model to analyze the spectral data for the 2007 September flare, we can extract a number of interesting physical quantities and functions such as the variable electric and magnetic field strengths, the variable electron energy distribution, the initial blob energy, the total radiated synchrotron energy, and the total event energy. We apply our model to the interpretation of the data for the W1 wave and the F2 sub-flare events in Sections (3.1) and (3.2), respectively.

3.1 F2 Sub-Flare

The light curve plotted in Figure 1 from Striani et al. (2013) depicts 12-hour binning of the flaring activity from MJD 54360-54390 (2007 September 27 - October 17), which comprises two waves and three sub-flares. The waves are plotted in black and labeled either W1 or W2. These waves have durations of several days, and hence they are clearly resolved using 12 hour time bins. On the other hand, the sub-flares F1, F2, and F3 are very short-duration episodes of enhanced emission, which are barely resolved by the 12 hour time bins. Hence no time-dependent spectra are available for the sub-flares, but instead only the peak spectra are available and therefore we will compare our model predictions with the peak spectra. The brightest sub-flare is the F2 event, and the peak spectrum for this feature is plotted in Figure 5 from Striani et al. (2013). We can use our model to interpret the F2 sub-flare by computing the peak spectrum using Equation (31) and then comparing the result with the corresponding AGILE spectrum taken from Striani et al. (2013). The theoretical parameters are varied until an acceptable qualitative fit is obtained.

In Figure 1, we plot the theoretical peak spectrum, $\mathscr{F}_{\nu}(\nu,t_{\rm pk})$, for the F2 event computed using Equation (32) along with the observational data taken from Figure 5 in Striani et al. (2013), based on the theoretical parameter values listed in Tables 1 and 2. The peak magnetic field for the F2 sub-flare is $B_{\rm pk}=713\,\mu$G, and the event exhibits strong electrostatic acceleration, with a peak electric to magnetic field ratio $E_{\rm pk}/B_{\rm pk}=1.23$. The peak electric field strength is $E_{\rm pk}=877\,\mu$G. The 12-hour duration of this sub-flare and the significant induced electrostatic fields at the peak imply that a substantial explosive reconnection event is taking place at the termination shock. The associated light curve above 100 MeV is compared with the observational data in Figure 2. It is important to note that according to Figure 1 from Striani et al. (2013), the rising side of the F2 sub-flare is “contaminated” by the decaying side of the F1 sub-flare. Since the F1 feature is not modeled here (due to a lack of spectral data), we only plot the theoretical light curve for F2. Hence the theoretical light curve slightly underestimates the observed flux on the rising side of F2, as expected.

3.2 W1 Wave

During the 2007 September flare from the Crab Nebula, the W1 wave event persisted over several days, during the time interval MJD 54368-54373 (2007 September 25-30). Striani et al. (2013) provide a single integrated spectrum for the W1 wave in their Figure 5, but a true time-dependent spectrum for this feature is not currently available in the literature. The corresponding light curve for the W1 event is depicted in both Figures 1 and 4 from Striani et al. (2013) using 12-hour and 24-hour time binning, respectively. The peak is clearly resolved in their Figure 1 but not in Figure 4 due to the coarse time bins. In this study, we will therefore test our model using the 12-hour binned light curve for W1 taken from Figure 1 of Striani et al. (2013).

Although our theoretical model is fully time-dependent, in the case of the W1 transient, we must perform a time average over the interval MJD 54368-54373 in order to make contact with the spectrum for this event plotted in Figure 5 from Striani et al. (2013). Hence the theoretical spectrum for the W1 event is obtained using the time integration

[TABLE]

where $t_{1}$ and $t_{2}$ are the starting and ending times of the W1 flare, and $\mathscr{F}_{\nu}(\nu,t)$ is computed using Equation (32). The resulting spectrum, $\bar{\mathscr{F}}_{\nu}(\nu)$, is compared with the observed time-averaged spectrum for W1 taken from Figure 5 in Striani et al. (2013), and the theoretical parameters are varied until adequate agreement is obtained. Our results for the theoretical W1 spectrum are plotted along with the observational data in Figure 3, and the corresponding light curve above 100 MeV is plotted and compared with the observational data in Figure 4. The W1 wave event is not as bright at its peak as the F2 sub-flare, and therefore the peak magnetic field for W1 is only $B_{\rm pk}=472\,\mu$G. Further parameters are provided in Tables 1 and 2.

4 Discussion and Conclusion

The remarkable series of bright $\gamma$-ray flares emitted by the Crab Nebula and observed by AGILE or Fermi between 2007-2013 has presented considerable challenges to classical particle acceleration models, due in part to the apparent violation of the classical synchrotron burn-off limit, implying particle acceleration on sub-Larmor timescales. We have previously demonstrated in Paper 1 that the time-dependent spectra observed during the intense 2011 April $\gamma$-ray flare can be described using a relatively simple analytical model. A time-independent version of the model was also compared with all of the available flare data by Kroon et al. (2016). While the earlier model comparisons have been promising, the 2007 September Crab Nebula $\gamma$-ray flare was a unique event since it comprised both short-duration sub-flares of the sort that have been seen during other transients, combined with long-duration wave events lasting about one week. Hence the analysis of these data represents the strongest test of the theoretical model we have developed. This has motivated us to test our model by making detailed comparisons with all of the available data for the 2007 September event, as presented by Striani et al. (2013). We discuss our primary findings below.

The 2011 April light curve is characterized by two distinct peaks, known as sub-flares (Bühler et al., 2012), which in our model are treated as emission from two individual blobs that pass separately through the termination shock, but are both included in the same overall flare event. The first 2011 sub-flare was found to have a peak magnetic field $B_{\rm pk}=706\,\mu$G and an electric/magnetic field ratio $E_{\rm pk}/B_{\rm pk}=1.84$ at the peak. These values are similar to those obtained here for the F2 sub-flare observed during the 2007 event, for which we obtain a peak magnetic field $B_{\rm pk}=713\,\mu$G and a peak field ratio $E_{\rm pk}/B_{\rm pk}=1.23$. In our analysis of the 2007 September flare, we find that the mean Lorentz factor for the incident electrons is $\bar{\gamma}_{0}\sim 10^{9}$ (see Table 2), which is significantly larger than the upstream Lorentz factor in the cold pulsar wind, $\Gamma\sim 10^{6}$ (Cerutti et al., 2014). This implies that the electrons in the incident plasma blob are produced as a result of impulsive reconnection on the upstream side of the termination shock (Montani & Bernardini, 2014). The values for the magnetization parameter, $\sigma$, obtained here are $\sigma=0.036$ for W1 and $\sigma=7.4\times 10^{-5}$ for F2, which agree with the expected values in the downstream region (Sironi & Spitkovsky, 2009). On the other hand, PIC simulations carried out by Lyutikov et al. (2017, 2018) and Cerutti et al. (2014) suggest that $\sigma\gg 1$. This is known as the “$\sigma$-problem” in pulsar astrophysics and is currently unresolved (Bühler & Blandford, 2014; Lyutikov et al., 2018).

Striani et al. (2013) posit that sub-flares are powered by emission from blobs with radii $R_{\rm b}\sim(1-5)\times 10^{15}\,$cm possessing relatively strong magnetic fields, on the order of $10^{3}\,\mu$G. Conversely, the waves are associated with larger plasma blobs with radii R_{\rm b}\lower 1.72218pt\hbox{;\buildrel>\over{\scriptstyle\sim};}10^{16}\,cm, possessing weaker magnetic fields that the sub-flare blobs. It is interesting to compare the estimates obtained by Striani et al. (2013) with the physical parameters resulting from the application of our model to the same event. During the sub-flare F2, we find a blob radius $R_{\rm b}=10^{15}$ cm and a peak magnetic field $B_{\rm pk}=713\,\mu$G. By contrast, during the W1 wave we obtain a blob radius $R_{\rm b}=10^{16}$ cm and a peak magnetic field $B_{\rm pk}=472\,\mu$G. These values agree reasonably well with the estimates obtained by Striani et al. (2013), which are $(1.5\pm 0.5)\times 10^{3}\,\mu$G for F2 and $(0.8\pm 0.2)\times 10^{3}\,\mu$G for W1. However, it is important to emphasize that the parameters obtained by Striani et al. (2013) were based on an analysis of the energy budget for the event, whereas our parameters are obtained via detail spectral comparisons. Hence the two determinations are separate but complementary, and the fact that they agree strongly supports the interpretation that the sub-flares and waves are powered by explosive reconnection.

The “sub-flare/wave” paradigm adopted by Striani et al. (2013) represents an interesting classification scheme, since the waves tend to be longer-lasting but less intense episodes of enhanced $\gamma$-ray activity, compared to the more powerful but brief “sub-flare” counterparts. This classification can be further explored by examining the energetics of the various temporal features observed during the 2007 flare. The energy integrals for each channel are defined in Equations (58)-(65) from Paper 1, and are denoted by $\mathscr{E}_{\rm inj}$ (initial electron energy in blob), $\mathscr{E}_{\rm part}(t)$ (total electron energy in blob), $\mathscr{E}_{\rm elec}(t)$ (total electrostatic acceleration energy), $\mathscr{E}_{\rm synch}(t)$ (total energy emitted in synchrotron radiation), and $\mathscr{E}_{\rm esc}(t)$ (total energy in escaped electrons). We have computed the time-dependent energy in each of these channels and plotted them as functions of time for W1 and F2 in Figure 5. The asymptotic values for each channel energy obtained as $t\to\infty$ are listed in Table 3.

Despite the difference in duration between W1 and F2, we find that they each produce a total energy of $\sim 2\times 10^{41}\,$erg in the synchrotron channel. The average power of the sub-flare F2 is about five times greater than that of the W1 wave due to its shorter duration. Figure 5 indicates that most of the energy advects downstream from the shock and goes on to power the synchrotron emission in the outer region of the nebula. Depending on the strength of the magnetic field in the nebula, this could result in an observable X-ray afterglow. In fact, Kroon et al. (2016) considered the afterglow associated with this theoretical model for each of the various Crab Nebula $\gamma$-ray flares. They found that an observable X-ray afterglow could persist for up to two weeks, depending on the magnetic field strength. The theoretical prediction for the 2007 September flare was plotted in Figure 5 from Kroon et al. (2016), but no X-ray observations were performed after this particular $\gamma$-ray transient, and therefore the prediction cannot be compared with data.

The standard observational estimate for the efficiency of the conversion of the pulsar’s spin-down power into synchrotron radiation is $\sim 30\%$ (Abdo et al., 2011). It is interesting to compare this estimate with the theoretical efficiency, defined as the ratio $\mathscr{E}_{\rm synch}/(\mathscr{E}_{\rm synch}+\mathscr{E}_{\rm esc})$, which can be evaluated using the asymptotic values for the channel energies listed in Table 3. The resulting theoretical efficiencies obtained for the F2 and W1 events are $\sim 10\%$ and $\sim 15\%$, respectively, which are similar to the results obtained in our analysis of the 2011 April flare in Paper 1. The discrepancy between our results and the spin-down estimates could be due to the fact that the spin-down estimate is a time average, and therefore it does not accurately represent the energetics during an individual transient. The asymptotic values of the electrostatic energy channels in the W1 and F2 events are both $\sim 5\times 10^{41}\,$erg, which is expected given the similar levels in the synchrotron energy channel.

It is interesting to compare the energy budget for the 2007 September event treated here with that from our previous analysis of the 2011 April flare in Paper 1. The total synchrotron energies in sub-flares 1 and 2 for the 2011 April event were found to be $4\times 10^{41}$ and $7\times 10^{41}$ erg, respectively, which exceed the synchrotron values derived here for W1 and F2 by a factor of $\sim 2$. Comparison of the energy in the electrostatic channel for each event can reveal useful information about the magnitude of the induced electrostatic fields primarily responsible for accelerating the electrons. For the 2007 September event, the maximum energy generated by the fields is seen to be about half of the initial energy in the blob for both W1 and F2. However, in the two sub-flares observed during the 2011 April event, we found in Paper 1 that the electrostatic energy considerably exceeds the initial energy of the blob, which implies a very energetic and efficient conversion of wind energy into synchrotron $\gamma$-rays. Hence the 2007 flare was a much more mild event. Another useful metric for establishing the importance of electrostatic acceleration is provided by the peak-field ratio $E_{\rm pk}/B_{\rm pk}$, which traces the explosiveness of magnetic reconnection during the transient. In the case of the W1 and F2 features, this ratio is of order unity (Table 2). However, during the more energetic and efficient sub-flares observed in 2011 April, we found that $E_{\rm pk}/B_{\rm pk}\sim 3$.

For the 2011 April sub-flares, our model indicates that the peak magnetic fields were in the range $B_{\rm pk}\sim 600-700\,\mu$G (see Paper 1). These results are similar to those obtained here for the 2007 flare, which implies that the same particle acceleration and emission mechanisms are operating throughout the range of $\gamma$-ray flares observed from the Crab Nebula. We conclude that our model is able to reproduce the observational data ($\gamma$-ray spectra and integrated light curves) for both the 2007 September and the 2011 April $\gamma$-ray flares from the Crab Nebula, despite the fact that the energy budgets were distinctly different for these two events. Analysis of the 2007 September $\gamma$-ray flare is especially important because this event comprised qualitatively different temporal features compared with any of the other Crab Nebula flares. The ability of our model to account for both the long-duration waves and the short sub-flares suggests that our theoretical framework for treating particle acceleration in the explosive reconnection environment surrounding the termination shock can help to provide a detailed understanding of the physics taking place in these extreme astrophysical sources.

The authors are grateful to the anonymous referee for providing several useful comments that helped to strengthen and improve the presentation. J.J.K. was supported at NRL by NASA under contract S-15633Y. J.D.F. was supported by the Chief of Naval Research.