A new approach to the GeV flare of PSR B1259-63/LS2883

TL;DR

This paper proposes a novel model where GeV gamma-ray flares in PSR B1259-63/LS2883 are caused by inverse Compton scattering of photons from an accretion disk formed by captured matter, explaining observed emissions.

Contribution

It introduces a new accretion disk-based inverse Compton scattering model to explain the GeV flares in PSR B1259-63/LS2883, linking disk formation to flare timing.

Findings

Model reproduces observed SEDs and light curves.

Accretion disk formation is linked to disk-crossing events.

Explains the timing and intensity of gamma-ray flares.

Abstract

PSR B1259-63/LS2883 is a binary system composed of a pulsar and a Be star. The Be star has an equatorial circumstellar disk (CD). The {\it Fermi} satellite discovered unexpected gamma-ray flares around 30 days after the last two periastron passages. The origin of the flares remain puzzling. In this work, we explore the possibility that, the GeV flares are consequences of inverse Compton-scattering of soft photons by the pulsar wind. The soft photons are from an accretion disk around the pulsar, which is composed by the matter from CD captured by the pulsar's gravity at disk-crossing before the periastron. At the other disk-crossing after the periastron, the density of the CD is not high enough so that accretion is prevented by the pulsar wind shock. This model can reproduce the observed SEDs and light curves satisfactorily.

Figures

Tables

Table 1. Table 1: Different mass transfer scenarios under variance on the CD parameters. v ¯ rel , ± subscript ¯ 𝑣 rel plus-or-minus \bar{v}_{\rm{rel},\pm} and T ± subscript 𝑇 plus-or-minus T_{\pm} are the averaged relative velocity and the duration of mass transfer at the “ ± plus-or-minus \pm ” crossing respectively. Configuration 1: θ n = 40 ∘ subscript 𝜃 n superscript 40 \theta_{\rm{n}}=40^{\circ} , ϕ n = 24 ∘ subscript italic-ϕ n superscript 24 \phi_{\rm{n}}=24^{\circ} ; Configuration 2: θ n = 45 ∘ subscript 𝜃 n superscript 45 \theta_{\rm{n}}=45^{\circ} , ϕ n = 19 ∘ subscript italic-ϕ n superscript 19 \phi_{\rm{n}}=19^{\circ} ; Configuration 3: θ n = 35 ∘ subscript 𝜃 n superscript 35 \theta_{\rm{n}}=35^{\circ} , ϕ n = 14 ∘ subscript italic-ϕ n superscript 14 \phi_{\rm{n}}=14^{\circ} . Three values of each calculated parameters correspond to n = 3 , n = 3.3 , n = 3.5 formulae-sequence 𝑛 3 formulae-sequence 𝑛 3.3 𝑛 3.5 n=3,n=3.3,n=3.5 . A bar in the value means mass transfer condition is not satisfied under the corresponding CD configuration. The bold values are used in the following parts of the paper.

		Configuration 1	Configuration 2	Configuration 3
Prograde	${\overline{v}}_{\mathrm{rel},+}$ ($\times {10}^7$ cm/s)	1.0	1.1	1.2
	${\overline{v}}_{\mathrm{rel},-}$	1.3	1.3	1.4
	$T_{+}$ (days)	4.9, –, –	7.7, –, –	8.7, –, –
	$T_{-}$	8.6, 5.9, –	9.5, 5.1, –	10.2, 3.5, –
	$M_{\mathrm{tot},+}$ ($\eta \times {10}^{24}$ g)	0.7, –, –	0.8, –, –	1.2, –, –
	$M_{\mathrm{tot},-}$	3.4, 1, –	2.8, 0.7, –	2.4, 0.4, –
	$R_{\mathrm{circ},+}$ ($\times {10}^{10}$ cm)	1.5, –, –	1.1, –, –	0.85, –, –
	$R_{\mathrm{circ},-}$	1.8, 2.2, –	1.2, 1.4, –	0.8, 1, –
Retrograde	${\overline{v}}_{\mathrm{rel},+}$	1.7	1.8	1.9
	${\overline{v}}_{\mathrm{rel},-}$	2.7	2.6	2.4
	$T_{+}$	–	–	–
	$T_{-}$	2.1, –, –	–	–
	$M_{\mathrm{tot},+}$	–	–	–
	$M_{\mathrm{tot},-}$	0.1, –, –	–	–
	$R_{\mathrm{circ},+}$	–	–	–
	$R_{\mathrm{circ},-}$	0.02, –, –	–	–

Equations

$${\mathbf{n}}_{\mathrm{cd}}=(\sin {\theta }_{\mathrm{n}}\cos {\varphi }_{\mathrm{n}},\sin {\theta }_{\mathrm{n}}\sin {\varphi }_{\mathrm{n}},\cos {\theta }_{\mathrm{n}}).$$

$$\begin{array}[]{l}v_{\rm{p},x}=e\sqrt{p\mu}\sin\phi\cos\phi/p-\sin\phi\sqrt{p\mu}/r\\ v_{\rm{p},y}=e\sqrt{p\mu}\sin^{2}\phi/p+\cos\phi\sqrt{p\mu}/r\\ \end{array}$$

$${\mathbf{v}}_{\mathrm{cd}}=\sqrt{\frac{G{M}_{\star }}{r}}{\mathbf{n}}_{\mathrm{cd}}\times \mathbf{r}/r.$$

$$r_{\mathrm{BH}}=\frac{2G{M}_{\mathrm{p}}}{v_{\mathrm{rel}}^2},$$

$$r_{\mathrm{s}}=\sqrt{\frac{L_{\mathrm{spin}}}{4\pi {\rho }_{\mathrm{cd}}{v}_{\mathrm{rel}}^{2}c}},$$

$${\rho }_{\mathrm{cd}}$$

$${\dot{M}}_{\mathrm{trans},\mathrm{BH}}=\pi r_{\mathrm{BH}}^2{\rho }_{\mathrm{cd}}{v}_{\mathrm{rel}}.$$

$$l(t)=\frac{(GM_{\rm{p}})^{2}}{v^{3}_{\rm{rel}}}\big{(}\frac{|\nabla v_{\rm{rel}}|}{v_{\rm{rel}}}+\frac{|\nabla\rho_{\rm{cd}}|}{\rho_{\rm{cd}}}\big{)}.$$

$$M_{\mathrm{tot}}=\eta \int {\dot{M}}_{\mathrm{trans},\mathrm{BH}}dt,$$

$$\bar{l}=\eta \int l(t){\dot{M}}_{\mathrm{trans},\mathrm{BH}}dt/M_{\mathrm{tot}}.$$

$$r_{\mathrm{A}}=5.1\times 10^8{\dot{M}}_{\mathrm{acc},16}^{-2/7}{m}_{\mathrm{p}}^{-1/7}{\mu }_{30}^{4/7}\text{cm}$$

$${\dot{M}}_{\mathrm{acc}}={\dot{M}}_{\mathrm{acc},0},(r_{\mathrm{in}}>r_{\mathrm{M}})$$

$$\dot{M}_{\rm{acc}}=\dot{M}_{\rm{acc,0}}\big{(}\frac{t}{\tau}\big{)}^{-\beta}\quad,\quad(r_{\rm{in}}=r_{\rm{M}})$$

$$r_{\mathrm{in}}=R_{\mathrm{circ}}-v_{r}t,(t<\tau )$$

$$r_{\mathrm{in}}=r_{\mathrm{M}},(t\ge \tau ),$$

$$M_{\mathrm{tot}}={\int }_0^{\infty }{\dot{M}}_{\mathrm{acc}}dt,$$

$${\dot{M}}_{\mathrm{acc},0}=\frac{(\beta -1)M_{\mathrm{tot}}}{\beta \tau }.$$

$$\tau$$

$$\frac{\dot{M}_{\rm{acc,0}}}{10^{16}\,\rm{g}/\rm{s}}=\big{[}2.7\alpha^{4/5}\frac{\beta-1}{\beta}M_{\rm{tot,22}}m_{\rm{p}}^{-1/4}R_{\rm{circ,10}}^{-5/4}\big{]}^{10/7}$$

$$T=1.4\times 10^4{\alpha }^{-1/5}{\dot{M}}_{\mathrm{acc},16}^{3/10}{m}_{\mathrm{p}}^{1/4}{r}_{\mathrm{in},10}^{-3/4}\text{K},$$

$$E_{\gamma }\approx 300(kT/1\text{eV}){\Gamma }_{0,4}^2\text{MeV},$$

$$E_{\gamma }\approx 5.11{\Gamma }_{0,4}\text{GeV},$$

$${\dot{N}}_{\gamma }=\int d{N}_{\mathrm{e}}(l){\int }_0^{\infty }\frac{2}{{\Gamma }^2{ϵ}_{\mathrm{ph}}}\pi r_0^{2}cn({ϵ}_{\mathrm{ph}},l)f(x)d{ϵ}_{\mathrm{ph}},$$

$$d{N}_{\mathrm{e}}(l)=\frac{L_{\mathrm{spin}}}{\Gamma m_{\mathrm{e}}{c}^2}\frac{dl}{c}.$$

$$n({ϵ}_{\mathrm{ph}},l)=\frac{4\pi l}{h^3{c}^3}{\int }_{R_{\mathrm{in}}}^{R_{\mathrm{out}}}\frac{{ϵ}_{\mathrm{ph}}^{2}RdR}{(R^2+l^2{)}^{3/2}(\exp \frac{{ϵ}_{\mathrm{ph}}}{kT(R)}-1)}.$$

$${\dot{N}}_{\gamma }=\frac{8{\pi }^2{r}_0^2}{h^3{c}^3}\frac{L_{\mathrm{spin}}}{{\Gamma }^3{m}_{\mathrm{e}}{c}^2}\int {ϵ}_{\mathrm{ph}}f(x){\int }_{R_{\mathrm{in}}}^{R_{\mathrm{out}}}\frac{dR}{\exp \frac{{ϵ}_{\mathrm{ph}}}{kT(R)}-1}d{ϵ}_{\mathrm{ph}},$$

$${\dot{N}}_{\gamma }=\frac{8{\pi }^2{r}_0^2}{h^3{c}^3}\frac{L_{\mathrm{spin}}{r}_{\mathrm{in}}}{{\Gamma }^3{m}_{\mathrm{e}}{c}^2}\int \frac{{ϵ}_{\mathrm{ph}}f(x)}{\exp \frac{{ϵ}_{\mathrm{ph}}}{kT(r_{\mathrm{in}})}-1}d{ϵ}_{\mathrm{ph}},$$

$$F_{\gamma }=\frac{{\dot{N}}_{\gamma }}{4\pi D^2},$$

$$F_{\nu }=\frac{4\pi h{\nu }^3\cos i}{c^2{D}^2}{\int }_{R_{\mathrm{in}}}^{R_{\mathrm{out}}}\frac{RdR}{\exp h\nu /kT-1},$$

$$F_{\nu ,\star }=\frac{2h{\nu }^3\pi R_{\star }^2}{c^2{D}^2}\frac{1}{\exp h\nu /k{T}_{\star }-1},$$

Full text

gamma rays: stars–pulsars: individual (PSR B$1259-63$)

1 Introduction

PSR B1259-63/LS 2883 is a binary system composed of a radio pulsar and a massive main sequence Be star (Johnston et al., 1992). The pulsar PSR B1259-63 is a non-recycled, spin-down powered radio pulsar, with a spin period of 47.76 ms, and a period derivative of 2.27$\times 10^{-15}$ (Shannon et al., 2014). The companion star is a Be star with the mass of $\sim 31\,M_{\odot}$, and rotating at a near-break-up rate (Negueruela et al., 2011). The binary orbit is highly eccentric, with the eccentricity $e=0.87$, orbital period $P=1237$ days and semimajor axis $a=7.2$ AU. The stellar wind of the Be star interacts with the pulsar wind fiercely, forming a termination shock front, which is believed to be the origin of un-pulsed, orbital modulated high energy radiations of X-rays (Hirayama et al., 1999; Chernyakova et al., 2006) and TeV gamma-rays (Kirk et al., 1999; Aharonian et al., 2005). The ejected matter at the equator of the Be star builds up a Keplerian decretion disk, which transports the angular momentum of the star outwardly (Rivinius et al., 2013). The obstruction of the circumstellar disk (CD) leads to eclipse of the pulse radio emission near the periastron (Johnston et al., 1992, 1996), and the interaction between the pulsar and the CD gives rise to the double peaks structure of the X-ray light curve (Kong et al., 2011).

During the periastron passage in 2010/2011, the first passage after the launch of the Fermi satellite, Tam et al. (2011) and Abdo et al. (2011) discovered the gamma-rays emission (100 MeV-100 GeV) from the system with the Fermi-LAT detector. Marginal gamma-rays were detected within period around the time of the periastron $t_{\rm{p}}$, before it faded at $\sim t_{\rm{p}}+15$ days. Then at $\sim t_{\rm{p}}+30$ days, an sudden re-emergence of the gamma-ray emission surprisingly occurred, with the flux fast rising up to $20-30$ times of its previous value over a few days, followed by a slow decay about two weeks before its disappear after $~{}t_{\rm{p}}+60$ days. The GeV flare reappeared during the following periastron passage in 2014 (Malyshev et al., 2014; Tam et al., 2014; Wood et al., 2014) at approximate the same orbital phase. The similarities of the flux, on-set phases, duration and sub-scale structures in the light curves between GeV flares in two periastron passages (Tam et al., 2015; Caliandro et al., 2015) suggest its repeating nature and the origin due to the interaction between the pulsar and the circumstellar matter. The cause of the GeV flare is puzzling: as the first-considered radiation mechanism, inverse Compton scattering (IC) of the cool pulsar wind electrons with seed photons from either the Be star or the CD gives the maximum of GeV flux at the periastron, where the external photons field is densest111Considered the anisotropy of IC, the GeV light curve peaks slightly before the periastron.. Whereas the flare occurs $\sim 30$ days after the periastron (Khangulyan et al., 2007; Takata & Taam, 2009; Pétri & Dubus, 2011; van Soelen et al., 2012) ; on the other hand, if we trust the geometry of the CD implied by the X-ray light curve (Chernyakova et al., 2006), then the flare appears when the pulsar has well left the CD. Therefore it is difficult to attribute the flare to the immediate interaction between the pulsar and the CD.

The mystery of the GeV flare attracts many explaining attempts: Khangulyan et al. (2012) still tried to explain the flare as a consequence of cool pulsar wind IC with the soft photons from the CD, but they also took into account of the time variance of the pulsar wind zone length (PWZ) towards the observer. PWZ is the place where the IC process take place. They argued that the rapid rising of the gamma-ray flux corresponds to the fast growth of the PWZ immediately after the pulsar’s exiting from the CD. The main difficulty of this model is insufficient soft photons from CD for IC. As mentioned by the authors, the required target photon luminosity is $\sim 40\%$ of LS 2883. Significant heating of the CD by the passage of the pulsar is supposed by the authors to reconcile this problem.

Using also IC, Dubus & Cerutti (2013) sought for seed photons in X-rays from the synchrotron radiation of the shock-heated electrons. Since in this model the seed photons have higher energy of a few keV, the bulk Lorentz factor of the pulsar wind should be $\sim 500$ to give the observed gamma-ray emission. The Lorentz factor is much less than the typical values of pre-shock pulsar wind (see van Soelen et al. 2012 and references therein). It is difficult for this model to explain the delay between the GeV flare and the X-ray lightcurve maximum.

Another approach to the GeV flare, proposed by Tam et al. (2011) and modeled in detail by Kong et al. (2012), is the Doopler-boosted synchrotron radiation in the shock tail: The stellar wind from the Be star collides with the pulsar wind to form a termination shock. Electrons (and positrons) in the pulsar wind are shock-heated and flow along the shock-tail at a mildly relativistic velocity. The hot electrons radiate away their energy through the synchrotron radiation, which is in X-ray in the flow co-moving frame and is Doopler-boosted to gamma-ray at the rest frame. This model is favored by the observational fact that the GeV flare appears when the pulsar is around the inferior conjunction, where the shock tail is close to the line of sight. However this model predicts that the light curves of the X-ray and gamma-ray peak at the same orbital phase, which violates the observation.

Mochol & Kirk (2013) explored the possibility that, the flare is the result of a transit between the driven reconnection and the electromagnetic precursor scenarios of the termination shock. The above mentioned scenarios are both ways in which the oscillating magnetic field of the rotating pulsar dissipates its energy. The radiation efficiency is low in the driven reconnection scenario, and is high in the electromagnetic precursor scenario. When the condition of transition is satisfied, the radiation efficiency jumps thus gives the observed flare. This model also predicts another GeV flare at a different orbital phase when the transition condition is satisfied again (the peak is not necessary in the GeV band as mentioned by the authors, depending on the wind conditions.) The prediction is not supported yet by observations.

A recent work of de la Cita et al. (2017) studied numerically the impact of stellar wind clumps or inhomogeneities on the high-energy non-thermal radiation in binary system like this. They applied their study to PSR B1259-63/LS2883, and explored the scenario that a dense matter clump impacts on the two-wind interaction region. Although enlightening, their simulation is not able to reproduce the GeV flare satisfactorily.

We explore the possibility in this paper that, when the pulsar passes across the CD, some of the matter in the CD is captured by the gravity of the neutron star. After the pulsar exits the CD, the captured matter forms an accretion disk. Soft photons emitted from the accretion disk are up-scattered by the pulsar wind, results in the observed gamma-ray flares. The gravity capture happens only when the termination shock front are inside the capture radius, inside which the velocity of the CD matter in the pulsar rest frame is less than the escape velocity from the gravity field of the pulsar. If the radius of the termination shock front is larger than the capture radius, the matter flow will be redirected before it enters the capture radius. We should also check whether the captured matter has enough specific angular momenta to form an accretion disk. The detailed calculations in section 2.1 show that the formation of an accretion disk by this scenario is likely.

This paper is organized as follows. In section 2 we describe the model: In section 2.1 we study how the CD matter is captured by the gravity of the pulsar; Then in section 2.2 we describe how an accretion disk is formed from the captured matter, and how the accretion disk evolves with time; In section 2.3 we study the IC process, in which the pulsar wind up-scatters the accretion disk soft photons to gamma-rays. In section 3 we compare our model calculations with observations. We conclude and discuss in section 4.

2 Model

2.1 Accretion from the stellar disc

To describe the motion of the pulsar and the CD matter, we construct the Cartesian coordinates frame as follows: the origin is on the barycenter of Be star and the pulsar, the $x$-axis is towards the periastron, $z$-axis is along the angular momentum of the orbit. The norm vector of the Be star’s CD is

[TABLE]

Where $\theta_{\rm{n}}$ and $\phi_{\rm{n}}$ are the polar angle and the azimuthal angle of the norm vector respectively. The pulsar intersects with the mid-plane of the CD at the true anomaly $\phi_{\rm{cd,\pm}}=\phi_{n}\pm\pi/2$. Based on the Kepler equations, the velocity of the pulsar as a function of the true anomaly $\phi$ is:

[TABLE]

where $r=p/(1+e\cos\phi)$ is the pulsar-star distance, $e$ is the eccentricity of the orbit, $\mu\equiv G(M_{\star}+M_{\rm{p}})$ is the gravity constant times the sum of the total mass of the Be star and the pulsar, $p\equiv a(1-e^{2})$ and $a$ is the semi-major axis. The CD is thought to be Keplerian, therefore we assume the velocity of the disk matter is pure tangential and ignore the outward velocity. As a result, at any point $\mathbf{r}$ the velocity of the disk matter is:

[TABLE]

Here we presume that the rotation of CD is prograde with the orbit of the pulsar. The retrograde case is discussed later. The velocity of the disk matter viewed at the pulsar rest frame is $\mathbf{v_{\rm{rel}}}=\mathbf{v_{\rm{cd}}}-\mathbf{v}_{\rm{p}}$. Any matter with the impact parameter less than the radius

[TABLE]

is thought to be captured by the gravity of the neutron star. $r_{\rm{BH}}$ is the Bondi-Hoyle value (Bondi & Hoyle, 1944) when $v_{\rm{rel}}$ is much larger than the sound speed. Mass transfer via gravity-capture is only possible when the shock front is within $r_{\rm{BH}}$, i.e., $r_{\rm{s}}<r_{\rm{BH}}$ as discussed in above section, where

[TABLE]

and $L_{\rm{spin}}$ is the spin down power of the pulsar, $\rho_{\rm{cd}}$ is the local density of the CD. We adopt the following parameters:

• •

$e=0.87$, $a=7.2\,AU$

• •

$M_{\star}=31\,M_{\odot}$ (Negueruela et al., 2011), $M_{\rm{p}}=1.4\,M_{\odot}$

• •

$\theta_{\rm{cd}}=45^{\circ}$, $\phi_{\rm{cd}}=19^{\circ}$ (Okazaki et al., 2011)

• •

$L_{\rm{spin}}=8\times 10^{35}$ ergs/s,

The density distribution of the decretion disk of the Be stellar has been modeled by previous researchers (Lee et al., 1991; Porter, 1999; Carciofi & Bjorkman, 2006) as:

[TABLE]

where $n$ is in the range $3\sim 3.5$, with 3.5 corresponds to an steady state isothermal outflow, $H$ is the scale height of the disk, $z$ is the cylindrical coordinate in the vertical direction, $\Delta\phi$ is the half-opening angle of the disk projected on the orbital plane, $R_{\star}$ is the radius of the Be star, which is $\sim 10\,R_{\odot}$ (Negueruela et al., 2011). The second part of above equation is the disk density profile projected to the orbital plane, and $\Delta\phi=18.5^{\circ}$ as modeled by Chernyakova et al. (2006). Takata et al. (2012) argued that a large base density ($\rho_{0}\sim 1\times 10^{-9}$ g/cm3) is needed for this system to account for the double peak structure of the X-ray light curve.

In the upper panel of figure 1, we plot the $r_{\rm{s}}$ and $r_{\rm{BH}}$ as function of $\phi$ under different $\rho_{\rm{cd}}$ profile index $n$, with $\rho_{0}=10^{-9}$ g/cm3. $\phi$ is the true anomaly and the zero point is at the periastron. We ignore the perturbation of the density profile of the disk by the pulsar. We can see that the scenario of mass transfer is sensitive to $n$: when $n=3$, mass transfer occurs in both CD crossing; When $n=3.3$, mass transfer only happens during a part of the CD crossing before the periastron (“–” crossing hereafter, the other CD crossing is denoted as “+”); When $n=3.5$, no mass transfer from the CD during both crossings. Changing the $n$ from 3 to 3.5 is equivalent to changing $\rho_{0}$ from $\sim$0.5 to $\sim$2 times the current value. The illustration of the orbit and the position of the disk (grey shade region) is plotted in the bottom panel of figure 1 with $n=3.3$, with the mass transfer region shaded in red.

When the condition that $r_{\rm{s}}<r_{\rm{BH}}$ is met, the Bondi-Hoyle like process transfers matter into the gravity-bounded sphere at a rate:

[TABLE]

Since $v_{\rm{rel}}$ and $\rho_{\rm{cd}}$ are not uniform and have gradient across the radius of the CD, the gravity-captured matter possess the specific angular momentum of (Illarionov & Sunyaev, 1975; Shapiro & Lightman, 1976):

[TABLE]

During the mass transfer process, the total accreted mass is:

[TABLE]

where $\eta$ accounts for the inefficiency of the Bondi-Hoyle accretion. According to numerical study by Blondin (2013), under the influence of the bow shock of the pulsar $\eta$ can be as small as $\sim 1\%$. Here we leave $\eta$ to be determined by fitting to the observation.

After ramming into the capture radius, the gases dissipate their energy and shift to the orbit of lowest energy for a given angular momentum, i.e., a torus of gas at the so called circular radius: $R_{\rm{circ}}\equiv{\bar{l}}^{2}/(GM_{\rm{p}})$, where

[TABLE]

After that, the matter losses angular momentum via viscosity and spread from $R_{\rm{circ}}$ to form an accretion disk.

In table 1 we list different mass transfer scenarios under variance of disk parameters. We see from this table that the mass transfer scenarios are not sensitive to slight changes of disk orientation and inclination (see the definition of three CD configurations in the caption of table 1). However, whether or not the mass transfer occurs depends strongly on the rotation direction of the CD: mass transfer is expected to occur when the CD is prograde, but not retrograde with the orbit. By study the spin-orbit coupling of this system, Shannon et al. (2014) suggested the misalignment of the spin axis of LS 2883 and orbit axis is less than $90^{\circ}$, which supports a prograde CD.

For the purpose of further modeling, we work with $\theta_{\rm{n}}=45^{\circ}\,,\phi_{\rm{n}}=19^{\circ}$ and the assumption that the CD is prograde and $n=3.3$. In this case, $M_{\rm{tot}}=\eta 7\times 10^{23}$ g and $R_{\rm{circ}}=1.4\times 10^{10}$ cm. $R_{\rm{circ}}$ is much larger than the radius of the light cylinder ($7.6\times 10^{7}$ cm), therefore the angular momentum of the accreted matter is enough to form an accretion disk. The mass transfer process lasts $\sim 5$ days, while the flare lasts more than 60 days after the disk passage. Therefore the accretion rate of the disk should be much less than the mass transfer rate. Thus, as an approximation, we consider the gravity-captured matter accumulated around the pulsar within $r_{\rm{BH}}$ during the mass transfer process. When the matter transfer ends, the accumulated matter begins to form an accretion disk.

2.2 Evolution of the accretion disk

As described above, after being captured by the gravity of the neutron star, complex physics processes (including shocks) will convert much of the kinetic energy of the gas into radiation and redistribute the remaining kinetic energy and angular momenta, before a torus of gas is formed at $R_{\rm{circ}}$. This phase is much alike the situation where a tidal disrupted star’s debris is accreted onto a black hole (Rees, 1988; Cannizzo et al., 1990).

The torus is the predecessor of an accretion disk. After the accretion disk is developed, the inner edge decreases until it is disrupted by the torque of the magnetic field of the pulsar. The inner most radius $r_{\rm{M}}\approx\,r_{\rm{A}}$ (Frank et al., 2002), where

[TABLE]

is the Alfven radius, $\dot{M}_{\rm{acc,16}}$ is the accretion rate of the accretion disk in units of $10^{16}$ g/s, $\mu_{30}$ is the magnetic dipole in units of $10^{30}$ G cm3, $m_{\rm{p}}\equiv M_{\rm{p}}/M_{\odot}$. After the accretion disk reaches the $r_{\rm{M}}$, the propeller effect ejects the accreting matter thus the mass of the accretion disk declines steadily. This is much alike the situation where a pulsar is accreting from a fossil disk. Therefore, we follow the strategy of Chatterjee et al. (2000), to assume the accretion rate in two phases:

[TABLE]

[TABLE]

where $t=0$ is the onset of disk formation and $t=\tau$ is the moment when the disk descends to $r_{\rm{M}}$, $\beta$ is the index to describe the decreasing of the accretion rate, $\beta=19/16$ for an electron scattering dominated disk opacity, and $\beta=1.25$ for a Kramer opacity. $\beta=5/3$ is a well-known value given by Rees (1988), when considering a tidal disrupted star’s debris accreting onto a black hole. $r_{\rm{in}}$ evolves with time as:

[TABLE]

[TABLE]

The constant accretion rate $\dot{M}_{\rm{acc,0}}$ is normalized to the total bounded mass:

[TABLE]

which gives:

[TABLE]

On the other hand, if we assume the accretion disk to be a Shakura-Sunyaev disk (Shakura & Sunyaev, 1973), then

[TABLE]

where $\alpha$ is the viscosity index.

Combining equations (14, 15), the accretion rate can be solved as:

[TABLE]

where $R_{\rm{circ,10}}$ is $R_{\rm{circ}}$ in units of $10^{10}$ cm.

Since the temperature of the accretion disk increase inwardly, the inner most region of the accretion disk dominates the radiation, where the temperature is:

[TABLE]

where $r_{\rm{in,10}}$ is the inner radius of the accretion disk in units of $10^{10}$ cm. With the time evolution of $\dot{M}_{\rm{acc}}$ and $r_{\rm{in}}$ from equations (14, 12), the temperature as a function of $t$ is shown in figures 2 and 3.

2.3 Inverse Compton scattering from the pulsar wind

The soft photons from the accretion disk are inverse-Compton scattered to high energy by electrons in the pulsar wind. If we assume a monochromatic electron energy of the pulsar wind, with the Lorentz factor $\Gamma$, then the characteristic energy of the scattered photon is (Blumenthal & Gould, 1970):

[TABLE]

in the Thomson regime ($kT\Gamma/0.511\,\text{MeV}\ll 1$), and

[TABLE]

in the Klein-Nishina regime ($kT\Gamma_{0}/0.511\,\text{MeV}\gtrsim 1$), where $\Gamma_{4}\equiv\Gamma/10^{4}$. For $kT\sim 10$ eV as we calculated above, a pulsar wind with $\Gamma\sim 10^{3}-10^{4}$ will produce inverse Compton scattered photons in the Fermi-LAT’s energy range (100 MeV–100 GeV). Assuming an isotropic soft photons field, the energy differential power distribution (photons per energy range, per unit time) of inverse Compton scattering in Thomson limit is (Blumenthal & Gould, 1970):

[TABLE]

where $dN_{\rm{e}}(l)$ is the number of electrons at distance $l$ along the line of sight, within range $dl$. In figure 4, we illustrate the IC process. The accretion disk is not necessarily face on as plotted; If we suppose nearly all the spin down energy of the pulsar is converted into the kinetic energy of electrons in the pulsar wind at $l$ (as adopted by Tavani & Arons (1997), $\sim 96\%$ of the pulsar wind energy is in the kinetic energy of electrons for this system. For a discussion on the role of pulsar wind in spin braking, see Tong 2016 and references therein), then we have:

[TABLE]

$\epsilon_{\rm{ph}}$ in equation (20) is the energy of the soft photons, $r_{0}$ is the classic radius of the electron, $n(\epsilon_{\rm{ph}},l)$ is energy differential number density of the soft photons at $l$, $f(x)=2x\ln x+x+1-2x^{2}$ for $0<x<1$, $f(x)=0$ for $x>1$ and $x\equiv E_{\gamma}/(4\Gamma^{2}\epsilon_{\rm{ph}})$.

For simplicity, we consider the case where the accretion disk is viewed face-on. The soft photon field at $l$ is:

[TABLE]

A change of the inclination will slightly reduce the density of photon field. This will eventually lead to larger fitted $\eta$. Combining equations (20,21,22), and work out the integrate over $l$,

[TABLE]

Since the surface brightness of a black body $\propto T^{4}$, and $T\propto R^{-3/4}$, the soft photon field from the accretion disk can be represented by a ring with width of $r_{\rm{in}}$ and temperature of $T(r_{\rm{in}})$. Therefore the integration over $R$ in Equation (23) can be simplified as:

[TABLE]

Assuming the pulsar wind is isotropic, the spectrum of the inverse Compton process is:

[TABLE]

With above equations, and the known temporal function of $T$ and $r_{\rm{in}}$ in equations (17) and (12), we know how the spectrum energy distribution (SED) evolves with time, as shown in figure 5. For each time, we show the SED corresponding for $(\beta,\eta)=(19/16,0.059),(1.25,0.047),(5/3,0.024)$ in black dashed, red dash-dotted and blue solid curves respectively. $\alpha$ is 0.4 as the best fitted value found in the next section. The numbers labeled besides the apexes of each group of curves indicates the corresponding days after the periastron.

3 Comparing with observations

In order to evaluate this model, we compare it against observations. Caliandro et al. (2015) presented the light curves and SEDs of the gamma-ray flares in 2010/2011 and 2014 passages. They divided the light curves of the flares into two phases based on their fluxes level: the peak phase and the tail phase. The peak phase of the first/second passage is from the 31 st day to 40 th/42 nd day after the periastron, and the tail phase is from the 40 th/42 nd to 71st day after the periastron. They obtained the SEDs by integrating the photons over the peak phases, tail phases and all the flares duration. To compare with them, we also average our calculated SEDs over the corresponding time ranges. In figures 6 and 7, sub-panels titled with “All”, “Peak” and “Tail” correspond to SEDs averaged over the whole durations, peak and tail phases. The black dashed, red dash-dotted and blue solid curves are SEDs calculated with different pairs of $(\beta,\eta)$ (see the captions of the figures for details). The green points with error bars are observed SED data from Caliandro et al. (2015), those points with downward arrows are observational upper limits. Integrating equation (25) over the Fermi-LAT energy range (100 MeV-100 GeV) at each instance gives the flux as function with time, i.e., the light curves. We compare the calculated light curves with the observed ones in the panels titled with “Light curve” of figures 6 and 7.

4 conclusion and discussion

In this paper, we study a new explanation of the observed post-periastron GeV flares of the binary system PSR B1259-63/LS2883. This phenomenon could be results as the pulsar wind inverse Compton-scattering of the soft photons, which are emitted from an accretion disk around the pulsar. The accretion disk composed by the matter from the circumstellar disk (CD) captured by the pulsar’s gravity at disk-crossing. The pulsar crosses the CD twice in each orbital period, one pre-periastron and another post-periastron. In the post-periastron crossing, the density of the CD is not enough so that the pulsar wind prevents the matter from being accreted. While in the other crossing the density is sufficient. That explain the GeV flares are only observed once in each orbit.

With certain parameters of the accretion disk, this model can reproduce the observed SEDs and light curves satisfactorily (see figures 6 and 7).

4.1 Sub-structures in the light curves and the disk instability

It is worth mentioning that the light curves in Caliandro et al. (2015) show some significant sub-structures other than a pure decay. The disk instability might produce the above-mentioned small-scaled outbursts (see Mineshige & Wheeler 1989 and references therein; see Lasota 2001 for a review). The recurrence time of those outbursts is of the viscous timescale near the inner edge of the disc, which is $\sim 0.1$ days estimated with the best fitted parameters. However the variance time scale of the sub-structures in the light curves is of several days. Therefore, if the disk instability accounts for the sub-structures in the light curves, a different $\alpha\lesssim 0.05$ is needed in the decaying phase than that in the rising phase. It is possible that the $\alpha$ of an accretion disk is not a constant and is time-dependent via other disk properties (Smak, 1984). In this work, we only consider the quasi-stable disk as a simplification and do not try to reproduce the sub-structures.

4.2 UV excess from the accretion disk

One critical evidence to prove or to disprove the formation of the accretion disk is the expected UV emission excess. The frequency differential energy flux at frequency $\nu$ from the accretion disk is:

[TABLE]

where $i$ is the inclination of the disk and $D$ is the distance to the system. We assume $i=45^{\circ}$ and use the simplification that the integration over the whole accretion disk can be approximated by a isothermal ring with width $r_{\rm{in}}$ and temperature $T(r_{\rm{in}})$. While the flux from the Be star is:

[TABLE]

where $T_{\star}=30200$ K is the temperature of the Be star, and the flux from the CD $F_{\nu,\rm{sd}}$ can be represented by an isothermal disk with a temperature of $0.6\,T_{\star}$ and an effective radius of the pseudophotosphere $R_{\rm{eff}}(\nu)$ (Carciofi & Bjorkman, 2006).

The SED in the optical band is shown in figure 8. An excess of the SED in UV band is a natural prediction. The energy flux ($\nu F_{\nu}$) from $2\times 10^{16}$ Hz (82.7 eV) to $3\times 10^{16}$ Hz (165 eV) is from $3.2\times 10^{-11}$ ergs/cm2/s to $\sim 1\times 10^{-12}$ ergs/cm2/s at the peak of the flare. It is higher than the low-energy tail of the soft X-ray emission from the shock, extrapolating to the same energy range from the figure 3 of Kong et al. (2012). Therefore, an extreme-UV observation covering $2\times 10^{16}-3\times 10^{17}$ Hz (10-15 nm) during the GeV flare period can be used to test the proposed model.

4.3 The dependence on the CD density profile of the model

We see from table 1 that the mass transfer scenario dependents on the density profile of the CD. If $n=3$, or equivalently $\rho_{0}$ halves from the current assumed value of $1\times 10^{9}$ g/cm3, the mass transfer condition will not be satisfied at disk-crossing; on the other hand, if $n=3.5$ or equivalently $\rho_{0}$ doubles from the current value, the mass transfer can occur at both disk-crossings. Thus people expect a smaller GeV peak $\sim 20$ days after the main flare with $\sim 1/4$ of the flux. The fluxes uncertainties of the current light curves are too large to confirm or to exclude it.

4.4 The accretion of angular momentum

In our treatment, the accretion of angular momentum is due to the velocity and density gradient along the CD. Early studies suggested the specific angular momentum captured in this way is equation (7). However Davies & Pringle (1980) and Livio et al. (1986) argued analytically and numerically respectively that, very little (if not none) angular momentum can be accreted by the center mass in this way. Later more sophisticated numerical calculations (Ruffert & Anzer, 1995; Ruffert, 1997) showed that up to 70% of the total available angular momentum can be accreted. In our paper, as long as the actual specific angular momentum is enough, so that the circular radius $R_{\rm{circ}}$ is well outside the light cylinder, the validity of our model is not damaged. In our model, $R_{\rm{circ}}$ influences the gamma-ray flux in the combination $\alpha^{4/5}R^{-5/4}_{\rm{circ}}$, thus any change of $R^{-5/4}_{\rm{circ}}$ will be absorbed by a change of $\alpha$ during the fitting.

In fact, as $\alpha$ varies in a reasonable range, the time delay between the disk-crossing and the flare can be adjusted accordingly in a range of a couple of weeks. As a result, this model has a weak predictive power over the position of the CD. As we mentioned in above paragraphs, that this model is potentially testable with observation of the UV excess. Besides, we expect that the spin period derivative of PSR B1259-63 increases during the GeV flare, due to the additional torque acting on the pulsar magnetic fields by the accretion disk.

SXY thanks Prof. J. Takata for helpful instructions. The authors thanks anonymous referee for his/her careful reviewing and helpful comments, which brings much improvement of the manuscript. This work is partially supported by a GRF grant under 17302315.