{'fraction_non_alphanumeric': 0.03657259733619032, 'fraction_numerical': 0.0430008035257737, 'mean_word_length': 4.386935991605457, 'pattern_counts': {'":': 0, '<': 2, '<?xml version=': 0, '>': 4, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 1}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
A Luminous Gamma-ray Binary in the Large Magellanic Cloud

Gamma-ray binaries consist of a neutron star or a black hole interacting with a normal star to produce gamma-ray emission that dominates the radiative output of the system. Only a handful of such systems have been previously discovered, all within our Galaxy. Here we report the discovery with the Fermi Large Area Telescope (LAT) of a luminous gamma-ray binary in the Large Magellanic Cloud from a search for periodic modulation in all sources in the third Fermi LAT catalog. This is the first such system to be found outside the Milky Way. The system has an orbital period of 10.3 days and is associated with a massive O5III star located in the supernova remnant DEM L241, previously identified as the candidate high-mass X-ray binary (HMXB) CXOU J053600.0-673507. X-ray and radio emission are also modulated on the 10.3 day period, but are in anti-phase with the gamma-ray modulation. Optical radial velocity measurements suggest that the system contains a neutron star. The source is significantly more luminous than similar sources in the Milky Way at radio, optical, X-ray and gamma-ray wavelengths. The detection of this extra-galactic system, but no new Galactic systems raises the possibility that the predicted number of gamma-ray binaries in our Galaxy has been overestimated, and that HMXBs may be born containing relatively slowly rotating neutron stars.


INTRODUCTION
Although hundreds of interacting binary systems are known X-ray emitters (Liu et al. 2006(Liu et al. , 2007, very few systems produce detectable gamma-ray emission. We here classify gamma-ray binaries as those systems where most of the electromagnetic output of the system is at gamma-ray energies. In order to generate gamma-rays, non-thermal emission mechanisms are required such as the particle acceleration in a shock between the wind from a rapidly rotating neutron star and its companion (Dubus 2006), where the Fermi mechanism (Marcowith et al. 2016) may accelerate particles to high energies, or in the highvelocity jets from an accreting black-hole "microquasar" (e.g. Mirabel & Rodríguez 1998). Gamma-ray emission has been detected from the microquasar Cygnus X-3 (e.g. Fermi LAT Collaboration et al. 2009;Tavani et al. 2009;Corbel et al. 2012) and possibly from Cygnus X-1 (Bodaghee et al. 2013). However these two sources can be viewed as "gamma-ray emitting X-ray binaries" since emission predominantly occurs at X-ray energies (Dubus 2015).
Some evolutionary models of HMXBs predict that these systems pass through such a brief gamma-ray binary phase shortly after the formation of a short rotation period neutron star in a supernova explosion (Meurs & van den Heuvel 1989) in a binary system with an O or B spectral type companion. A population of gamma-ray binaries is thus expected, the observable number depending on factors including the duration of the gamma-ray binary phase and the gamma-ray luminosity. From their binary population synthesis study Meurs & van den Heuvel (1989) predicted about 30 binaries containing neutron stars during their pulsar phase which could thus be gamma-ray binaries. Following the launch of the Fermi Gamma-ray Space Telescope mission in 2008, its Large Area Telescope (LAT; Atwood et al. 2009) was used to confirm GeV gamma-ray emission from the systems PSR B1259-63, LS I +61 • 303, and LS 5039 which had been suspected to be members of this class of object (e.g. Dubus 2013, and references therein). All three of these systems display modulation of their gamma-ray fluxes on their orbital periods. This suggests that detection of periodic emission from a gamma-ray source can be a powerful way to find new binaries. From an earlier search for periodic modulation in cataloged Fermi LAT sources a 16.5 day period was found from 1FGL J1018.6-5856 which multi-wavelength observations confirmed to be a gamma-ray binary (Corbet et al. 2011;Fermi LAT Collaboration et al. 2012). These initial results suggested that the predicted population of ∼30 Galactic gamma-ray binaries might indeed exist. How-ever, since the discovery of 1FGL J1018.6-5856 in 2011 no additional gamma-ray binaries had been found with Fermi, possibly calling into question the number of such sources in the Galaxy. On the other hand, we note that the gamma-ray binary HESS J0632+057, although detectable at higher TeV energies, is at most only a faint source at the GeV energies accessible to Fermi (Malyshev & Chernyakova 2016), and the pulsar PSR J2032+4127 might become a gamma-ray binary at the periastron passage of its highly eccentric 20-30 year long orbit around its Be star companion (Lyne et al. 2015).
We describe here our program to search for additional gamma-ray binaries from the detection of periodic modulation in the Fermi LAT light curves of sources in the third Fermi source catalog. This has resulted in the discovery of periodic modulation with a 10.3 day period from the direction of the Large Magellanic Cloud (LMC). The LMC is a neighbor galaxy of the Milky Way located ∼50 kpc from the Earth (Macri et al. 2006;Pietrzyński et al. 2013;de Grijs et al. 2014). This source was then localized to a position that suggested identification with the "P3" point-like source in a LAT survey of the LMC by Ackermann et al. (2016).
We next identified the counterpart as a previously proposed high-mass X-ray binary (HMXB) CXOU J053600.0-673507 . A previous observation of the supernova remnant (SNR) DEM L241 surrounding CXOU J053600.0-673507 with ATCA (Bozzetto et al. 2012) had revealed a point-like source which had been interpreted, following an earlier suggestion ) based on XMM-Newton data, as a pulsar-wind nebula.
To confirm the identification of CXOU J053600.0-673507 as the counterpart of the periodic gamma-ray source we obtained X-ray and radio observations of CXOU J053600.0-673507, using the Swift satellite and the Australia Compact Telescope Array (ATCA) respectively, and found that the X-ray and radio fluxes from this source are also modulated on the 10.3 day period. CXOU J053600.0-673507 has a previously identified O5III(f) counterpart (Crampton et al. 1985;Seward et al. 2012). This V = 13.5 optical counterpart had previously been reported to show up to 30 km s −1 velocity differences from day to day Crampton et al. 1985). We therefore obtained optical radial velocity measurements with the Southern Astrophysical Research (SOAR) 4.1 m and South African Astronomical Observatory (SAAO) 1.9 m telescopes which allowed us to determine that the source probably contains a neutron star. This is thus the first gamma-ray binary to be found outside the Milky Way and also the first gamma-ray binary to be found with the Fermi LAT since 1FGL J1018.6-5856.
We describe LAT observations and our program to search for modulated gamma-ray sources in Section 2.1. The Swift X-ray observations of CXOU J053600.0-673507 are presented in Section 2.2 and the ATCA observations in Section 2.3. Optical spectroscopy to obtain radial velocity measurements is described in Sections 2.4.1 and 2.4.2 and optical photometry obtained with Optical Gravitational Lensing Experiment (OGLE) is described in Section 2.4.3. The discovery of periodic gamma-ray emission from the direction of the LMC is presented in Section 3.1, the identification of the counterpart as CXOU J053600.0-673507 in DEM L241 from the detection of modulated X-ray and radio emission is given in Sections 3.2 and 3.4 respectively. Constraints on the system from optical radial velocity measurements are given in Section 3.5.1 and OGLE upper limits on photometric modulation in Section 3.6. The nature of LMC P3/CXOU J053600.0-673507 and the implications of the discovery of a gamma-ray binary outside the Milky Way, but the lack of new sources in our Galaxy, for the overall population of gamma-ray binaries are discussed in Section 4 with an overall conclusion in Section 5. Unless otherwise stated, uncertainties are given at the 1σ level. For luminosity calculations we adopt a distance of 50 kpc.

Gamma-ray Observations and Analysis
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 to > 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 somewhat shorter data set obtained 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 calculated power spectra of these to investigate the presence of periodic modulation. For the strongest peak in each power spectrum the False Alarm Probability (FAP, Scargle 1982), 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 energy range of 100 MeV to 500 GeV, were created using a modified version of aperture photometry 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 degree radius and make3FGLxml. Photon probabilities were calculated using gtsrcprob and then summed for a 3 degree 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 from 0.05 days (1.2 hrs) 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 −1 and ∼1/2578 days −1 respectively. Since background is not fitted for each bin, artifact signals can be seen at Fermi's orbital period, the survey period at twice this, one day, the Moon's 27.3 day sidereal period, the 53 day precession period of the Fermi satellite, and one quarter of a year 11 .
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 (Fermi LAT Collaboration et al. 2012;Kerr 2011). 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.

X-ray Observations and Analysis
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 23.6 ′ ×23.6 ′ . 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 2 at 1.5 keV. At 8.1 keV, the effective area of the XRT is ∼20 cm 2 . The count rate for a source of 1 mCrab is 0.7 counts s −1 (Hill et al. 2004).
XRT observations of CXOU J053600.0-673507 took place from 2015 November 21 to 2016 January 19 (MJD 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.
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 h 35 m 47 s . 4, Dec = −67 • 31 ′ 55 ′′ . 4, which was excised from the background extraction. Because of the presence of the SNR 11 http://fermi.gsfc.nasa.gov/ssc/data/analysis/LAT caveats temporal.html we restrict our X-ray analysis to energies above 2 keV where SNR emission is reduced .
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 −2 counts s −1 (full energy range, source plus background), which is significantly less than the count rate of 0.5 counts s −1 where pileup becomes important.
Individual spectra were not useful for analysis, as each spectrum was found to have 22-123 counts in the 0.3-10 keV energy band. A cumulative spectrum was therefore extracted, which has a total of 784 counts. Photons in the 2-10 keV energy band were considered to reduce contamination from the DEM L241 supernova remnant ). 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 and 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 χ 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 5 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.  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 22m-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 Xray 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 together with the phase calibrator PKS 0530-727.

Radio Observations and Analysis
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.

Optical Observations and Analysis
Long-slit spectra of the optical counterpart of CXOU J053600.0-673507 were acquired using the Goodman spectrograph (Clemens et al. 2004) on the 4.1 m SOAR telescope. All observations used a 1.03 ′′ slit and a 2100 l mm −1 grating, yielding a resolution of 0.9Å over a wavelength range ∼ 4300-4970Å. At each epoch 20 minutes of exposure time was obtained. The spectra were optimally extracted and calibrated in the usual manner.
Qualitatively the spectra appear identical to those presented of the system in previous work (Crampton et al. 1985;Seward et al. 2012), with absorption lines of H and He and the Bowen blend in emission. We determined barycentric radial velocities through cross-correlation with a field O star taken with a similar setup around the region of the He II line at 4542Å.
For this paper, the first observations were obtained on 2015 December 10 (UT) and the last on 2016 April 16 (Table 3) comprising 13 separate epochs listed as Barycentric Julian Dates on the TDB system (Eastman et al. 2010). Nonetheless, due to the long period of the system the phase coverage is still not optimal.

SAAO Observations
Spectra were obtained using the newly upgraded spectrograph at the Cassegrain focus of the 1.9 m Radcliffe telescope at SAAO. The spectrograph's collimator, optics and detector were all replaced during the second half of 2015. All associated software was also rewritten during the upgrade.
Blue spectra were obtained using two gratings: The first has a broad spectral range of 800Å and a resolution of 1Å. The second has a second order spectrum in the blue, which allowed for an improved resolution of 0.5 A while sacrificing some spectral range and throughput. Observations using the two gratings were performed as part of science verification tests.
Three observations were obtained during November and December 2015 after the initial testing of the upgraded spectrograph. Exposure times of 1200 s were used for both gratings, achieving a good signal-to-noise ratio in each spectrum.

OGLE Photometry
Optical I-band data from OGLE (Udalski et al. 2015) Phase IV were used to investigate the long-term and periodic behavior of the optical counterpart in this system. The source is identified as LMC518.09.16278 within the OGLE IV phase of the project. The data cover a period of several years starting at MJD 55,260 and continuing to MJD 57,339. The coverage was almost nightly for the first 2-3 years, but more recently the frequency of observations has been reduced to approximately once every 3-4 nights.

Gamma-ray Results
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 −9 ). 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 −4 ) and for the period not to coincide with a known artifact. Some sources were found with modulation above this threshold but which could be interpreted as 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 so 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 −5 ) from 3FGL J0526.6-6825e. This is the LMC, that 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 = 84.0 • , declination = -67.55 • (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 −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 epoch of maximum flux for sinusoidal modulation of MJD 57,410.25 ± 0.34 ( Figure  2). We also performed a likelihood fit to phase-resolved LAT spectra which revealed a large modulation amplitude with a profile that is more complex than purely sinusoidal (Figure 3). For the phase-resolved likelihood analysis we again used the LMC model (Ackermann et al. 2016) for a 10 degree 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 the knowledge of the precise location of P3. Such an analysis is beyond the scope of this paper.

X-ray Results
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).
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 −13 erg cm −2 s −1 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 find that the neutral hydrogen column density for the fully covered absorption could not be accurately con- strained. This is not surprising as only energies above 2 keV are included in the analysis. Therefore, we froze N H to 1.9×10 21 atoms cm −2 , which was found from Chandra data . This value of the fully covered absorber is comparable with the Galactic HI value given in the Leiden/Argentine/Bonn survey (Kalberla et al. 2005), which is 1.62×10 21 atoms cm −2 . 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 and such cutoffs are often seen at higher energies (Coburn et al. 2002).

Comparison with Previous X-ray Observations
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 to 5 keV from (3.71 ± 0.10) to (4.70 ± 0.12) × 10 −13 erg cm −2 s −1 . 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 cycleto-cycle variability if confirmed by additional observa-

Radio Results
A point source was found in the LAT error region of the gamma-ray source together with the detection of a point source 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 mod- ulation. 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 which will result 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.

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.
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 d. 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 2 , the semi-amplitude of the radial velocity due to the orbital motion of the O star, = 5.5 ± 2.7 km s −1 , 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 −1 and K 2 = 10.7±2.4 km s −1 . This fit is nonetheless poor (χ 2 /ν = 70/10; rms = 5.2 km s −1 ), 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 −1 ; van der Marel et al. 2002). We determined the uncertainties in the time of conjunction, K 2 , 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" (see e.g. Strader et al. 2015), i.e. f (M ) = P K 3 2 /(2πG) = (M 1 sin i) 3 /(M 1 +M 2 ) 2 , where M 1 and M 2 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 +1.1 −0.6 ) × 10 −3 M ⊙ . Assuming the O star has a mass between 25 and 42 M ⊙ Martins et al. 2005) and a neutron star mass of 1.4M ⊙ , the inclination is 50 +16 −13 • . For a neutron star of 2.0M ⊙ the inclination is 35 +11 −9 • . We cannot rule out that the compact object is a black hole, but it would require a low inclination: 14 +4 −3 • for 5M ⊙ and 8 ± 2 • for 10M ⊙ .
As a check, we also fit the radial velocities with both the period and phase free. The best fit period is 10.1 d, 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 χ 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 −1 , 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, because of small altitudeangle dependent effects on wavelength calibration we do not use the SAAO spectra here for determination of orbital parameters. The altitude angle effects are hoped to be removed in the future.

Previous SAAO Optical Spectroscopy
Radial velocity measurements of the optical counterpart of CXOU J053600.0-673507 were obtained by Seward et al. (2012) in 2011 October and 2012 April, and these, along with additional measurements between 2012 Nov 2013 Jan are presented by . These observations confirm radial velocity changes, but with significantly larger uncertainties of ∼5 to > 10 km s −1 than our SOAR measurements. Folding these measurements on the 10.3 day orbital period gives ambiguous results due to the larger measurement uncertainties.

OGLE Photometric Results
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.4M ⊙ 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.

4.1.
Properties and Nature of LMC P3/CXOU J053600.0-673507 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 gammaray orbital modulations is similar to that seen in the other two systems with O star companions, LS 5039 and 1FGL J1018.6-5856 (Dubus 2015;An et al. 2015;Hadasch et al. 2012). 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 ⊙ companion at periastron. 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 20M ⊙ ) 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 ∼ 10 −13 ph cm −2 s −1 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 1.5 from O5V to O5III (Martins et al. 2005) and the orbital separations  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Ė ≥ L γ ≈ 4.3 × 10 36 erg s −1 , to account for the gamma-ray luminosity. For comparison, this is a factor 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 ∼ 10 5 years if it is associated with the surrounding SNR . Assuming no magnetic field decay and a spin-down power given byĖ = µ 2 Ω 4 /c 3 (Spitkovsky 2006), with µ the pulsar magnetic moment and Ω its angular frequency, the magnetic field cannot be higher than 4 × 10 11 G and the current spin period must be shorter than 39 ms in order to be consistent withĖ ≥ 4.3 × 10 36 erg s −1 10 5 years after the birth of the pulsar.

Population of Gamma-ray Binaries
The lack of discovery of additional Milky Way gammaray binaries since 1FGL J1018.6-5856, but 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 refer-ences 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 (∼ 10 −3 year −1 ) and the expected lifetime of gamma-ray emission (∼ 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.

CONCLUSION
We have discovered a luminous gamma-ray binary with a 10.3 day period located in an SNR in the LMC using the Fermi LAT. The radio and X-ray counterparts also exhibit flux modulation on this period. The source properties, including radial velocity measurements of the O5 III (f) counterpart, suggest that the system contains a rapidly rotating neutron star. The system may eventually evolve into an X-ray binary. Further multiwavelength observations have the potential to enable a better modeling of the system. Deep radio and X-ray observations should be made to search for the spin period of the neutron star. The discovery of this source suggests that the number of Galactic gamma-ray binaries may have been overestimated, however searches for additional systems should still continue.
We thank C. C. Cheung for support and important comments during the production of this paper. We also thank an anonymous referee for useful comments. Based on observations obtained at the Southern Astrophysical Research (SOAR) telescope, which is a joint project of the Ministério da Ciência, Tecnologia, e Inovação (MCTI) da República Federativa do Brasil, the U.S. National Optical Astronomy Observatory (NOAO), the University of North Carolina at Chapel Hill (UNC), and Michigan State University (MSU). This work was partially supported by NASA Fermi grant NNX15AU83G. The OGLE project has received funding from the National Science Centre, Poland, grant MAESTRO 2014/14/A/ST9/00121 to AU. The Australia Telescope Compact Array is part of the Australia Telescope National Facility which is funded by the Australian Government for operation as a National Facility managed by CSIRO. JBC was supported by an appointment to the NASA Postdoctoral Program at the Goddard Space Flight Center administered by Universities Space Research Association through a contract with NASA. J. Strader acknowledges support from the Packard Foundation. We thank the Swift team for undertaking observations. The Fermi LAT Collaboration acknowledges generous ongoing support from a num-ber of agencies and institutes that have supported both the development and the operation of the LAT as well as scientific data analysis. These include the National Aeronautics and Space Administration and the Department of Energy in the United States, the Commissariat a l'Energie Atomique and the Centre National de la Recherche Scientifique / Institut National de Physique Nucléaire et de Physique des Particules in France, the Agenzia Spaziale Italiana and the Istituto Nazionale di Fisica Nucleare in Italy, the Ministry of Education, Culture, Sports, Science and Technology (MEXT), High Energy Accelerator Research Organization (KEK) and Japan Aerospace Exploration Agency (JAXA) in Japan, and the K. A. Wallenberg Foundation, the Swedish Research Council and the Swedish National Space Board in Sweden. Additional support for science analysis during the operations phase is gratefully acknowledged from the Istituto Nazionale di Astrofisica in Italy and the Centre National d'Études Spatiales in France.   Note. -The stated errors combine the statistical error, determined from RMS values in the region surrounding CXOU J053600.0-673507, and a systematic error conservatively taken to be 5% in the flux density scale between epochs. Note. -Times are Barycentric Julian Dates (TDB) -2400000.5