Consequences of Energetic Magnetar-like Outbursts of Nearby Neutron Stars: ¹⁴C Events and the Cosmic Electron Spectrum

Observations have revealed the existence of strong magnetic fields on neutron stars, with the dipole component ranging from 10¹¹–10¹³ G for pulsars and 10¹⁴–10¹⁵ G for magnetars. Neutron stars might also have multipole components much larger than its dipole field. The strong multipole components can be a consequence of dynamos inside the young fast spinning neutron star (Thompson & Duncan 1993) and they are buried under the surface by mass accretion (Romani 1990; Cumming et al. 2001), especially the prompt accretion of supernova fallback material (Chevalier 1989; Geppert et al. 1999). Magnetars are a special type of neutron stars showing violent transient radiative phenomena, including giant flares and soft gamma-ray bursts (GRBs) with luminosity up to 10⁴⁷ erg s⁻¹ (Kaspi & Beloborodov 2017). These outbursts are characterized by their X-ray luminosity exceeding the spin-down power of the neutron star, and it is believed they are powered by the high magnetic energy stored inside the star (Thompson & Duncan 1995). Such bursts and outbursts are not limited to magnetars. Recent observations have found that pulsars with high magnetic field also exhibit magnetar-like activities, e.g., PSR J1846–0285 (Gavriil et al. 2008) and PSR 1119–6217 (Göğüs et al. 2016; Archibald et al. 2017). These events can be powered by the high magnetic energy of the multipole components buried inside the neutron star.

The cosmogenic nuclide ¹⁴C is produced in atmosphere either by high-energy particles (e.g., protons) or γ-ray photons from high-energy phenomena. Through the global carbon cycle, some of the ¹⁴CO₂ produced in the atmosphere can be retained in annual tree rings (Damon et al. 1995; Usoskin et al. 2006). Recently, four events of rapid increase of the ¹⁴C content occurring in AD 994, AD 775, BC 660, and BC 3371 were found using tree rings (Miyake et al. 2012, 2013; Park et al. 2017; Wang et al. 2017; Büntgen et al. 2018). However, the physical origin of these events is mysterious. Several models have been proposed, such as gamma-ray emissions from nearby supernovae (Miyake et al. 2012) or GRBs (Hambaryan & Neuhäuser 2013; Pavlov et al. 2013), and large solar proton events (SPEs; Melott & Thomas 2012; Miyake et al. 2012; Usoskin et al. 2013; Wang et al. 2017). Some recent works (Mekhaldi et al. 2015; Büntgen et al. 2018) argued that large SPEs caused the ¹⁴C rapid-increase events. However, there are several problems with the interpretation of these ¹⁴C events as SPEs (Neuhäuser & Hambaryan 2014). First, there are no definite historic records of strong aurorae or sunspots around AD 775 and AD 994 (Chai & Zou 2015; Stephenson 2015). Second, if the AD 775 event would have been an SPE only a few times larger than the Carrington event in 1859 (Usoskin et al. 2013), it would have been the solar origin, there must also be a rapid ¹⁴C increase in the Carrington event, which is not the case (Miyake et al. 2012; Neuhäuser & Hambaryan 2014). Third, Cliver et al. (2014) compared energetics and spectrum of the 1956 SPE with the AD 775 case, and found the inferred solar fluence (>30 MeV) value for this event is inconsistent with the occurrence probability distribution for SPEs (Usoskin 2017), which casts strong doubts on the solar interpretation for the AD 775 event (Cliver et al. 2014). Fourth, whether the Sun can produce such large proton events is still debated (Schrijver et al. 2012). Recent studies shows that the occurrence frequency of superflares on solar-type stars is once in about three thousands years, which is too low (Notsu et al. 2019). Fifth, charged particles from SPEs are affected by the geomagnetic field, and will produce ¹⁴C symmetrically in both hemispheres. However, the latest measurements of Büntgen et al. (2018) showed that the ¹⁴C increase starts at different times for the AD 775 event in the southern and northern hemispheres. The short-duration outburst would shine on only one side of the Earth, which would affect the time and intensity of ¹⁴C increase at different places on Earth (Güttler et al. 2015).

The corresponding radiations caused a ¹⁴C rapid increase, which not only can deplete ozone, but can also damage communication, electronic, and power systems on the ground, and electronic equipment in satellites. Therefore, a study of the physical origin of ¹⁴C events is crucial. Here, we propose that the high-energy emissions from pulsar outbursts can cause rapid ¹⁴C increase. High-energy electrons are a very important probe of nearby cosmic-ray sources. Recent discoveries of an excess of electrons stimulated a number of works to discuss their possible origins. Outbursts of pulsars, which also generate high-energy electrons, may be responsible for the irregularities in the spectrum of the electrons observed by various space-borne experiments such as AMS-02 (Aguilar et al. 2014), Fermi-LAT (Abdollahi et al. 2017a), DAMPE (DAMPE Collaboration et al. 2017), and CALET (Adriani et al. 2018).

In Section 2, we describe the physical model of pulsar outbursts. Section 3 presents the gamma-ray emissions from outbursts and ¹⁴C events. Section 4 discusses the cosmic electron spectrum contributed by the outbursts. We search the possible neutron stars in our neighborhood in Section 5. Conclusions are given in Section 6. Throughout the paper, we adopt the convention Q_a = Q/10^a using cgs units.

In a pulsar with dipole field B = 10¹³ G, a magnetic multipole can be buried in the core with much higher magnetic field. Assuming the magnetic multipole has magnetic field strength B_m = 10¹⁶ G with length scale L = 10⁵ cm, the magnetic energy associated with this magnetic multipole is

$$\begin{eqnarray}&&{E}_{m}\sim \delta {B}^{2}{L}^{3}\sim {10}^{47}{B}_{m,16}^{2}{L}_{5}^{3}\,\mathrm{erg}.\end{eqnarray} \tag{ 1 }$$

Such magnetic energy is expelled from the core through ambipolar diffusion on the timescale (Beloborodov & Li 2016)

$$\begin{eqnarray}&&{t}_{\mathrm{amb}}=30{L}_{5}^{2/5}{B}_{m,16}^{2/5}/{B}_{13}^{8/5}{(\rho /{\rho }_{\mathrm{nuc}})}^{6/5}\,\mathrm{Myr}.\end{eqnarray} \tag{ 2 }$$

The dissipation of core magnetic energy can induce sudden global changes of magnetospheric structure, which powers the outbursts (Thompson & Duncan 1995). The origin of the sudden magnetospheric changes can be a result of crustal failure at the core–crust boundary when the high magnetic field escaping from the core breaks the crust (Thompson & Duncan 1995). Another possible mechanism to trigger the outburst is a slow buildup of magnetic energy in the magnetosphere followed by a sudden release (Lyutikov 2003; Gill & Heyl 2010; Parfrey et al. 2013). When the high magnetic field propagates outward in the crust, its high magnetic stress can break the crust and initiate plastic failure (Beloborodov & Levin 2014; Li et al. 2016). The plastic failure causes a surface displacement magnetic footpoint and leads to twisting of the surface magnetic field. If the twisting angle exceeds a critical value, the magnetosphere becomes unstable and produces a flare (Parfrey et al. 2013). In that case, the magnetosphere undergoes tearing instability and magnetic reconnection takes place with the launching of Alfvén waves and the ejection of plasmoids similar to the solar coronal mass ejection. The magnetic reconnection occurs within a very short time ∼100 R/c ∼ 10⁻² s (Uzdensky 2011), where R is the radius of neutron star.

Thompson & Duncan (2001) proposed that the extraordinarily high peak luminosity, and rapid variability of the initial spike emission from giant flares of magnetars imply that the spike emission must originate from a relativistically expanding fireball with an initial Lorentz factor of at least several tens, in order to avoid the pair-production problem. The observed light curve of the spikes of SGR 1806−20 (Hurley et al. 2005) is well explained by emission from relativistically expanding fireballs (Yamazaki et al. 2005). Meanwhile, the radio afterglow also can be fitted by a relativistic outflow with a Lorentz factor of a few tens (Wang et al. 2005). Therefore, similar to giant flares of SGRs (Hurley et al. 2005), we assume a relativistic outflow with Lorentz factor Γ₀ ∼ 10 is launched during the outburst of pulsars. The ejecta propagates into the interstellar medium (ISM) and builds up a relativistic forward shock and a reverse shock structure, including shock-accelerated electrons with associated synchrotron emission. Synchrotron radiation of the shocks can provide high-energy photons for rapid ¹⁴C increase and the accelerated electrons can responsible for the irregularities in the spectrum of cosmic electrons.

We first calculate the high-energy emissions from the outbursts of pulsars. From theory and numerical simulations (Blandford & Eichler 1987; Sironi et al. 2013), Fermi acceleration of charged particles in collisionless shocks leads to power-law distributions of energetic particles, described by ${dN}/d{\gamma }_{e}\propto {\gamma }_{e}^{-p}$, where γ_e is the electron Lorentz factor and p = 2.1–2.2. In the case of relativistic collisionless shocks, the Bohm approximation is not valid for maximum energy electrons (Kirk & Reville 2010; Lemoine 2013). The residence timescale downstream is given by ${t}_{r}={{Nt}}_{\lambda }={({R}_{{\rm{L}}}/\lambda )}^{2}(\lambda /c)\,=({R}_{{\rm{L}}}/\lambda )({R}_{{\rm{L}}}/c)$, where R_L is the Larmor radius, $N={({R}_{{\rm{L}}}/\lambda )}^{2}$ is the scattering time that the electrons experience before returning to the upstream, λ is the size of coherence cell, and t_λ is the time spent in crossing the coherence cell. The Lorentz factor of the maximal energy electrons can be derived by equating the synchrotron cooling time ${t}_{\mathrm{syn}}=6\pi {m}_{e}c/({\sigma }_{{\rm{T}}}{\gamma }_{e}{B}^{2})$ to the residence time t_r, which reads (Kirk & Reville 2010)

$$\begin{eqnarray}&&{\gamma }_{e,\max }={\left(6\pi \lambda {e}^{2}/{\sigma }_{{\rm{T}}}{m}_{e}{c}^{2}\right)}^{1/3}=1.2\times {10}^{7}{n}_{0}^{-1/6}{\lambda }_{1}^{1/3},\end{eqnarray} \tag{ 3 }$$

where $\lambda \equiv 10{\lambda }_{1}c/{\omega }_{{pi}}$ with plasma frequency ω_pi and n is the number density of the ISM.

When the relativistic outflow encounters the ISM, the external shock is generated, including forward and reverse shocks (Sari et al. 1998). The standard afterglow model of GRBs (Sari et al. 1998) is used. We use the standard afterglow model of GRBs (Mészáros 2002; Wang et al. 2015; Zhang 2018). _e and _B are the shock energy equipartition parameters for the shock-accelerated electrons and the magnetic fields, respectively. For the forward shock emission, the cooling frequency ν_c, the typical synchrotron frequency ν_m, and the maximum spectral flux F_ν,max are (Sari et al. 1998)

$$\begin{eqnarray}&&{\nu }_{c}=1.0\times {10}^{20}\,\mathrm{Hz}\,{E}_{\mathrm{iso},47}^{-1/2}{\epsilon }_{B,-2}^{-3/2}{n}_{0}^{-1}{t}_{\mathrm{obs},1}^{-1/2}{(1+Y)}^{-2},\end{eqnarray} \tag{ 4 }$$

$$\begin{eqnarray}&&{\nu }_{m}=8.4\times {10}^{15}\,\mathrm{Hz}\,{C}_{p}^{2}{E}_{\mathrm{iso},47}^{\displaystyle \frac{1}{2}}{\epsilon }_{B,-2}^{\displaystyle \frac{1}{2}}{\epsilon }_{e,-1}^{2}{t}_{\mathrm{obs},1}^{-\displaystyle \frac{3}{2}},\end{eqnarray} \tag{ 5 }$$

$$\begin{eqnarray}&&{F}_{\nu ,\max }=4.9\times {10}^{7}\,\mathrm{Jy}\,{E}_{\mathrm{iso},47}{\epsilon }_{B,-2}^{1/2}{n}_{0}^{1/2}{\left(\displaystyle \frac{d}{150\mathrm{pc}}\right)}^{-2},\end{eqnarray} \tag{ 6 }$$

where ${C}_{p}=3(p-2)/(p-1)$ and Y is the inverse-Compton parameter, which can be calculated from $Y\,\simeq [-1+\sqrt{1+4x{\epsilon }_{e}/{\epsilon }_{B}}]/2$ (Sari & Esin 2001), where $x\,=\min \{1,2.67{({\gamma }_{m}/{\gamma }_{c})}^{(p-2)}\}$ is the radiation coefficient of the shocked electrons, and γ_c is the electron cooling Lorentz factor (Sari et al. 1998)

$$\begin{eqnarray}&&{\gamma }_{c}\approx \displaystyle \frac{7.7\times {10}^{8}}{(1+Y)}\displaystyle \frac{1}{{\rm{\Gamma }}{B}^{{\prime} 2}{t}_{\mathrm{obs}}}.\end{eqnarray} \tag{ 7 }$$

We estimate the observed flux at 20 MeV as

$$\begin{eqnarray}\begin{array}{rcl}{f}_{{\nu }_{\mathrm{obs}}} & = & {F}_{\nu ,\max }\,{\nu }_{c}^{1/2}{\nu }_{m}^{(p-1)/2}{\nu }_{\mathrm{obs}}^{-p/2}=4.0\times {10}^{4}\,\mathrm{erg}\,{\mathrm{cm}}^{-2}\\ & & \times {\mathrm{MeV}}^{-1}\,{\epsilon }_{e,-1}^{p-1}{\epsilon }_{B,-2}^{\displaystyle \frac{p-2}{4}}{E}_{47}^{\displaystyle \frac{p+2}{\ 4}}{C}_{p}^{(p-1)}{\left(\displaystyle \frac{d}{150\mathrm{pc}}\right)}^{-2}{t}_{\mathrm{obs},1}^{\displaystyle \frac{2-3p}{4}}\\ & & \times {\left(1+Y\right)}^{-1}{\left(\displaystyle \frac{h{\nu }_{\mathrm{obs}}}{20\mathrm{MeV}}\right)}^{\displaystyle \frac{-(p-1)}{2}}.\end{array}\end{eqnarray} \tag{ 8 }$$

The value of Y is in order of 1 for typical parameters.

With typical parameters, the synchrotron radiation from forward shock can reach the GeV range. Meanwhile, the synchrotron-self-Compton radiation and inverse-Compton radiation will also contribute high-energy photons. For the chosen parameters, these are much weaker than the synchrotron component. Therefore, we ignore them. For typical parameters, the gamma-ray emission from a pulsar outburst is described by a power-law spectrum with an index of −1.1 (Sari et al. 1998). From the X-ray and radio afterglow observations of GW170817/GRB 170817, the spectral index is about −0.6 (Ruan et al. 2018; Troja et al. 2018), which indicates that the cooling break has not passed through the X-ray band. If the break passes through the X-ray band, the spectral index becomes steeper to −1.1, which is consistent with the standard afterglow model.

Using the GEANT4 simulation code with QGSP-BERT-HP, we computed the production yield of ¹⁴C due to gamma-rays in the atmosphere. High-energy photons with energy from 1 to 300 MeV are considered. In order to calculate the energy required by the increase of ¹⁴C, we apply GEANT4 code (version 10.5) to simulate the production rate of ¹⁴C. This simulation code uses the Monte Carlo method to trace particles involved in the interactions. The American standard atmosphere is used in our analysis. Considering the plane-parallel atmosphere model, we construct an atmosphere 87 km deep with 1 km resolution. The physics listed in QGSP_BERT_HP are used to simulate the physics process. For the gamma-ray photons, a power-law spectrum with the index α = −1.1 (Sari et al. 1998) is adopted with the minimum energy E_min = 1 MeV and the maximum energy E_Max = 100 MeV. 3000,000 photons satisfying the given spectrum are used in simulation. From the simulation, we collect 6894 ¹⁴C atoms and the average production rate of ¹⁴C is 60.67 atoms erg⁻¹. For the AD 775 event, the required ¹⁴C in atmosphere is Q = 1.3 × 10⁸ atoms cm⁻². The required energy is 7 × 10²⁴ erg in the atmosphere. The incident energy necessary for the increase of ¹⁴C content in the atmosphere for AD 775 event is about ${E}_{\mathrm{atm}}\sim 7\times {10}^{24}$ erg. Yields of ¹⁰Be and ³⁶Cl are about 0.36 and 0.12 atoms erg⁻¹, respectively. The transport and deposition of ¹⁰Be would be affected by gamma-ray radiation from outbursts (Thomas et al. 2005), and is not fully understood (Heikkilä et al. 2009). Hambaryan & Neuhäuser (2013) found that the ratio between ¹⁴C and ¹⁰Be can be explained in the gamma-ray photon case.

Below, we calculate the gamma-ray energy received in the atmosphere from pulsar bursts. The fluence F received in the atmosphere between 1 and 300 MeV from 1 to 10⁴ s,

$$\begin{eqnarray}&&F={\int }_{1\,\mathrm{MeV}}^{300\,\mathrm{MeV}}{\int }_{1\,{\rm{s}}}^{{10}^{4}\,{\rm{s}}}{f}_{{\nu }_{\mathrm{obs}}}{{dt}}_{\mathrm{obs}}d{\nu }_{\mathrm{obs}},\end{eqnarray} \tag{ 9 }$$

can be calculated numerically, where ${f}_{{\nu }_{\mathrm{obs}}}$ is defined in Equation (8). The energy received by the atmosphere can be calculated from ${E}_{\gamma }=\pi {R}_{\oplus }^{2}\times F$ (with Earth radius R_⊕). In addition, the relativistic outflow from pulsar outbursts usually has a small opening angle θ (Yamazaki et al. 2005; Lyutikov 2006). For example, Yamazaki et al. (2005) derived that the opening half-angle of the jet of the SGR 1806−20 giant outburst is 11^◦ from light curve (Yamazaki et al. 2005). From E_atm = E_γ, in order to produce these ¹⁴C events, the distance of the pulsar should be less than ∼200 pc, if θ ∼ 16° is adopted.

There are four ¹⁴C events in about 5000 yr. The rate of ¹⁴C events is about 1000 yr per event. Below, we estimate the pulsar outburst rate. The characteristic size of multipole field is L ∼ 10⁵cm and hence the maximum number of multiple field patches on the neutron star surface is N ∼ 4πR²/L² ∼ 10³. It takes about 10⁶ yr for the core field to diffuse out (Beloborodov & Li 2016). Therefore, assuming each patch outburst randomly connects over 10⁶ yr, the mean outburst interval is 10⁶/10³ ∼ 10³ yr, which is consistent with the rate of ¹⁴C events. Assuming there are about 10⁴ young pulsars in the Milky Way galaxy, a typical rate of pulsar outbursts is ten per year, which could be detectable by Fermi/GBM. Obviously, there is lack of detection of such pulsar outbursts. The possible reason is the beaming effect of pulsar outbursts. Interestingly, a handful of flares from non-SGR sources have been reported (Tavani et al. 2011; Tam et al. 2018), but these amount to three over the past 25 yr of high-energy observations.

Below, we compare the physical properties of the SGR 1806−20 giant outburst with our proposed pulsar outburst. For the SGR 1806−20 giant outburst, the photon spectrum of the hard spike for giant outburst is uncertain. It can be described by a blackbody spectrum (Hurley et al. 2005) or a non-thermal spectrum (Palmer et al. 2005). For pulsar outbursts, we use the standard afterglow model of GRBs. The high-energy photon spectrum is power-law form with index about −1.1 (Sari et al. 1998). The isotropic gamma-ray energy of the SGR 1806−20 giant outburst is about 3.7 × 10⁴⁶ erg in the spike. It must be noted that the spike’s intensity drove all X-ray and gamma-ray detectors into saturation (Hurley et al. 2005; Palmer et al. 2005). Therefore, the above gamma-ray energy may be a lower limit. The peak gamma-ray luminosity of the SGR 1806−20 giant outburst is 2 × 10⁴⁷ erg s⁻¹ (Hurley et al. 2005). From the radio afterglow fitting, the total kinetic energy in the ejecta of SGR 1806−20 giant outburst is about 10⁴⁵ erg (Taylor et al. 2005). The total energy of pulsar outburst in our model is about 10⁴⁷ erg (Equation (1)). The rate of giant outbursts is very uncertain. Lazzati et al. (2005) derived an upper limit for the rate of giant outbursts: <1/130 yr⁻¹ per galaxy. Popov & Stern (2006) found that the rate of giant outburst is less than 10⁻³ yr⁻¹ for the Milky Way galaxy. In order to explain the rate of ¹⁴C rapid-increase events, the rate of pulsar outbursts is about 10⁻³ yr⁻¹. Yamazaki et al. (2005) derived that the half opening angle of the jet of SGR 1806−20 giant outburst is 11^◦ from the light curve. The half opening angle of pulsar outburst used in this paper is 16^◦.

3.1. Transport and Deposition of Atmospheric ¹⁰Be

The electrons produced during the outbursts of pulsars can also reach the Earth and affect the spectrum of cosmic-rays. However, due to scattering on the interstellar turbulence it takes much longer for electrons to arrive in the solar system as compared to gamma-rays. Therefore, the spectrum observed in the Earth vicinity may be formed by multiple eruptions. Given the propagation time is much longer than the period between two successive eruptions, one can ignore variability of the sources and assume that they are almost stationary. One should note, however, that this assumption is only valid for calculations of the spectrum of particles far enough from the source. In the close vicinity of the source, its temporal variability strongly affects the spectrum and spatial distribution of particles.

There are two main irregularities in the spectrum of cosmic-ray electrons: the positron excess observed by PAMELA (Adriani et al. 2009) and AMS-02 (Aguilar et al. 2013) and spectral break in the spectrum observed by AMS-02 (Aguilar et al. 2014), CALET (Adriani et al. 2018), and DAMPE (DAMPE Collaboration et al. 2017). We note that data from AMS-02 (Aguilar et al. 2014) and CALET (Adriani et al. 2018) are different from data obtained by DAMPE (DAMPE Collaboration et al. 2017) as well as data obtained by Fermi-LAT (Abdollahi et al. 2017a). The reason for the discrepancy is unknown. In addition, recent observations of the electron spectrum by VERITAS found no irregularities (Archer et al. 2018); however, their systematic uncertainties are higher than the effect in question. Keeping this in mind, we rely on data obtained by CALET and we estimate the energy requirements necessary to make an observable contribution to the spectrum of electrons.

These two features, namely the positron excess and the spectral break, should not be confused. They share a similar concept implying that there is an additional local source of positrons and electrons. But properties of sources required for positron excess and are different from the ones required for spectral break (Aguilar et al. 2019).

The positron excess was originally attributed to nearby pulsars (see, e.g., Hooper et al. 2009; Yüksel et al. 2009, and references therein). Therefore, below we are going to check if pulsars can also explain the spectral break. Similarly to the papers mentioned above we take Geminga as an example; however, we do not exclude that other pulsars can contribute to the spectrum as well.

In the case of uniform background conditions, the time-dependent spectrum of particles arising from a distant source can be calculated from the diffusion equation

$$\begin{eqnarray}&&\displaystyle \frac{\partial f}{\partial t}-{\rm{\nabla }}\left(D{\rm{\nabla }}f\right)+\displaystyle \frac{\partial }{\partial E}\left(\displaystyle \frac{{dE}}{{dt}}f\right)=Q(E,{\boldsymbol{r}},t),\end{eqnarray} \tag{ 10 }$$

where dE/dt is the rate of energy losses, D is the diffusion coefficient, and Q describes spatial, temporal, and spectral properties of the source.

If the medium around the source is uniform and the source is small enough to be considered point-like, one can use a spherically symmetric geometry. After introducing new variables

$$\begin{eqnarray}&&\tau (E,{E}_{0})={\int }_{{E}_{0}}^{E}\displaystyle \frac{{dE}}{{dE}/{dt}}\quad {\rm{and}}\quad \lambda ={\int }_{{E}_{0}}^{E}\displaystyle \frac{D(E)}{{dE}/{dt}}{dE},\end{eqnarray} \tag{ 11 }$$

the Green function can be expressed as

$$\begin{eqnarray}&&{G}_{k}({\boldsymbol{r}},E,t;{E}_{0},{t}_{k})=\displaystyle \frac{| {dE}/{dt}({E}_{0})| }{| {dE}/{dt}(E)| }\displaystyle \frac{\delta (t-{t}_{0}-\tau )}{{\left(4\pi \lambda \right)}^{3/2}}\exp \left[-\displaystyle \frac{{{\boldsymbol{r}}}^{2}}{4\lambda }\right],\end{eqnarray} \tag{ 12 }$$

and the solution of the equation is

$$\begin{eqnarray}&&f({\boldsymbol{r}},E,t)={\int }_{0}^{t}{{dt}}_{0}{\int }_{E}^{\infty }{{dE}}_{0}Q({E}_{0},{t}_{0}){G}_{k}({\boldsymbol{r}},E,t;{E}_{0},{t}_{0}).\end{eqnarray} \tag{ 13 }$$

If a stationary source ($\tfrac{\partial Q}{\partial t}=0$) is located close enough to the observer for energy losses to be irrelevant (r² ≪ λ), but far enough for diffusion approximation to be valid, spatial distribution of the particles can be approximated by

$$\begin{eqnarray}&&f(r)\propto {r}^{-1}.\end{eqnarray} \tag{ 14 }$$

In this case the differential flux of particles through a sphere of radius r₀ is

$$\begin{eqnarray}&&S(E)=4\pi {r}_{0}^{2}D\displaystyle \frac{\partial f}{\partial r}=4\pi {{Dr}}_{0}f({r}_{0},E).\end{eqnarray} \tag{ 15 }$$

There are two ways of estimating the magnitude of f, and therefore calculating the contribution of the source into the cosmic-rays spectrum. The most natural way is to limit the total power of source Q. However, in the past the total energy output of pulsar may have been significantly higher. Therefore, the pulsar may have been producing many more electrons in the past.

The second way, which takes into account past history, is to use gamma-ray data. Recent observations by HAWC (Abeysekara et al. 2017) and Milagro (Abdo et al. 2009) revealed an extended gamma-ray emission around two pulsars, Geminga and B0656+14. This emission is most likely leptonic and therefore its existence confirms that pulsars may act like sources of electrons and positrons.

The spatial extent of the emission, however, is very narrow (≈10–20 pc) and is of an order of mean free path of TeV electrons if the diffusion coefficient is of the order of the intragalactic one. In order to reproduce the observed halo, it is necessary to confine electrons and positrons within this distance, and confinement radius usually significantly exceeds the mean free path. Therefore, it is necessary to assume that the value of diffusion coefficient in the neighborhood of the pulsars is at least 100 times smaller than in the Galaxy and that it is similar to that predicted for standard Bohmian diffusion (Hooper et al. 2017).

The problem is that for D ≈ 10²⁷ cm² s⁻¹ it requires about 6 Myr for particles to cross 200 pc, which significantly exceeds the age of Geminga. To solve this problem one may assume that there is an additional convective transport of electrons (Hooper et al. 2017) or that the area of low diffusion coefficient is of finite size (Fang et al. 2018).

The two-diffusion model can be approximated in the following way. Let us assume that the diffusion coefficient within radius r_D around the pulsar is D₁, while the diffusion coefficient in ISM is D₀. According to Abeysekara et al. (2017), D₁/D₀ ≈ 0.01. Assuming that energy losses are not essential, one can obtain that

$$\begin{eqnarray}f(r)\propto \left\{\begin{array}{ll}{f}_{0}\tfrac{{D}_{0}}{{D}_{1}}\left(\tfrac{{r}_{D}}{r}-1+\tfrac{{D}_{1}}{{D}_{0}}\right), & {\rm{for}}\ r\lt {r}_{d}\\ {f}_{0}\tfrac{{r}_{D}}{r}. & {\rm{for}}\ r\gt {r}_{d}.\end{array}\right.\end{eqnarray} \tag{ 16 }$$

On the other hand, column density of the gamma-ray emitting electrons is

$$\begin{eqnarray}&&{N}_{e}=\int f(r){dl}=\displaystyle \frac{{D}_{0}}{{D}_{1}}{f}_{0}{r}_{D}\left\{\mathrm{ln}\left(\displaystyle \frac{{r}_{\max }}{{r}_{\min }}\right)-1+\displaystyle \frac{{D}_{1}}{{D}_{0}}\right\},\end{eqnarray} \tag{ 17 }$$

where radius of confinement of electrons responsible for gamma-ray emission observed by HAWC is ${r}_{\max }\approx 10$ pc and r_min is of an order of mean free path. Therefore, ${r}_{\min }\,\approx 3{D}_{1}/c\approx 0.03\mbox{--}0.1$ pc.

Surface brightness of the gamma-ray emission can be obtained by the integration of the column density of the electrons with inverse-Compton cross-section

$$\begin{eqnarray}&&{I}_{\gamma }=\displaystyle \frac{1}{4\pi }\int {dE}\int {{dE}}_{\mathrm{soft}}{N}_{e}(E)c\rho ({E}_{\mathrm{soft}})\displaystyle \frac{d{\sigma }_{\mathrm{IC}}(E,{E}_{\gamma },{E}_{\mathrm{soft}})}{{{dE}}_{\gamma }},\end{eqnarray} \tag{ 18 }$$

where E_soft and ρ(E_soft) are the energy and density of soft photons, respectively.

Assuming that the spectrum of electrons is a power law with an index of −2, like we already did before, one can estimate from gamma-ray data that the flux of 1 TeV electrons at the boundary r = r_D is equal to

$$\begin{eqnarray}\begin{array}{rcl}{E}^{3}{j}_{0}(1\,{\rm{TeV}},{r}_{D}) & = & 145\,{{\rm{m}}}^{-2}\,{{\rm{sr}}}^{-1}\,{{\rm{s}}}^{-1}\,{{\rm{GeV}}}^{2}\\ & & \cdot \,{\left(\displaystyle \frac{{D}_{0}/{D}_{1}}{100}\right)}^{-1}{\left(\displaystyle \frac{{r}_{D}}{10{\rm{pc}}}\right)}^{-1}.\end{array}\end{eqnarray} \tag{ 19 }$$

At the position of Earth the flux of electrons is ${j}_{0}({r}_{0}){r}_{0}\,\approx {j}_{0}({r}_{D}){r}_{D}$, and therefore

$$\begin{eqnarray}\begin{array}{rcl}{E}^{3}{j}_{0}(1\,{\rm{TeV}},{r}_{0}) & \approx & 9\,{{\rm{m}}}^{-2}\,{{\rm{sr}}}^{-1}\,{{\rm{s}}}^{-1}\,{{\rm{GeV}}}^{2}\\ & & \cdot \,{\left(\displaystyle \frac{{D}_{0}/{D}_{1}}{100}\right)}^{-1}{\left(\displaystyle \frac{{r}_{0}}{200{\rm{pc}}}\right)}^{-1},\end{array}\end{eqnarray} \tag{ 20 }$$

which is lower than the value observed by CALET by a factor of 3–4. One should note, however, that we completely ignored energy losses. Therefore, more accurate calculations of the electrons flux will most likely produce lower values. Therefore, Geminga alone cannot explain the spectral break observed by CALET. Yet, it potentially has the power to explain positron excess observed by AMS-02, but one needs to assume that it produces electrons and positrons in a 1:1 ratio, and that was successfully demonstrated by Hooper et al. (2017) and Fang et al. (2018).

One should notice that TeV gamma-ray emission requires high-energy electrons. However, the electrons produced during the outburst have relatively low energy (see Equation (3)) and therefore it is unlikely that they will be visible. In this case we can restrict the amount of electrons using the total power of the source.

As discussed above, the spectrum of electrons from the relativistic outflow is a power law and is ∝E⁻². If we know at least one point of differential spectrum f(E₀, r₀) and energy losses of electrons with energy E₀ are negligible, one can estimate the total power of the source as

$$\begin{eqnarray}&&\begin{array}{l}W=4\pi {{Dr}}_{0}\int {dE}\,{Ef}({E}_{0},{r}_{0}){\left(\displaystyle \frac{E}{{E}_{0}}\right)}^{-2}\\ =\,4\pi {{Dr}}_{0}{E}_{0}^{2}f({E}_{0},{r}_{0})\mathrm{ln}\left(\displaystyle \frac{{E}_{\max }}{{E}_{\min }}\right),\end{array}\end{eqnarray} \tag{ 21 }$$

where D is the diffusion coefficient.

According to the data of CALET (Adriani et al. 2018), there is a feature in the spectrum of electrons with a magnitude of about 30 m⁻² sr⁻¹ s⁻¹ GeV² at the energy of about E₀ = 1.2 TeV. We can obtain

$$\begin{eqnarray}&&W=2.3\times {10}^{34}\,\mathrm{erg}\,{{\rm{s}}}^{-1}\,\left(\displaystyle \frac{D}{3\times {10}^{29}\,{{\rm{cm}}}^{2}\,{{\rm{s}}}^{-1}}\right)\left(\displaystyle \frac{{r}_{0}}{150\,{\rm{pc}}}\right),\end{eqnarray} \tag{ 22 }$$

where we assumed that E_max = 10³E_min. Note that in the case of Geminga (r₀ ≈ 200) pc, the total power required $\dot{W}\approx 3\times {10}^{34}$ erg s⁻¹ is comparable to the total spin-down power of the pulsar. Therefore, Geminga can contribute only partially to the observed spectral break. The positron excess on the other hand requires less energy and can be completely supplied by spin-down power alone. This point agrees with the conclusions we made above, which were based on gamma-ray data.

The magnetic outbursts may produce up to 3 × 10³⁵ erg s⁻¹ (assuming 10% acceleration efficiency), which is more than enough to reproduce the break. This conclusion is also in agreement with the fact that positron excess and spectral break are of a different nature.

The expected contribution of the nearby pulsar to the total flux of electrons and positrons obtained near the Earth is shown in Figure 1. The spectrum was calculated using Equation (13) and takes into account the finite age of the pulsar (was assumed of order of 4 × 10⁵ yr) and energy losses. Therefore, total power required to reproduce the observed irregularity equals to W = 3.6 × 10³⁴ erg s⁻¹ for distance of r₀ = 150 pc.

Zoom In Zoom Out Reset image size

If we rely on data obtained by DAMPE, we need a slightly higher power of electrons source. The possible fit to DAMPE data is shown in Figure 1 by a dashed–dotted line. The total power required is W = 5.5 × 10³⁴ erg s⁻¹ for distance of r₀ = 150 pc.

Since the source locates relatively close to the Earth, it should produce a noticeable anisotropy in the arrival directions of electrons. It also may be useful for differentiation between the multi-burst model and the stationary model, since spectral shape in both cases is effectively the same.

The dipole anisotropy in the arrival direction of electrons can be estimated as

$$\begin{eqnarray}&&a=\displaystyle \frac{3D}{c}\displaystyle \frac{| \partial f/\partial r| }{{f}_{\mathrm{tot}}},\end{eqnarray} \tag{ 23 }$$

where f_tot is total spectrum of electrons. If the source is stationary and energy losses are negligible, the expression can be simplified to

$$\begin{eqnarray}&&a=\displaystyle \frac{3D}{{{cr}}_{0}}\displaystyle \frac{{f}_{0}(E,{r}_{0})}{{f}_{\mathrm{tot}}(E,{r}_{0})}.\end{eqnarray} \tag{ 24 }$$

Therefore, a(E) can be easily estimated if we know the energy dependence of diffusion coefficient D(E), total spectrum of electrons ${f}_{\mathrm{tot}}(E,{r}_{0})$, and contribution of local source f₀(E, r₀). In our calculations we assumed that $D={D}_{4}{\left(\tfrac{E}{4\mathrm{GeV}}\right)}^{0.3}$ (Ackermann et al. 2012).

In the case of a non-stationary source that injects electrons with a series of bursts, a solution of Equation (13) is required. In addition to the injections responsible for the events that happened during AD 994, AD 775, BC 660, and BC 3371, we also assumed that the source was stationary until 6000 yr ago to take into account all possible injections that happened in the past. The total power of each injection was assumed to be 2 × 10⁴⁵ erg.

The expected dipole anisotropy with experimental upper limits obtained by Fermi-LAT (Abdollahi et al. 2017b) is shown in Figure 2. One can see that for the burst-like source expected anisotropy is higher and also the dependence on the energy is steeper. This difference may be the way to differentiate the two given models.

Zoom In Zoom Out Reset image size

We also search for the pulsars in our neighborhood. There are two criteria to select candidates. First, the distance of pulsars must be less than 200 pc. Second, due to the outburst, X-ray emissions of pulsars should be detected as residual emission of crust heating. We start by using the Australia Telescope National Facility Pulsar Catalog (version 1.59; Manchester et al. 2005) to compile a list of candidates at distance ≤200 pc from the Sun. Although the detection of a TeV gamma-ray halo around PSR B0656+14 provides good reason to consider it a source of energetic electrons and positrons (Abeysekara et al. 2017), the distance of PSR B0656+14 is ${288}_{-27}^{+33}$ pc (Brisken et al. 2003); therefore, we exclude it. Among these candidates, we exclude the millisecond pulsars (i.e., P ≤ 20 ms) as there is no evidence that outburst from these old objects can attain the required level of energy release. For the remaining candidates, if any of these sources has exhibited such magnetar-like outbursts, X-rays should be powered by crust heating. With these selection criteria, we narrow down to four nearby X-ray emitting pulsars as promising candidates, namely PSR B1055−52, PSR B0834+06, Geminga, and RX J1856.5−3754. Due to the pulsar beaming effect, it must be noted that our census of the young nearby pulsars is not complete. It may be improved via HAWC TeV gamma-ray survey observations in the near future, i.e., via further TeV halo discoveries. Sudoh et al. (2019) have predicted that HAWC can eventually detect up to ∼80 TeV halos. On the other hand, with the upcoming Galactic plane survey with Cerenkov Telescope Array, ∼30–160 new TeV halos can possibly be uncovered. Also, it has been speculated that a few tens of currently unidentified sources in TeV catalogs can possibly be TeV halos (Sudoh et al. 2019). From their observational properties, PSR B0834+06, Geminga, and RX J1856.5−3754 could have had outbursts in the past. Especially, the spin-down luminosity of the X-ray dim isolated neutron star (XDINS) RX J1856.5–3754 at ∼108–134 pc as inferred from period P, and the spin-down rate $\dot{P}$ is $\dot{E}\sim 3\times {10}^{30}$ erg s⁻¹, which is smaller, by an order of magnitude, than its X-ray luminosity ${L}_{{\rm{x}}}\sim 3\times {10}^{31}$ erg s⁻¹, estimated at a distance of $d\sim 120$ pc (Burwitz et al. 2003). This clearly indicates that the X-ray emission from RX J1856.5−3754 is not powered by its rotational energy. On the other hand, the energy stored in its strong magnetic field is sufficient to power the observed X-rays. Comparing the observed surface temperature of RX J1856.5−3754 with various neutron star-cooling curves (Aguilera et al. 2008), it is far higher than the expected temperature at its characteristic age, even with the effect of Joule heating taken into account. Although the discrepancy can possibly be reconciled if its true age, 4 × 10⁵ yr, is closer to its kinematic estimate (Mignani et al. 2017), it may also indicate this neutron star might have been heated by a certain event, such as an outburst, in the past.

Recently, on top of the well-studied thermal radiation from a two-temperature blackbody, an X-ray excess in the keV regime has been discovered from RX J1856.5−3754 (Yoneyama et al. 2017). While the nature of this excess remains uncertain, a possible scenario of resonant cyclotron scattering has been proposed (Yoneyama et al. 2017). This suggests the additional component arises from the scattering between the thermal surface photons and plasma flowing along the close magnetic field lines, which have also been observed in the X-ray spectra of magnetars (Rea et al. 2008). If this explanation is correct, it will further indicate the presence of magnetic activities on RX J1856.5−3754 (Beloborodov 2013). X-ray polarization measurement in the near future can provide a test for this scenario (Weisskopf 2018).

PSR B1055−52 is one of the earliest detected γ-ray pulsars (Fierro et al. 1993). Its dispersion measure suggests a distance of only d ∼ 90 pc. It has a spin period of P = 0.197 s and a spin-down rate of $\dot{P}=5.8\times {10}^{-15}$ s s⁻¹. These imply a characteristic age of $\tau \sim 5\times {10}^{5}$ yr and a dipolar surface magnetic field strength of ${B}_{s}\sim 1.1\times {10}^{12}$ G. The X-ray emission of PSR B1055−52 consists of two thermal components and a non-thermal magnetospheric component (De Luca et al. 2005). The high temperature (T ∼ 2 × 10⁶ K) thermal component originates from a hot spot with a radius of ∼460 m, which is consistent with the conventional polar cap size as defined by the footprints of dipolar magnetic field lines on the surface. Hence, there is no hint of the presence of stronger multipolar magnetic field on the surface. We do not find any compelling feature from PSR B1055−52 that makes it a possible contributor to the episodic cosmic-ray increment.

The distance of PSR B0834+06 is estimated to be d ∼ 190 pc, as derived from its dispersion measure. Its spin parameters P = 1.27 s and $\dot{P}=6.8\times {10}^{-15}$ s s⁻¹ imply a characteristic age of τ ∼ 3 × 10⁶ yr and a dipolar surface magnetic field strength of B_s ∼ 3 × 10¹² G. The X-ray emission from PSR B0834+06 is found to be thermal dominant and can be described by a hot spot on the stellar surface with a temperature of T ∼ 2 × 10⁶ K (Gil et al. 2008). The size of the hot spot is much smaller than the polar cap defined by the dipolar field lines. This suggests the surface field might be multipolar and could be ∼50 times stronger than the dipolar estimate given above. A field strength of the order of ∼10¹⁴ G makes outburst activities possible in PSR B0834+06, though there is so far no observational indication for any transient behavior from this pulsar.

Geminga (PSR J0633+1746) is the second strongest persistent γ-ray source in the sky. The distance of this radio-quiet γ-ray pulsar is in the range of d ∼ 190–370 pc, as determined by the parallax measurement of its optical counterpart (Faherty et al. 2007). Its spin period and the proper-motion corrected spin-down rate (P = 0.237 s, $\dot{P}=1.1\times {10}^{-14}$ s s⁻¹) suggest a characteristic age of τ ∼ 3 × 10⁵ yr and a dipolar surface field strength of B_s ∼ 1.6 × 10¹² G. X-ray observations have identified a small hot spot on its surface (Caraveo et al. 2004). The area of this hot spot is ∼25 times smaller than the dipolar polar cap. Although the discrepancy can be due to a geometrical viewing effect, it can also possibly be a result of multipolar surface field. In such a scenario, the surface field strength can attain a value of ∼4 × 10¹³ G, which is comparable to the field of low magnetic field magnetars (Turolla & Esposito 2013). Some previous studies of the distribution of thermal emission on the surface of Geminga have even suggested the thermal X-rays can possibly be a “mark” left by tectonic activity in the past (Halpern & Ruderman 1993; Page et al. 1995).

Geminga is also powering a complex pulsar wind nebula (PWN), which can be resolved in several spatial components using X-ray (Hui et al. 2017; Posselt et al. 2017). In a multi-epoch analysis, Hui et al. (2017) have discovered fast X-ray variabilities in various components of Geminga PWN. In the region of the post-shock flow (see Figure 2 in Hui et al. 2017), clumps of X-rays have been found to appear in several epochs. This resembles the magnetic reconnection in the Sun that gives rise to the eruptive coronal mass ejection (see Figure 1 in Lin et al. 2005). Energetic charged particles can possibly be ejected by Geminga PWN in a similar way.

Different from the above three candidates, RX J1856.5−3754 is not a rotation-powered pulsar. It belongs to a class of XDINSs (see Haberl 2013 for a detailed review). Its distance has been bracketed in a range of ∼108–134 pc (Mignani et al. 2017). The X-ray emission from RX J1856.5−3754 is purely thermal (Sartore et al. 2012). While the thermal X-rays of the hot spot from the rotation-powered pulsars likely comes from the bombardment of the back-flow current from the accelerating regions in their magnetosphere (Cheng & Zhang 1999), the thermal radiation from RX J1856.5−3754 should have a different origin. A low amplitude X-ray pulsation at a period of P = 7.06 s (Tiengo & Mereghetti 2007) and its subsequent spin-down rate $\dot{P}=3.0\times {10}^{-14}$ s s⁻¹ (van Kerkwijk & Kaplan 2008) have been found. These imply a surface dipolar field strength of ${B}_{s}\sim 1.5\times {10}^{13}$ G and a characteristic age of $\tau \sim 4\times {10}^{6}$ yr. Recently, optical polarization measurements have confirmed the presence of a strong magnetic field on RX J1856.5–3754 (Mignani et al. 2017). Its spin-down luminosity as inferred from P and $\dot{P}$ is $\dot{E}\sim 3\times {10}^{30}$ erg s⁻¹, which is smaller than its X-ray luminosity L_x ∼ 3 × 10³¹ erg s⁻¹ (Burwitz et al. 2003) by an order of magnitude. On the other hand, the energy stored in its strong magnetic field is sufficient to power the observed X-rays. Comparing the observed surface temperature of RX J1856.5−3754 with various neutron star-cooling curves (Aguilera et al. 2008), it is far higher than the expected temperature at its characteristic age even with the effect of Joule heating taken into account. Although the discrepancy can possibly be reconciled if its true age is closer to its kinematic estimate (4 × 10⁵ yr; Mignani et al. 2017), it may also indicate this neutron star may have been heated by a certain event, such as an outburst, in the past. The heat can come from internal heating by plastic flow (Li et al. 2016) or external heating by Alfvén waves launched in the magnetosphere (Li & Beloborodov 2015). Such heating is transient and produces delayed thermal afterglow. For its impact on the long term cooling of the neutron star, we can estimate the increased thermal luminosity by

$$\begin{eqnarray}&&{ \mathcal L }\sim E/{t}_{\mathrm{amb}}\sim {10}^{32}{L}_{5}^{13/5}\delta {B}_{16}^{8/5}{B}_{13}^{8/5}{\left(\rho /{\rho }_{\mathrm{nuc}}\right)}^{6/5}\,\mathrm{erg}\,{{\rm{s}}}^{-1}.\end{eqnarray} \tag{ 25 }$$

The currently known TeV halos around pulsars are resulted from the interactions between the continuous outflow of electrons/positrons with the ambient photon field. While the case that we consider here was resulted from the episodic injection by RX J1856.5−3754 in the past, the diffusion can make the flux of its possible relic halo below the detection limit of the existing TeV observing facilities.

In this paper, we propose that the outbursts of nearby pulsars (d < 200 pc) can cause the ¹⁴C events, and contribute a significant fraction of cosmic electron spectrum in the trillion electron volts energy range. The relativistic outflow produced by outbursts encounters the ISM, and the external shock is generated, including forward and reverse shocks, which can accelerate electrons. The synchrotron emission from accelerated electrons generates high-energy photons, which interact with the atmosphere producing the cosmogenic nuclide ¹⁴C.

Meanwhile, these same relativistic electrons contribute a significant fraction of the cosmic electron spectrum in the trillion electron volts energy range. They may be responsible for the spectral break recently observed in the spectrum. We also review the properties of cataloged pulsars in the solar neighborhood, and find that PSR B0834+06, Geminga, and RX J1856.5−3754 could have had outburst activities in the past.

Our model can be tested with future experimental facts and astronomical observations. First, the short time and the narrow beam of the outburst would have an influence on the time and intensity of the ¹⁴C increase at different locations on Earth. Second, our model predicts a hard energy dependence of dipole anisotropy in the arrival direction of electrons. Third, a non-thermal component is expected in soft X-rays band. Since these outburst events can significantly affect our environment, monitoring the nearby pulsar outbursts is important in the future.

We thank the referee for detailed and very constructive suggestions that have allowed us to improve our manuscript. We also thank Xiang Li and Lei Feng for helpful discussion on GEANT4 code, and Prof. Kevin Mackeown for language editing. F.Y.W. is supported by the National Natural Science Foundation of China (grant U1831207). D.O.C. is supported in part by the grant RFBR 18-02-00075 and by the foundation for the advancement of theoretical physics and mathematics “BASIS.” C.Y.H. is supported by the National Research Foundation of Korea grant 2016R1A5A1013277. K.S.C. is supported by a GRF grant under 17302315.