The likely Fermi detection of the Be X-ray binary GRO J1008-57: hint for gamma-ray production in neutron-star X-ray binaries

TL;DR

This study reports a likely gamma-ray detection from the Be X-ray binary GRO J1008-57 using Fermi LAT data, suggesting episodic high-energy emission possibly linked to giant X-ray outbursts, and highlights the need for further observations.

Contribution

First evidence of gamma-ray emission from GRO J1008-57, indicating episodic high-energy activity in neutron-star X-ray binaries and challenging existing emission models.

Findings

Likely gamma-ray detection in orbital phase 0.8--0.9

Detection linked to a giant X-ray outburst around 2012-2013

Emission luminosity estimated at 10^{34} erg s^{-1}

Abstract

In our search for gamma-ray emission from Be X-ray binaries from analysis of the data obtained with the Large Area Telescope (LAT) on board the Fermi Gamma-Ray Space Telescope, we find likely detection of GRO J1008-57. The binary has an orbital period of 249.48 days, and it is only significantly detected in its orbital phase 0.8--0.9 (>4 sigma). Further light curve analysis indicates that the detection is probably largely due to an emitting event in one orbital cycle around year 2012--2013, following a giant X-ray outburst of the source. This property of having occasional gamma-ray emitting events is similar to that seen in another high-mass X-ray binary 4U 1036-56. However, models considering possible gamma-ray emission from an accreting neutron star have difficulty in explaining the observed 10^{34} erg s^{-1} luminosity of the source, unless the distance was largely over-estimated. Further observational studies are required, in order to more clearly establish the high-energy emission properties of GRO J1008-57 or similar high-mass X-ray binaries and find clues for understanding how gamma-ray emission is produced from them.

Tables

Table 1. Table 1 Fermi  LAT flux measurements of GRO J1008 − - 57

		Phase I		Phase II
$E$	Band	$F/{10}^{-12}$	TS	$F/{10}^{-12}$	TS
(GeV)	(GeV)	(erg cm-2 s-1)		(erg cm-2 s-1)
0.15	0.1–0.2	9.1	0	6.8	0
0.36	0.2–0.5	5.2	0	3.5	0
0.84	0.5–1.3	5.7$\pm$2.8$\pm$1.0	9	1.1	0
1.97	1.3–3.0	5.1	0	2.2	0
4.62	3.0–7.1	2.4$\pm$1.2$\pm$0.1	6	2.4	2
10.83	7.1–16.6	4	1	2.0	1
25.37	16.6–38.8	2.58$\pm$1.28$\pm$0.01	7	2.68$\pm$1.30$\pm$0.03	10
59.46	38.8–91.0	5.2	0	3.1	0
139.36	91.0–213.3	10.8	0	5.3	0
326.60	213.3–500.0	22.2	0	30.8	2

Table 2. Table 2 Three TS > > 9 detections in 0.1 orbital phase bins

Central Time	Orbital Phase	TS
(MJD)
55135.728	0.8–0.9	9.7
55559.844	0.5–0.6	9.1
56383.128	0.8–0.9	17.5

Full text

1. Introduction

X-ray binaries (XRBs) contain a compact primary, either a black hole or a neutron star, and when the companion is an early-type (O/B) massive star, such XRBs are classified as high-mass X-ray binaries (HMXBs). Generally the compact star in an HMXB system accretes from matter carried in the stellar wind of the companion and in most cases the accretion rates are low (e.g., Davidson & Ostriker 1973; Lamers et al. 1976; Walter et al. 2015). High X-ray luminosities (1035–1040 erg s*-1*) can be observed when the compact star is interacting with the dense part of the stellar wind of a Be star companion. Be XRBs account for a majority of the known HMXBs ($\sim$50%; Walter et al. 2015), in most of which the compact stars are magnetized neutron stars (Liu et al., 2006; Lutovinov & Tsygankov, 2009).

There are three Be XRBs, PSR B1259$-$63/LS 2883 (Aharonian et al., 2009; Tam et al., 2011; Abdo et al., 2011), LS I $+$61°303 (Abdo et al., 2009; Hadasch et al., 2012), and HESS 0632$+$057 (Hinton et al., 2009; Bongiorno et al., 2011; Li et al., 2017), have been seen to emit photons in GeV–TeV $\gamma$-ray band. They are also classified as $\gamma$-ray binaries as their spectral energy distributions peak in $\gamma$-ray energies (Dubus, 2013). $\gamma$-ray emission from them is orbitally modulated, and interestingly a superorbital periodicity is seen in LS I $+$61°303 (Ackermann et al. 2013; Ahnen et al. 2016 and references therein). PSR B1259$-$63/LS 2883 is the only $\gamma$-ray binary with a known $\sim$47.8 ms radio pulsar (Johnston et al., 1992). The pulsar moves around the Be-type companion with a long orbital period of $\sim$3.4 years and a high eccentricity of 0.87 (Johnston et al., 1994), and during the periastron passage a strong flare particularly at high-energies of X-ray to $\gamma$-ray is induced due to the close interaction between the winds of the two stars (e.g., Chernyakova et al. 2014; Tam et al. 2015; see also Xing et al. 2016). Recently another pulsar system PSR 2032$+$4127/MT91 213 has also been found to be a candidate $\gamma$-ray binary system with a pulsar moving around a Be-type companion in a $\sim$50 year long orbit (Lyne et al., 2015; Ho et al., 2017).

Currently there are only five confirmed $\gamma$-ray binaries in the Galaxy (including the above three). While the limited number may provide constraints on binary evolution and high-energy physical processes (e.g., Meurs & van den Heuvel 1989; Dubus 2013), searches for more members of this class or related sources with $\gamma$-ray emission will help improve our understanding. Be XRBs are good candidates for such searches because of the existence of the circumstellar disks, providing an environment for different possible interactions with the neutron star primaries. For example, Romero et al. (2001) has reported the likely detection of the variable $\gamma$-ray emission from a Be XRB system A 0535+26. The $\gamma$-ray emission was suggested to originate in hadronic processes, in which hadrons accelerated from the magnetosphere of the neutron star could impact the surrounding accretion disk and produce $\gamma$-ray photons via the neutral pion decay process (e.g., Cheng & Ruderman 1991). However, no $\gamma$-ray emission from this source was detected during a giant X-ray outburst (Acciari et al., 2011). It was also predicted that A 0535+26 would have a high neutrino yield (Anchordoqui et al., 2003), which is not found with IceCube. Motivated by these and taking advantage of the all-sky monitoring capability of the Large Area Telescope (LAT) onboard the Fermi Gamma-ray Space Telescope (Fermi), we conducted search for $\gamma$-ray emission from Be XRBs. The targets were selected from those given in Walter et al. (2015). Considering the X-ray and $\gamma$-ray emission from this type of sources are generally orbitally modulated, we only included the sources with known orbital parameters, allowing to conduct orbital-phase resolved search. Such search may be more sensitive due to orbitally-dependent physical processes (see, e.g., Xing et al. 2016).

In our search, likely detections of GRO J1008$-$57 in its certain orbital phase ranges were found. This Be XRB, discovered in 1993 (Stollberg et al., 1993), contains a spin period $P=93.5$ s transient X-ray pulsar, whose pulsed emission was detectable during the source’s X-ray outbursts (e.g., Coe et al. 2007; Kühnel et al. 2013). The binary has an orbit with a period of 249.48 days and an eccentricity of 0.68, determined from long-term timing of the pulsar’s pulsed emission (Coe et al., 2007; Kühnel et al., 2013). The magnetic field of the pulsar is known to be the highest among the Be XRBs, likely as high as $\sim$8$\times$1012 G given the suggested cyclotron line at $\sim$88 keV (Shrader et al., 1999). GRO J1008$-$57 exhibits type-I X-ray outbursts periodically at each periastron passage (e.g., Tsygankov et al. 2017), due to the interaction between the neutron star and the circumstellar disk around the Be companion (e.g., Reig 2011), and occasionally type-II outbursts. The latter type of the outbursts in Be XRBs, which are much stronger than type-I, is believed to cause major changes in the structure of the circumstellar disk, sometimes even leading to the disappearance of the disk (Reig, 2011).

In this paper, we report our results obtained for GRO J1008$-$57 from our analysis of the Fermi LAT data.

2. Fermi LAT Data Analysis and Results

2.1. LAT Data and Source Model

LAT scans the whole sky every three hours in the energy band of 0.1–500 GeV (Atwood et al., 2009). In the analysis, we selected LAT events from the Fermi Pass 8 database in the time period from 2008-08-04 15:43:36 (UTC) to 2017-10-19 00:56:35 (UTC). For GRO J1008$-$57, a $\mathrm{20^{o}\times 20^{o}}$ region centered at its position was selected. We followed the recommendations of the LAT team111http://fermi.gsfc.nasa.gov/ssc/data/analysis/scitools/ by including the events with zenith angles less than 90 degrees (to prevent the Earth’s limb contamination) and excluding the events with quality flags of ‘bad’.

All sources within a 20 degree region centered at the target in the Fermi LAT 4-year catalog (Acero et al., 2015) were included to make the source model. Spectral forms of these sources are provided in the catalog. Spectral parameters of the sources within 5 degrees from the target were set as free parameters, and parameters of other sources were fixed at their catalog values. GRO J1008$-$57 was included in the source model as a point source with power-law emission. In addition, the background Galactic and extragalactic diffuse emission were added in the source model using the spectral model gll_iem_v06.fits and file iso_P8R2_SOURCE_V6_v06.txt respectively. The normalizations of the diffuse components were free parameters in the analysis.

2.2. Source Search in Whole Data

Using the LAT science tools software package v11r5p3, we first performed standard binned likelihood analysis to the LAT data for GRO J1008$-$57. Since below 300 MeV, the instrument response function of the LAT has relatively large uncertainties and Galactic background emission is also strong, we chose events above 300 MeV for the likelihood analysis. The obtained Test Statistic (TS) value at the position of GRO J1008$-$57 was $\simeq$11. The TS value at a position indicates the fit improvement for including a source, and is approximately the square of the detection significance of the source (Abdo et al., 2010). Therefore we found possible detection with $>$3$\sigma$ significance.

We checked whether the $\gamma$-ray emission was due to contamination by nearby sources which were not included in the Fermi LAT 4-year catalog. A preliminary LAT 8-year point source list was released in early 2018, which contains nearly 2500 new sources, although it is not encouraged to use this list directly222https://fermi.gsfc.nasa.gov/ssc/data/access/lat/fl8y/. In addition, the extended source templates and the Galactic and extragalactic diffuse emission models have not been updated accordingly. We thus only added nearby new sources to the source model, and re-performed the maximum likelihood analysis. In this analysis, we only found TS$\simeq$7 at the position of GRO J1008$-$57, which is low for a possible detection.

2.3. Orbital-phase Resolved Search

Because Be XRBs often show enhanced emission at certain orbital phase ranges, particularly indicated by those $\gamma$-ray binaries, we searched for possible $\gamma$-ray emission in 10 orbital phase ranges of GRO J1008$-$57 (i.e., 0.0–0.1, …, 0.9–1.0, with phase zero at the periastron). Likelihood analysis to the data in each of the phase bins was performed. We found during phase 0.8–0.9 and 0.5–0.7 (defined as Phase I and II, respectively), the source was possibly detected with TS$>$9. In Figure 1, its 0.3–500 GeV orbital light curve and TS curve in 10 orbital-phase bins are shown. For the data points not in Phase I and II, their TS values were smaller than 5 ($<$2$\sigma$ significance), and the 95% (at 2$\sigma$ level) flux upper limits were derived.

2.3.1 Phase I

The $\gamma$-ray emission was significantly detected with a TS value of $\sim$20. Photon index $\Gamma=2.2\pm$0.2 and 0.3–500 GeV flux $F_{0.3-500}=13\pm 4\times 10^{-12}$ erg s*-1* cm*-2* were obtained during Phase I from the maximum likelihood analysis. We extracted the $\gamma$-ray spectrum by performing maximum likelihood analysis of the LAT data in 10 evenly divided energy bands in logarithm from 0.1–500 GeV. In the extraction, the spectral normalizations of the sources within 5 degrees from GRO J1008$-$57 were set as free parameters, while the other parameters of the sources were fixed at the values obtained from the above maximum likelihood analysis. A point source with power-law emission was assumed for GRO J1008$-$57, and the $\Gamma$ value was fixed to 2. We kept only spectral data points with the flux values 2 times greater than the uncertainties, and otherwise derived 95% (at 2$\sigma$ level) flux upper limits. The obtained flux and TS values of the spectral data points are given in Table 1. We also evaluated the systematic uncertainties of the spectral points induced by the uncertainties of the LAT effective area and the Galactic diffuse model. The former are 10% at 100 MeV, 5% at 500 MeV, and 20% at 10 GeV (Rando & for the Fermi LAT Collaboration, 2009), and the latter is a dominant one in the full LAT energy band, which was estimated by varying the normalization by $\pm$6%. The obtained systematic uncertainties are also given in Table 1.

We also calculated the residual TS map during Phase I in the energy range of 0.5–39 GeV. This energy range covers the three spectral data points with TS$>$4 (Table 1). The TS map is shown in the left panel of Figure 2. There is $\gamma$-ray emission present at the position of GRO J1008$-$57 with a TS of $\sim$22. We ran gtfindsrc to determine the position of the $\gamma$-ray emission, and obtained R.A.=152$\fdg$5, Decl.=$-$58$\fdg$3, (equinox J2000.0), with 1$\sigma$ nominal uncertainty of 0$\fdg$1. GRO J1008$-$57 is 0$\fdg$07 from this position and within the 1$\sigma$ error circle, indicating the likely association. As a comparison, the 0.5–39 GeV TS map excluding the data in Phase I was also calculated (shown in the right panel of Figure 2). No $\gamma$-ray emission is present as TS$\sim$0 at the position of GRO J1008$-$57. The comparison of the two TS maps support the detection of GRO J1008$-$57 in Phase I.

2.3.2 Phase II

There are also marginal detections with $>$2$\sigma$ significance during phase 0.5–0.7. We searched for possible detection during this wide phase range, and found that the $\gamma$-ray emission can be detected with a TS of 10 ($>$3$\sigma$ significance). $\Gamma=1.5\pm$0.3 and $F_{0.3-500}=9\pm 7\times 10^{-12}$ erg s*-1* cm*-2* were obtained from the maximum likelihood analysis. We extracted the $\gamma$-ray spectrum during Phase II, and the obtained flux and TS values of the spectral data points are provided in Table 1. There is actually only one spectral data point possibly detected, TS$\simeq$10 in energy band of 16.6–38.8 GeV. Comparing to the results obtained in Phase I, $\Gamma$ is lower in Phase II, although the uncertainties are too large for drawing a conclusion.

We calculated the residual TS map in 16.6–38.8 GeV band during Phase II, which is shown in the middle panel of Figure 3. The TS maps in this energy band during Phase I and the remaining phase ranges were also calculated, which are shown in the left and right panel of Figure 3 respectively. While the residual at the position of GRO J1008$-$57 in Phase II has TS$\sim$10, the map is noisy with other residuals. The situation is worse in Phase I, as a residual with TS$\sim$9 is slightly off the position of GRO J1008$-$57 and there are several other residuals with similar TS values. In the remaining phase ranges, TS$\sim$0 at the source position. Based on the TS maps, we concluded that the result of the detection of GRO J1008$-$57 in Phase II is uncertain.

2.3.3 Possible Emission Contamination

GRO J1008$-$57 is located in a complex region with a few sources detected relatively nearby and the $\gamma$-ray emission could only be detected during some phase ranges. We thus checked whether the nearby catalog sources are variable, whose flux variations might affect our results. There are six nearby sources listed in the LAT 4-year catalog (see Figure 2). The brightest one (3FGL J1018.9$-$5856) is a $\gamma$-ray binary, two of them (3FGL J1016.3$-$5858 and 3FGL J1020.0$-$5749) are $\gamma$-ray pulsars (PSR J1016$-$5857 and PSR J1019$-$5749, respectively), and the other three are unidentified sources. We first checked the 0.1–300 GeV variability indices (Acero et al. 2015) of them, and found none of them has variability index greater than the threshold to be identified as a variable source. There are also eight new sources listed in the LAT 8-year point source catalog (see Figure 2). For them, we calculated the 0.3–500 GeV variability indices for two of them (FL8Y J1010.8$-$5751 and FL8Y J1013.6$-$5833, the two closest to GRO J1008$-$57) following the method used in Acero et al. (2015). Both of these two sources are unidentified sources, and have variability indices (80.3 for J1010.8$-$5751 and 82.4 for J1013.6$-$5833) smaller than the threshold (149.7 for 112 degrees of freedom) considered for a variable source. We then checked the 0.3–500 GeV fluxes of these nearby sources during the 10 orbital phase bins of GRO J1008$-$57. For each of the sources, the fluxes during each bins have $<$2$\sigma$ deviations from those derived from the likelihood analysis of the whole data, again indicating no significant variations. We also checked the SIMBAD Astronomical Database and found that except GRO J1008$-$57, there are only a few normal stars known within the 2$\sigma$ error circles of the position obtained during Phase I. Therefore no evidence was found to indicate any possible emission contamination from nearby sources.

We also checked possible contamination due to the two nearby pulsars. Ideally, we could gate off their pulsed emission for this analysis. However no ephemerides covering the whole $>$9 years of LAT observations are available for the two pulsars. We instead extracted the 0.5–39 GeV TS map during Phase I with these two pulsars kept in the map. The TS map (Figure 4) shows that excess $\gamma$-ray emission at the position of GRO J1008$-$57 is present with TS$\sim 28$ (comparing to TS$\sim 22$ when the two pulsars were removed; cf. Figure 2), and it is clearly resolved from that of the two pulsars; in other words, the excess is not due to non-clean removal of the pulsars.

In addition, we evaluated how much the uncertainty of the Galactic diffuse emission model, the dominant one among the systematic uncertainties, affected our detection result. By increasing the normalization of the model by 6%, we extracted the 0.5–39 GeV TS map during Phase I and found that the excess emission at the position of GRO J1008$-$57  is still present but with TS reduced to $\sim$18. The TS value corresponds to $\sim$4$\sigma$ detection significance, indicating that the detection result was not changed by considering the systematic uncertainty.

3. Discussion

Having analyzed nine years of the Fermi LAT Pass 8 data, we searched for $\gamma$-ray emission from Be X-ray binaries with known orbital parameters. The search was conducted in the whole data and orbital-phase resolved data. Only for GRO J1008$-$57, possible residual emission was found in the whole data, but without sufficiently high significance. In the orbital-phase resolved analysis, the GRO J1008$-$57 region was found to have excess $\gamma$-ray emission significantly detected during the orbital phase range of 0.8–0.9. The $\gamma$-ray emission may be described with a power law with $\Gamma\sim 2$, although the spectra suffer large uncertainties. The likely orbital dependence of the excess $\gamma$-ray emission supports its association with GRO J1008$-$57. (In addition, given that the Fermi LAT collaboration very recently published an updated source catalog (4FGL catalog) and updated the diffuse background models, we checked whether our results would be affected when the new source catalog and diffuse models were used. Nearly the same results were obtained. We conclude that the results presented here are valid without significantly changes.)

Since GRO J1008$-$57 is considered as a transient X-ray pulsar, with pulsed X-ray emission detected during the outbursts, accretion onto the neutron star certainly occurs at least in outbursts. When in quiescence, the X-ray luminosity of GRO J1008$-$57 is generally $>$1034 erg s*-1*, and in one Swift monitoring of the binary over a full orbital cycle, the luminosity stayed at a level of $\sim 10^{35}$ erg s*-1*, higher than that for the onset of the propeller phase (Tsygankov et al., 2017). It is thus likely that there is always sufficiently strong mass accretion in the accretion disk surrounding the neutron star, not allowing a possible switch for the neutron star from being accretion powered to rotation powered (such as in the accretion-powered millisecond pulsar binary SAX J1808.4$-$3658; for searches for $\gamma$-ray emission see Xing & Wang 2013; de Oña Wilhelmi et al. 2016). The scenario for the $\gamma$-ray production from PSR B1259$-$63/LS 2883 probably does not work for GRO J1008$-$57.

In addition to neutron star $\gamma$-ray binaries, another small group of so-called transitional millisecond pulsar (MSP) binary systems, with PSR J1023+0038 as a prototype (Archibald et al., 2009), are also known to produce much enhanced $\gamma$-ray emission during their active state when an accretion disk is present around the pulsar (Wang et al., 2009; Stappers et al., 2014; Takata et al., 2014). Although the radiation mechanism for the enhanced $\gamma$-ray emission in the active state is not certain, the scenarios of an active (or partially active; see Xing et al. 2018) radio pulsar interacting with the accretion disk (Takata et al., 2014) or a propellering neutron star (Papitto et al., 2014; Papitto & Torres, 2015) have been proposed. Possible evidence for the former scenario is the coherent X-ray and optical pulsations detected in emission from PSR J1023+0038 in the active state (Archibald et al. 2015; Ambrosino et al. 2017; note no radio pulsations are seen). However, the distance to GRO J1008$-$57 was estimated to be 5.8$\pm$0.5 kpc (Riquelme et al. 2012; note its $V\simeq 15$ mag, but no distance information is provided in Gaia Archive), based on which the $\gamma$-ray luminosity due to the detection in Phase I is $L_{\gamma}\simeq 3.6\times 10^{34}$ erg s*-1* in 0.5–39 GeV (or $\simeq 6.7\times 10^{34}$ erg s*-1* in 0.1–300 GeV). This high luminosity value cannot be explained with the scenarios. Even assuming a pulsar wind might be turned on in GRO J1008$-$57, the long-term spin-down rate of the neutron star would only provide a rotational power of $\sim 10^{31}$ erg s*-1* (Kühnel et al., 2013), too low to produce the observed $\gamma$-ray emission (J. Takata, private communication). Similarly in the propellering neutron star scenario, the predicted $\gamma$-ray luminosity of GRO J1008$-$57 would be $\sim 10^{30}$ erg s*-1* (cf. Section 3.1 in Papitto & Torres 2015).

Among previous searches for $\gamma$-ray emission from neutron star X-ray binaries (e.g., Xing & Wang 2013; de Oña Wilhelmi et al. 2016), Li et al. (2012) have provided a very similar case to GRO J1008$-$57 by justifying that the $\gamma$-ray transients GRO J1036$-$55 and AGL J1037$-$5708 are associated with the HMXB 4U 1036$-$56. The likely association indicates that the HMXB only had occasional (detectable) $\gamma$-ray emission during several days. Following this feature, we performed the likelihood analysis to the LAT data in each time intervals of 0.1 orbital phase of GRO J1008$-$57  (i.e., 24.9 days). We found that there were three possible detection events with TS$>$9, which are listed in Table 2. Two were in phase 0.8–0.9 and one in 0.5–0.6, and their occurrences are shown in Figure 5, comparing to the X-ray outburst activity. It can be seen that the detection with the highest TS value of 17.5 was in the phase of 0.8–0.9, after a bright type-II outburst (occurred in November 2012 at an orbital phase of $\sim$0.3) during the same orbital cycle. The other two possible events only had TS$\simeq$9, too marginal to be considered. We made the TS maps for the TS=17.5 detection and showed them in Figure 6, indicating clear excess emission at the position of GRO J1008$-$57. A fine light curve (5 day binned) during this orbital phase range of 0.8–0.9 was constructed, but no clear pattern, such as a flaring event, could be identified (which actually excludes possibilities of contamination by a Solar flare or a $\gamma$-ray burst; both would appear as bright and sharp flaring events in daily LAT light curves). Therefore similar to that seen in 4U 1036$-$56, the occurences of the $\gamma$-ray emitting events did not coincide with any X-ray outbursts.

In Bednarek (2009), a model for $\gamma$-ray emission from neutron star HMXBs was proposed, in which a turbulent and magnetized region near an accreting neutron star would be a site accelerating particles to relativistic energies. Li et al. (2012) applied this model to 4U 1036-56, but the energy, considered to be transferred from the rotating neutron star to the high-energy particles in the model, is not large enough for the observed $\gamma$-ray emission from the HMXB, and thus they suggested that the available accretion energy should be considered. For GRO J1008$-$57, the same situation presents. The available power $L_{\rm rot}$ from the rotating neutron star is low, $L_{\rm rot}\sim 3\times 10^{30}B^{8/7}_{12}\dot{M}^{3/7}_{16}P^{-2}_{100}$ erg s*-1* (Bednarek, 2009; Li et al., 2012), where $B_{12}$ is the magnetic field of the neutron star in units of $10^{12}$ G, $\dot{M}_{16}$ is the mass accretion rate in units of $10^{16}$ g s*-1*, and $P_{100}$ is the spin period in units of 100 sec. The accretion power at the magnetospheric radius (i.e., at which the accretion disk is truncated by the magnetic field of the neutron star) is $L_{\rm acc}\sim 4.6\times 10^{33}B^{-4/7}_{12}\dot{M}^{9/7}_{16}$ erg s*-1* (Li et al., 2012); when $\dot{M}_{16}\gtrsim 12$, $L_{\rm acc}\gtrsim L_{\gamma}$. Note that both 4U 1036$-$56 and GRO J1008$-$57 require high mass accretion rates in order to match the observed $\gamma$-ray luminosities, if the face values of the estimated source distances are used.

We conclude that the $\gamma$-ray emission likely from GRO J1008$-$57 was detected with the Fermi LAT. Similar to that in 4U 1036$-$56, the emission was not persistent but occasionally detectable. The detailed physical process that gives rise to such high-energy emission is not clear, as the current models considering an accreting neutron star have difficulty in explaining the observed relatively high luminosities. Further observational studies of these HMXBs at high-energies should help establish their emission properties more clearly and may provide a clue for identifying the emission mechanism.

This research made use of the High Performance Computing Resource in the Core Facility for Advanced Research Computing at Shanghai Astronomical Observatory, and also made use of the MAXI data provided by RIKEN, JAXA and the MAXI team. This research was supported by the National Program on Key Research and Development Project (Grant No. 2016YFA0400804) and the National Natural Science Foundation of China (U1738131, 11633007).