# High-energy emissions from the gamma-ray binary 1FGL J1018.6-5856

A. M. Chen11affiliation: Institute of Astrophysics, Central China Normal
University, Wuhan 430079, China 22affiliation: Department of Physics, The
University of Hong Kong, Pokfulam Road, Hong Kong , C. W. Ng22affiliation:
Department of Physics, The University of Hong Kong, Pokfulam Road, Hong Kong ,
J. Takata33affiliation: School of physics, Huazhong University of Science and
Technology, Wuhan 430074, China , Y. W. Yu11affiliation: Institute of
Astrophysics, Central China Normal University, Wuhan 430079, China ,K. S.
Cheng22affiliation: Department of Physics, The University of Hong Kong,
Pokfulam Road, Hong Kong<EMAIL_ADDRESS><EMAIL_ADDRESS>
###### Abstract

The gamma-ray binary 1FGL J1018.6-5856 is known to be composed of a massive
star with a compact object with an orbital period of 16.54 days. We study
mechanisms of high-energy emissions (X-ray, GeV, and TeV $\gamma$-rays) from
this binary system. The paper is composed of two parts. In the first part, we
re-analyze the GeV emission from 1FGL J1018.6-5856 using 7 yrs of the
observational data of Fermi/LAT. The data analysis of phase-resolved spectrum
suggests that the 0.1-10 GeV emissions may contain two different components:
(1) a simple power-law with exponential cutoff component in the LOW state; (2)
an additional component extended to $\sim$10 GeV in the HIGH state. In the
second part, we develop a model to explain the multi-wavelength emissions of
this binary. Assuming that 1FGL J1018.6-5856 includes a pulsar, the GeV
spectrum in the LOW state can be well explained by the two-layer outer gap
model of pulsar magnetosphere. As the pulsar is orbiting around the star, the
pulsar is considered to be wrapped by two kinds of termination shock: a bow
shock due to the interaction between pulsar wind and stellar wind, and a back
shock due to the effect of orbital motion. We proposed that X-rays are mainly
produced by the synchrotron radiation at the apex of the bow shock, and TeV
$\gamma$-rays are mainly produced by the inverse-Compton scattering at the
back shock. The sharp peak in X-ray and TeV light curves are produced by in
the shock tail due to the Doppler boosting effect at inferior-conjunction
while another broad sinusoidal modulations are possibly produced by the apex
of the bow shock and the back shock due to variations of shock conditions.

binaries: close — gamma rays: stars — X-rays: binaries — raiation mechanisms:
non-thermal

##  1 Introduction

The gamma-ray binary is a class of high-mass X-ray binary systems containing a
compact object and a massive star and emit high-energy $\gamma$-rays beyond
1MeV (see Dubus 2013 for a review). The compact object is widely thought to be
a rotation-powered pulsar, although it still remains to be proven for most
binary systems. Until now, 7 high-mass gamma-ray binaries have been detected,
namely, PSR B1259-63/LS2883 (Aharonian et al. 2005), LS 5039 (Aharonian et al.
2006), LS I+61∘303 (Albert et al. 2006), H.E.S.S. J0632+057 (Hinton et al.
2009), 1FGL J1018.6-5856 (Ackermann et al. 2012), PSR J2032+4127/MT91
213111Strictly speaking, PSR J2032+4127/MT91 213 is a pulsar/Be star binary
system resemble to PSR B1259-63/LS2883, but its TeV counterpart is still not
detected. (Lyne et al. 2015) and LMC P3 (Corbet et al. 2016, the first gamma-
ray binary to be found outside the Milky Way). Only PSR B1259-63/LS 2883 and
PSR J2032+4127/MT91 213 have been confirmed to contain a pulsar so far.

The $\gamma$-ray source 1FGL J1018.6-5856 (J1018.6 hereafter for short) was
identified as a gamma-ray binary by Ackermann et al. (2012) based on the blind
search for periodic sources of $\gamma$-ray data taken by Fermi/LAT. The
follow-up observations in radio, optical and X-ray also confirmed the binary
nature of this source with a period of 16.6 days. The optical spectroscopy
indicates that the optical counterpart is a massive main-sequence star with a
spectral type of O6V((f)). The distance derived from the interstellar
absorption lines of the companion star is $d_{\rm L}=5\pm 2\rm{kpc}$. Later
Abramowski et al. (2012, 2015) also discover the TeV counterpart of this
system using the HESS telescope. The TeV light curves show similar orbital
modulations as seen in X-ray band. An et al. (2013, 2015) gave detail reports
on NuSTAR, XMM-Newton, and Swift-XRT observations of J1018.6 and constrained
the orbital period to be $16.544\pm 0.008$ days. They found that the X-ray
light curves have a regular spike with small variability and a broad
sinusoidal modulation. The behavior of J1018.6 at different orbital phases is
similar in X-ray, GeV, and TeV bands, showing in all cases a maximum flux near
phase 0 222The $\gamma$-ray maximum is denoted as phase 0 in Ackerman et al.
(2012), and we use the same denotes in the paper.(see Fig.3 in Abramowski et
al. 2015 for the TeV, GeV and X-ray light curves of J1018.6). Strader et al
(2015) presented medium-resolution optical spectroscopy of the O star
secondary of J1018.6 using the SOAR telescope. They find that the O star has a
low radial velocity with semi-amplitude of 11-12 $\rm{km\cdot s^{-1}}$ which
favors a neutron star primary. They also find that the high energy (X-ray and
$\gamma$-ray) flux maxima occur when the star is behind the compact object
along the line of sight.

Unfortunately, the orbital parameters and accurate conjunction phases of
J1018.6 are still unknown yet, and only limited conclusions can be drawn on
the relation between the compact object and the companion star (Abramowski et
al. 2015). However, J1018.6 shares many properties with another well studied
gamma-ray binary LS 5039\. Both two binary systems are composed of an almost
identical massive star (O6.5(f) for LS 5039, and O6V(f) for J1018.6) and a
compact object orbiting on a timescale of days ($P_{\rm orb}=3.9$ days for LS
5039, and $P_{\rm orb}=16.5$ days for J1018.6). They are fairly steady
$\gamma$-ray sources and their periodic modulations have not shown large
changes (Ackerman et al. 2012). The spectrum of J1018.6 also resemble that of
LS 5039 (i.e., a power-law in X-ray band and a power-law with exponential cut-
off in GeV band). Both of them show a maximum flux when the star is behind the
compact object in X-ray and TeV light curves. With the facts given above, we
can give rough constraints on the system parameters and use the emission model
which has been applied to LS 5039 (e.g., Takata et al. 2014) to explain the
high energy emissions from J1018.6.

The multi-wavelength emissions of LS 5039 has been well studied under the
pulsar model. In the rotation-powered pulsar scenario, the interactions
between a pulsar wind and a stellar wind will form a bow shock (Tavani & Arons
1997; Dubus 2006a). The effect of the orbital motion of pulsar also produce a
termination shock in the opposite direction of companion star (Bosch-Ramon et
al. 2012; Zabalza et al. 2013; Takata et al. 2014). The shocks accelerate
electrons and/or positrons from the pulsar wind, and then emit X-rays and TeV
$\gamma$-rays through synchrotron radiation and inverse-Compton process,
respectively. The emissions from the shocked regions are expected to show
orbital variations due to a combination of the following effects: anisotropic
inverse-Compton scattering, the $\gamma$-ray absorption, the Doppler-boosting
effect, etc.

Besides the shock emissions, the pulsar itself will also contribute to the
high energy emissions of the system in different ways (Sierpowska-Bartosik
$\&$ Torres 2007; Kapala et al. 2010; Torres 2011; Takata et al. 2014). The
first component is due to pulsed emission originating in outer-gaps of the
pulsar which peaks around 1 GeV. This magnetospheric emission is not expected
to vary along its orbit in a binary. Alternatively, the high-energy emission
could also arise from the un-shocked pulsar wind via up-scattering of ambient
photons provided by the companion star (Kirk, Ball $\&$ Skjaeraasen 1999;
Cerutti, Dubus $\&$ Henri 2008). A narrow component appears at the high-energy
spectra peaking at an energy $\gamma_{\rm w}m_{\rm e}c^{2}$ (in the Klein-
Nishina regime), with $\gamma_{\rm w}$ being the Lorentz factor of the wind.
This modulated anisotropic inverse Compton emission is expect to peak when the
star is in front of the pulsar along the line of sight, with a stronger effect
for more edge-on inclinations.

This work purports to explain the high energy emissions from J1018.6 under the
rotation-powered pulsar scenario mentioned above. The main purpose of this
study are composed of two parts: (1) to present results of the observational
analysis using 7 yr Fermi/LAT data, which provide us with more detailed
information on the GeV emission from J1018.6, and (2) to develop an emission
model to discuss the radiation processes of the X-ray, GeV, and TeV emissions.
The paper is organized as follow: in Section 2, we report on the results of
our analysis of the Fermi/LAT data. We describe the emission from the pulsar
magnetosphere and shocked wind regions in Section 3 and Section 4. In Section
5 we compare the model predictions with the results of the multi-wavelength
observations. Finally, a brief discussion and summary are given in Sections 6.

##  2 Data Analysis

We analyzed the gamma-ray orbital light curve of J1018.6 using the data from
the Fermi Large Area Telescope. Photon events with energy ranged from 100 MeV
to 100 GeV and time ranged from 2008-08-09 to 2015-08-08 were selected, from
the “Pass 8 Source” event class with the corresponding instrumental response
function of “P8R2_SOURCE_V6”. The region of interest (ROI) is a
$20^{\circ}\times 20^{\circ}$ square centered at the epoch J2000 position of
the source:
$(\textrm{R.A.},\textrm{Dec})=(10^{\textrm{h}}18^{\textrm{m}}55.18^{\textrm{s}},-58^{\circ}56^{\prime}44.2^{\prime\prime})$.
To reduce the contamination from the Earth’s albedo, time intervals were
excluded when the ROI is observed at a zenith angle greater than $90^{\circ}$.
The Fermi Science Tools package version v10r0p5333Available at
http://fermi.gsfc.nasa.gov/ssc/data/analysis/software/ was used in the data
analyses in this study. The gtlike tool was used to perform maximum binned
likelihood analysis to obtain the spectral model for all the 3FGL catalog
sources (gll_psc_v16.fit) (Acero et al. 2015) that are within $25^{\circ}$
from the center of ROI, the galactic diffuse emission (gll_iem_v06) and the
isotropic diffuse emission (iso_P8R2_SOURCE_V6_v06), available from the Fermi
Science Support Center (FSSC)444http://fermi.gsfc.nasa.gov/ssc/. The spectral
parameters for sources that are non-variable and located $5^{\circ}$ away from
the center are fixed to their catalog values. Four extended sources within the
region: HESS J1303-631, Puppis A, Vela Jr and Vela X, were modeled by the
extended source templates provided by the FSSC. With the spectral indices
fixed to the global fit and leaving only the normalization parameters free,
this model is then used to calculate the orbital light curve. The orbital
phase is assigned to each of the photons using the TEMPO2 package (Hobbs et
al. 2006) with the Fermi pulg-in (Ray et al. 2011). In the timing model, we
assumed the orbital period to be 16.544 days (An et al. 2015) and the binary
model to be a main-sequence/pulsar binary (Wex, 1998). The whole data set is
then separated into 10 phase bins. Using the gtlike tool (binned likelihood
analysis) and the model obtained in the global fit, we derived the energy flux
of J1018.6 in each orbital phase (Fig. 1).

We performed spectral analysis in the selected orbital phase intervals to
investigate if the spectrum of 1FGL J1018.6-5856 is varying throughout the
orbital period. We defined the phase interval between 0.0 - 0.3 and 0.8 - 1.0
as the “HIGH” state and the interval between 0.3 - 0.8 as the “LOW” state. We
used the same data set described above and sub-selected these two states.
Also, we used the same spectral model as the initial guess. The spectral form
of J1018.6 is modeled by a power-law with simple exponential cutoff:

$\dfrac{\textrm{d}N}{\textrm{d}E}=N_{0}\left(\dfrac{E}{E_{0}}\right)^{-\Gamma}\textrm{exp}\left(-\dfrac{E}{E_{C}}\right),$
(1)

where $N_{0}$ is the normalization constant, $E_{0}$ is the scale factor of
energy, $\Gamma$ is the spectral power-law index and $E_{C}$ is the cut-off
energy. Using binned likelihood analysis, the best-fit parameters for HIGH
state are $N_{0}=(1.20\pm 0.02)\times 10^{-9}$, $\Gamma=-1.82\pm 0.02$ and
$E_{C}=3000\pm 190$. For LOW state, $N_{0}=(1.19\pm 0.03)\times 10^{-9}$,
$\Gamma=-1.94\pm 0.01$ and $E_{C}=3277\pm 175$. Then, we obtained the spectral
points by performing the model fit in each energy band, using a simple power-
law form for J1018.6:

$\dfrac{\textrm{d}N}{\textrm{d}E}=N_{0}\left(\dfrac{E}{E_{0}}\right)^{-\Gamma},$
(2)

where $\Gamma$ is assumed to be $-2.1$. Upper limit is derived when the
detection significance is less than $3\sigma$. The phase-resolved spectra are
shown in Fig. 2.

##  3 Magnetospheric Emissions

The power-law with exponential cut-off spectrum in GeV band of J1018.6 is
typical of the $\gamma$-ray pulsar observed by Fermi/LAT. It is indicated that
the GeV emissions are mainly contributed by the pulsar magnetospheric
radiations.

Currently, because there is no existing published result on the timing
parameters of the pulsar in J1018.6, the properties of the pulsar remains
unknown. Yet, in order to explain the complete emission spectrum of J1018.6,
the magnetospheric contribution cannot be ignored. We used the standard outer
gap model (Cheng et al. 1986a & 1986b) to simulate the curvature spectrum from
the charge accelerations in the outer gap, which extends from the null charge
surface to the light cylinder. The separation of the oppositely charged
particles induces an electric potential in the space between them, leading to
the growth of the outer gap. On the other hand, the curvature photons can
undergo photon-photon pair creation with the softer photons from the pulsar
surface. The accumulation of charges will reduce the electric potential,
resulting in the depletion of the outer gap. These two instantaneously
occurring processes can be resembled by the two-layer structure (Wang et al.
2010 & 2011) which defines contrasting charge densities for the primary
acceleration and the screening regions. In this study, we simulate the
magnetospheric emission from a pulsar that has a magnetic inclination angle
$\alpha=40^{\circ}$ and a viewing angle at $\beta=120^{\circ}$ (or
$60^{\circ})$. The period and the surface magnetic field strength are assumed
to be $P=0.05$ s and $B=10^{12}$ G respectively. The pulsar has an average
charge density of $0.7$ Goldreich-Julian charge density and the primary
acceleration region occupies $70\%$ of the gap thickness.

The observed gamma-ray magnetospheric spectrum of J1018.6 (LOW state in Fig.
2) is flat in the energy range below 1 GeV but at the same time has detection
beyond 10 GeV. These features indicate that the pulsar could have an outer gap
which gap size is varying in time. This is a result of the dynamical
development of the local charge density and electric field strength, which
depends on the instant pair creation rate and the current flow (Takata et al.
2016). When the gap size increases, the charges are accelerated by a greater
electric potential, and thus the maximum energy of photon radiated can be
higher, and vice versa. In the simulation, we simplify the dynamics by
superposing the curvature spectrum of the pulsar in different gap fractions,
together with a weighting factor that describes the fraction of time that the
gap fraction remains at the corresponding value. A similar non-stationary
superposition method has been used by Chai et al. (2016) to explain the sub-
exponentiality in the cutoff spectrum for the Vela and Geminga pulsars. In
their model, they assumed a power-law distribution of the curvature spectrum
in different cutoff energies derived from the very fine phase-resolved
spectrum and superpose the results to obtain the sub-exponential cutoff
spectrum. In this study, we adopted a similar power-law relation between the
gap fraction ($f$) and the fraction of time ($w_{f}$):

$w_{f}=\frac{f^{\lambda}}{\sum{f^{\lambda}}},$ (3)

where $\lambda$ is assumed to be $-2.487$ from fitting the resultant
superposed spectrum to the $Fermi$-observed spectral energy distributions at
the LOW state.

##  4 Shock emissions

The interactions between the pulsar wind and the stellar wind will form a bow
shock with a hollow-like tail, and the effect of the orbital motion of pulsar
will also produce a back shock in the opposite direction of companion star.
The shocks accelerate the electrons and positrons of the pulsar wind, and then
emit the X-rays and TeV $\gamma$-rays through synchrotron radiation and
inverse-Compton process, respectively.

###  4.1 Orbital Phase

Considering an energetic pulsar in orbit around the massive companion star,
the orbital separation $d$ is given by (e.g., Dubus 2006b)

$d=\frac{a(1-e^{2})}{1+e\cos\theta},$ (4)

with the semi-major axis $a=(GMP_{\textrm{orb}}^{2}/4\pi^{2})^{1/3}$, $M$ the
total binary mass, $e$ the eccentricity, $\theta$ the true anomaly. The
orbital phase (mean anomaly) is defined as

$\phi=\frac{\varphi-e\sin\varphi}{2\pi}+\phi_{0},$ (5)

with the eccentric anomaly

$\varphi=2\arctan\left(\sqrt{\frac{1-e}{1+e}}\tan\frac{\theta}{2}\right),$ (6)

and $\phi_{0}$ is the orbital phase at periastron which equals 0 for most
cases.

The system parameters and orbital phases of J1018.6 are still unknown yet.
However, J1018.6 shares many similarities with another well studied gamma-ray
binary LS 5036, a limited conclusions can be drawn by comparing the properties
between J1018.6 and LS 5039. The period of J1018.6 is about 4 times longer
than LS 5039, which means that a semi-major axis is larger by a factor of 2.5
than in LS 5039. The low amplitude of the flux modulation can be interpreted
as a sign of a low eccentricity orbit and low inclination. In our calculation,
we adopted a modest value of $e\simeq 0.35$ and $i\simeq 40^{\circ}$ which is
similar to the case of LS 5039. According to the optical spectroscopy analysis
done by Strader et al (2005), finding that the high-energy flux maxima occur
when the star is behind the compact object along the line of sight, we
estimate the value of periastron angle with $\omega\simeq 90^{\circ}$. We
illustrate the orbital phase of J1018.6 applied in our model in Fig. 3. We
should point out that the accurate determination of system parameters and
orbital phases still need further observations.

###  4.2 Shock Geometry

The location of the terminal bow shock is determined by the dynamical pressure
balance between the pulsar wind and the stellar wind (Dubus 2006a):

$\frac{L_{\textrm{sd}}}{4\pi r_{\rm s}^{2}c}=\rho_{\rm w}v_{\rm w}^{2},$ (7)

with $\rho_{\rm w}=\dot{M}/4\pi R^{2}v_{\rm w}$ is the wind density. Then the
distance from the shock contact discontinuity in the shock apex to the pulsar
is determined by

$r_{\textrm{s}}=d\frac{\eta^{1/2}}{1+\eta^{1/2}},$ (8)

with $\eta=L_{\textrm{sd}}/c\dot{M}v_{\textrm{w}}$. $L_{\textrm{sd}}$ is the
spin down luminosity of the pulsar, $c$ is the speed of light, $\dot{M}$ is
the mass loss rate of the massive star and $v_{\textrm{w}}$ is the velocity of
the stellar wind. We assume that the value of $\eta$ is constant along the
whole orbit.

Away from the apex, the shock surface becomes a hollow cone. The half opening
angle of the shock contact discontinuity can be well described by (Eichler &
Usov 1993)

$\theta_{\textrm{shock}}=2.1(1-\overline{\eta}^{2/5}/4)\overline{\eta}^{1/3},$
(9)

where $\overline{\eta}=\min(\eta,\eta^{-1})$.

The effect of the orbital motion of pulsar will also produce a termination
back shock in the opposite direction of companion star. The distance from the
back shock to pulsar can be determined by the balance between the ram pressure
of the stellar wind due to Coriolis force and the ram pressure of the pulsar
(Bosch-Ramon & Barkov 2011):

$\frac{L_{\textrm{sd}}}{4\pi x_{\rm s}^{2}c}=\rho_{\rm w}(2\Omega x_{\rm
s})^{2},$ (10)

with $\Omega=2\pi/P_{\rm orb}$ is the angular velocity of the binary system.
Then the distance from the shock contact discontinuity in the shock apex to
the pulsar is determined by (Bosch-Ramon & Barkov 2011; Zabalza et al. 2013):

$x_{\rm s}\sim\sqrt{\frac{L_{\rm sd}v_{\rm w}}{\dot{M}c(2\Omega)^{2}}},$ (11)

with $\Omega\sim 10^{-6}/s$ for J1018.6 yields $x_{\rm s}\sim 2-4d$.

###  4.3 Shock Radiations

Although some gamma-ray binaries show correlated X-ray and TeV $\gamma$-ray
flux variations in light curves, however, for J1018.6, it is difficult to
model both X-ray and TeV $\gamma$-ray emissions using a one zone population of
electrons. Furthermore, the phase-averaged TeV $\gamma$-ray spectrum of
J1018.6 extends up to $\sim 20\rm{TeV}$ with $E_{\rm{break}}\geq 10\rm{TeV}$,
posing constraints on the magnetic field with $0.001\leq B\leq 0.1\rm G$, with
such a low magnetic field can not explain the observed X-ray flux level
(Abramowski et al 2015). A similar situation also appears in the case of LS
5039. Therefore, with the facts given above point toward to the presence of at
least two particle populations or two radiation zones to explain the X-ray and
TeV $\gamma$-ray emission. In this paper, we use a similar emission model
proposed by Takata et al. (2014): we assumed that the X-ray emissions are
mainly produced by the apex region of the bow shock due to a higher magnetic
field and the TeV $\gamma$-rays are mainly produced by the back shock with a
lower magnetic field and small absorption. The intense stellar photon field at
the shock apex, and corresponding high opacities owing to pair production will
reduce the TeV emissions at this region greatly. The synchrotron radiation
from the back shock can be ignored because the magnetic field is small compare
to the shock apex. In our calculations, we approximate that the observed
emissions are mainly produced in three regions: the bow shock (with apex and
tail), and the back shock as illustrated in Fig. 4. The radiation in the shock
apex and the back shock are isotropic since the bulk motion is non-
relativistic. As the particle flow propagating away from the shock apex, the
bulk Lorentz will increase gradually to relativistic in the shock tail
(Bogovalov et al. 2008, 2012). So the emissions from the shock tail should be
beamed. When the line of sight is near the beaming direction, we will receive
the Doppler-boosted radiations (e.g., Kong et al. 2012). The sharp peaks
around phase 0 in X-ray and TeV light curves can be produced by the Doppler
boosting effect in the tail of the bow shock.

####  4.3.1 Magnetic Field

In the free pulsar wind region, the magnetic energy will be gradually
converted into the particle energy from the light cylinder to the terminal
shock. We assume that the variation of the magnetization of pulsar wind
$\sigma$ follows a simple power-law form (e.g., Kong et al. 2011, 2012)

$\sigma=\sigma_{\textrm{L}}\left(\frac{r}{r_{\textrm{L}}}\right)^{-\alpha},$
(12)

with $\sigma_{\textrm{L}}$ being the magnetization at the light cylinder
$r_{\textrm{L}}$. According to the standard pulsar wind theory, the laboratory
frame magnetic field $B_{\rm p}$ at $r_{\textrm{s}}$ of pre-shocked region is
given by (Kennel & Coroniti 1984a,b):

$\frac{B_{\rm p}^{2}}{4\pi}\left(\frac{1+\sigma}{\sigma}\right)=\frac{L_{\rm
sd}}{4\pi r_{\rm s}^{2}c}.$ (13)

Then the downstream magnetic field for $\sigma\ll 1$ at shock region can be
calculated using perpendicular MHD shock jump conditions (Kennel & Coroniti
1984a,b; Dubus et al. 2015)

$B=\frac{1}{\hat{\gamma}-1}\left[\frac{L_{\rm sd}}{r_{\rm
s}^{2}c}\left(\frac{\sigma}{1+\sigma}\right)\right]^{1/2},$ (14)

with $\hat{\gamma}=4/3$ is the adiabatic index of the flow. We can estimate
the magnetic field at the shock regions using the equations above: with
$\sigma(r_{\rm s})\sim 0.1$ and $\sigma(x_{\rm s})\sim 0.008$, we have
$B(r_{\rm s})\sim 5\rm G$ and $B(x_{\rm s})\sim 0.2\rm G$ for the bow shock
and the back shock region, respectively.

####  4.3.2 Electron Distribution

It is usually assumed that the electrons in the shocked pulsar wind can be
accelerated to a power-law distribution:

$Q(\gamma_{\rm
e})=K\gamma_{\textrm{e}}^{-p}\qquad\textrm{with}\quad\gamma_{\textrm{e},\min}\leq\gamma_{\textrm{e}}\leq\gamma_{\textrm{e},\max},$
(15)

and $K$ is the normalization of particle distribution at the shock. The
minimum Lorentz factor of the shocked pulsar wind is assumed to be the average
Lorentz factor of cold pulsar wind at terminal shock with
$\gamma_{\textrm{e},\min}\approx\gamma_{\textrm{w}}$. The maximum Lorentz
factor is determined by the balance between the accelerate timescale
$\tau_{\textrm{acc}}=2\pi\zeta R_{\textrm{L}}/c$ and the synchrotron loss
timescale
$\tau_{\textrm{syn}}=\gamma_{\textrm{e}}m_{\textrm{e}}c^{2}/P_{\textrm{syn}}$,
so $\gamma_{\textrm{e},\max}=(3e/\zeta\sigma_{\textrm{T}}B)^{1/2}$. The
electron numbers can be calculated from

$\dot{N}=\frac{\Delta\Omega}{4\pi}\frac{L_{\rm sd}}{\gamma_{\rm
w}(1+\sigma)m_{\rm e}c^{2}},$ (16)

with $\Delta\Omega/4\pi$ is the fraction of pulsar wind electrons injected
into the shock front555Strictly speaking, the power in injected electrons is a
fraction of the spin-down luminosity $L_{\rm inj}\simeq\epsilon L_{\rm sd}$.
We leave the uncertain coefficient in $\Delta\Omega/4\pi$..

The electrons will lose their energies through radiative cooling or adiabatic
cooling processes. The evolved spectrum can be obtained from the continuity
equation of the electron distribution (Ginzburg & Syrovatshii 1964):

$\frac{\partial n(\gamma_{\rm e},t)}{\partial
t}+\frac{\partial\dot{\gamma}_{\rm e}n(\gamma_{\rm e},t)}{\partial\gamma_{\rm
e}}=Q(\gamma_{\rm e}),$ (17)

where $\dot{\gamma_{\rm e}}$ is the total energy loss rate of the electrons
and $Q(\gamma_{\rm e})$ is the injection rate. Assuming $\partial
n(\gamma_{\rm e},t)/\partial t=0$ , the spectrum at a steady state can be
acquired:

$n(\gamma_{\rm e})=\frac{1}{|\dot{\gamma_{\rm e}}|}\int^{\gamma_{\rm
e,\max}}_{\gamma_{\rm e}}Q(\gamma_{\rm e}^{\prime}){\rm d}\gamma_{\rm
e}^{\prime}.$ (18)

For a steady-state, non-thermal source where electrons are rapidly cooling by
emitting synchrotron radiation and inverse-Compton scattering, the total
energy loss rate of the electrons is (Moderski et al. 2005)

$\dot{\gamma}_{\rm e}=\dot{\gamma}_{\rm syn}+\dot{\gamma}_{\rm
IC}=\frac{4c\sigma_{\rm T}u_{\rm B}}{3m_{\rm e}c^{2}}\gamma_{\rm
e}^{2}(1+qF_{\rm KN}),$ (19)

where $q=u_{\star}/u_{B}$, with $u_{B}=B^{2}/8\pi$ being the magnetic energy
density, $u_{\star}=L_{\rm star}/4\pi cR^{2}$ being the seed photon energy
density. For a Planckian (black-body) distributed seed photons from the
companion star can be well approximated as a mono-energetic photon
distribution with energy $\epsilon_{0}=2.8kT_{\star}$. So $F_{\rm
KN}\simeq(1+b)^{-1.5}$ and $b=4\gamma_{\rm e}(2.8kT_{\star}/m_{\rm e}c^{2})$
(Moderski et al. 2005). The radiative cooling timescale $\tau_{\rm c}$ for an
electron with Lorentz factor $\gamma_{\rm e}$ is

$\tau_{\rm c}(\gamma_{\rm e})=\frac{3m_{\rm e}c}{4\sigma_{\rm T}\gamma_{\rm
e}u_{B}}/(1+qF_{\rm KN}).$ (20)

For a source with the dynamical timescale $\tau_{\rm dyn}$, the electron
number at a given Lorentz factor that can accumulate in the source could be
simply calculated by (Moderski et al. 2005; Kong et al. 2012)

$N(\gamma_{\rm e})\simeq\dot{N}Q(\gamma_{\rm e})\min[\tau_{\rm c}(\gamma_{\rm
e}),\tau_{\rm dyn}],$ (21)

with $\tau_{\rm dyn}=3r_{\rm s}/c$ is the typical dynamical timescale for the
post-shock of $\gamma$-ray binary.

In Fig. 5, we show the cooling timescales of electrons in the bow shock and
the back shock. We can see that the radiative loss timescales are much shorter
than the dynamical timescales of the flow in the bow shock region. The the
radiative cooling of electrons with $\gamma_{\rm e}<\gamma_{\rm KN}=m_{\rm
e}c^{2}/4\epsilon_{0}$ and $\gamma_{\rm KN}<\gamma_{\rm e}<\gamma_{\rm
s}=\gamma_{\rm KN}(q^{2/3}-1)$ are mainly dominated by the inverse-Compton
process in the Thomson regime and Klein-Nishina regime, respectively. At
higher energies, the Compton rate decreases due to the KN effects and
eventually synchrotron cooling always dominated for $\gamma_{\rm
e}>\gamma_{\rm s}$. The IC spectra have high-energy break at $h\nu(\gamma_{\rm
s})$ caused by the rapidly cooling of synchrotron radiation. The value of
$\gamma_{\rm s}$ is very sensitive to the strength of magnetic field, which
means that the strength of magnetic field in the shock regions will have a
strong impact on the TeV spectrum.

####  4.3.3 Radiation process

The synchrotron radiation power at frequency $\nu$ from a single electron with
Lorentz factor $\gamma_{\rm e}$ is given by (Rybicki & Lightman 1979)

$P^{\rm SYN}_{\nu}(\gamma_{\rm e})=\frac{\sqrt{3}e^{3}B}{m_{\rm
e}c^{2}}F(\frac{\nu}{\nu_{\rm c}}),$ (22)

where $\nu_{\rm c}=3\gamma_{\rm e}^{2}eB/(4\pi m_{\rm e}c)$. The function
$F(x)$ is defined as

$F(x)=x\int_{x}^{+\infty}K_{5/3}(k){\rm d}k,$ (23)

with $K_{5/3}$ being the modified Bessel function of order of 5/3.

The external inverse-Compton radiation power at frequency $\nu$ from a single
electron with Lorentz factor $\gamma_{\rm e}$ is given by (Aharonian & Atoyan
1981)

$\frac{{\rm d}P^{\rm EIC}_{\nu}(\gamma_{e})}{{\rm d}\Omega}=\frac{3\sigma_{\rm
T}}{4\pi}\int^{\infty}_{\nu_{\rm s,min}}{\rm d}\nu_{\rm s}\frac{\nu
f_{\nu_{\rm s}}^{\rm STAR}}{4\gamma^{2}_{e}\nu^{2}_{\rm s}}h(\xi,b_{\theta}),$
(24)
$h(\xi,b_{\theta})=1+\frac{\xi^{2}}{2(1-\xi)}-\frac{2\xi}{b_{\theta}(1-\xi)}+\frac{2\xi^{2}}{b^{2}_{\theta}(1-\xi)^{2}},$
(25)

where $f_{\nu_{\rm s}}^{\rm STAR}=\pi(R_{\star}/R_{\rm
s})^{2}2h\nu_{s}^{3}/c^{2}[{\rm exp}(h\nu_{s}/kT_{\rm{}_{\star}})-1]$ is the
flux density of the massive star photons, $\xi=h\nu/(\gamma_{\rm e}m_{\rm
e}c^{2})$, $b_{\theta}=2(1-{\rm cos}\theta_{\rm SC})\gamma_{\rm e}h\nu_{\rm
s}/(m_{\rm e}c^{2})$, $h\nu_{\rm s}\ll h\nu\leq\gamma_{\rm e}m_{\rm
e}c^{2}b_{\theta}/(1+b_{\theta})$, $\theta_{\rm SC}$ is the angle between the
injecting photons and the scattered photons, and is varied along with the
orbital phase.

####  4.3.4 $\gamma$-ray absorption

The intense radiation field provided by the companion star will absorb TeV
photons through pair production. The differential absorption opacity for a
$\gamma$-ray of energy $E$ is given by (Gould & Schreder 1967; Dubus 2006b):

$\textrm{d}\tau_{\gamma\gamma}=(1-\cos\theta_{\gamma\gamma})n_{\epsilon}\sigma_{\gamma\gamma}\textrm{d}l\textrm{d}\epsilon,$
(26)

where $\theta_{\gamma\gamma}$ is the collision angle and $n_{\epsilon}$ is the
photon density from the massive star. The $\gamma\gamma$ interaction cross
section is given by(Gould & Schreder 1967):

$\sigma_{\gamma\gamma}=\frac{\pi r_{\rm
e}^{2}}{2}(1-\beta^{2})\left[(3-\beta^{4})\ln\frac{1+\beta}{1-\beta}-2\beta(2-\beta^{2})\right],$
(27)

where $r_{\rm e}=e^{2}/m_{\rm e}c^{2}$ is the classical radius of electrons,
and $\beta=\sqrt{1-1/s}$, $s=(\epsilon
E/2m_{e}^{2}c^{4})(1-\cos\theta_{\gamma\gamma}).$ Absorption can only happen
when $s>1$ which give the threshold energy
$\epsilon_{\min}=2m_{e}^{2}c^{4}/E(1-\cos\theta_{\gamma\gamma})$ for
interaction with a $\gamma$-ray photon of energy $E$. The total opacity for a
$\gamma$-ray energy $E$ can be obtained by integrating Eq. (26) over the
travel distance and seed photon energy distribution. A more detailed
description of $\gamma$-ray absorption in massive X-ray binaries can be found
in Dubus (2006b).

The absorption spectra at inferior-conjunction and super-conjunction of the
bow shock and the back shock are shown in Fig. 6. The orbital variation of
$\gamma$-ray absorption at 1 TeV of the bow shock and the back shock in Fig.
7. We can see that the maximum absorption occur around 0.5 TeV at super-
conjunction and shift to 1 TeV at inferior conjunction. The TeV photons
produced at apex of the bow shock will suffer strong absorption due to a
denser photon field. We assume that the TeV emissions of J1018.6 are mainly
produced at the back shock region and ignore the contributions from the apex
of the back shock.

####  4.3.5 Doppler boosting effect

According to the relativistic hydrodynamical simulations by Bogovalov et
al.(2008), the shocked wind flows away from the apex to the tail which can
reach highly relativistic speeds. It means that the high energy emissions from
the shock tail should be beamed due to the Doppler boosting effect. The
relativistic boost is given by:

$D_{\rm obs}=\frac{1}{\Gamma(1-\beta\cos\theta_{\rm obs})},$ (28)

with $\theta_{\rm obs}$ being the angle between the direction of the flow and
the line of sight, $\Gamma$ being the Lorentz factor of the flow, and
$\beta=\sqrt{1-\Gamma^{-2}}$. The outgoing energy will be modified by
$\nu=D_{\rm obs}\nu^{{}^{\prime}}$ and the outgoing flux will be
$F_{\nu}=D_{\rm obs}^{3}F_{\nu^{{}^{\prime}}}^{{}^{\prime}}$. When the line of
sight is near the beaming direction, we will receive the Doppler-boosted
radiations.

The relativistic aberration will also affect the inverse-Compton process
(Dubus et al. 2010). The energy density from the star in the flow frame will
be reduced by a factor of $D_{\star}^{2}$, with the relativistic boost
involved is given by:

$D_{\star}=\frac{1}{\Gamma(1-\beta\cos\theta_{\star})},$ (29)

$\theta_{\star}$is the angle of the interaction photons and shock flows in the
comoving frame.

Finally, the observed total flux densities from all the shock accelerated
electrons at frequency $\nu$ with can be calculated by

$F_{\nu}=\frac{1}{4\pi d_{\rm L}^{2}}\cdot D_{\rm obs}^{3}\cdot\int{\rm
d}\gamma_{\rm e}N(\gamma_{\rm e})\cdot\left[P^{\rm
SYN}_{\nu}(\gamma_{e})+\frac{{\rm d}P^{\rm EIC}_{\nu}(\gamma_{e})}{{\rm
d}\Omega}\right]\cdot e^{-\tau_{\gamma\gamma}}.$ (30)

##  5 Results

In this section, we present calculated results using the model described above
and compare with observations. The detailed model parameters assumed in the
model fitting are given in Table. 1. Since the system parameters are still
unknown yet, we adopt a modest value of eccentricity and inclination with
$e\simeq 0.35$ and $i\simeq 40^{\circ}$ which is similar to LS 5039. The
observed fluxes are explained by a spin-down power $L_{\rm sd}\sim 4\times
10^{36}\rm{ergs^{-1}}$. If the pulsar has a typical magnitude of the surface
magnetic field $B_{\rm s}\sim 10^{12}\rm G$, the dipole radiation model
predicts the pulsar spin-down predicts the rotation period of $P\sim 0.1\rm
s$. The Lorentz factor of the free pulsar wind is estimated by $\gamma_{\rm
w}\sim\sigma_{\rm L}\gamma_{\rm L}\sim 4\times 10^{5}$. The best fit value of
the injection spectral slope is $p=1.5$, which is smaller than the prediction
of standard Fermi first order shock acceleration model. This may hints that
the electrons are accelerated via reconnection rather than Fermi process
(Dubus et al. 2015).

In Fig. 8, we calculate the spectrum and compare with observations. We
proposed that the X-ray emissions are mainly produced by the synchrotron
radiation in the bow shock with $B\sim 5\rm G$ at apex. The synchrotron
cooling timescales in the bow shock for the particles with Lorentz factor
$\gamma_{\rm e}>10^{6}$ is shorter than the IC cooling timescale in Klein-
Nishina regime. The strong magnetic field and $\gamma\gamma$ absorption in
this region will greatly reduce the TeV photons. We proposed that observed TeV
$\gamma$-rays are mainly produced in the back shock region with a lower
magnetic field of $B\sim 0.2\rm G$, then the IC cooling dominates the
synchrotron cooling for the electrons a Lorentz factor $\gamma_{\rm e}\sim
10^{7}$, and hence the spectrum of IC emissions can extend to 10 TeV. It is
shown that J1018.6 exhibits extremely powerful emissions at GeV band with a
spectral shape consistent with a power-law and an exponential cut-off at
energies of a few GeV (see Fig. 2) which is typical spectrum of $\gamma$-ray
pulsars observed by Fermi/LAT. Here, we used the two-layer outer-gap model to
fit the spectrum around GeV band. However, the flux variation along the
orbital phase (see Fig. 1) indicates an additional component in GeV band. We
leave this component discuss in Sec. 6.

In Fig. 9 and Fig. 10, we show our calculated light curves of X-ray and TeV
$\gamma$-ray and compare them with observations. The sharp peak around phase 0
in X-ray and TeV light curves are mainly emitted in the tail of the bow shock
due to the Doppler-boosting effect at inferior-conjunction. When the line of
sight cross the beaming direction, we will receive the Doppler-boosting
radiations. There is another broad sinusoidal modulation in X-ray light curves
while the modulation in TeV light curves are not so obvious as in X-ray. This
broad peak can be produced by the apex of the bow shock due to variations of
shock conditions. In previous model, it is usually believed that synchrotron
originated X-ray emissions from the shock apex would show a peak around
periastron due to a higher magnetic field. However, we show that it can be
also possibly peak at the aparastron. This is because in the shocked regions
with strong magnetic fields and with denser seed photon fields, the radiative
cooling timescales of electrons would be less than the dynamical timescales of
these regions. For the energy loss of the electrons produced the X-ray
emissions are mainly dominated by the IC cooling in Klein-Nishina regime,
which means that $f_{\nu}\propto N_{\rm e}\propto\tau_{\rm c}\propto
1/u_{\star}\propto R^{2}\propto d^{2}$. The orbital distance varies with the
orbital phase as sinusoidal function, which means that the X-ray flux will
also vary in the same way. For the case of the back shock, the dynamical
timescales will have a strong effect on the electron distribution, and the
flux variation in the TeV light curves is not so obvious as in X-ray band. The
IC cooling timescale at aparastron is longer than that in periastron, which
means that of electron remains much longer at higher energy. The hardening in
the electron spectrum will compensates the effect of magnetic field, and thus
show a flux maxima at aparastron. As we can see that the cooling process of
electrons will have a strong effect on the spectrum and light curves.

##  6 Conclusion and Discussion

In this paper, we discussed the mechanisms of high energy emissions from the
gamma-ray binary 1FGL J1018.6-5856. We presented the data analysis using 7 yr
data of Fermi/LAT. The GeV spectrum with a power-law and exponential cutoff at
the LOW state can be well explained the two-layer outer-gap model of pulsar
magnetosphere. For the spectrum at HIGH state, an enhancement at 0.1-10 GeV is
seen exclusively. As seen in the $\gamma$-ray orbital light curves and phase-
resolved spectrum at HIGH state, a pure magnetospheric emission can not
explain the GeV emission of J1018.6 well since the magnetospheric emission is
not expected to vary as the pulsar moves along its orbit in a binary. It is
indicated that there should exist another additional components around GeV
band emission.

Alternatively, the emission from the un-shocked pulsar wind up-scattering
stellar photons can also contribute to the GeV emissions (e.g. Kapala et al.
2010). However, the Compton flux are expected to peak at super conjunction
which is contrast to observations. Besides, an IC component around 1 GeV
require the Lorentz factor of cold pulsar wind with $\gamma_{\rm w}\approx
10^{4}$ which can not explain the spectrum in 10-100GeV. The inverse-Compton
origin of this component seems can be ruled out. Kong et al. (2012) proposed
that Doppler-boosting of synchrotron emission in the bow shock tail can
explain the GeV flare and spectrum of PSR B1259-63/LS 2883\. However, a
synchrotron spectrum cut off at the energy of 10 GeV will require a very small
magnetic fields or a very high Lorentz factor of the bulk flow in the shock
tail which seems inconsistent with observations. The variation of GeV light
curves is difficult to be re-produced in via inverse-Compton scattering and it
also does not favor the origin of Doppler boosting synchrotron radiation. A
more detailed study of the Fermi/LAT data will provide us a better
understanding of the GeV emissions.

We also investigated the emission mechanisms of X-ray and TeV $\gamma$-rays
from J1018.6. In previous studies, it is generally believed that the observed
non-thermal broad emissions from gamma-ray binaries are produced by the
energetic particles at the shock apex close to the wind standoff region(e.g.,
Dubus 2006a). However, the intense stellar photon field at this location and
corresponding high opacities owing to $\gamma\gamma$ absorption will reduce
the TeV fluxes greatly. Furthermore, the phase-averaged $\gamma$-ray spectrum
of J1018.6 extends up to $\sim 20\rm{TeV}$ with $E_{\rm{break}}\geq
10\rm{TeV}$, posing constraints on the magnetic field ($0.001\leq B\leq 0.1\rm
G$). Such a low magnetic field can not explain the observed X-ray flux. We
proposed that the X-ray emissions are mainly produced by the synchrotron
radiation in the shock apex with $B\sim 5\rm G$, while TeV $\gamma$-rays are
mainly produced in the back shock region with a lower magnetic field of $B\sim
0.2\rm G$ and small absorption of TeV photons. The emissions from the tail of
the bow shock will also contribute to the observed flux when the star is
behind the pulsar (at inferior conjunction). The sharp peaks around phase 0 in
X-ray and TeV light curves are produced in this region due to the Doppler
boosting effect while another broad peak around phase 0.35 are produced by the
shock apex and back shock due to the shock variations.

In our model, we assume that the binary has a symmetrical shock geometry and
simply divided the emission region in to 3 parts (see Fig. 4). However, the
realistic shock geometries in gamma-ray binaries are much more complicated
than our assumptions. According to the relativistic hydrodynamical simulations
done by Bosch-Ramon et al. (2012, 2015), the orbital motion will also deflect
the shock flows, and strong instabilities lead to the development of
turbulence and wind mixing in the flows will also disrupted the shock
structure and affect the shocked flow evolution. Detail numerical simulations
on gamma-ray binaries will provide us more accurate estimates of shock
regions, flow evolutions and emission models to explore the high-energy
emissions from these systems.

When we have finished the draft of this paper, we noted that An & Romani
(2017) appeared. According to their analysis, they find that the orbital
modulation in the low-energy $\gamma$-ray band is similar to that in the X-ray
band, suggesting a common spectral component. The orbital light curve above a
GeV changes significantly, suggesting that the GeV emission contains
significant flux from a pulsar magnetosphere. They find that a simple one-zone
model is inadequate to explain the shock emission, but that beamed
Synchrotron-self Compton radiation from adiabatically accelerated plasma in
the shocked pulsar wind can reproduce the complex multiband light curves,
including the variable X-ray spike coincident with the gamma-ray maximum.
Although the physical situation we consider is similar, there are notable
differences between our model and those of An & Romani (2017). Assuming that
J1018.6 includes a pulsar, the GeV spectrum in the LOW state can be well
explained by the two-layer outer gap model of pulsar magnetosphere. As for the
shock radiations, we proposed that X-rays are mainly produced by the
synchrotron radiation at the bow shock, and TeV $\gamma$-rays are mainly
produced by the external inverse-Compton scattering at the back shock. The
sharp peak in X-ray and TeV light curves are produced by in the shock tail due
to the Doppler boosting effect at inferior-conjunction, while another broad
sinusoidal modulations are produced by the apex of the bow shock and the back
shock due to variations of shock conditions.

Table 1: Model parameters of 1FGL J1018.6-5856  System parameters |  |  |  |
---|---|---|---|---
eccentricity of orbit |  | $e$ |  | 0.35
inclination angle |  | $i$ |  | $45^{\circ}$
periastron angle |  | $\omega$ |  | $90^{\circ}$
distance |  | $d_{\rm L}$ |  | 5.4 kpc
orbital period |  | $P_{\textrm{orb}}$ |  | 16.54 days
Pulsar and pulsar wind |  |  |  |
spin down power |  | $L_{\textrm{sd}}$ |  | $4\times 10^{36}{\rm erg~{}s^{-1}}$
rotation period |  | $P$ |  | 0.05 s
surface magnetic field |  | $B_{\textrm{s}}$ |  |  $1\times 10^{12}$G
Lorentz factor of free pulsar wind |  | $\gamma_{\textrm{w}}$ |  | $4\times 10^{5}$
magnetization at the light cylinder |  | $\sigma_{\textrm{L}}$ |  | $5\times 10^{3}$
magnetization evolution index |  | $\alpha$ |  | $1.2$
Star and stellar wind |  |  |  |
mass |  | $M_{\star}$ |  | $31M_{\odot}$
radius |  | $R_{\star}$ |  | $10R_{\odot}$
temperature |  | $T_{\star}$ |  | $38900\textrm{K}$
mass loss rate |  | $\dot{M}_{\textrm{w}}$ |  | $2\times 10^{-8}M_{\odot}/\textrm{yr}$
velocity of stellar wind |  | $v_{\textrm{w}}$ |  | $2\times 10^{8}\textrm{cm/s}$
Shocked region |  |  |  |
electron distribution index |  | $p$ |  | 1.5
acceleration efficiency |  | $\zeta$ |  | 1.0
bulk Lorentz factor of the flow |  | $\Gamma$ |  | 2.0

Figure 1: The orbital light curve of 1FGL J1018.6-5856 in 0.1 - 100 GeV
obtained from binned likelihood analysis. The red dash-line indicates the mean
energy flux.

Figure 2: The orbital phase resolved spectra of 1FGL J1018.6-5856. Red and
blue curves represent the LOW and HIGH states respectively.

Figure 3: Schematic diagram of the orbital phases of J1018.6 in our model. The
accurate determination of the orbital phase still need further observations.

Figure 4: Geometry of the termination shock. The interaction between the
pulsar wind and the stellar wind forms a bow shock with a hollow cone-like
tail. As the particle flow propagating away from the apex, the bulk Lorentz
will increase gradually to mild relativistic in the tail. The effect of
orbital motion of the pulsar will also produce a termination shock in the
opposite direction of the companion star (the back shock). The synchrotron
radiation in the apex of the bow shock and the IC process in the back shock
produce the X-rays and TeV $\gamma$-rays, respectively. The emissions from the
tail of the bow shock will also contribute to the observed flux when the star
is behind the pulsar (at inferior conjunction).

Figure 5: Cooling timescales with respect to the electron energies.

Figure 6: Spectrum of $\gamma$-ray absorption at inferior-conjunction
($\phi=0$) and super-conjunction ($\phi=0.7$) of the bow shock (solid lines)
and the back shock (dashed lines).

Figure 7: Orbital variation of absorption for $E$=1 TeV at the bow shock
(solid lines) and the back shock (dashed lines).

Figure 8: Calculated spectra as compared with observations. The X-ray data are
taken from Dubus (2013) and the TeV data are taken from Abramowski et al.
(2015) without error bars. The GeV data show the analysis result of the
present work (Section 2). The solid lines show the results of the emission
model discussed in Section 3 and Section 4. The X-ray emissions are mainly
produced by the synchrotron radiation at the shock apex, and the TeV emissions
are mainly produced by the inverse-Compton process at the back shock. The GeV
emissions are mainly contributed by the outer-gap of the pulsar magnetosphere.

Figure 9: Calculated X-ray light curve as compared with observation. The X-ray
data are taken from Ackerman et al. (2012) without error bars. The dotted line
and the dashed line corresponding to the emission from the shock apex and the
shock tail, respectively, and the solid line correspond to the total flux.

Figure 10: Calculated TeV light curve as compared with observation. The TeV
data are taken from Abramowski et al. (2015) without error bars. The dotted
line and the dashed line corresponding to the emission from the back shock and
the shock tail, respectively, and the solid line correspond to the total flux.

## References

  * (1) Acero, F., Ackermann, M., Ajello, M., et al. 2015, ApJS, 218, 23
  * H. E. S. S. Collaboration et al. (2012) Abramowski, A., Acero, F., et al. 2012, A&A, 541, A5
  * H. E. S. S. Collaboration et al. (2015) Abramowski, A., Aharonian, F., et al. 2015, A&A, 577, A131
  * Fermi LAT Collaboration et al. (2012) Ackermann, M., Ajello, M., et al. 2012, Science, 335, 189
  * Aharonian et al. (2005) Aharonian, F., Akhperjanian, A. G., Aye, K.-M., et al. 2005, A&A, 442, 1
  * Aharonian et al. (2006) Aharonian, F., Akhperjanian, A. G., Bazer-Bachi, A. R., et al. 2006, A&A, 460, 743
  * Aharonian & Atoyan (1981) Aharonian, F. A., & Atoyan, A. M. 1981, Ap&SS, 79, 321
  * Albert et al. (2006) Albert, J., Aliu, E., Anderhub, H., et al. 2006, Science, 312, 1771
  * An et al. (2013) An, H., Dufour, F., Kaspi, V. M., & Harrison, F. A. 2013, ApJ, 775, 135
  * An et al. (2015) An, H., Bellm, E., Bhalerao, V., et al. 2015, ApJ, 806, 166
  * An & Romani (2017) An, H., & Romani, R. W. 2017, arXiv:1703.04214
  * Bogovalov et al. (2008) Bogovalov, S. V., Khangulyan, D. V., Koldoba, A. V., Ustyugova, G. V., & Aharonian, F. A. 2008, MNRAS, 387, 63
  * Bogovalov et al. (2012) Bogovalov, S. V., Khangulyan, D., Koldoba, A. V., Ustyugova, G. V., & Aharonian, F. A. 2012, MNRAS, 419, 3426
  * Bosch-Ramon & Barkov (2011) Bosch-Ramon, V., & Barkov, M. V. 2011, A&A, 535, A20
  * Bosch-Ramon et al. (2012) Bosch-Ramon, V., Barkov, M. V., Khangulyan, D., & Perucho, M. 2012, A&A, 544, A59
  * Bosch-Ramon et al. (2015) Bosch-Ramon, V., Barkov, M. V., & Perucho, M. 2015, A&A, 577, A89
  * Cerutti et al. (2008) Cerutti, B., Dubus, G., & Henri, G. 2008, A&A, 488, 37
  * (18) Chai, Y., Cheng, K. S., & Takata, J. 2016, JASS, 33, 75
  * Cheng et al. (1986a) Cheng, K. S., Ho, C., & Ruderman, M. 1986a, ApJ, 300, 500
  * Cheng et al. (1986b) Cheng, K. S., Ho, C., & Ruderman, M. 1986b, ApJ, 300, 522
  * Corbet et al. (2016) Corbet, R. H. D., Chomiuk, L., Coe, M. J., et al. 2016, ApJ, 829, 105
  * Dubus (2006a) Dubus, G. 2006a, A&A, 456, 801
  * Dubus (2006b) Dubus, G. 2006b, A&A, 451, 9
  * Dubus et al. (2008) Dubus, G., Cerutti, B., & Henri, G. 2008, A&A, 477, 691
  * Dubus et al. (2010) Dubus, G., Cerutti, B., & Henri, G. 2010, A&A, 516, A18
  * Dubus et al. (2015) Dubus, G., Lamberts, A., & Fromang, S. 2015, A&A, 581, A27
  * Dubus (2013) Dubus, G. 2013, A&A Rev., 21, 64
  * Eichler & Usov (1993) Eichler, D., & Usov, V. 1993, ApJ, 402, 271
  * Ginzburg & Syrovatskii (1964) Ginzburg, V. L., & Syrovatskii, S. I. 1964, The Origin of Cosmic Rays, New York: Macmillan, 1964,
  * Gould & Schréder (1967) Gould, R. J., & Schréder, G. P. 1967, Physical Review, 155, 1404
  * Hinton et al. (2009) Hinton, J. A., Skilton, J. L., Funk, S., et al. 2009, ApJ, 690, L101
  * (32) Hobbs, G. B., Edwards, R. T., & Manchester, R. N. 2006,MNRAS, 369, 655
  * Kapala et al. (2010) Kapala, M., Bulik, T., Rudak, B., Dubus, G., & Lyczek, M. 2010, 25th Texas Symposium on Relativistic Astrophysics, 193
  * Kennel & Coroniti (1984a) Kennel, C. F., & Coroniti, F. V. 1984a, ApJ, 283, 694
  * Kennel & Coroniti (1984b) Kennel, C. F., & Coroniti, F. V. 1984b, ApJ, 283, 710
  * Kirk et al. (1999) Kirk, J. G., Ball, L., & Skjæraasen, O. 1999, Astroparticle Physics, 10, 31
  * Kong et al. (2012) Kong, S. W., Cheng, K. S., & Huang, Y. F. 2012, ApJ, 753, 127
  * Kong et al. (2011) Kong, S. W., Yu, Y. W., Huang, Y. F., & Cheng, K. S. 2011, MNRAS, 416, 1067
  * Lyne et al. (2015) Lyne, A. G., Stappers, B. W., Keith, M. J., et al. 2015, MNRAS, 451, 581
  * Moderski et al. (2005) Moderski, R., Sikora, M., Coppi, P. S., & Aharonian, F. 2005, MNRAS, 363, 954
  * (41) Ray, P. S., Kerr, M., Parent, D., et al. 2011, ApJS, 194, 17
  * Rybicki & Lightman (1979) Rybicki, G. B., & Lightman, A. P. 1979, New York, Wiley-Interscience, 1979. 393 p.,
  * Sierpowska-Bartosik & Torres (2007) Sierpowska-Bartosik, A., & Torres, D. F. 2007, ApJ, 671, L145
  * Strader et al. (2015) Strader, J., Chomiuk, L., Cheung, C. C., Salinas, R., & Peacock, M. 2015, ApJ, 813, L26
  * Takata et al. (2014) Takata, J., Leung, G. C. K., Tam, P. H. T., et al. 2014, ApJ, 790, 18
  * (46) Takata, J., Ng, C. W., & Cheng, K. S., 2016, MNRAS, 455, 4249T
  * Tavani & Arons (1997) Tavani, M., & Arons, J. 1997, ApJ, 477, 439
  * Torres (2011) Torres, D. F. 2011, Astrophysics and Space Science Proceedings, 21, 531
  * Wang et al. (2010) Wang, Y., Takata, J., & Cheng, K. S. 2010, ApJ, 720, 178
  * Wang et al. (2011) Wang, Y., Takata, J., & Cheng, K. S. 2011, MNRAS, 414, 2664
  * (51) Wex, N. 1998, MNRAS, 298, 67
  * Zabalza et al. (2013) Zabalza, V., Bosch-Ramon, V., Aharonian, F., & Khangulyan, D. 2013, A&A, 551, A17