LS 5039 is one of a handful of ray binary systems in the milky Way comprising a pulsar and a massive Otype companion star with an orbital period of 3.9 days. Recently we conducted a data analysis using approximately 16 yr of Fermi Large Area Telescope observations spanning from 2008 August 4 to 2024 July 8. In our timing analysis we discovered two new periodic signals with frequencies higher and lower than the known orbital period. The higherfrequency signal has a period of 3.63819 days with a 7.1 confidence level while the lowerfrequency signal has a period of 4.21654 days with a 6.3 confidence level. Additionally in data from the High Energy Stereoscopic System of Cherenkov telescopes two potential signals are found with periods similar to the two newly discovered ones. Considering that these two signals fall within the same frequency interval as the orbital period we suggest the possibility of a third body orbiting the barycenter of the LS 5039 binary system with the new periodic signals arising from specific frequency combinations of the two orbital periods.
['Gamma-rays', 'Gamma-ray sources', 'Periodic variable stars', '637', '633', '1213', 'Astrophysics - High Energy Astrophysical Phenomena']
Investigating the Potential of LS 5039 as a Triple System Using FermiLAT Data
{'Investigating the Potential of LS 5039 as a Triple System Using Fermi-LAT Data': "Lujun Zeng, 1 Mengqing Zhang, 1 Chongyang Ren, 2, 3 Pengfei Zhang, 1 and Jingzhi Yan 2, 3  \n- 1 Department of Astronomy, School of Physics and Astronomy, Key Laboratory of Astroparticle Physics of Yunnan Province, Yunnan University, Kunming 650091, People's Republic of China; zhangpengfei@ynu.edu.cn\n- 2 Key Laboratory of Dark Matter and Space Astronomy, Purple Mountain Observatory, Chinese Academy of Sciences, Nanjing 210034, People's Republic of China; jzyan@pmo.ac.cn\n- 3 School of Astronomy and Space Science, University of Science and Technology of China, Hefei, Anhui 230026, People's Republic of China", 'ABSTRACT': 'LS 5039 is one of a handful of γ -ray binary systems in the Milky Way, comprising a pulsar and a massive O-type companion star with an orbital period of 3.9 day. Recently, we conducted a data analysis using approximately 16 year of Fermi-LAT observations, spanning from 2008 August 4 to 2024 July 8. In our timing analysis, we discovered two new periodic signals with frequencies higher and lower than the known orbital period. The higher-frequency signal has a period of 3.63819 day with a 7 . 1 σ confidence level, while the lower-frequency signal has a period of 4.21654 day with a 6 . 3 σ confidence level. Additionally, in data from the High Energy Stereoscopic System of Cherenkov Telescopes, two potential signals with periods similar to the two newly discovered ones. Considering that these two signals fall within the same frequency interval as the orbital period, we suggest the possibility of a third body orbiting the barycenter of the LS 5039 binary system, with the new periodic signals arising from specific frequency combinations of the two orbital periods.  \nKeywords: Gamma-rays(637); Gamma-ray sources(633); Periodic variable stars(1213)', '1. INTRODUCTION': "To date, only a few γ -ray binary systems have been identified in our Galaxy (Aharonian et al. 2006; Albert et al. 2006; Hinton et al. 2009; Fermi LAT Collaboration et al. 2012a; Li et al. 2017; Corbet et al. 2019), with LS 5039 being one of them. LS 5039 is a well-studied high-mass X-ray binary, consisting of a rapidly rotating neutron star with a spin period of 9 s (Yoneda et al. 2020) and a massive O-type companion star (Ma'ız Apell'aniz et al. 2016) with a mass of 23 M ⊙ (Casares et al. 2005a). The two objects orbit each other with a period of 3.9 day in a moderately eccentric orbit with an eccentricity of 0.35 and an orbital inclination of i = 25 · , at a distance of 2.5 kpc from Earth within our Galaxy (Casares et al. 2005a). The orbital phase zero (i.e. ϕ 0 ) is defined at its periastron passage, with the reference epoch T 0 set at HJD 2451943.09 (i.e. MJD 51942.59; Casares et al. 2005a). Notably, LS 5039 is one of only three known γ -ray binary systems in which the compact object is a neutron star, the other two being LS I + 61 · 303 (Weng et al. 2022) and PSR B1259 -63/SS 2883 (Johnston et al. 1992).  \nLS 5039's electromagnetic emissions span a broad spectrum, ranging from radio frequencies to TeV γ -ray energies (Mold'on et al. 2012; Casares et al. 2005a; Takahashi et al. 2009; Abdo et al. 2009; Aharonian et al. 2006) and potentially reaching PeV (Cao et al. 2021; Bykov et al. 2024). Its spectral energy distribution, which is dominated by MeV-GeV γ -rays, characterizes LS 5039 as a γ -ray binary. The 3.9 day orbital modulation of LS 5039 has been identified by multi-wavelength observations and is well-supported by extensive literature. Notably, at energies above 0.1 TeV, Aharonian et al. (2006) presented its orbital modulation with a period of 3 . 9078 ± 0 . 0015 day based on the dataset from the High Energy Stereoscopic System (HESS) of Cherenkov Telescopes. In GeV, LS 5039 was detected by Abdo et al. (2009) using the 11 month Fermi -LAT observations, who also reported that its γ -ray emissions are modulated by a 3 . 903 ± 0 . 005 day orbital period. Recently, Zhang (2024) discovered a new periodic signal with a period of 26.3 day in the γ -ray binary LS I + 61 · 303 using Fermi -LAT observations. Motivated by these findings, we conducted a new timing analysis of LS 5039, hoping  \nTable 1. Results of likelihood  \nNotes. Best-fit parameters of the likelihood, F ph in unit of × 10 -7 photons cm -2 s -1 , E b and E c in unit of GeV.  \nto uncover phenomena beyond its known orbital period in GeV gamma-rays.  \nIn the fourth catalog Data Release 4 (4FGL-DR4; Ballet et al. 2023), LS 5039 is associated with the γ -ray counterpart 4FGL J1826.2 -1450, which has been monitored by Fermi-LAT for approximately 16 year. In our timing analysis, in addition to the detection of the orbital period of 3.90609 day, we identified two new periodic signals with frequencies higher and lower than the orbital period. These signals have periods of 3.63819 day and 4.21654 day, with confidence levels of 7.1 σ and 6.3 σ , respectively. Moreover, two similar periodic signals were observed in data from the HESS of Cherenkov Telescopes (Aharonian et al. 2006). The details of our data analysis and findings are presented below. In them, the numbers in parentheses indicate the uncertainties associated with the values.", '2.1. Data Reduction': "For our data analysis, we utilized the Pass 8 Front+Back events (evclass = 128 and evtype = 3), focusing on the energy range of 0.1500.0 GeV within a 20 · × 20 · region of interest centered on 4FGL J1826.2 -1450 (R.A. = 276.5637, Decl. = -14.8496), covering the period from MJD 54682.687 to 60499.299. Events with zenith angles < 90 · were selected, and only high-quality events from good time intervals were retained using the expression 'DATA QUAL > 0 && LAT CONFIG==1'. The Galactic and isotropic diffuse emissions were modeled using the templates of gll iem v07.fits and iso P8R3 SOURCE V2 v1.txt, respectively. The data was reduced and analyzed using the Fermitools version 2.2.0.  \nA model file was created based on 4FGL-DR4, and a binned maximum likelihood analysis was performed to update the model parameters. For 4FGL J1826.2 -1450, the 4FGL-DR4 provides a log-parabola (LP) spec-  \ntral shape, dN/dE = N 0 ( E/E b ) -[ α + β log( E/E b )] . Additionally, we also employed a power-law with an exponential cutoff (PLEC) model, dN/dE = N 0 ( E/E 0 ) -Γ exp[ -( E/E c ) b ], to fit the events around 4FGL J1826.2 -1450. All parameters for these two models are summarized in Table 1. As shown in it, we conclude that the LP model, which yields a higher TS value of 29,003.2, provides a better representation of the γ -ray emissions from LS 5039 compared to the PLEC model, which has a TS value of 28,912.7. The LP model results were then saved as the bestfit model. Based on it, we generated a TS map for 4FGL J1826.2 -1450's γ -ray emissions (see left panel of Figure 1) and a residual TS map excluding γ -rays from all sources listed in the 4FGL-DR4 (right panel of Figure 1). These maps indicate that the γ -ray emissions from 4FGL J1826.2 -1450 is well-described by the LP model.", '2.2. Light Curve Construction and Timing Analysis': "A modified version of the aperture photometry (AP) method was employed to construct the light curve for 4FGL J1826.2 -1450 in the 0.1-500.0 GeV. To optimize the signal-to-noise ratio, an aperture radius of 3 · . 16 was selected, taking into account the performance characteristics of the LAT instrument, with a criterion of θ < max(6.68 -1.76log 10 ( E MeV ) , 1 . 3) · , 1.3), following that performed in Abdo et al. (2010). The light curve was binned at 500 s intervals for each time bin. To mitigate contamination from γ -rays originating from the Sun and Moon, we excluded events occurring when the target was within 5 · of either the Sun or the Moon, using gtmktime . To address significant exposure time variations across different time bins, the exposure time for each bin was calculated using gtexposure . Using the best-fit model, event probabilities for 4FGL J1826.2 -1450 were determined by employed gtsrcprob and used as weights for constructing the AP light curve, rather than simply counting photon numbers (Kerr 2011; Fermi LAT Collaboration et al. 2012b; Corbet et al. 2019). Additionally, the times in the light curve were corrected to the barycenter using gtbary .  \nIn our timing analysis, a method of the generalised Lomb-Scargle periodogram (LSP; Lomb 1976; Scargle 1982; Zechmeister & Kurster 2009) was employed to derive the power spectrum for the AP light curve. Compare to other period search algorithms, the LSP offers more accurate frequency determination for unevenly spaced data, is less susceptible to aliasing, provides better frequency resolution, improves the detection of weak signals, and facilitates the assessment of the confidence level for periodic signals. This same methodology was  \nFigure 1. TS maps in the 0.1-500.0 GeV range, covering a 3 · × 3 · region centered on 4FGL J1826.2 -1450, with γ -ray sources from the 4FGL-DR4 marked by crosses. Left panel: TS map representing γ -ray emissions from 4FGL J1826.2 -1450, generated using the best-fit model with the target excluded. Right panel: Residual TS map created using the same model, except with 4FGL J1826.2 -1450 included.  \n<!-- image -->  \nFigure 2. LSP power spectra (black histogram) derived from the 0.1-500.0 GeV AP light curve of 4FGL J1826.2 -1450. Left panel: The three periodicities, P H , P orb , and P L , are in close proximity. The gray dashed-dotted line marks the location of the 53.4 day precession period of the Fermi -LAT orbit. Right panel: LSP power spectrum zoomed in the left. The green dotted and pink dashed lines indicate the signals of P H and P L , respectively. The horizontal red dashed and blue dashed-dotted lines stand for 5 σ and 7 σ confidence levels. HESS results that drawn from the Figures 2 and 3 of Aharonian et al. (2006) are shown in the inset plot.  \n<!-- image -->  \nFigure 3. LSP power spectra for 4FGL J1827.4 -1445 (black histogram) and 4FGL J1829.4 -1500c (orange histogram). The vertical green and pink arrows show the positions of the two signals of P H and P L . The gray dasheddotted line is the precession period of the Fermi-LAT orbit. And the red dashed line is 5 σ confidence level.  \n<!-- image -->  \nalso employed by Zhang (2024). The power spectrum covers a frequency range from f max = 1/0.06 day -1 to f min = 1/5816.612 day -1 (e.g., the inverse of the length of Fermi -LAT observations, 1 /T obs ). We show the power spectrum in left panel of Figure 2, in it a prominent main peak appears at 3.90609(3) day, corresponding to the orbital period ( P orb ) of LS 5039. In addition to P orb , two peaks can be observed on either side of it. One period appears at the higher frequency with period of 3.63819(14) day ( P H ) and another at the lower frequency at 4.21654(20) day ( P L ). For the three periods, their uncertainties were calculated by δP = 3 8 P 2 T obs √ p n , where P is a period, T obs is the length of Fermi -LAT observations, and p n is the normalized height of the peak power, as that provided in Horne & Baliunas (1986). These periods, P orb , P H , and P L , correspond to frequencies of f orb = 0 . 256010(2) day -1 , f H = 0 . 274862(10) day -1 , and f L = 0 . 237161(11) day -1 , respectively. Moreover, the frequency difference between P H and P orb ( f H -f orb = 0 . 018851) is nearly equal to the difference between P orb and P L ( f orb -f L = 0 . 018849). In left panel of Figure 2, we marked the P orb , P H , and P L with orange solid, green dotted, and pink dashed lines, respectively. And the first, second, third, and fourth harmonics of P orb also appear, we mark them with four purple arrows. For clarity we zoom in the power spectrum in the gray shaded region in right panel of Figure 2.  \nThe heights of the peaks of P orb , P H and P L are ∼ 865 . 88, ∼ 36 . 38, and ∼ 33 . 10, respectively, compared to the mean power level, based on the normalization  \nprovided in Horne & Baliunas (1986). For P H and P L , the probabilities ( p lsp ) of a periodic signal reaching a power level by chance fluctuation assuming the Gaussian white noise were calculated based on method provided by Lomb (1976) and Scargle (1982). And their values are 1 . 20 × 10 -17 and 4 . 21 × 10 -15 , respectively. In our timing analysis, the number of independent frequencies (i.e., the trial factor) was calculated by N = ( f max -f min ) /δf = 96943, where δf is frequency resolution determined by the length of the Fermi -LAT observations. Then the False Alarm Probability (FAP) is calculated by the form of FAP=1 -(1 -p lsp ) N ∼ p lsp × N , after taking into account the trial number N . And their FAPs are derived to be 1 . 16 × 10 -12 and 4 . 09 × 10 -10 for P H and P L , which corresponding to ∼ 7 . 1 σ and ∼ 6 . 3 σ confidence levels. In right panel of Figure 2, we show 5 σ and 7 σ confidence levels with red dashed and blue dashed-dotted lines, respectively. Interestingly, the HESS data reported by Aharonian et al. (2006) also present two potential signals with periods of 3.42(5) and 4.56(10) day, which are similar to the two signals reported here. The errors associated with these periods were derived from the fullwidths at half-maximum of the power peaks. For easy reference, we also include a schematic of the results for LS 5039 in the inset of right panel of Figure 2, adapted from Aharonian et al. (2006).  \nConsidering the broad point spread function of the Fermi -LAT, periodic signals from neighboring sources may also appear in the power spectrum of the target. From the residual TS map (right panel of Figure 1), we can see that there are no new γ -ray sources beyond those listed in the 4FGL-DR4. Therefore, we constructed power spectra for the two γ -ray sources closest to the target, 4FGL J1827.4 -1445 and 4FGL J1829.4 -1500, using the same process, and present them in Figure 3 with black and orange histograms, respectively. We observe a prominent peak corresponding to the P orb of LS 5039, along with a weaker peak corresponding to its first harmonic. However, no signals corresponding to P H and P L are detected.  \nTo test whether these two signals are energydependent, we constructed light curves in three energy bands, 0.1-0.3, 0.3-1.0, and 1.0-500.0 GeV, using the same aperture photometry method. Their power spectra are shown in Figure 4 with black, red, and green histograms, respectively. We observe that the peak values of these two new signals, as well as the orbital period signal, decrease as the energy increases across the three energy ranges.  \nFigure 4. LS 5039's energy-dependent LSP power spectra in 0.1-0.3, 0.3-1.0, and 1.0-500.0 GeV. The black, red, and green histograms stand for the spectra constructed with 0.10.3, 0.3-1.0, and 1.0-500.0 GeV AP light curves, respectively. As in Figure 2, the red dashed and blue dashed-dotted lines are 5 σ and 7 σ confidence levels.  \n<!-- image -->  \nFigure 5. Phase-resolved light curves derived from the 0.1500.0 GeV Fermi-LAT events, folded over the periods P H (top panel) and P L (bottom panel). Both panels share the same y-axis labels.  \n<!-- image -->  \nTo examine the periodic profiles of the two signals, we derived phase-resolved light curves based on the periods of P H and P L by dividing the 0.1-500.0 GeV FermiLAT events into 20 phase intervals. For each phase interval, we performed a likelihood analysis to determine the fluxes, using phase zero at periastron and setting the reference epoch T 0 at MJD 51,942.59, as described by Casares et al. (2005b). We present the phase-resolved light curves for P H and P L in the upper and lower panels of Figure 5, respectively. These light curves clearly demonstrate periodic modulation profiles that differ from the modulation profile of the orbital period. Moreover, there are also significant differences between the profiles of P H and P L .", '3. SUMMARY AND DISCUSSION': "Here, we conducted a timing analysis using the LSP method on approximately 16 years of Fermi -LAT observations. Our analysis identified two new periodic signals with frequencies higher and lower than the orbital period (3.90609 day). These signals have periods of 3.63819 day and 4.21654 day, with confidence levels of 7.1 σ and 6.3 σ , respectively. Moreover, these new periodic signals share the same frequency interval as the orbital period. All three signals, including the two new periodic signals and the orbital period, show a trend of decreasing strength with increasing energy. This trend may be attributed to statistical factors, as there are much more γ -rays in lower energies than that in higher energies. Additionally, in previous studies, the HESS data appeared two potential signals with similar periods (Aharonian et al. 2006) to those we discovered in our Fermi-LAT analysis, which independently support the authenticity of our new findings. We know that in astrophysics, there are many causes that can produce periodic signals, but it is quite rare for two signals to have the same frequency interval. Considering that the P H and P L signals share the same interval as the orbital period in the frequency domain, we will discuss two scenarios that could generate these signals in the following discussion.  \nDuring timing analysis, peaks with the same frequency interval in the power spectrum can originate from the sidelobes of the main peak (i.e., the orbital period; Bloomfield 1976), caused by the finite length of observations and the use of windowing functions. These sidelobes are evenly distributed on both sides of the main peak, with their peak values gradually decreasing, which makes the first sidelobe the most likely candidates. In this scenario, we can obtain the locations of the first sidelobe using the formula: P sidelobe = P orb ± 7 5 P 2 orb /T obs in our timing analysis, where P sidelobe is the periods corresponding to the first sidelobe, and T obs is the length  \nof Fermi-LAT observations during the timing analysis. The values for the first sidelobes are calculated to be 3.902420 or 3.909765 day, which differ significantly from the P H and P L values. Therefore, we believe that the P H and P L signals cannot originate from the sidelobes of the orbital period of LS 5039. Additionally, if the two signals originated from the sidelobes, they should also appear in the power spectra of the two neighboring γ -ray sources. However, the power spectra of these two sources are significantly influenced by the orbital period of LS 5039, and they do not show the P H and P L signals. This further indicates that the P H and P L signals are related to LS 5039 and not associated with the sidelobes.  \nWe then introduce an alternative scenario that could also produce two periodic signals with the same frequency interval. We propose that there may be a third body in the LS 5039 system, suggesting that LS 5039 could be a three-body system, as that reported in Ransom et al. (2014); Rebassa-Mansergas et al. (2022). The third body orbits the barycenter of the original binary system, with a hypothesized orbital period of P hpy . And P H and P L signals likely originate from certain frequency combinations involving the P orb and P hyp . To simplify this process, we omit the phase shift parameters and assume that the P hyp modulation has a sinusoidal function. And the γ -ray flux ( F γ ) is modulated by a form of  \nF γ = A sin(2 πtf hyp ) + C. (1)  \nSince the P orb signal is much stronger than the P H and P L in the power spectrum (Figure 2), hence we assume the amplitude ( A ) and constant ( C ) of P hyp are strongly modulated by the main periodic signal ( P orb ) and they have the forms of  \nA = A 1 sin(2 πtf orb ) + C 1 (2)  \nand  \nThen we simplify this process again and assume A 1 = A 2 = 1 and C 1 = C 2 = 0. Then,  \nF γ = 0 . 5 { cos[2 πt ( f orb -f hyp )] -cos[2 πt ( f orb + f hyp )] } +sin(2 πtf orb ) . (4)  \nConsequently, we will see two peaks (with the corresponding frequencies to f H = f orb + f hyp and f L = f orb -f hyp ) around the orbital period in the power spectrum. The hypothesized period can be calculated as P hyp = 1 f hyp = 2 f H -f L = 53 . 049(21) day, with the uncertainty derived from the propagation of errors in f H and f L .  \nAdditionally, we note that Fermi -LAT operations can introduce artificial periodic signals (Atwood et al. 2009). One such signal, which is close to the P hyp , corresponds to the Fermi -LAT's orbital precession period of 53.4 day 1 , we indicate it in Figures 2 (left) and 3 by gray dashed-dotted lines. However, there are no peaks at this precession period in the power spectra of 4FGL J1826.2 -1450 and the two neighboring γ -ray sources, and neither P H or P L appears in the power spectra of the neighboring sources. Therefore, we conclude that our light curves are not influenced by this precession period, furthermore the 53.4 day period does not fall within the uncertainties of P hyp . Moreover, the HESS dataset also reveal two potential signals (Aharonian et al. 2006) with periods similar to those reported here. Thus, we infer that P H or P L are not related to the Fermi -LAT's orbital precession period.  \nRegardless, this discovery enhances our knowledge of LS 5039 and offers new insights into the dynamics and evolution of massive binary systems. The cause of these two new periodic signals remains uncertain. Further observations, particularly multi-wavelength ones, are encouraged to confirm the authenticity of the two periodic signals and to reveal their origins.  \n- 10\n- 11\n- Wewould like to thank the anonymous referee for helpful suggestions and Jian Li for useful comments. This work is supported in part by the National Natural Science Foundation of China under grant Grants Nos. 12233006, U2031205, and 12163006, the Basic Research Program of Yunnan Province No. 202201AT070137, and the joint foundation of Department of Science and Technology of Yunnan Province and Yunnan University No. 202201BF070001-020. P.Z. acknowledges the support by the Xingdian Talent Support Plan - Youth Project. 1 2 3 4 5 6 7 8 9  \nC = A 2 sin(2 πtf orb ) + C 2 . (3)"}