A LUMINOUS GAMMA-RAY BINARY IN THE LARGE MAGELLANIC CLOUD

All gamma-ray observations were obtained with the LAT on board the Fermi satellite (Atwood et al. 2009). The LAT is a pair-conversion telescope, sensitive to gamma-ray photons with energies between ∼20 MeV and > 300 GeV. The LAT data used here were obtained between 2008 August 5 and 2016 March 24 (MJD 54,683 to 57,471). The initial search for gamma-ray binaries was performed with a data set covering a somewhat briefer time span, with the same start date up to 2015 August 27 (MJD 57,261). Analysis was performed using version v10r0p5 of the Fermi Science Tools, with Pass 8 "Source" class, front plus back data, for an energy range of 100 MeV to 300 GeV.

The third Fermi LAT catalog (3FGL, Acero et al. 2015) contains 3033 sources. To search for gamma-ray binaries we created light curves of all 3FGL sources and then calculated their power spectra to investigate the presence of periodic modulation. For the strongest peak in each power spectrum, the False Alarm Probability (FAP, Scargle 1982), i.e., the estimated probability of a signal reaching a power level by chance, under the assumption of white noise, was calculated. This FAP takes into account the number of independent frequencies searched, but does not include the effect of searching for periodicity in multiple sources. The light curves, covering an initial energy range from 100 MeV to 500 GeV, were created using a modified version of aperture photometry. In this version, where the probability that a photon comes from a source of interest is summed, rather than simply the number of photons (Fermi LAT Collaboration et al. 2012, and references therein). To estimate the probabilities, models were created for each source using the 3FGL catalog, using sources within a 10° radius and make3FGLxml. Photon probabilities were calculated using gtsrcprob and then summed for a 3° radius aperture centered on each source. Time bins of 500s were used for all sources.

Power spectra of weighted-photon aperture photometry LAT light curves were calculated, weighting each data point's contribution by its relative exposure, after first subtracting the mean count rate. This is required because of the large exposure changes from time bin to time bin, which are particularly apparent because of the short time bins (Fermi LAT Collaboration et al. 2009). For each source, the calculated power spectrum covered a period range of 0.05 days (1.2 hr) to the length of the light curve, i.e., ∼2788 days for the full data set and ∼2578 days for the initial search. The power spectra were oversampled by a factor of 5, compared to the nominal frequency resolutions of ∼1/2788 days⁻¹ and ∼1/2578 days⁻¹, respectively. Since the background is not fitted for each bin, artifact signals can be seen at Fermi's orbital period, at the survey period at twice that orbital period, at one day, at the Moon's 27.3 day sidereal period, at the 53 day precession period of the Fermi satellite, and at one quarter of a year.¹²

After detecting likely modulation from the region of the LMC, light curves were then generated by using a model for the LMC (Ackermann et al. 2016). Weighting each photon by its probability of coming from a source of interest has been found to increase the signal-to-noise level of the light curve (Kerr 2011; Fermi LAT Collaboration et al. 2012). However, it can reduce the apparent modulation of the light curve compared to its actual variability. When a source is brighter than its average level, the probability that a photon came from the source is underestimated; conversely the probability will be overestimated when the source is fainter than its mean level. This effect is especially pronounced for a faint source in the presence of brighter emission from neighboring sources or overall background level. For this reason, likelihood analysis provides a much more reliable estimate of the amplitude of source variability.

The Swift X-ray Telescope (XRT; Burrows et al. 2005) is a Wolter I X-ray imaging telescope, with a focal length of 3.5 m, fitted with a CCD chip covering a region of 236 × 236. The energy resolution FWHM in the XRT at the time of launch was ∼140 eV at 5.9 keV. Sensitive to X-rays ranging from 0.3 to 10 keV, the effective area of the XRT is ∼125 cm² at 1.5 keV. At 8.1 keV, the effective area of the XRT is ∼20 cm². The count rate for a source of 1 mCrab is 0.7 counts s⁻¹ (Hill et al. 2004).

XRT observations of CXOU J053600.0–673507 took place from 2015 November 21 to 2016 January 19 (MJD 57,347–57,406), with exposures ranging from ∼1.1 ks to ∼4.9 ks. Table 1 gives the observation log. The data were reduced and analyzed using the Swift XRT product generator (Evans et al. 2007) and the standard criteria given in the Swift XRT Data Reduction Guide (Capalbi et al. 2005). These procedures are described below.

Table 1.  Swift XRT Observation Log of CXOU J053600.0-673507

ObsID	Start Time (UT)	End Time (UT)	Phasea	Exposureb	Count Ratec	Fluxd
00034169001	2015 Nov 21 09:22:58	2015 Nov 21 11:57:36	0.898–0.908	2.5	<0.39e	<5.2e
00034169002	2015 Nov 24 18:51:58	2015 Nov 24 23:02:51	0.227–0.244	2.6	0.3 ± 0.1	${3.4}_{-1.3}^{+1.7}$
00034169003	2015 Nov 27 20:15:59	2015 Nov 27 22:59:00	0.524–0.535	2.4	0.7 ± 0.2	${9.2}_{-2.3}^{+2.8}$
00034169004	2015 Nov 30 13:41:58	2015 Nov 30 16:17:56	0.789–0.799	2.4	0.2 ± 0.1	${3.0}_{-1.3}^{+1.8}$
00034169005	2015 Dec 03 11:58:57	2015 Dec 03 14:27:02	0.073–0.083	1.4	<0.52e	<6.98e
00034169006	2015 Dec 08 09:58:58	2015 Dec 08 12:26:48	0.550–0.560	1.6	${0.8}_{-0.2}^{+0.3}$	${10.9}_{-3.0}^{+3.7}$
00034169007	2016 Jan 09 07:47:58	2016 Jan 09 09:03:12	0.648–0.653	1.1	${1.0}_{-0.4}^{+0.5}$	${13.0}_{-4.8}^{+6.3}$
00034169008	2016 Jan 11 04:35:58	2016 Jan 11 20:27:05	0.829–0.894	4.0	0.3 ± 0.1	${4.6}_{-1.3}^{+1.6}$
00034169009	2016 Jan 13 02:40:42	2016 Jan 13 08:20:50	0.016–0.039	4.1	<0.27e	<3.62e
00034169010	2016 Jan 15 04:16:57	2016 Jan 15 14:35:20	0.216–0.258	4.7	0.5 ± 0.1	6.8 ± 1.6
00034169011	2016 Jan 17 07:13:58	2016 Jan 17 13:41:43	0.423–0.449	4.9	0.8 ± 0.2	10.5 ± 2.0
00034169012	2016 Jan 19 01:10:57	2016 Jan 19 19:35:58	0.592–0.667	3.6	0.5 ± 0.1	6.7 ± 1.8

Notes.

^aPhase zero is defined as the epoch of maximum flux in the Fermi LAT. ^bThe net exposure time spread over several snapshots. Units are ks. ^cCount Rate is in the 2–10 keV energy band. Units are 10⁻² counts s⁻¹. Errors are at the 1σ level. ^dUnabsorbed Swift XRT flux in the 2.0–10.0 keV bandpass converted with PIMMS. Units are 10⁻¹³ erg cm⁻² s⁻¹. ^e3σ upper limits.

Download table as:  ASCIITypeset image

CXOU J053600.0–673507 was observed in Photon Counting (PC; Hill et al. 2004) mode, with a readout time of 2.5 s, adopting the standard grade filtering (0–12 for PC). Data were reduced and screened using xrtgrblc and xrtgrblcspec in HEAsoft v.6.19. The data were reprocessed with the XRTDAS data pipeline package xrtpipeline, using the standard filtering procedure to apply the newest calibration and default screening criteria. The source spectra were extracted from count-dependent regions generated by xrtgrblc. An unrelated field source was found at (J2000) R.A. = 05^h35^m474, decl. = −67°31'554, but was excised from the background extraction. Because of the presence of the SNR, we restrict our X-ray analysis to energies above 2 keV where SNR emission is reduced (Bamba et al. 2006).

The ancillary response files, accounting for vignetting, point-spread function correction, and different extraction regions, were generated and corrected for exposure using the FTOOLS packages xrtmkarf and xrtexpomap, respectively. We find the count rates in the observations to be ∼(1.1–3.4) × 10⁻² counts s⁻¹ (full energy range, source plus background), which is significantly less than the count rate of 0.5 counts s⁻¹ where pileup becomes important.

Individual spectra were not useful for analysis, because each spectrum was found to have 22–123 counts in the 0.3–10 keV energy band. Therefore, a cumulative spectrum was extracted, with a total of 784 counts. Photons in the 2–10 keV energy band were considered to reduce contamination from the DEM L241 SNR (Bamba et al. 2006). This gave a total of 119 counts for this energy band, and the total exposure is ∼34.9 ks. Data in the spectral files produced by xselect were further processed using grppha, which is designed to define the binning, quality flags, and systematic errors of the spectra. We used the quality flag to further eliminate bad data from the PHA files. Initially, we grouped the bins to ensure a minimum of 20 counts to fit the spectra, using ${\chi }^{2}$ statistics. However, insufficient bins were produced in the resulting spectrum. Due to the small number of counts, we therefore used the "C" statistic (Cash 1979) for the spectral analysis. The cumulative spectrum was grouped to have five counts per bin. The spectra were analyzed using XSPEC v12.9.0k. We made use of the XSPEC convolution model cflux to calculate the fluxes and associated errors of CXOU J053600.0–673507.

Radio observations were obtained using ATCA (Wilson et al. 2011). Dedicated follow-up observations were made between 2015 November 29 and 2016 March 12 (MJD 57,355 to 57,459, see Table 2), with observations centered at 5.5 and 9.0 GHz, with 2 GHz bandwidths for both bands. The ATCA, which consists of six 22 m diameter antennas, was in several different array configurations over this period: 1.5A (minimum baseline 153 m, maximum baseline 4.5 km), 750C (46 m, 5.0 km), EW352 (31 m, 4.4 km), and 6B (214 m, 6.0 km). Observations were reduced following standard procedures in Miriad (Sault et al. 1995), with the flux density scale set by observations of calibrators PKS 1934–638 and/or PKS 0823–500. Initially, mosaiced observations were made covering both the nominal position of the gamma-ray source P3, and the proposed X-ray counterpart CXOU J053600.0–673507. Observations differed in length, hour-angle coverage, angular resolution, and sensitivity to extended radio emission in the vicinity, resulting in a heterogeneous data set. The final observations of the series were conducted as targeted observations rather than mosaics of the region, with both a source inside the LAT error region and CXOU J053600.0–673507 radio sources observed along with the phase calibrator PKS 0530–727.

We also examined the Australia Telescope On-line Archive (atoa.atnf.csiro.au) for other observations of this region. Mosaiced observations were made in 2002, at 4.8 GHz with a 128 MHz bandwidth, with simultaneous, but spatially undersampled, observations at 8.4 GHz. Investigation of this data did not reveal any counterpart stronger than several mJy, with the image dominated by the nearby bright H ii region PKS 0535−676, ∼1.5 Jy at 4.8 GHz.

2.4.1. SOAR

2.4.2. SAAO Observations

2.4.3. OGLE Photometry

We examined the power spectra of the light curves of all 3FGL sources to search for evidence of periodic modulation that could be the sign of a binary system. The binary systems LS I+61°303, LS 5039, and 1FGL J1018.6–5856 were all detected at extremely high levels of statistical significance (ratio of peak power to mean power >100; FAP < 10⁻⁹). In addition, these binary signals were also often strongly seen in the power spectra of nearby sources. For candidate new binaries, our threshold for further investigation was for a source to have a peak power ≥18 × mean power level (FAP < 5 × 10⁻⁴), in addition for the period not to coincide with a known artifact. Some sources were found with modulation above this threshold, which could be interpreted as being due to non-periodic modulation from an active galactic nucleus (AGN). Such sources could be identified if the source was already associated with an AGN, if the source was located in a region where the presence of an early-type star would be unlikely—such as far from the Galactic plane, or the peak appeared to be part of underlying low-frequency noise, and thus would be due to non-periodic variability. A five sigma detection of a binary nominally arises with a relative peak height of >25. However, to ensure that a signal does not have some other cause—such as an unknown artifact—full confidence in the discovery of a new binary requires the identification of a counterpart at other wavelengths, and the detection of the same period in the counterpart.

From the power spectra of the 3FGL light curves, we noted a peak near 10.3 days (height ∼22 × mean power with an FAP of ∼10⁻⁵) from 3FGL J0526.6–6825e. This is the LMC, which is treated as a single object in the 3FGL catalog. A "difference" image, made by subtracting an image of 10.3 day minimum phase only from an image of 10.3 day maximum phase only, showed the modulation to be localized, and somewhat offset from the catalog position at roughly right ascension = 840, declination = −6755 (J2000). The modulation center was near "P3," an unassociated source recently found in LAT observations of the LMC by Ackermann et al. (2016), who noted several potential counterparts, including the SNR DEM L241 just outside the localization error region.

The spectral-spatial model for gamma-ray emission from the LMC (Ackermann et al. 2016) was used to create probability-weighted aperture photometry light curves, including a position centered on DEM L241. The 10.3 day peak in the power spectrum increased to ∼44× mean power (Figure 1) with an FAP of less than 10⁻¹⁰, strongly suggesting this was the location of the modulated gamma-ray source. The period was found to be 10.301 ± 0.002 days, with an epoch of maximum flux for sinusoidal modulation of MJD 57,410.25±0.34 (Figure 1). We also performed a likelihood fit to phase-resolved LAT spectra that revealed a large modulation amplitude with a profile that is more complex than purely sinusoidal (Figure 2). For the phase-resolved likelihood analysis, we again used the LMC model (Ackermann et al. 2016) for a 10° radius, centered on P3, for an energy range of 100 MeV to 300 GeV. The parameters of all sources—apart from P3—were frozen to their previously determined best fit values. Leaving the power law index of P3 free for the phase-resolved likelihood fits resulted in unreasonably steep spectra with power-law indices of ≥5 for phases of lower flux. We, therefore, also froze the power-law index to the value of 2.77 from Ackermann et al. (2016). To obtain a more robust measurement of spectral variability, it may be necessary to rederive a model for this complex region using Pass 8 data and knowledge of the precise location of P3. Such an analysis is beyond the scope of this paper.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

The XRT observations of CXOU J053600.0–673507 monitor the source for more than three orbital periods. We binned the light curves to a resolution of one bin per observation to investigate the X-ray orbital modulation of the system. The Swift observations revealed strong, approximately sinusoidal, modulation on the 10.3 day gamma-ray period (Figure 3). However, X-ray minimum occurs near the phase of gamma-ray maximum.

To fit the cumulative XRT spectrum, we used several models that are used to describe systems that host a neutron star: a power law (see Figure 4), a power law with a high-energy cutoff (highecut × power in XSPEC), and a cutoff powerlaw (cutoffpl in XSPEC). All models were modified by an absorber that fully covers the source using appropriate cross-sections (Balucinska-Church & McCammon 1992) and abundances (Wilms et al. 2000).

Zoom In Zoom Out Reset image size

The model that provides a good fit (C statistic of 20.15 for 19 degrees of freedom) to the data is a power law with photon index Γ = 1.3 ± 0.3, modified by a fully covered absorber (tbabs in XSPEC). We find the unabsorbed flux of CXOU J053600.0–673507 to be (3.2 ± 0.4) × 10⁻¹³ erg cm⁻² s⁻¹ in the 2.0–7.5 keV band. Our photon index is thus consistent with the values of 1.51–1.62, found by Bamba et al. (2006), and 1.28 ± 0.08, reported by Seward et al. (2012).

We found that the neutral hydrogen column density for the fully covered absorption could not be accurately constrained. This is not surprising, as only energies above 2 keV are included in the analysis. Therefore, we froze ${N}_{{\rm{H}}}$ to 1.9 × 10²¹ atoms cm⁻², which was found from Chandra data (Seward et al. 2012). This value of the fully covered absorber is comparable to the Galactic HI value of 1.62 × 10²¹ atoms cm⁻², as given in the Leiden/Argentine/Bonn survey (Kalberla et al. 2005). Therefore, we assume that the fully covered absorber is interstellar in origin. While a good fit does not require a high-energy cutoff, which is typically found in accreting pulsars, our spectra are limited to energies below 10 keV; such cutoffs are often seen at higher energies (Coburn et al. 2002).

In Table 4, we summarize our current and also previous X-ray observations of CXOU J053600.0–673507. We note that the two Chandra observations from Seward et al. (2012) showed an increase in the flux between 0.5 and 5 keV from (3.71 ± 0.10) to (4.70 ± 0.12) × 10⁻¹³ erg cm⁻² s⁻¹. The observation times correspond to phases of ∼0.17 and ∼0.26 and the flux increase is thus consistent with the behavior we observe at these phases with the Swift XRT. However, the XMM observation of Bamba et al. (2006) was obtained at a phase of ∼0.6, which corresponds to orbital maximum but did not show a particularly high flux. This could be indicative of cycle-to-cycle variability if confirmed by additional observations.

Table 4.  Summary of X-Ray Observations

Mission/Instrument	Time	Γ	Flux	${L}_{{\rm{x}}}$ a	Reference
	(MJD)		(10−13 erg cm−2 s−1)	(1035 erg s−1)
XMM-Newton/EPIC	53368	${1.57}_{-0.06}^{+0.05}$	6.4 ± 0.4b	2.32 ± 0.14	Bamba et al. (2006)
Chandra/ACIS	55599	1.28 ± 0.08	3.71 ± 0.10c	2.52 ± 0.07	Seward et al. (2012)
Chandra/ACIS	55600	1.28 ± 0.08	4.70 ± 0.12c	3.19 ± 0.08	Seward et al. (2012)
Swift/XRT	57347–57406	1.3 ± 0.3	<3.62–13.0d	<1.09–3.9	Present work

Notes. Summary of the X-ray spectral parameters derived from the XMM-Newton (Bamba et al. 2006), Chandra (Seward et al. 2012), and Swift observations (this work).

^aLuminosity converted to the 0.3–10.0 keV bandpass with PIMMS assuming derived spectral parameters for a distance of 50 kpc. ^bAbsorbed XMM-Newton flux in the 0.5–10.0 keV bandpass. ^cAbsorbed Chandra flux in the 0.5–5.0 keV bandpass. ^dUnabsorbed Swift XRT flux in the 0.3–10.0 keV bandpass. Upper limits are 3σ.

Download table as:  ASCIITypeset image

A point source was found in the LAT error region of the gamma-ray source; and another was detected at the position of CXOU J053600.0–673507. The radio flux densities folded on the 10.3 day period (Figure 3) show modulation of the emission on the gamma-ray period, but again out of phase with the gamma-ray modulation. We note that the folded 5.5 GHz light curve is more "noisy" than that at 9 GHz. In particular, at phase ∼0.68 there is a low flux 5.5 GHz point without a corresponding decrease at 9 GHz. The 9 GHz light curve is expected to be less susceptible to systematic effects due to the smaller beam size at this wavelength that results in reduced contamination from other LMC sources detected in sidelobes of the telescope response. However, we note that the large scatter in 5.5 GHz fluxes seen around the phase range ∼0.7–0.8 might also suggest cycle-to-cycle variability.

3.5.1. SOAR and SAAO Optical Spectroscopic Results

The overall optical continuum and line strengths in the SOAR and SAAO spectra showed no obvious changes from previous observations. Radial velocities were determined from the SOAR spectra and these are shown in Figure 5, folded on the 10.3 day period. We find a clear orbital variation consistent with binary motion.

Zoom In Zoom Out Reset image size

We performed circular Keplerian fits to the radial velocities with the period fixed to the value determined from analysis of the LAT data: 10.301 days. Two sets of fits were performed: one with the phase fixed, such that the peak gamma-ray flux occurs when the compact object is behind the O star; the other with the phase free. For the former fit, we find K₂, the semi-amplitude of the radial velocity due to the orbital motion of the O star, $\,=\,5.5\pm 2.7$ km s⁻¹, i.e., marginal evidence of radial velocity variations with this period and phase.

With the phase free, the best-fit parameters are: systemic velocity 295.8 ± 2.0 km s⁻¹ and ${K}_{2}=10.7\pm 2.4$ km s⁻¹. This fit is nonetheless poor (${\chi }^{2}/\nu$ = 70/10; rms = 5.2 km s⁻¹), suggesting that residual variations due to, e.g., a time-variable wind are present. The systemic velocity is consistent with that expected for the LMC disk at this location (277 ± 20 km s⁻¹; van der Marel et al. 2002). We determined the uncertainties in the time of conjunction, K₂, and systemic velocity using a standard bootstrap.

Constraints on the mass of the compact object in the system as a function of inclination angle (i) were obtained by calculating the "mass function" (e.g., Strader et al. 2015), i.e., $f{(M)={{PK}}_{2}^{3}/{(2\pi G)=({M}_{1}\sin i)}^{3}/({M}_{1}+{M}_{2})}^{2}$, where M₁ and M₂ are the masses of the compact object and the O star, respectively, and G is the gravitational constant. The mass function is $f(M)=({1.3}_{-0.6}^{+1.1})\times {10}^{-3}\,{M}_{\odot }$.

Assuming the O star has a mass between 25 and 42 M_⊙ (Seward et al. 2012; Martins et al. 2005) and a neutron star mass of $1.4\,{M}_{\odot }$, the inclination is ${50}_{-13}^{+16}^\circ$. For a neutron star of $2.0\,{M}_{\odot }$, the inclination is ${35}_{-9}^{+11}^\circ$. We cannot rule out that the compact object is a black hole, but it would require a low inclination: ${14}_{-3}^{+4}^\circ$ for $5\,{M}_{\odot }$ and 8 ± 2° for $10\,{M}_{\odot }$.

As a check, we also fit the radial velocities with both the period and phase free. The best-fit period is 10.1 days, slightly smaller than, but consistent with, the more precise gamma-ray period. We also investigated eccentric orbit fits to the radial velocities. We find that for eccentricities less than the estimated upper limit of 0.7, there is no improvement in the reduced ${\chi }^{2}$ of the fits. Additional radial velocity measurements will be required to better constrain or measure the eccentricity of the system.

Very low inclinations can, in principle, be constrained through a comparison of the projected rotational velocity to the breakup velocity. We do not have the appropriate comparison data to make a precise estimate of the projected rotational velocity; we roughly estimate a value of about 80 km s⁻¹, which would suggest the inclination cannot be lower than about 6^◦, but we emphasize that this estimate is uncertain.

In sum, the radial velocity data are consistent with a neutron star and a wide range of inclinations, with a black hole allowable if the system is very close to face on. Observations at a wider range of phases are necessary to fully constrain the system parameters.

For the SAAO observations, although the radial velocities inferred are generally consistent with those obtained from the SOAR observations, we do not use the SAAO spectra here for determination of orbital parameters because of small altitude-angle dependent effects on wavelength calibration. The altitude angle effects will hopefully be removed in the future.

3.5.2. Previous SAAO Optical Spectroscopy

The OGLE light curve folded on the orbital period of CXOU J053600.0–673507 is shown in Figure 6. There is no statistically significant evidence for any orbital modulation in the folded OGLE data. This does not strongly constrain the system parameters, as the expected tidal distortion of the O star due to a 1.4 M_⊙ companion would be small. We note that for the gamma-ray binary LS 5039, which has a 3.9 day period, any photometric variability is below 2 millimagnitudes (Sarty et al. 2011). Thus, the lack of detection of photometric variability for CXOU J053600.0–673507, which has a longer period and a giant rather than main-sequence primary, is not surprising.

Zoom In Zoom Out Reset image size

The modulation of the X-ray and radio fluxes of CXOU J053600.0–673507 on the gamma-ray period confirms that it is the counterpart of LMC P3. The detection of radio emission is unusual for an HMXB, but common for gamma-ray binaries (Dubus 2015). The gamma-ray to X-ray luminosity ratio and identification with a massive star firmly classify P3 as a gamma-ray binary. The anti-correlation of the X-ray and gamma-ray orbital modulations is similar to that seen in the other two systems with O star companions, LS 5039 and 1FGL J1018.6–5856 (Hadasch et al. 2012; An et al. 2015; Dubus 2015). With a 10.3 day orbit and an O5III companion, Kepler's third law requires that the orbital eccentricity is less than 0.7 for the compact object to avoid hitting the surface of its 15 ${R}_{\odot }$ companion at periastron. The fact that these three systems show regular modulations when the three gamma-ray binaries with Be-star companions (LS I+61°303, HESS J0632+057, PSR B1259–63) show significant orbit-to-orbit variability (Hadasch et al. 2012; Caliandro et al. 2015; Aliu et al. 2014) is likely to be related to the lack of a circumstellar disk around the O star.

The compact object in the system, which was presumably formed in the supernova explosion that also created DEM L241, appears to be a neutron star. Thus, the progenitor of the SNR, which has been proposed to have had a mass $\gtrsim 20\,{M}_{\odot }$ (Bamba et al. 2006) was still below the threshold for the formation of a black hole as considered by Seward et al. (2012).

The process resulting in gamma-ray emission must be steady, given the lack of evidence for long-term X-ray and gamma-ray variability. As with the other gamma-ray binaries, the level of emission and modulations are more likely imprinted by the interaction of a pulsar wind with the stellar wind of the O-star companion. The modulations are thought to be due to a combination of anisotropic inverse Compton scattering and relativistic Doppler boosting, with maximum high-energy gamma-ray emission at superior conjunction (Dubus 2015). However, we note that gamma-ray maximum occurs somewhat after superior conjunction (Figure 5), raising the possibility that the system may have an eccentric orbit. More precise orbital parameters will be needed to test this. The analogy with LS 5039 and 1FGL J1018.6–5856 also suggests that P3 may be a TeV source with a modulation correlated with the X-ray modulation. Scaling from the TeV/X-ray ratio of LS 5039, the maximum flux at 1 TeV is expected to be $\sim {10}^{-13}\,\mathrm{ph}\,{\mathrm{cm}}^{-2}\,{{\rm{s}}}^{-1}\,{\mathrm{TeV}}^{-1}$, only slightly fainter than the currently known TeV gamma-ray sources in the LMC (H.E.S.S. Collaboration et al. 2015), raising the possibility that it may be detectable with H.E.S.S.

LMC P3 is the most luminous gamma-ray binary observed yet. It is at least four times more luminous in GeV gamma-rays and 10 times more luminous in radio and X-rays than LS 5039 (Dubus 2015; Marcote et al. 2015) and 1FGL J1018.6–5856 (Fermi LAT Collaboration et al. 2012). The luminosity of the companion, a factor of 1.5 from O5V to O5III (Martins et al. 2005), and the orbital separations (0.1–0.4 AU) are comparable in all three systems. Hence, the higher luminosity is more likely to be due to the injected power in non-thermal particles rather than to higher radiation or matter densities.

The pulsar spin-down power must be $\dot{E}\geqslant {L}_{\gamma }\approx 4.3\times {10}^{36}$ erg s⁻¹, to account for the gamma-ray luminosity. For comparison, this is a factor of 5 greater than the spin-down power of PSR B1259–63, which apparently converts nearly all this power to high-energy gamma-rays over parts of its orbit (Caliandro et al. 2015). The pulsar age is also constrained to $\sim {10}^{5}$ years if it is associated with the surrounding SNR (Seward et al. 2012). Assuming no magnetic field decay and a spin-down power given by $\dot{E}={\mu }^{2}{{\rm{\Omega }}}^{4}/{c}^{3}$ (Spitkovsky 2006), with μ the pulsar magnetic moment and Ω its angular frequency, the magnetic field cannot be higher than $4\times {10}^{11}$ G, and the current spin period must be shorter than 39 ms in order to be consistent with $\dot{E}\geqslant 4.3\times {10}^{36}$ erg s⁻¹ 10⁵ years after the birth of the pulsar.

The lack of discovery of additional Milky Way gamma-ray binaries since 1FGL J1018.6–5856, but with the detection of an LMC source, suggests that we may have detected at least the majority of persistent higher-luminosity Galactic gamma-ray binaries, even if we caution that some gamma-ray binaries may be missing from the 3FGL catalog because of their low duty cycle (e.g., Caliandro et al. 2015) or because they are GeV-faint, such as HESS J0632+057 (Malyshev & Chernyakova 2016, and references therein). One possibility is that we detect as gamma-ray binaries only those systems with the fastest rotating neutron stars at birth, whereas most are born as slower rotators. For example, adopting a normal distribution for birth spin periods (Faucher-Giguère & Kaspi 2006) with a mean of 300 ± 150 ms, these neutron stars would be very faint in the rotation-powered stage, or would go directly to the HMXB stage, with no gamma-ray emission. In the latter case, earlier predictions of an extensive Galactic population of such objects based on the birth rate of HMXBs ($\sim {10}^{-3}\,{\mathrm{yr}}^{-1}$) and the expected lifetime of gamma-ray emission ($\sim {10}^{5}$ years) would have been optimistic (Meurs & van den Heuvel 1989).

Although LMC P3 is considerably brighter than its Galactic counterparts, we note that even this source would not be detectable at the ∼780 kpc distance of the Andromeda galaxy. Thus, unless significantly brighter sources exist, we are currently limited to gamma-ray detection of binaries within the Milky Way and its satellites.