<figure><img src="image_1.png" /></figure>

Discovery of the Galactic High-mass Gamma-Ray Binary 4FGL J1405.1–6119

R. H. D. Corbet<sup>1,2</sup>, L. Chomiuk<sup>3</sup>, M. J. Coe<sup>4</sup>, J. B. Coley<sup>5,6</sup>, G. Dubus<sup>7</sup>, P. G. Edwards<sup>8</sup>, P. Martin<sup>9</sup>, V. A. McBride<sup>10</sup>, J. Stevens<sup>8</sup>, J. Strader<sup>3</sup>, and L. J. Townsend<sup>11</sup>

<sup>1</sup> University of Maryland, Baltimore County, and X-ray Astrophysics Laboratory, Code 662 NASA Goddard Space Flight Center, Greenbelt Rd., MD 20771, USA

<sup>2</sup> Maryland Institute College of Art, 1300 W Mt. Royal Avenue, Baltimore, MD 21217, USA

<sup>3</sup> Department of Physics and Astronomy, Michigan State University, East Lansing, MI 48824, USA

<sup>4</sup> School of Physics and Astronomy, University of Southampton, Southampton SO17 1BJ, UK

<sup>5</sup> Department of Physics and Astronomy, Howard University, Washington, DC 20059, USA

<sup>6</sup> CREŞŞT/Code 661 Astroparticle Physics Laboratory, NASA Goddard Space Flight Center, Greenbelt Road, MD 20771, USA

<sup>7</sup> Institut de Planétologie et d'Astrophysique de Grenoble, Univ. Grenoble Alpes, CNRS, F-38000 Grenoble, France

<sup>8</sup> Commonwealth Scientific and Industrial Research Organisation Astronomy and Space Science, P.O. Box 76, Epping, New South Wales 1710, Australia

<sup>9</sup> Institut de Recherche en Astrophysique et Planétologie, Université de Toulouse, CNRS, F-31028 Toulouse cedex 4, France

<sup>10</sup> South African Astronomical Observatory, P.O. Box 9, Observatory, 7935, South Africa

<sup>11</sup> Department of Astronomy, University of Cape Town, Private Bag X3, Rondebosch, 7701, South Africa

Received 2019 June 14; revised 2019 August 22; accepted 2019 August 22; published 2019 October 14

Abstract

We report the identification from multiwavelength observations of the *Fermi* Large Area Telescope (LAT) source 4FGL J1405.1–6119 (=3FGL J1405.4–6119) as a high-mass gamma-ray binary. Observations with the LAT show that gamma-ray emission from the system is modulated at a period of 13.7135 ± 0.0019 days, with the presence of two maxima per orbit with different spectral properties. X-ray observations using the *Neil Gehrels Swift Observatory* X-ray Telescope show that X-ray emission is also modulated at this period, but with a single maximum that is closer to the secondary lower-energy gamma-ray maximum. A radio source, coincident with the X-ray source is also found from Australia Telescope Compact Array observations, and the radio emission is modulated on the gamma-ray period with similar phasing to the X-ray emission. A large degree of interstellar obscuration severely hampers optical observations, but a near-infrared counterpart is found. Near-infrared spectroscopy indicates an O6 III spectral classification. This is the third gamma-ray binary to be discovered with the *Fermi*-LAT from periodic modulation of the gamma-ray emission, the other two sources also have early O star, rather than Be star, counterparts. We consider at what distances we can detect such modulated gamma-ray emission with the LAT, and examine constraints on the gamma-ray binary population of the Milky Way.

*Unified Astronomy Thesaurus concepts: Gamma-ray sources (633); X-ray binary stars (1811); High mass X-ray binary stars (733); Neutron stars (1108)*

1. Introduction

High-mass gamma-ray binaries (HMGBs) are very rare objects. They consist of an OB star in orbit with a compact object, where interactions between the two objects cause emission whose spectral energy distribution (SED) peaks above 1 MeV (e.g., Dubus 2015; Dubus et al. 2017, and references therein). The principal emission mechanism is thought to be interactions between the wind of a rapidly rotating neutron star and the wind from the OB companion. However, while the presence of a rapidly rotating neutron star is suspected, it has only been directly detected in the PSR B1259–63 and PSR J2032+4127 binary systems where the gamma-ray binary phase is confined to orbital phases near periastron passage of these highly eccentric systems with long orbital periods (e.g., Li et al. 2017; Abeysekara et al. 2018; Johnson et al. 2018, and references therein). Gamma-ray binaries are expected to be precursors to high-mass X-ray binaries, which they can evolve into after the neutron star has spun down sufficiently to allow accretion to occur. Systems that contain jet-producing accreting black holes may also be gamma-ray sources. However, in this case the SED will peak at X-ray energies and these systems can be classified as “gamma-ray emitting X-ray binaries.” Such systems would include Cygnus X-3 (e.g., Corbel et al. 2012; Zdziarski et al. 2018) and Cygnus X-1 (e.g., Bodaghee et al. 2013; Zanin et al. 2016; Zdziarski et al. 2017).

Generally, HMGBs may be divided into those systems that contain a Be star, and those that contain an O star. Be star systems differ from the O star systems in that Be stars posses circumstellar disks, that may be transient. The known O-star HMGBs are LS 5039, 1FGL J1018.6–5856, and LMC P3, while the known Be-star gamma-ray binaries are LS I+61° 303, HESS J0632+057, PSR B1259–63, and PSR J2032+4127 (e.g., Dubus et al. 2017, and references therein). Note that HESS J0632+057, while detectable at ~TeV energies, appears only weakly detected at ~GeV energies (Malyshev & Chernyakova 2016; Li et al. 2017), and for PSR J2032+4127, while the pulsations away from periastron were discovered at ~GeV energies with the Large Area Telescope (LAT; Abdo et al. 2009), binary related activity was primarily detected at TeV rather than GeV energies (Abeysekara et al. 2018).

The number of detectable Galactic gamma-ray binaries will depend on factors that include the distribution of initial neutron star spin periods, and the lifetimes of gamma-ray binaries. Initial estimates of the possible population of gamma-ray binaries came from Meurs & van den Heuvel (1989) who predicted ~30 systems in the appropriate evolutionary stage. More recently, based on light-curve modeling and the currently known number of systems, Dubus et al. (2017) estimated that ~100 Galactic systems may be detectable. Due to the large point-spread function (PSF) of gamma-ray telescopes such as the *Fermi*-LAT (e.g., The Fermi-LAT Collaboration 2019),
there are very many known gamma-ray sources whose nature is not yet known, and so are potentially binary systems.

Our program to discover gamma-ray binaries initially searches for signs of periodic modulation in gamma-ray light curves obtained with the Fermi-LAT. We then search for counterparts at other wavelengths, and for confirmation of modulation on the same period in the counterparts. From this program we previously identified 1FGL J1018.6-5856 and LMC P3 as high-mass gamma-ray binaries (Corbet et al. 2011, 2016; The Fermi-LAT Collaboration 2012). Here we present the discovery of 4FGL J1405.1-6119 as an additional gamma-ray binary.

We present previous observations and analysis of 4FGL J1405.1-6119 in Section 2. We describe LAT observations and our program to search for modulated gamma-ray sources in Section 3.1. The Swift X-ray observations of the region are presented in Section 3.2, and the Australia Telescope Compact Array (ATCA) observations in Section 3.3. Near-infrared spectroscopy is described in Section 3.4. The discovery of periodic gamma-ray emission from the direction of 4FGL J1405.1-6119 is presented in Section 4.1, the identification of the counterpart as CXOGSG J140514.4-611827 from the detection of modulated X-ray and radio emission is given in Sections 4.2 and 4.3 respectively. The nature of 4FGL J1405.1-6119 and the implications for the overall population of gamma-ray binaries in the Galaxy are discussed in Section 5 with an overall conclusion in Section 6. Unless otherwise stated, uncertainties are given at the 1σ level.

## 2. Previous Observations of 4FGL J1405.1-6119

4FGL J1405.1-6119 is in the fourth LAT catalog (The Fermi-LAT Collaboration 2019) and counterparts were also present in previous LAT catalogs as 1FGL J1405.1-6123c (Abdo et al. 2010), 2FGL J1405.5-6121 (Nolan et al. 2012), and 3FGL J1405.4-6119 (Acero et al. 2015), and was present in the LAT eight year source list as FL8Y J1405.3-6119. In the third catalog of Fermi sources detected in the 10 GeV-2 TeV energy range, it is identified as 3FHL J1405.1-6118 (Ajello et al. 2017).

Lee et al. (2012) examined sources in the second Fermi-LAT catalog to search for pulsars and identified 2FGL J1405.5-6121 as a candidate based on variability and spectral criteria. Saz Parkinson et al. (2016) also undertook a classification of sources in the 3FGL catalog into pulsars and active galactic nuclei, which are the two main categories of identified LAT sources. From this analysis they found that 3FGL J1405.4-6119 was most likely to be a pulsar. Saz Parkinson et al. (2016) also noted the presence of X-ray Telescope (XRT) and Chandra sources in the Fermi error region, with the Chandra source located at R.A. = 14h05m14s47.47 decl. = -61°18'27"7 (±0"19 statistical, ±0"8 systematic). Saz Parkinson et al. (2016) also noted a possible coincidence with the supernova remnant (SNR) G311.5-00.3 (there listed as “G311.5+0.3”). 3FGL J1405.4-6119 was also listed as a strong pulsar candidate by Wu et al. (2018). Clark et al. (2017) included 3FGL J1405.4-6119 in a search for pulsations from LAT pulsar candidates and obtained a 95% upper limit on the pulsed fraction of 0.58.

Using IBIS on board the International Gamma-Ray Astro-physics Laboratory (INTEGRAL) covering 2-200 keV, Landi et al. (2017) reported activity from IGR J14059-6116 between 2003 December 7 to 2009 December 14 (MJD 52,980.45-

55,179.04) and proposed that this IBIS source was the counterpart of 3FGL J1405.4-6119. Landi et al. (2017) also reported on XRT observations which showed the detection of one source within the IBIS error region at R.A. = 14h05m13s93 decl. = -61°18'29"62 with a 5"2 (90% confidence) error radius. The XRT source was also positionally coincident with the near-infrared source 2MASS J14051441-6118282/ allWISE J140514.40-611827.7. Based on the infrared colors of this source, Landi et al. (2017) considered it unlikely that 3FGL J1405.4-6119 was a blazar.

## 3. Observations and Analysis

### 3.1. Gamma-Ray Observations and Analysis

The Fermi-LAT (Atwood et al. 2009) is a pair conversion telescope sensitive to gamma-ray photons with energies between ~20 MeV to >300 GeV. The LAT data used in this paper were obtained between 2008 August 5 and 2019 March 13 (MJD 54,683 to 58,555). During this time, Fermi was primarily operated in a sky survey mode where the LAT pointing position is alternately rocked away from the zenith to the orbit north for one spacecraft orbit, then toward the orbit south for one orbit. In this way, the entire sky is observed every two spacecraft orbits, approximately every three hours. For LAT analysis we used the *fermitools* version 1.0.1. We used the updated Pass 8 LAT data files ("P8R3," Bruel et al. 2018) and the weekly photon files provided by the Fermi Science Support Center which include precomputed diffuse response columns.

In our continuing search for new gamma-ray binaries, we create light curves for all sources in Fermi-LAT catalogs and then calculate power spectra of these to investigate the presence of periodic modulation. The third LAT catalog contained 3033 sources (Acero et al. 2015). The most recent LAT catalog is 4FGL (The Fermi-LAT Collaboration 2019) which contains 5099 sources. In addition, the eight year source list ("FL8Y"<fnref n="12" />) which was a precursor to the fourth LAT catalog contains 5523 sources. In our search for gamma-ray binaries, we analyzed all 3FGL sources and FL8Y sources. At the time of writing a search of all the sources in the 4FGL catalog was still in progress. However, after identifying 3FGL J1405.4-6119/ FL8Y 1405.3-6119 as a binary (Corbet et al. 2019), we then used the parameters of this source from the 4FGL catalog for a refined analysis of 4FGL J1405.1-6119.

Light curves covering an energy range of 100 MeV to 500 GeV were created using a variant of aperture photometry where, instead of simply summing the number of photons within an aperture, we estimate the probability that each photon comes from a source of interest and sum these probabilities (e.g., Kerr 2011; The Fermi-LAT Collaboration et al. 2012). To estimate the probability of a photon coming from a source, models were created for each source using the 3FGL catalog and FL8Y source list using sources within a 10° radius and make3FGLxml and makeFL8Yxml respectively. Photon probabilities were calculated using gtsrcprob and then summed for a 3° radius aperture centered on each source. We note that although in general the use of “probability photometry” increases the signal-to-noise of the light curves, it affects the photometric properties as probabilities are based on a constant source brightness. Thus, when a source is brighter

<footnote n="12">¹² https://fermi.gsfc.nasa.gov/ssc/data/access/lat/fl8y/</footnote>
than the model predicts, the probability of a photon coming from the source is underestimated, and, when the source is fainter than the model, the probability is overestimated. This results in a decrease of the apparent modulated amplitude. In addition, there is an energy-dependent effect as the smaller PSF of the LAT at higher energies results in higher-energy photons having higher weights. Time bins of 500 s were used for all sources.

Power spectra of these probability-weighted 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 beneficial because of the substantial exposure changes from time bin to time bin (The Fermi-LAT Collaboration et al. 2009). For each source the calculated power spectrum covered a period range from 0.05 days (1.2 hr) to the length of the light curve, i.e., ~3510 days, giving ~70,200 independent frequencies. The power spectra were oversampled by a factor of 5 compared to the nominal resolution, which we take to be the inverse of the length of the light curve (e.g., VanderPlas 2018, and references therein), i.e., ~1/3510 day⁻¹. 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. In addition, possible statistical problems with the FAP have been noted (e.g., Koen 1990; Baluev 2008; Süveges 2014). In our photometric analyses the background is not fitted for each time bin, and artifact signals can be seen at several periods including Fermi's ~90 minute 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.<fnref n="13" /> In addition, because of the broad PSF of the LAT, particularly at lower energies, variability in nearby sources can also cause apparent modulation in a light curve. Because of these potential artifacts, when apparent evidence of periodic modulation is found in a LAT light curve, it is highly desirable to be able to confirm the modulation using observations at other wavelengths.

## 3.2. X-Ray Observations and Analysis

The Swift-XRT (Burrows et al. 2005) is a Wolter I X-ray imaging telescope sensitive to X-rays ranging from 0.3 to 10 keV. Although some previous XRT observations of the region had previously been made dating back to 2011 October, these were infrequent, had short observation durations, and had been made over a long interval. We therefore requested additional Swift-XRT Target of Opportunity (TOO) observations to cover two orbital cycles more intensively and with longer durations. The XRT TOO observations of 4FGL J1405.1–6119 (CXOGSG J140514.4–611827) took place from 2018 May 16 to June 17 (MJD 58,254–MJD 58,286) with exposures ranging from ~3.2 to ~4.0 ks. For completeness, we additionally analyzed the 15 short archival XRT observations, which were performed between 2011 October 4 and 2017 December 13 (MJD 55,838 and MJD 58,100). The exposures of these archival XRT observations ranged from ~90 s to ~4.6 ks.

4FGL J1405.1–6119 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). We reduced and analyzed the data using the Swift-XRT product generator (Evans et al. 2007) and the HEAsoft v.6.20 package and calibration files dated 2017 May 1, following the procedures defined in the XRT Data Reduction Guide (Capalbi et al. 2005). The data were reprocessed with the XRTDAS standard data pipeline package xrtpipeline using the standard filtering procedure to apply the newest calibration and default screening criteria.

We find the background-subtracted count rates of our TOO observations to be between 1.0±1.8 $_{-0.9}^{+1.9}$ × 10⁻³ and 7 ± 2 × 10⁻³ counts s⁻¹. Since our observations of 4FGL J1405.1–6119 were not affected by pile-up, we extracted the source spectra from count-dependent circular regions generated by the Swift-XRT product generator. The ancillary response files, accounting for vignetting, PSF correction, and different extraction regions, were generated and corrected for exposure using the FTOOL packages xrtmkarf and xrtexpomap, respectively. Due to the large neutral hydrogen column density, the source was not detected at energies below 2 keV. We, therefore, restrict our spectral analysis to energies above 2 keV.

Individual spectra were not useful for analysis, as each spectrum was found to have between 0 and 21 counts. A cumulative spectrum was therefore extracted, which has a total of ~165 counts, and the total exposure is ~34.9 ks. We further processed the spectral data produced by xselect using the FTOOL grppha, which defined the binning and quality flags of the spectra. We used the quality flag to further eliminate bad data. 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 5 counts per bin.

## 3.3. Radio Observations and Analysis

Radio observations were obtained using the ATCA (Wilson et al. 2011). Dedicated follow-up observations were made between 2017 November 17 and 2018 December 8 (MJD 58,074 to 58,460, 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, with the more compact arrays somewhat more sensitive to the bright extended emission in the vicinity. Details of the array configurations are given in Table 2.

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, observations were made covering the nominal position of the 3FGL source, and the positions of radio sources cataloged by Schinzel et al. (2017) and candidate X-ray sources. Observations differed in length, hour-angle coverage, angular resolution, and sensitivity to extended radio emission in the vicinity, resulting in a heterogeneous data set. Schinzel et al. (2017), who observed with the ATCA in the relatively compact H168 and H214 array configurations, detected four sources within the 3σ error region around 4FGL J1405.1–6119, all of which were brighter than 15 mJy, and are resolved out on longer ATCA baselines. Schinzel et al. (2017)

<footnote n="13">¹³ http://fermi.gsfc.nasa.gov/ssc/data/analysis/LAT_caveats_temporal.html</footnote>
also list a further 60 objects in the field detected outside the $3\sigma$ region, which appear to be associated with the H II region cataloged by Caswell & Haynes (1987) at (l, b) = (311.627, +0.27). Observations of the candidate X-ray sources in more extended ATCA array configurations revealed the fainter (~2 mJy), variable radio counterpart. The final observations of the series were conducted as targeted observations, with the phase calibrator PKS 1420–679 used throughout.

### 3.4. Near-infrared Observations and Analysis

We observed the candidate near-IR counterpart to 4FGL J1405.1–6119 (2MASS J14051441–6118282, allWISE J140514.40–611827.7; Landi et al. 2017) with FLAMINGOS-2 (Eikenberry et al. 2004) on Gemini South on 2018 June 7 through Program ID GS-2018A-Q-412. We used the R3K grism and a 4 pixel slit, giving a resolution of $R \sim 1800$ across the K-band. The target was observed with 10 separate 120 s exposures in a standard ABBA nod, giving 1200 s of exposure time total. We also observed the nearby A0V star HD 119942 as a telluric standard immediately before our science observations.

The data were reduced using standard near-IR procedures, with optimal extraction of the spectra from subtracted, flat-fielded nodded pairs. The final combined spectrum has a signal-to-noise of about 140 per resolution element in the continuum. After interpolating over the intrinsic $Br\gamma$ feature in the telluric standard, we corrected the 4FGL J1405.1–6119 spectrum for telluric absorption by scaling the telluric spectrum to minimize the residuals in the regions of strong telluric lines.

## 4. Results

### 4.1. Gamma-Ray Results

The power spectra of the light curves of all 3033 3FGL and 5523 FL8Y sources were examined for indications of periodic modulation that could arise from previously unknown binary systems. For candidate new binaries our usual threshold for further investigation has been for a source to have a peak power $\ge 18 \times$ mean power level (FAP $< 5 \times 10^{-4}$) and for the period not to coincide with a known artifact. We also exclude very long periods that are suggestive of arising from red-noise type behavior, such as from active galactic nuclei.

From our power spectra, initially from 3FGL sources, and persisting into the FL8Y list and the 4FGL catalog (The Fermi-LAT Collaboration 2019), we noted two peaks in the power spectra of 3FGL J1405.4–6119 and FL8Y J1405.3–6119 which, although lower than our usual threshold for further investigation, were consistent at a $1.7\sigma$ level with being harmonics. From the power spectrum of 4FGL J1405.1–6119 (Figure 1, bottom panel) the two peaks are at $6.8586 \pm 0.0013$ days and $13.7235 \pm 0.0046$ days. The heights of the peaks are ~15.1 and ~17.4 compared to the mean power level, respectively, with FAP values of 0.02 and 0.002, for 77436 trials. For a two value trial, i.e., performing a test on whether there is modulation at either twice or half the period of stronger peak, then the peak height of the 6.86 day period implies an FAP of $6 \times 10^{-7}$.

Folding the LAT light curve on the 13.7 day period (Figure 1, bottom panel inset) shows a profile with a single sharp peak, which occurs near MJD 56,498.7. The presence of two peaks in the power spectrum can be ascribed to the sharpness of this peak which deviates from a simple sine

<figure><img src="image_2.png" /></figure>

**Figure 1.** Power spectra of LAT light curves of 4FGL J1405.1–6119. The bottom panel shows the power spectrum of the probability-weighted LAT light curve ($E > 100$ MeV) while the upper panel shows the power spectrum from conventional aperture photometry without probability weighting ($E > 200$ MeV). In both panels the vertical dashed lines indicate the harmonically related peaks at $13.7135 \pm 0.0019$ (red) and $6.85675 \pm 0.00096$ (blue) days. The peak at 1 day in the lower panel is an artifact commonly seen in LAT power spectra. The dotted-dashed horizontal gray lines show the white-noise significance levels. The insets show the light curves folded on a $13.7135$ day period. For clarity, two cycles are shown, and the count rates have been normalized by dividing by the mean rate.

modulation. Since the probability-weighted light curve can potentially suffer from photometric distortions (Section 3.1) we also created a light curve using conventional aperture photometry without the weighting. For this we employed a one degree radius aperture. From the power spectrum of this (Figure 1, top panel) we find instead a single peak at $6.85675 \pm 0.00096$ days and only a very small peak near $13.7$ days. The height of the peak for the power spectrum of the unweighted $>200$ MeV light curve is $26.5$ with an associated FAP of $3 \times 10^{-7}$, allowing only for the number of independent frequencies. Through the remainder of this paper we use a period of $13.7135 \pm 0.0019$ days, as the period implied by the harmonic in the power spectrum of the unweighted light curve. Period search techniques based on light-curve folding such as $\chi^2$ maximization (e.g., Leahy et al. 1983) or phase-dispersion minimization (Stellingwerf 1978) can naturally also produce signals at “subharmonics” of the intrinsic modulation period, i.e., multiples of the intrinsic modulation period. However, this effect does not occur in direct Fourier-based analyses which quantify the sine-wave components of the modulation. For this reason, we identify the longer, ~13.7 day period as the intrinsic modulation period of 4FGL J1405.1–6119. This is confirmed
<figure><img src="image_3.png" /><figcaption>**Figure 2.** Conventional aperture photometry of LAT observations of 4FGL J1405.1–6119 folded on the 13.7 day period separated by energy. The folded light curves are not background subtracted. Phase zero corresponds to MJD 56,498.7.</figcaption></figure>

by the single peak exhibited in the X-ray and radio light curves when folded on the longer period (Sections 4.2 and 4.3 respectively). The stronger peak near 13.7 days in the power spectrum of the conventional photometry light curve compared to the probability-weighted light curve is because the modulation in the conventional light curve is more nearly sinusoidal on this period. For the probability-weighted light curve, the modulation profile is more sharply peaked (i.e., less sinusoidal), and thus the power is spread over more than one Fourier component, principally the fundamental and the first harmonic. Such a change in profile shape between unweighted and weighted light curves is not necessarily expected to be a general feature of modulation in HMGBs. We note that if only the power spectrum of the conventional aperture photometry LAT light curve was available (Figure 1, bottom panel) then the shorter harmonic period might have been incorrectly determined to be the period of the system.

The unweighted light curve folded on 13.7 days (Figure 1, top panel) now shows, in addition to the primary peak at phase 0, an additional prominent peak near phase 0.5. We then examined the modulation as a function of energy (Figure 2). We find that below 200 MeV no modulation is detected. Above 200 MeV we find that there is a difference in the two peaks, where one is more prominent at lower energies, while the other becomes more prominent at higher energies, and coincides with the phasing of the single peak in the probability-weighted light curve. Thus, we ascribe the difference between the power spectra for the weighted and unweighted light curves as an energy-dependent effect, where the high-energy photons are effectively over-weighted because of the smaller PSF.

<figure><img src="image_4.png" /><figcaption>**Figure 3.** Gamma-ray flux of 4FGL J1405.1–6119 obtained from a phase-resolved likelihood analysis of the LAT data for energies between 200 MeV and 500 GeV. Spectral parameters were frozen to the values from the 4FGL catalog.</figcaption></figure>

To investigate the stability of the period, we divided the conventional aperture photometry light curve into three equal-length sections and calculated the power spectrum of each section separately. The periods of the orbital harmonic derived in this way (6.855 ± 0.006, 6.860 ± 0.004, and 6.853 ± 0.005 days) are all consistent with each other, and also the more precise period derived from the entire light curve.

Neither the power spectrum of the probability-weighted light curve nor that of the unweighted light curve (Figure 1) show any evidence of additional periodicities, or long-period/low-frequency variability, beyond the two peaks related to the ~13.7 day period, and an artifact at one day that is often seen in our power spectra of LAT light curves.

The spectral model derived for 4FGL J1405.1–6119 in the 4FGL catalog is a log normal function (LogParabola), as it was in the 3FGL catalog. For example,

$$ \frac{dN}{dE} = K \left( \frac{E}{E_0} \right)^{-\alpha - \beta \log_e (E/E_0)} \quad (1) $$

This model is used in the LAT catalogs for all sources with significantly curved spectra. Additional information on the 3FGL and 4FGL catalog results, including plots of spectral fits, are available from the Fermi Science Support Center.<fnref n="14" />,<fnref n="15" /> In Table 3 we show the spectral parameters in the 4FGL catalog of 4FGL J1405.1–6119 together with the other persistent HMGBs. In all cases the model employed in the catalog was LogParabola, and we note that the spectral parameters of all sources are broadly comparable.

We made a phase-resolved maximum likelihood analysis to also examine the modulation on 13.7 days. For this we generate a model based on the 4FGL source catalog, but held all parameters fixed, apart from the flux of 4FGL J1405.1–6119. We employed a 10° radius “region of interest.” The phase-resolved fluxes derived from this likelihood analysis are plotted in Figure 3. We find that the flux increase near phase zero is clearly detected, but only a modest increase near a phase of

<footnote n="14">¹⁴ https://fermi.gsfc.nasa.gov/ssc/data/access/lat/4yr catalogue/</footnote>

<footnote n="15">¹⁵ https://fermi.gsfc.nasa.gov/ssc/data/access/lat/8yr catalogue/</footnote>
<figure><img src="image_5.png" /><figcaption>**Figure 4.** Swift-XRT flux of the counterpart of 4FGL J1405.1−6119 folded on the 13.7 day period. The archival observations are shown as black circles, while our Target of Opportunity observations are shown as red diamonds. The green arrows are upper limits ($3\sigma$) from the archival observations.</figcaption></figure>

0.5 is observed. The overall folded light curve also appears somewhat noisier with a lower amplitude than the folded aperture photometry light curve. We investigated phase-resolved likelihood analysis for energy ranges between 200–1000, and 1000–500,000 MeV. However, the fits to the lower-energy range resulted in low Test Score (TS) numbers and so could not be used to investigate the secondary peak. We also made fits with the spectral parameters of FL8Y J1405.3−6119 left free, but these were found to result in spectral parameters which varied in an implausible way between phase bins.

While the likelihood analysis does not strongly show the double-peaked orbital profile shown by the conventional aperture photometry, we note that there is no known artifact that would produce double-peaked modulation in the aperture photometry. In addition, the aperture photometry has the benefit that it is model independent. For these reasons, along with the modest secondary peak that is seen in the likelihood analysis, we believe that the double-peaked structure more strongly found in the aperture photometry is indeed showing the intrinsic behavior of the source, while the secondary apparently softer peak is more difficult to fit due to the higher background at lower energies.

## 4.2. X-Ray Results

### 4.2.1. X-Ray Flux Variations

The XRT TOO observations of 4FGL J1405.1−6119 cover more than two orbital periods. The XRT light curve folded on the 13.7 days period is shown in Figure 4. From this, a strong, approximately sinusoidal, modulation on the 13.7 days gamma-ray period can be seen. However, we note that X-ray minimum occurs near the phase of primary gamma-ray maximum. Although the archival observations have shorter durations, and thus larger error bars, than our TOO observations, they also indicate modulation on the 13.7 day period.

To investigate the energy-dependent phase shift between the Fermi gamma-ray and Swift X-ray light curves, we calculated the maximum value of the cross-correlation function between the two folded light curves after applying

phase shifts to the X-ray light curve. For this we only used XRT data from the TOO observations. Due to the sparse coverage of the X-ray light curve, we did not bin this but instead linearly interpolated to cover observation gaps. From this we find the X-ray light curve leads the LAT light curve by a phase of 0.586. We also used a simple sine-wave fit to the entire XRT light curve to determine maximum of the X-ray light curve and find this occurs at a phase of $0.59 \pm 0.03$. This is thus consistent with the results of the cross-correlation analysis. However, this result should be treated with caution because we only have good coverage of two cycles, the individual flux measurements have relatively large errors, and the shape of the modulation is hence not yet well defined.

### 4.2.2. X-Ray Spectrum

The cumulative spectrum of 4FGL J1405.1−6119 were analyzed using `XSPEC v12.9.0k`. We made use of the `XSPEC` convolution model `cflux` to calculate the fluxes and associated errors of 4FGL J1405.1−6119. To fit the cumulative XRT spectrum, we used several models that are used to describe systems that host a neutron star: a power law (power), a power law with a high-energy cutoff (highecut in `XSPEC`), and a cutoff power law (cutoffpl in `XSPEC`). All models were modified by an absorber that fully covers the source using appropriate cross sections (Verner et al. 1996) and abundances (Wilms et al. 2000).

We initially allowed the neutral hydrogen column density and photon index to be free parameters and performed spectral fits on the cumulative spectrum. Due to the short exposure times of the archival observations (see Table 1), we chose to include only the TOO observations that were performed between 2018 May 16 and 2018 June 17 (MJD 58,254 and MJD 58,286). The model that provides the best fit (C statistic of 10.19 for 18 degrees of freedom) to the data is a power law. While a good fit does not require a high-energy cutoff, which is typically found in accreting pulsars (Coburn et al. 2002), our spectra are limited to energies below 10 keV. We note that while high-energy cutoffs are typically seen at higher energies, none were found in LS 5039 (Takahashi et al. 2009) or 1FGL J1018.6−5856 (An et al. 2015) up to 40 keV. We find the neutral hydrogen column density for the fully covered absorption and the photon index to be $1.1_{-0.8}^{+0.9} \times 10^{23} \text{ cm}^{-2}$ and $2 \pm 1$, respectively. We also attempted fitting the cumulative spectrum of the TOO observations along with the archival ones and found larger uncertainties in the model parameters. The cumulative XRT spectrum, and best-fitting power-law model are shown in Figure 5.

Since neither the neutral hydrogen column density for the fully covered absorber nor the photon index could be accurately constrained, we converted the mean derived value of optical reddening $E(B-V)$ of 10.2 (Section 4.4) into the predicted neutral hydrogen column density and chose to fit the cumulative spectrum with the $N_H$ frozen to it. From Equation (1) in Güver & Oumlzel (2009), we calculate the neutral hydrogen column density along the line of sight to be $6.9 \times 10^{22} \text{ cm}^{-2}$. Using the power model, we find the photon index to be $1.5 \pm 0.4$ and the unabsorbed X-ray flux in the 2–8 keV band to be $8 \pm 1 \times 10^{-13} \text{ erg cm}^{-2} \text{ s}^{-1}$. For comparison, Takahashi et al. (2009) found a power-law index for LS 5039 that varied between 1.45 and 1.61, for 1FGL J1018.6 An et al. (2015) found indices between ~1.4 and 1.7, and for LMC
<table><caption>**Table 1**

Swift-XRT Observation Log of CXOU J053600.0–673507(4FGL J1405.1–6119)</caption><thead><tr><th>ObsID</th><th>Start Time (UT)</th><th>End Time (UT)</th><th>Phase<sup>a</sup></th><th>Exposure<sup>b</sup></th><th>Count Rate<sup>c</sup></th><th>Flux<sup>d</sup></th></tr></thead><tbody><tr><td>00041805001</td><td>2011 Oct 4 02:55:00</td><td>2011 Oct 4 23:05:29</td><td>0.831–0.890</td><td>3.1</td><td>4<sup>+2</sup><sub>-1</sub></td><td>9<sup>+4</sup><sub>-3</sub></td></tr><tr><td>00041805002</td><td>2011 Oct 8 20:56:00</td><td>2011 Oct 8 21:48:29</td><td>0.177–0.178</td><td>1.0</td><td>1<sup>+2</sup><sub>-1</sub></td><td>3<sup>+5</sup><sub>-2</sub></td></tr><tr><td>00042313001</td><td>2011 Oct 10 14:52:59</td><td>2011 Oct 10 15:44:59</td><td>0.304–0.305</td><td>0.5</td><td>&lt;4.9</td><td>&lt;11</td></tr><tr><td>00041805003</td><td>2011 Nov 1 22:55:00</td><td>2011 Nov 1 23:49:51</td><td>0.933–0.934</td><td>1.3</td><td>2<sup>+1</sup><sub>-1</sub></td><td>4<sup>+4</sup><sub>-3</sub></td></tr><tr><td>00042320001</td><td>2012 Mar 10 10:58:00</td><td>2012 Mar 10 11:49:44</td><td>0.376–0.377</td><td>0.6</td><td>&lt;4.5</td><td>&lt;10</td></tr><tr><td>00041805005</td><td>2012 Sep 21 09:30:59</td><td>2012 Sep 22 23:05:25</td><td>0.592–0.704</td><td>4.6</td><td>5 ± 1</td><td>13 ± 3</td></tr><tr><td>00084747001</td><td>2015 Feb 19 18:29:00</td><td>2015 Feb 19 18:39:00</td><td>0.861–0.861</td><td>0.6</td><td>&lt;14.4</td><td>&lt;33</td></tr><tr><td>00084747003</td><td>2016 Dec 4 01:34:00</td><td>2016 Dec 4 01:39:00</td><td>0.500–0.500</td><td>0.3</td><td>8<sup>+14</sup><sub>-7</sub></td><td>20<sup>+33</sup><sub>-16</sub></td></tr><tr><td>00084747004</td><td>2016 Dec 7 02:53:00</td><td>2016 Dec 7 20:41:00</td><td>0.722–0.777</td><td>1.1</td><td>3<sup>+3</sup><sub>-2</sub></td><td>7<sup>+7</sup><sub>-4</sub></td></tr><tr><td>00084747005</td><td>2016 Dec 15 03:54:00</td><td>2016 Dec 15 03:59:00</td><td>0.309–0.309</td><td>0.2</td><td>8<sup>+13</sup><sub>-6</sub></td><td>18<sup>+29</sup><sub>-14</sub></td></tr><tr><td>00084747006</td><td>2017 Sep 22 02:07:00</td><td>2017 Sep 22 02:25:00</td><td>0.794–0.795</td><td>1.0</td><td>4<sup>+3</sup><sub>-2</sub></td><td>10<sup>+8</sup><sub>-5</sub></td></tr><tr><td>00084747007</td><td>2017 Dec 7 05:30:00</td><td>2017 Dec 7 05:47:00</td><td>0.346–0.347</td><td>0.9</td><td>10<sup>+4</sup><sub>-3</sub></td><td>22<sup>+10</sup><sub>-8</sub></td></tr><tr><td>00084747008</td><td>2017 Dec 9 00:35:00</td><td>2017 Dec 9 00:48:00</td><td>0.477–0.478</td><td>0.7</td><td>10<sup>+5</sup><sub>-4</sub></td><td>24<sup>+10</sup><sub>-9</sub></td></tr><tr><td>00084747009</td><td>2017 Dec 13 17:50:00</td><td>2017 Dec 13 17:55:00</td><td>0.821–0.822</td><td>0.2</td><td>&lt;11.5</td><td>&lt;27</td></tr><tr><td>00084747010</td><td>2018 May 16 10:02:00</td><td>2018 May 16 13:35:00</td><td>0.027–0.038</td><td>3.9</td><td>1.0<sup>+1.8</sup><sub>-0.9</sub></td><td>2<sup>+4</sup><sub>-2</sub></td></tr><tr><td>00084747011</td><td>2018 May 19 08:09:00</td><td>2018 May 19 11:48:00</td><td>0.241–0.251</td><td>3.9</td><td>1.5<sup>+1.0</sup><sub>-0.7</sub></td><td>4 ± 2</td></tr><tr><td>00084747012</td><td>2018 May 22 11:01:00</td><td>2018 May 22 14:46:00</td><td>0.468–0.479</td><td>3.9</td><td>6 ± 2</td><td>14<sup>+5</sup><sub>-4</sub></td></tr><tr><td>00084747013</td><td>2018 May 25 07:43:00</td><td>2018 May 25 11:05:00</td><td>0.677–0.687</td><td>3.7</td><td>5 ± 1</td><td>12 ± 3</td></tr><tr><td>00084747014</td><td>2018 May 28 04:14:00</td><td>2018 May 28 07:36:00</td><td>0.885–0.895</td><td>3.8</td><td>3 ± 1</td><td>8<sup>+3</sup><sub>-2</sub></td></tr><tr><td>00084747015</td><td>2018 May 31 08:51:00</td><td>2018 May 31 18:32:00</td><td>0.118–0.147</td><td>3.2</td><td>1.2<sup>+1.4</sup><sub>-0.9</sub></td><td>3<sup>+3</sup><sub>-2</sub></td></tr><tr><td>00084747016</td><td>2018 Jun 3 03:30:00</td><td>2018 Jun 3 06:53:00</td><td>0.320–0.330</td><td>3.9</td><td>4<sup>+2</sup><sub>-1</sub></td><td>10<sup>+4</sup><sub>-3</sub></td></tr><tr><td>00084747017</td><td>2018 Jun 6 00:21:00</td><td>2018 Jun 6 05:09:00</td><td>0.529–0.544</td><td>4.0</td><td>5 ± 1</td><td>12 ± 3</td></tr><tr><td>00084747018</td><td>2018 Jun 9 09:29:00</td><td>2018 Jun 9 20:54:00</td><td>0.776–0.810</td><td>3.9</td><td>5 ± 1</td><td>11 ± 3</td></tr><tr><td>00084747019</td><td>2018 Jun 17 15:25:00</td><td>2018 Jun 17 20:25:00</td><td>0.377–0.392</td><td>3.8</td><td>7 ± 2</td><td>16 ± 4</td></tr></tbody></table>

**Notes.**

<sup>a</sup> Phase zero is defined as the epoch of maximum flux in the Fermi-LAT (MJD 56498.7).

<sup>b</sup> The net exposure time spread over several snapshots. Units are ks.

<sup>c</sup> Count Rate is in the 2–10 keV energy band. Units are 10<sup>-3</sup> counts s<sup>-1</sup>. Errors are at the 1σ level.

<sup>d</sup> Unabsorbed Swift-XRT flux in the 1.0–10.0 keV bandpass converted with PIMMS. Units are 10<sup>-13</sup> erg cm<sup>-2</sup> s<sup>-1</sup>.

<sup>e</sup> 3σ upper limits.

<figure><img src="image_6.png" /><figcaption>**Figure 5.** Cumulative Swift-XRT spectrum of 4FGL J1405.1–6119 (=CXOGSG J140514.4–611827). The best-fit power-law model is shown as a histogram.</figcaption></figure>

P3 indices between ~1.3 and 1.6 have been reported (Bamba et al. 2006; Seward et al. 2012; Corbet et al. 2016). It is important to note that while this model provides a good fit to the XRT spectrum (C statistic of 10.39 for 19 degrees of freedom), this interpretation should be treated with caution because we assumed that the fully covered $N_H$ is entirely interstellar in origin. While this is the simplest interpretation of

the data, we cannot exclude the possibility that the fully covered $N_H$ is due to a combination of intrinsic and interstellar absorbers. However, for LS 5039 and 1FGL J1018.6–5856 the measured $N_H$ values are consistent with being only due to interstellar absorption (An et al. 2015; Takahashi et al. 2009).

### 4.3. Radio Results

A point radio source was found with ATCA in the LAT error region of the gamma-ray source which coincides with the position of CXOGSG J140514.4–611827 and also the proposed near-IR counterpart at R.A. = 14h05m14s42 ± 0.02, decl. = -61°18′28″33 ± 0.0′03 (J2000). The offset of this position from the 3FGL, FL8Y, and 4FGL positions and corresponding 95% semimajor error radii are 1.43 (1.88), 1.23 (1.31), and 1.50 (1.25)′, respectively.

The radio flux densities folded on the 13.7 day period (Figure 6) show modulation of the emission on the gamma-ray period. There is a gap in the coverage between phases ~0.1 and ~0.3, which hampers our determination of the phase of maximum flux, but it must be either within this gap or near phase ~0.3. We do have coverage at phase 0, the maximum of the gamma-ray flux, and find a minimum radio flux. Thus, similar to the X-ray flux, the radio flux is offset from the primary gamma-ray maximum. However, there is an indication that the radio flux maximum may occur somewhat before the X-ray maximum. We note a low intensity 9 GHz measurement at a phase of 0.7, while the flux at 5.5 GHz is not exceptionally
<table><caption>**Table 2**
Australia Telescope Compact Array Radio Measurements</caption><thead><tr><th>Time (MJD)</th><th>Half Duration (days)</th><th>Configuration</th><th>Flux Density 5.5 GHz (mJy)</th><th>Error 5.5 GHz (mJy)</th><th>Flux Density 9 GHz (mJy)</th><th>Error 9 GHz (mJy)</th></tr></thead><tbody><tr><td>58074.051677</td><td>0.054630</td><td>1.5C</td><td>1.43</td><td>0.22</td><td>0.59</td><td>0.14</td></tr><tr><td>58227.832580</td><td>0.034491</td><td>H168</td><td>0.35</td><td>0.09</td><td>0.36</td><td>0.09</td></tr><tr><td>58258.472279</td><td>0.143056</td><td>6D</td><td>1.81</td><td>0.06</td><td>1.28</td><td>0.06</td></tr><tr><td>58259.465913</td><td>0.148727</td><td>6D</td><td>1.49</td><td>0.11</td><td>1.25</td><td>0.08</td></tr><tr><td>58266.432869</td><td>0.129572</td><td>6D</td><td>1.36</td><td>0.04</td><td>0.86</td><td>0.03</td></tr><tr><td>58268.246122</td><td>0.064699</td><td>6D</td><td>0.85</td><td>0.22</td><td>0.51</td><td>0.08</td></tr><tr><td>58277.279860</td><td>0.024132</td><td>6D</td><td>1.22</td><td>0.08</td><td>0.34</td><td>0.06</td></tr><tr><td>58278.287325</td><td>0.023843</td><td>6D</td><td>1.32</td><td>0.07</td><td>0.70</td><td>0.06</td></tr><tr><td>58409.903240</td><td>0.052257</td><td>6A</td><td>1.60</td><td>0.17</td><td>1.12</td><td>0.09</td></tr><tr><td>58411.935300</td><td>0.049479</td><td>6A</td><td>1.03</td><td>0.40</td><td>1.15</td><td>0.17</td></tr><tr><td>58427.816723</td><td>0.054514</td><td>6B</td><td>1.04</td><td>0.14</td><td>0.63</td><td>0.12</td></tr><tr><td>58441.821527</td><td>0.003414</td><td>6B</td><td>1.31</td><td>0.26</td><td>0.77</td><td>0.11</td></tr><tr><td>58445.044617</td><td>0.121181</td><td>6B</td><td>0.92</td><td>0.15</td><td>0.57</td><td>0.16</td></tr><tr><td>58446.951851</td><td>0.091377</td><td>6B</td><td>0.82</td><td>0.05</td><td>0.53</td><td>0.05</td></tr><tr><td>58460.747800</td><td>0.117766</td><td>H168</td><td>0.55</td><td>0.07</td><td>0.54</td><td>0.05</td></tr></tbody></table>

**Note.** The stated errors combine the statistical error, determined from rms values in the region surrounding CXOGSG J140514.4–611827, and a systematic error conservatively taken to be 5% in the flux density scale between epochs. The ATCA array configurations are the standard names for the physical locations of the antennas: see https://www.narrabri.atnf.csiro.au/operations/array_configurations/configurations.html? for full details of the antenna spacings in each array configuration.

<figure><img src="image_7.png" /><figcaption>**Figure 6.** ATCA observations of the counterpart of 4FGL J1405.1–6119 folded on the 13.7 days period.</figcaption></figure>

low. The 9 GHz measurement appears to be reliable and not due to an instrumental effect.

### 4.4. Near-infrared Results: Spectral Classification and Distance

The near-infrared spectrum obtained with FLAMINGOS-2 is shown in Figure 7. To classify 4FGL J1405.1–6119, we use the near-IR library of early-type stars from Hanson et al. (2005). These spectra were obtained at higher resolution than our spectrum (R ~ 12000 rather than ~1800), so we smooth them to match our resolution before comparison.

The most prominent line in the spectrum is N III in emission (rest frame 2.1155 µm), and we also observe weak C IV emission at 2.078 µm. Absorption lines clearly present are Brγ (H) at 2.1661 µm and He II at 2.188 µm, and a weaker He I

<figure><img src="image_8.png" /><figcaption>**Figure 7.** Near-infrared spectrum of the counterpart of 4FGL J1405.1–6119 (2MASS J14051441–6118282, allWISE J140514.40–611827.7) obtained with FLAMINGOS-2 on Gemini South. Overplotted as a red dashed line is the spectrum of the spectral comparison star HD 190864.</figcaption></figure>

absorption line is present at 2.113 µm, nearly blended with the N III emission line.

In terms of the lines present and absent and their relative strengths, the best match from the Hanson et al. (2005) catalog is HD 190864, which is classified as an O6.5 III star. The object spectrum overplotted with this star is in Figure 7. By comparison, an O5 III star in the catalog, HD 15558, has much weaker He II than 4FGL J1405.1–6119. Similarly, the O5/O5.5 supergiants in the catalog have weaker He II and stronger C IV emission than 4FGL J1405.1–6119. Considering
<table><caption>**Table 3**
4FGL Spectral Parameters of Persistent High-mass Gamma-Ray Binaries</caption><thead><tr><th>Source</th><th><i>E</i><sub>0</sub> (MeV)</th><th>α (lp Index)</th><th>β</th></tr></thead><tbody><tr><td>4FGL J1405.1–6119</td><td>1910.99</td><td>2.792 ± 0.076</td><td>0.302 ± 0.048</td></tr><tr><td>4FGL J0240.5+6113<sup>a</sup></td><td>1177.76</td><td>2.396 ± 0.008</td><td>0.152 ± 0.005</td></tr><tr><td>4FGL J0535.2-6736<sup>b</sup></td><td>752.28</td><td>2.609 ± 0.090</td><td>0.153 ± 0.061</td></tr><tr><td>4FGL J1018.9-5856<sup>c</sup></td><td>1772.77</td><td>2.678 ± 0.023</td><td>0.236 ± 0.014</td></tr><tr><td>4FGL J1826.2-1450<sup>d</sup></td><td>1031.04</td><td>2.581 ± 0.023</td><td>0.124 ± 0.015</td></tr></tbody></table>

**Notes.** Spectral parameters are taken from the 4FGL LAT source catalog (The Fermi-LAT Collaboration 2019). For all sources the LogParabola model was used, see Section 4.1.

<sup>a</sup> = LS I +61°303 (Be star).

<sup>b</sup> = LMC P3 (O star).

<sup>c</sup> = 1FGL J1018.6–5856 (O star).

<sup>d</sup> = LS 5039 (O star).

<table><caption>**Table 4**
Distance Determinations</caption><thead><tr><th>Color</th><th>R&L<br>Distance (kpc) <i>E</i>(<i>B</i> − <i>V</i>)</th><th>S&F<br>Distance (kpc) <i>E</i>(<i>B</i> − <i>V</i>)</th></tr></thead><tbody><tr><td><i>H</i> − <i>K</i> (2MASS)</td><td>8.1 (8.76)</td><td>6.6 (11.56)</td></tr><tr><td><i>H</i> − <i>K</i> (VVV)</td><td>7.9 (9.12)</td><td>6.4 (12.0)</td></tr><tr><td><i>J</i> − <i>H</i> (VVV)</td><td>8.9 (8.65)</td><td>8.0 (11.0)</td></tr></tbody></table>

**Note.** Distances in kiloparsec determined using photometric measurements from 2MASS (Landi et al. 2017) and VVV (Minniti et al. 2010), the reddening laws of Rieke & Lebofsky (1985, R&L) and Schlafly & Finkbeiner (2011, S&F), and the absolute magnitudes of an O6.5 III star from Martins & Plez (2006).

O dwarfs, the absorption lines generally appear to be broader than observed in 4FGL J1405.1–6119, consistent with their higher rotational velocities and/or gravities.

We conclude that, on the basis of our spectrum, the best classification of 4FGL J1405.1–6119 is as an O6.5 III star, but emphasize that this classification should be revisited with a higher signal-to-noise spectrum in the future.

By cross-correlation between the standards and our spectrum in the region of the He II line, we find a barycentric radial velocity of 207 ± 16 km s⁻¹. We also note that the emission lines appear to lie at a different velocity, so the systemic velocity of this star is not necessarily well-determined.

In order to estimate the distance to 4FGL J1405.1–6119 we used the tabulated absolute magnitudes from Martins & Plez (2006) for an O6.5 III star (i.e., $M_J = -5.03$, $M_H = -4.92$, and $M_K = -4.82$), the 2MASS photometric measurements ($H = 14.369 \pm 0.068$, $K = 12.769 \pm 0.044$ Landi et al. 2017), and the reddening laws determined by Rieke & Lebofsky (1985), for $R_V = 3.1$, and also Schlafly & Finkbeiner (2011). The infrared counterpart is also present in the Vista Variables in the Via Lactea (VVV) catalog (Minniti et al. 2010) with measured values of $J = 17.270 \pm 0.028$, $H = 14.524 \pm 0.005$, $K = 12.854 \pm 0.003$. The distances derived from these photometric measurements and reddening curves are given in Table 4. The mean of these distances is 7.7 kpc, with a standard deviation of 1 kpc. For the remainder of this paper we use 7.7 kpc for estimates of source luminosity, but caution that there is considerable uncertainty on this. The mean derived value of $E(B-V)$ is 10.2 with a standard deviation of 1.5. We note that the reddening curve of Rieke & Lebofsky (1985)

<figure><img src="image_9.png" /></figure>

**Figure 8.** Maximum gamma-ray photon luminosity plotted against distance for HMGBs that show detectable periodic modulation in their LAT light curves. Sources with Be star primaries are plotted as open circles, sources with O star primaries as filled circles. The green diamonds show the approximate maximum distance at which modulation could be detected in each source, see the text for details. The green line is a fit to the maximum detectable distances with a slope fixed equivalent to distance⁻². The use of a different distance for an individual source would cause its location in this diagram to move parallel to the green dashed line.

gives somewhat higher distances and lower reddening than using the reddening curve from Schlafly & Finkbeiner (2011). In addition, the 2MASS and VVV measurements are formally inconsistent. If this is not due to, for example, different photometric bandpasses, it might indicate variability in the near-IR.

# 5. Discussion

## 5.1. Properties and Nature of 4FGL J1405.1–6119

The periodic modulation on 13.7 days found at gamma-ray, X-ray, and radio wavelengths, together with the identification of a stellar counterpart, clearly show that 4FGL J1405.1–6119 is a high-mass gamma-ray binary, and the 13.7 day period is expected to be the orbital period of the system. 4FGL J1405.1–6119 is the third gamma-ray binary, after 1FGL J1018.6–5856 and LMC P3, to be found from the initial discovery of periodic modulation of the LAT light curve. In all three cases, the primary star in the system is an early O star, rather than a Be star. The O6.5 III classification for 4FGL J1405.1–6119 indicates a mass for the primary star ~25–35 $M_\odot$ (Mahy et al. 2015).

For a distance of ~7.7 kpc, the implied maximum gamma-ray luminosity would be comparable to, and possibly higher than, that of LS 5039, while approximately half that of 1FGL J1018.6–5856 and a tenth that of LMC P3. The unabsorbed mean X-ray luminosity of 4FGL J1405.1–6119 is ~5.6 × 10³³ erg s⁻¹ (d/7.7 kpc)². For LMC P3 the unabsorbed luminosity measured with the Swift-XRT (2.0–7.5 keV) was also considerably larger at ~9.6 × 10³⁴ erg s⁻¹ (d/50 kpc)² (Corbet et al. 2016). For 1FGL J1018.6–5856 there are flares that reach unabsorbed luminosities (0.5–10 keV) of ~10³⁴ erg s⁻¹ (d/5.4 kpc)², while between these the luminosity can decline to ~2.6 × 10³³ erg s⁻¹ (d/5.4 kpc)² (An et al. 2013) and so on average it is more similar to 4FGL J1405.1
−6119. For LS 5039, the X-ray luminosity (1–10 keV) is also comparable at ~6 × 10³³ erg s⁻¹ (d/2.5 kpc)² (Bosch-Ramon et al. 2007; Takahashi et al. 2009; Rea et al. 2011).

Because of the similarities with the other systems, we also hypothesize that gamma-ray emission from 4FGL J1405.1−6119 is driven by the interaction between the wind from a rapidly rotating neutron star and the wind from the O star companion. Orbital modulation of observed flux in HMGBs can be driven by orbital phase-dependent changes in both viewing angle and varying source distance for systems with significant eccentricity. The overall pattern of variability in 4FGL J1405.1−6119 is similar to that seen in LMC P3 where the X-ray and radio modulations are approximately in phase with each other, but close to 0.5 out of phase from the primary gamma-ray peak. However, the energy-dependent gamma-ray modulation seen in 4FGL J1405.1−6119 with a secondary softer peak approximately 0.5 out of phase from the primary peak is unlike LMC P3. For 1FGL J1018.6−5856 the X-ray emission has a broad component that is out of phase from the overall gamma-ray modulation, although there is an additional X-ray component that exhibits sharp “flares” that is in phase with the gamma-ray modulation (The Fermi-LAT Collaboration et al. 2012; An et al. 2015). An & Romani (2017) have also reported for 1FGL J1018.6−5856 that the gamma-ray light curve below 200 MeV includes a component that is modulated with a similar phasing to that of the broad X-ray component. Anti-phasing between the modulation of the gamma-ray and X-ray emission in the O-star HMGBs may be explained if the gamma-ray variability is primarily due to inverse Compton scattering, which gives greater flux near superior conjunction, and X-ray modulation is due to Doppler boosting, which will give higher observed flux near inferior conjunction (e.g., Dubus et al. 2015).

To investigate the driving mechanisms behind these modulations in detail it will be important to determine a radial velocity curve for the system. In addition to the orbital variability of the radio flux, there also appears to be shorter-term variability which does not occur simultaneously at both frequencies. The origin of this radio variability is unclear, but it may be similar to what has been seen in LMC P3 (Corbet et al. 2016).

We note that 4FGL J1405.1−6119 is located ~40′ from the previously suggested SNR counterpart G311.5−00.3 (Saz Parkinson et al. 2016), which has a 5′ radio diameter (Green 2014), and so the two sources must be distinct. Thus, for 4FGL J1405.1−6119 we cannot directly estimate the age of the system from an associated SNR. Similarly, for 1FGL J1018.6−5856, Marcote et al. (2018) also find that their motion and distance measurements are not consistent with a previously proposed possible SNR association. This is unlike LMC P3 for which there does appear to be an association with an SNR (Seward et al. 2012; Corbet et al. 2016).

The Be star HMGB LS I +61°303 has shown, in addition to its ~26.5 days orbital period, a superorbital period near 1667 days that was originally detected from radio observations (e.g., Gregory 2002, and references therein). This superorbital modulation has also been seen in the LAT light curve of this source (e.g., Hadasch et al. 2012; Xing et al. 2017; Jaron et al. 2018). Such superorbital variability has not yet been shown to be a general property of HMGBs and, as noted in Section 4.1, we do not see such modulation in 4FGL J1405.1−6119. However, since superorbital modulation has been reported in

a number of high-mass X-ray binaries accreting from both Be star envelopes and the winds from OB star companions (Rajoelimanana et al. 2011; Corbet & Krimm 2013; Corbet et al. 2018), it is possible that superorbital modulation may still be found in other HMGBs beyond LS I +61°303.

## 5.2. Distance Limits on Detectability of Galactic Gamma-Ray Binaries

The Galactic population of gamma-ray binaries has been discussed by Dubus et al. (2017) and these authors found very large uncertainties on the possible population as 101⁻⁸⁹·⁵². We previously suggested that the discovery of a binary in the LMC but the lack of detection, at that time, of additional Galactic systems might suggest we had discovered all detectable sources. Here we investigate this in more detail by estimating the approximate maximum distance at which we could detect the known binaries from their periodic variability in LAT light curves.

In Figure 8, using Table 3 of Dubus et al. (2017) together with the parameters of 4FGL J1405.1−6119, we plot photon luminosity of those sources with detectable modulation against distance. We note that there is a general trend for the most distant known sources to be the most luminous. This is suggestive of a flux detection threshold.

As a detection threshold in our power spectra, we consider that a peak should have a minimum height of 20 times the mean power level to be detectable. We then determined the actual heights of the orbital peaks in the FL8Y LAT light curves using light curves covering from 2008 August 5 and 2018 July 5 (MJD 54,683 to 58,304). The ratio of the maximum detection distance to the actual distance was then taken as (peak height/20)¹/⁴. One square root to convert from power to amplitude, and one square root to allow for inverse square law decrease in flux. We then multiplied the source distances by the resulting factor, and these are plotted as green diamonds in Figure 8. The green line is a straight line fit (in log-log space) with a slope fixed to be equivalent to r⁻². The space to the left of this line is assumed to be detectable, while the space to the right would not be detectable. Note that changed distance measurements would move source locations along lines parallel to the green line, and thus the location of the green line itself would not change. In addition, due to the fourth root dependence, the location of the line does not strongly depend on the choice of detection threshold relative power level.

Although our binary detection limit is derived in a simple way, and does not account for different modulation patterns or differing backgrounds, it does give generally consistent results for the Galactic binaries. We note, however, that LMC P3 does lie inside the region where it would nominally not be detectable. Thus, the detectability of LMC P3 may be due to the considerably lower background due to it being located away from the Galactic plane. The Be star system HESS J0632 +057, which is at most weakly detected with the LAT, and so does not exhibit periodic modulation, is at a distance of 0.8 kpc and a luminosity of 2.9 × 10³⁵ ph s⁻¹ (Dubus et al. 2017) and hence is located to the right of the nominal detection threshold line.

With extended LAT observations, the detection line will slowly move to the right. We find empirically that the strengths of the peaks of the known periods in the power spectra relative to the mean power levels increase approximately linearly with
time (i.e., the mean amplitude noise decreases as the square root of time). Thus, the maximum detection distance would increase as $\propto t^{-4}$. The maximum detection volume, for sources near the plane, would increase as $\propto t^{-2}$. Convolved with this, for the discovery rate of new binaries from their modulation, is the poorly constrained luminosity function for gamma-ray binaries. We note that the Be star systems are relatively nearby, and lower luminosity, perhaps suggestive that there remains a more extensive population of Be stars remaining to be discovered.

Our searches for binaries based on detection of periodic modulation do require us to distinguish between periodic and nonperiodic variability. If we require $\gtrsim 10$ cycles to enable clear discrimination between a coherent signal and red noise, then periods $\lesssim 1$ yr would be required. Thus, discovery of Be star binaries with periods such as those of PSR B1259-63 and PSR J2032+4127 would be difficult or impossible. However, Be X-ray binaries generally exhibit binary periods between ~20 and ~400 days (e.g., Corbet et al. 2017, and references therein). Brief flaring behavior such as that exhibited by PSR B1259-63 (e.g., Johnson et al. 2018) would also be more readily distinguished from red-noise behavior with fewer cycles.

From the modest number of gamma-ray binaries known so far, there are two factors that may facilitate the discovery of O-star systems from the detection of modulated gamma-ray emission. The orbital periods of the O-star systems are shorter than those of the known Be systems. In addition, the Be star systems typically have lower luminosities and so the O-star systems are visible at greater distances.

Our detection threshold distance for sources with similar luminosities and gamma-ray modulation properties to LS 5039 and 4FGL J1405.1-6119 raises the possibility that the HMGB population estimate of Dubus et al. (2017) may have been an underestimate. Dubus et al. (2017) estimated that ~80% of LS 5039-like systems in the Galaxy would be detected. However, discovery from gamma-ray modulation alone within the Galactic plane appears currently restricted to ~7 kpc. Thus, multiwavelength observations of unidentified gamma-ray sources to aid variability searches will continue to be important. For example, the identification of an OB counterpart to an X-ray source within a LAT error ellipse could also facilitate the discovery of a gamma-ray binary.

## 6. Conclusion

We have identified 4FGL J1405.1-6119 as a gamma-ray binary with gamma-ray flux modulated on a period near 13.7 days. This periodicity is also seen at X-ray and radio wavelengths. As with 1FGL J1018.6-5856 and LMC P3, which were also detected from modulated gamma-ray emission with the LAT, 4FGL J1405.1-6119 contains an O, rather than a Be, star primary. The system is heavily obscured in the optical. Future infrared radial velocity measurements would be valuable for determining the system geometry and establishing the system orientation at maximum flux in the different wavebands. A determination of the system eccentricity would also be important. It is suspected that 4FGL J1405.1-6119 contains a rapidly rotating neutron star, although we do not yet have a direct detection of this. The Galactic population of gamma-ray binaries is unclear, but there may remain a lower gamma-ray luminosity population remaining to be discovered,

particularly since for the known systems, the Be star systems have lower luminosities.

We thank the anonymous referee for useful comments. This work was partially supported by NASA Fermi grant NNX15AU83G. 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. J.S. acknowledges support from a Packard Fellowship. This work was partially supported by NASA grant 80NSSC17K0507 and NSF grant AST-1714825. This paper was partially based on observations obtained at the Gemini Observatory, which is operated by the Association of Universities for Research in Astronomy, Inc., under a cooperative agreement with the NSF on behalf of the Gemini partnership: the National Science Foundation (United States), the National Research Council (Canada), CONICYT (Chile), Ministerio de Ciencia, Tecnología e Innovación Productiva (Argentina), and Ministério da Ciência, Tecnologia e Inovação (Brazil). We thank the Swift team for undertaking observations. The Fermi-LAT Collaboration acknowledges generous ongoing support from a number 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 Commissionariat à 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.

## ORCID iDs

R. H. D. Corbet https://orcid.org/0000-0002-3396-651X

L. Chomiuk https://orcid.org/0000-0002-8400-3705

J. B. Coley https://orcid.org/0000-0001-7532-8359

G. Dubus https://orcid.org/0000-0002-5130-2514

J. Strader https://orcid.org/0000-0002-1468-9668

## References

* Abdo, A. A., Ackermann, M., Ajello, M., et al. 2009, *Sci*, **325**, 840

* Abdo, A. A., Ackermann, M., Ajello, M., et al. 2010, *ApJS*, **188**, 405

* Abeysekara, A. U., Benbow, W., Bird, R., et al. 2018, *ApJL*, **867**, L19

* Acero, F., Ackermann, M., Ajello, M., et al. 2015, *ApJS*, **218**, 23

* Ajello, M., Atwood, W. B., Baldini, L., et al. 2017, *ApJS*, **232**, 18

* An, H., Bellm, E., Bhalerao, V., et al. 2015, *ApJ*, **806**, 166

* An, H., Dufour, F., Kaspi, V. M., et al. 2013, *ApJ*, **775**, 135

* An, H., & Romani, R. W. 2017, *ApJ*, **838**, 145

* Atwood, W. B., Abdo, A. A., Ackermann, M., et al. 2009, *ApJ*, **697**, 1071

* Baluev, R. V. 2008, *MNRAS*, **385**, 1279

* Bamba, A., Ueno, M., Nakajima, H., Mori, K., & Koyama, K. 2006, *A&A*, **450**, 585

* Bodaghee, A., Tomsick, J. A., Pottschmidt, K., et al. 2013, *ApJ*, **775**, 98

* Bosch-Ramon, V., Motch, C., Ribó, M., et al. 2007, *A&A*, **473**, 545

* Bruel, P., Burnett, T. H., Digel, S. W., et al. 2018, *arXiv:1810.11394*

* Burrows, D. N., Hill, J. E., Nousek, J. A., et al. 2005, *SSRv*, **120**, 165
Capalbi, M., Perri, M., Saija, B., Tamburelli, F., & Angelini, L. 2005, http://swift.gsfc.nasa.gov/analysis/xrt_swguide_v1_2.pdf

Cash, W. 1979, ApJ, 228, 939

Caswell, J. L., & Haynes, R. F. 1987, A&A, 171, 261

Clark, C. J., Wu, J., Pletsch, H. J., et al. 2017, ApJ, 834, 106

Coburn, W., Heindl, W. A., Rothschild, R. E., et al. 2002, ApJ, 580, 394

Corbel, S., Dubus, G., Tomsick, J. A., et al. 2012, MNRAS, 421, 2947

Corbet, R., Chomiuk, L., Coe, M., et al. 2019, AAS/High Energy Astrophysics Division, 17, 112.04

Corbet, R. H. D., Cheung, C. C., Kerr, M., et al. 2011, ATel, 3221, 1

Corbet, R. H. D., Chomiuk, L., Coe, M. J., et al. 2016, ApJ, 829, 105

Corbet, R. H. D., Coley, J. B., & Krimm, H. A. 2017, ApJ, 846, 161

Corbet, R. H. D., Coley, J. B., Krimm, H. A., et al. 2018, ATel, 11918, 1

Corbet, R. H. D., & Krimm, H. A. 2013, ApJ, 778, 45

Dubus, G. 2015, CRPhy, 16, 661

Dubus, G., Guillard, N., Petrucci, P.-O., & Martin, P. 2017, A&A, 608, A59

Dubus, G., Lamberts, A., & Fromang, S. 2015, A&A, 581, A27

Eikenberry, S. S., Elston, R., Raines, S. N., et al. 2004, Proc. SPIE, 5492, 1196

Evans, P. A., Beardmore, A. P., Page, K. L., et al. 2007, A&A, 469, 379

Green, D. A. 2014, BASI, 42, 47

Gregory, P. C. 2002, ApJ, 575, 427

Güver, T., & Oumlzel, F. 2009, MNRAS, 400, 2050

Hadasch, D., Torres, D. F., Tanaka, T., et al. 2012, ApJ, 749, 54

Hanson, M. M., Kudritzki, R.-P., Kenworthy, M. A., Puls, J., & Tokunaga, A. T. 2005, ApJS, 161, 154

Hill, J. E., Burrows, D. N., Nousek, J. A., et al. 2004, Proc. SPIE, 5165, 217

Jaron, F., Massi, M., Kiehlmann, S., et al. 2018, MNRAS, 478, 440

Johnson, T. J., Wood, K. S., Kerr, M., et al. 2018, ApJ, 863, 27

Kerr, M. 2011, ApJ, 732, 38

Koen, C. 1990, ApJ, 348, 700

Landi, R., Bassani, L., Bazzano, A., et al. 2017, MNRAS, 470, 1107

Leahy, D. A., Elsner, R. F., & Weisskopf, M. C. 1983, ApJ, 272, 256

Lee, K. J., Guillemot, L., Yue, Y. L., Kramer, M., & Champion, D. J. 2012, MNRAS, 424, 2832

Li, J., Torres, D. F., Cheng, K.-S., et al. 2017, ApJ, 846, 169

Mahy, L., Rauw, G., De Becker, M., Eenens, P., & Flores, C. A. 2015, A&A, 577, A23

Malyshev, D., & Chernyakova, M. 2016, MNRAS, 463, 3074

Marcote, B., Ribó, M., Paredes, J. M., Mao, M. Y., & Edwards, P. G. 2018, A&A, 619, A26

Martins, F., & Plez, B. 2006, A&A, 457, 637

Meurs, E. J. A., & van den Heuvel, E. P. J. 1989, A&A, 226, 88

Minniti, D., Lucas, P. W., Emerson, J. P., et al. 2010, NewA, 15, 433

Nolan, P. L., Abdo, A. A., Ackermann, M., et al. 2012, ApJS, 199, 31

Rajoelimanana, A. F., Charles, P. A., & Udalski, A. 2011, MNRAS, 413, 1600

Rea, N., Torres, D. F., Caliandro, G. A., et al. 2011, MNRAS, 416, 1514

Rieke, G. H., & Lebofsky, M. J. 1985, ApJ, 288, 618

Sault, R. J., Teuben, P. J., & Wright, M. C. H. 1995, adass IV, 77, 433

Saz Parkinson, P. M., Xu, H., Yu, P. L. H., et al. 2016, ApJ, 820, 8

Scargle, J. D. 1982, ApJ, 263, 835

Schinzel, F. K., Petrov, L., Taylor, G. B., & Edwards, P. G. 2017, ApJ, 838, 139

Schlaflay, E. F., & Finkbeiner, D. P. 2011, ApJ, 737, 103

Seward, F. D., Charles, P. A., Foster, D. L., et al. 2012, ApJ, 759, 123

Stellingwerf, R. F. 1978, ApJ, 224, 953

Süveges, M. 2014, MNRAS, 440, 2099

Takahashi, T., Kishishita, T., Uchiyama, Y., et al. 2009, ApJ, 697, 592

The Fermi-Lat Collaboration, Abdo, A. A., Ackermann, M., et al. 2009, Sci, 326, 5959

The Fermi-LAT Collaboration, Abdollahi, S., Acero, M., et al. 2019, arXiv:1902.10045

The Fermi-LAT Collaboration, Ackermann, M., Ajello, M., et al. 2012, Sci, 335, 6065

VanderPlas, J. T. 2018, ApJS, 236, 16

Verner, D. A., Ferland, G. J., Korista, K. T., & Yakovlev, D. G. 1996, ApJ, 465, 487

Wilms, J., Allen, A., & McCray, R. 2000, ApJ, 542, 914

Wilson, W. E., Ferris, R. H., Axtens, P., et al. 2011, MNRAS, 416, 832

Wu, J., Clark, C. J., Pletsch, H. J., et al. 2018, ApJ, 854, 99

Xing, Y., Wang, Z., & Takata, J. 2017, ApJ, 851, 92

Zanin, R., Fernández-Barral, A., de Oña Wilhelmi, E., et al. 2016, A&A, 596, A55

Zdziarski, A. A., Malyshev, D., Chernyakova, M., & Pooley, G. G. 2017, MNRAS, 471, 3657

Zdziarski, A. A., Malyshev, D., Dubus, G., et al. 2018, MNRAS, 479, 4399