A Chandra Search for Periodic X-ray Sources in the Bulge of M31

We present a systematic search for periodic X-ray sources in the bulge of M31, using ~ 2 Ms of archival Chandra observations spanning a temporal baseline of 16 years. Utilizing the Gregory-Loredo algorithm that is designed for photon-counting, phase-folded light curves, we detect seven periodic X-ray sources, among which four are newly discovered. Three of these sources are novae, the identified periods of which range between 1.3-2.0 hour and is most likely the orbital period. The other four sources are low-mass X-ray binaries, the identified periods of which range between 0.13-19.3 hour and are also likely orbital due to a clear eclipsing/dipping behavior in the light curve. We address implications on the X-ray binary population of the M31 bulge. Our study demonstrates the potential of using archival X-ray observations to systematically identify periodic X-ray sources in external galaxies, which would provide valuable information about the underlying exotic stellar populations.


INTRODUCTION
The Andromeda galaxy (= M31, at a distance of 785 kpc; Mc-Connachie et al. 2005) is a massive spiral with a substantial stellar bulge.Thanks to its proximity to the Milky Way, M31 provides a unique opportunity for studying galaxies similar in physical properties with our Galaxy.In particular, the X-ray emissions of M31 have been observed since the era of the Einstein satellite (van Speybroeck et al. 1979) and have subsequently been studied with improved spatial resolution, spectral resolution and sensitivity through the utilization of many X-ray missions, especially Chandra (Kong et al. 2002;Kaaret 2002) and XMM-Newton (Shirey et al. 2001;Osborne et al. 2001;Pietsch et al. 2005).
M31 holds hundreds of bright X-ray point sources, mostly accretion-powered X-ray binaries (XRBs) involving a stellar-mass black hole (BH), a neutron star (NS) or a white dwarf (WD), which can be studied in detail thanks to low Galactic foreground absorption.These sources have been the subject of systematic investigations over the past decades, e.g., by Kong et al. (2002); Pietsch et al. (2005); Voss & Gilfanov (2006); Hofmann et al. (2013); Barnard et al. (2014), among others.In the bulge of M31, where the stellar population is predominantly old, low-mass X-ray binaries (LMXBs) make up the majority, if not all, of the XRBs (Hofmann et al. 2013).Another important component is cataclysmic variables (CVs), which experience mass transfer from a main sequence star or a red giant star to a WD (Robinson 1976).The thermonuclear runaway in the accreted matter triggers classical nova events.The classical novae in the central region of M31 have been monitored by XMM-Newton, Chandra and Swift from 2006 to 2011 (Henze et al. 2010(Henze et al. , 2011(Henze et al. , 2014)).Additionally, ultraluminous X-ray sources (ULXs), super-soft sources (SSSs, Kahabka & van den Heuvel 1997), and quasi-soft sources (QSSs) ★ E-mail: zhangjiachang@smail.nju.edu.cn† E-mail: baotong@smail.nju.edu.cn‡ E-mail: lizy@nju.edu.cn are known to exist among the X-ray populations of M31.These diverse stellar populations serve as a natural and promising laboratory to advance our understanding of the underlying physical processes related to accretion and binary evolution.
Many X-ray sources in M31 have exhibited strong flux variability, a characteristic of accretion-powered systems.A particularly interesting feature is periodic variability, which arise due to either binary orbital motion or the spin of the accretor.Over the years, a number of periodic X-ray signals have been identified in M31 by either dedicated searches or serendipitous discoveries.Table 1 summarizes the previously known X-ray periodic sources in M31.These include four novae with periodic signals, three SSSs with pulsation periods, four LMXBs with orbital or super-orbital periods, among which one also shows a pulsation period.These periodic variations provide valuable insights into the dynamics and characteristics of the X-ray binary populations in M31.Yet, compared to the fraction of the Galactic XRBs with known orbital periods (145 in 348 for LMXBs, Avakyan et al. 2023;112 in 169 for high-mass X-ray binaries [HMXBs], Neumann et al. 2023), far less is currently known for the M31 XRB populations.
In this work, we conduct a systematic search for periodic X-ray signals from a large number of X-ray sources in the bulge of M31, based on extensive Chandra monitoring observations taken over the past two decades.These observations provide a rare but promising dataset for probing periodic X-ray signals from an extragalactic population of XRBs.In Section 2, we describe the preparation of the X-ray data and the construction of a list of X-ray sources in the M31 bulge.In Section 3, we detail the period searching method, procedures and results.The X-ray spectra of the periodic X-ray sources are analyzed in Section 4, followed by an examination of the physical nature of the periodic sources based on their temporal and spectral properties in Section 5.In Section 6, we address the implications and limitations of our results.A brief summary of our study is given in Section 7.

X-RAY DATA PREPARATION
The bulge of M31 has been extensively observed by Chandra using its Advanced CCD Imaging Spectrometer (ACIS) and High Resolution Camera (HRC).We require that the aim-point of a selected observation was ≲ 1 ′ from the center of M31, which helps to minimize variations in the point-spread function (PSF) across different observations at any given position, ensuring consistent and accurate photometry for the X-ray point sources.This results in a total of 116 ACIS observations taken between October 13, 1999 and November 28, 2015 with a total exposure time of ∼ 996.9 ks, and 64 HRC observations taken between November 30, 1999 and June 1, 2012 with a total exposure of ∼ 1032.4 ks.These observations together span a temporal baseline of over 16 yrs, but with a rather irregular cadence.
The information of the ACIS and HRC observations utilized in this work is given in Table B1 and Table B2.
We downloaded and uniformly reprocessed the raw data for each observation, following the standard procedure 1 and using CIAO v.4.14 and calibration data files CALDB v.4.9.8.The ACIS and HRC observations were treated separately.The CIAO tool acis_process_event and hrc_process_events were used to reprocess the level 1 event files to create level 2 events, correcting charge transfer inefficiency and filtering cosmic-ray afterglows.For the ACIS observations, we included only data from the I0, I1, I2 and I3 chips if the aim-point was placed on the I3 chip, and only data from the S2 and S3 chips if the aim-point was placed on the S3 chip.The analysis of the light curve for each observation revealed that the instrumental background remained quiescent for the majority of the time intervals.Hence all science exposures were preserved for the subsequent timing analysis, ensuring an uninterrupted light curve within each observation.Additionally, the CIAO tool axbary was employed to correct the photon arrival time in each registered event to the Solar system barycenter (i.e.Temps Dynamique Barycentrique time).
To align the astrometry among individual ACIS or HRC observations, we performed centroid matching for commonly detected point sources using the CIAO tool reproject_aspect.The observation with the longest exposure time (i.e., ObsID 14196 for ACIS and ObsID 1912 for HRC) was selected as the reference frame for this alignment process.The stacked counts map over the energy range of 0.5-8 keV for ACIS and the PI range of 50-300 for HRC was then created respectively.The stacked HRC counts map of the inner 15 ′ × 15 ′ (∼ 3.4 kpc × 3.4 kpc) of M31 is shown in Figure 1.Exposure and PSF maps were generated, with the PSF maps being computed according to a 1 http://cxc.harvard.edu/ciaospecified enclosed count radius (ECR).A fiducial source spectrum was employed for both the exposure and PSF maps.This spectrum is characterized by an absorbed power-law model with a photon index of 1.7 and a hydrogen column density consistent with the Galactic foreground value ( H = 7.0 × 10 20 cm −2 ), representative of the bulge X-ray sources (Hofmann et al. 2013).
Point source detection was conducted following the procedures outlined in Zhu et al. (2018).Briefly, the CIAO tool wavdetect was applied to the stacked energy = 0.5-8 keV counts map of ACIS and stacked PI = 50-300 counts map of HRC.We restricted the region of interest within a projected distance of 500 ′′ from the M31 center, which not only ensures an optimal PSF, but also minimizes the fractional contamination by cosmic X-ray background (CXB) sources and potential HMXBs related to the star-forming disk of M31.The algorithm was fed with the stacked exposure map and the 50% ECR PSF map, setting a false-positive probability threshold of 10 −6 .The centroid of each source provided by wavdetect was further refined using a maximum-likelihood method that iterates over the recorded positions of the individual counts detected within the 90% ECR.Our detection resulted in a raw list of 456 ACIS sources and 403 HRC sources, which form the basis of the period searching (Section 3).A cross-matching of ACIS-detected and HRC-detected sources suggests a total of 647 positionally independent sources.Voss & Gilfanov (2006) published an X-ray source catalog based on a large fraction of the ACIS observations used here, whereas Hofmann et al. (2013) also published a source catalog using the same set of 64 HRC observations.We find that almost all sources within the same region of interest in these two catalogs are recovered in our source lists.We defer a detailed description and analysis of the updated source catalog of the M31 bulge to future work.
For each source in each ObsID, we adopted a default 90% ECR to extract the source events, again focusing on the energy range of 0.5-8 keV for ACIS and PI range of 50-300 for HRC.Source events from all ACIS or all HRC observations together form the time series of a given source to be fed to the period searching algorithm (Section 3.1).Crowding of X-ray sources is not a general concern for the M31 bulge.For a handful of sources which have a close neighbor source, we have also tested an extraction region of 75% ECR, but found no significant difference in the result.Local background counts were extracted from a concentric annulus with inner-to-outer radii of 2-4 times the 90% ECR, masking pixels falling with the 90% ECR of neighboring sources, if any.
We examined possible flux variability by quantifying the net photon flux in the individual observations.The total source counts, background counts, source area, background area, and local exposure  2.
were fed to the CIAO tool aprates to calculate the net photon flux and associated error, which accounts for the Poisson statistics in the low-count regime.If the 3 lower limit of the photon flux hits zero in a certain observation, the source is considered non-detected and we derive a 3 upper limit of the source photon flux for this particular observation.

Method
We apply the Gregory-Loredo (GL) algorithm (Gregory & Loredo 1992) to detect potential periodic signals, which is a Bayesian-based, phase-folding method.The GL algorithm is effective in identifying periodic signals from X-ray data that is often characterized by a moderate number of photon events and/or an irregular observing cadence, as is the case of M31 here.This algorithm has been successfully employed to detect periodic X-ray sources in a Galactic bulge field (Bao & Li 2020), the Chandra Deep Field South (Bao & Li 2022), as well as a number of Galactic globular clusters (Bao et al. 2023(Bao et al. , 2024)), all based on deep Chandra observations.Readers are referred to these papers for details of the GL algorithm and specific implementation for Chandra data, as well as potential caveats.A concise description of the working principle of the GL algorithm is given in Appendix A. Main steps of the period searching procedure are as follows.
By design, the GL algorithm folds the time series of a given source at trial frequencies (periods).The chosen resolution and frequency range is a compromise of efficiency and computational resource.We primarily explore period ranges of (100, 500), (500, 1000), (1000, 5000), (5000, 10000) and (10000, 20000) seconds, with a frequency resolution of 10 −5 , 10 −6 , 10 −7 , 10 −8 and 10 −9 Hz, respectively.Given the timespan of ∼ 5 × 10 8 sec between the first and last Chandra observations, the chosen frequency resolutions are optimal for an efficient search of periodic signals.Since the GL algorithm only determines the most probable period, once a tentative period is identified, a second search is performed excluding a narrow interval around the identified period, which ensures that a possible second period within the same period searching range will not be missed.
In the second step, we rerun the algorithm on the time series from each individual observation.If a signal of roughly the same period were detected in more than one observations, we rerun the GL algorithm on the combined time series from those observations to refine the period.We do not attempt to combine the ACIS and HRC observations due to their substantially different instrumental response.
Although we are also interested in signals of lower frequencies (longer periods), they are hard to identify due to potentially significant flux variations among observations that may cause numerous spurious signals.Therefore, we search for periods between 20000-40000 sec using only a series of closely-packed long observations (ObsID 13825,13826,13827,13828,14195,15267,14196 for ACIS) and (ObsID 5925,6177,5926,6202,5927,5928 for HRC) at a frequency resolution of 10 −9 Hz.
We note in passing that the ACIS and HRC data in principle allow for probing shorter periods (< 100 sec), e.g., due to spin modulation of neutron stars.For ACIS, a natural lower limit is set by the CCD frame time, which is typically 3.2 sec; For HRC, the actual limit inversely depends on the rate of registered counts, which could be also as low as ∼20 ms for the present case.However, essentially all the X-ray sources in M31 have a moderate or low count rate, which introduces a substantial Poisson noise in the phase-folded light curve, challenging the identification of a true periodic signal.Nevertheless, this has been attempted by us, which resulted in several periods on the order of a few tens of sec.However, all such short periodic signals were only found in a single observation and had a relatively low value (∼ 0.9) of probability threshold ( GL ).Given the large number of single-epoch light curves fed to the GL algorithm, we cannot rule out the possibility that these short-period signals are simply false detections due to noise.Moreover, these periods are not only too low to be associated with an orbital modulation, but also too short-lived to be related to a spin modulation, which ought to be persistent.Therefore, we do not consider them genuine periodic signals.
We adopt a probability threshold  GL of 90% for selecting the tentative periods returned by the GL algorithm.Some of these tentative periodic signals, however, could be spurious due to one of the following reasons.Therefore we apply further filtering to the tentative periodic signals: (i) Eliminate periodic signals caused by the dither pattern of the ACIS/HRC detectors.For observations with ACIS, the nominal peri-ods are 1000 s in yaw and 707 s in pitch, while for observations with the HRC, dithering follows nominal periods of 1087 s in yaw and 768 s in pitch2 .When part of the source region falls onto a bad pixel or goes outside the CCD edge, a spurious signal modulated by the dithering could occur.Any signals detected at the two dither periods or their harmonics to within 1 per cent are considered spurious and excluded.
(ii) Exclude fake signals produced by strong, aperiodic source variability.For sources exhibiting strong (up to a factor of 10 variation in the net photon flux) long-term (i.e.inter-observation) variations, we repeat the period search in two subsets of the light curve: one covering only the high state (defined as the observation[s] with the highest photon flux) and the other excluding the high state.If a periodic signal arising from the full time series could not be reproduced in either subset, we exclude it as a false detection.Moreover, for each tentative periodic source, we inspect its light curve to identify short-term (i.e.intra-observation) flares, and subsequently repeat the period search after discarding any observation(s) in which significant short-term flares are present.Only when a periodic signal survives this procedure it is considered a genuine signal.
In Section 3.3, we further address the possibility of false detection of periodic signals due to an intrinsic red noise of the X-ray sources.

Resultant periodic signals
The above procedures result in a total of seven periodic signals/sources in the M31 bulge, among which four are newly discovered.Table 2 lists these periodic sources in the order of increasing period, giving information on the source ID, position, period, the GL probability and the relevant ObsID, as well as tentative classification (see details in Section 5).The detected periods range from 463.24 sec to 69417.50 sec.It is noteworthy that the longest period, found for Seq.7, actually exceeds our primary period searching range.The GL algorithm originally reported a period half of this value (∼ 34707 sec).After further examination of multiple, adjacent observations, we determine the true period, which is ∼ 69417.50sec, from the eclipsing behavior of this source (see details in Section 3.3).
Among the seven sources, three (Seq.4,Seq.6 and Seq.7 in Table 2) have a persistent periodic signal (i.e., detected over all ACIS and all HRC observations).Seq.1, while itself being a persistent source, has its periodic signal detected in only 5 HRC observations and 3 ACIS observations (marked in the top middle panel of Figure 2).For these four sources, the periods derived from the ACIS and HRC observations differ by ∼0.01%-0.3%,which can be understood as the typical uncertainty of periodic signals determined by the GL algorithm (Bao et al. 2023).The phase-folded light curve, interobservation light curve and cumulative ACIS spectrum (Section 4) of these four sources are shown in Figure 2. To phase-fold the light curve, we adopt either the ACIS-detected or HRC-detected period, whichever has the more source counts.A phase consistency between the ACIS and HRC light curves is clearly evident in Seq.4,Seq.6 and Seq.7.For Seq.1, due to a ∼0.3% difference between the ACIS-and HRC-derived periods, a significant phase shift is seen.The remaining three sources (Seq.2,Seq.3 and Seq.5) have their periodic signal detected in only 2-5 observations (all are HRC) when each source was caught in an outburst.The phase-folded light curves of these three sources are shown in Figure 3.We defer more detailed discussions of the periodic signals to Section 5.The periodic signal of Seq.1 is only detected in five HRC and three ACIS observations, which are marked with a square.For this source, a spectrum combining the three ACIS observation is also fitted (shown in red).The fitted spectral models are denoted.(1) Source sequence number assigned in the order of increasing period.( 2) and (3) Right Ascension and Declination (J2000) of the source centroid.(4)The ID of observations based on which the periodic signal is detected.( 5)-( 8) The period and the  GL determined by the GL algorithm, the number of total counts in the 0.5-8 keV band and the number of estimated background counts of ACIS observations.( 9)-( 12) The period and the  GL determined by the GL algorithm, the number of total counts in the PI range of 50-300 and the number of estimated background counts of HRC observations.( 13

The effect of red noise
Accretion-powered systems, such as CVs and LMXBs, are known to exhibit aperiodic variability across a wide range of time scales.The so-called red noise, which is a significant component of the aperiodic variation, can potentially introduce false periodic signals, particularly at lower frequencies (Warner 1989).It is therefore instructive to estimate the possibility of false alarms among the reported signals by the GL algorithm.Following Bao et al. (2023), we adopt a power-law model to describe the source power spectrum in order to account for the presence of red noise: Here  is the normalization factor,  is the spectral index, and  P represents the Poisson noise dictated by the mean photon flux of the source.To mitigate the potential effect of interrupted observations in the Fourier analysis, we utilize the longest ACIS or HRC observation for each periodic source to characterize their power spectrum.This approach ensures a continuous and extended observation period, avoiding the influence of interruptions on the analysis.The power spectrum of a given source is fitted with Eq. 1 using the Markov Chain Monte Carlo approach (with the python emcee package, Foreman-Mackey et al. 2013) to determine the best-fit parameters and errors.
Next, following the procedure proposed by Timmer & Koenig (1995), we simulate 1000 time series using the best-fit power spectrum model of each source, which are fed to the GL algorithm.The histogram of the resultant  GL for each source is shown in Figure 4.The fraction of false signals (i.e., those with  GL above the default threshold of 0.9) represents the false alarm probability (FAP) for the source.
The majority of detected periodic sources have a FAP < 1%, which can be taken as strong evidence for an intrinsic periodicity.The two sources with the longest reported periods, Seq.6 and Seq.7, have their FAP much greater than 1%, which suggests a substantial false alarm probability due to red noise.However, this does not necessarily exclude their possibility of being true periodic variations, since the GL algorithm can still uncover intrinsic periodic signals in the presence of substantial red noise.To confirm that Seq.6 and Seq.7 are genuine periodic signals, we select for both sources eight adjacent and long observations (ObsID 14197, 14198, 13825, 13826, 13827, 13828, 14195, 14196) and align the light curves from each individual observation according to the phase determined by the identified period.The aligned light curves are shown in Figure 5, which clearly reveal an eclipsing behavior in both sources.Furthermore, the consistency in their eclipse phases across these observations serves as compelling evidence for the accuracy of the identified period.Therefore, we conclude that these two periodic signals are genuine, despite their large FAP inferred from the simulated red noise-based light curves.

X-RAY SPECTRAL ANALYSIS
It is instructive to investigate the X-ray spectral properties of the periodic sources, which may provide insights into their physical nature.For this purpose, we focus on the ACIS data, but discard the HRC data due to its poor energy resolution.However, three periodic sources, Seq.2, Seq.3 and Seq.5, are previously known to be novae (see Section 5 for details).Their detection by Chandra occurs solely during the outburst phase, which were all through HRC observations (Table 2), and all three sources fall below the detection limit in the Figure 4.The distribution of  GL from 1000 simulated light curves based on the modelled red noise of each periodic source (coloured histograms).The red dashed line marks  GL = 0.90, the adopted threshold for a periodic signal to be identified.
ACIS observations, precluding the availability of a meaningful spectrum.Hence, we estimate their peak 0.5-8 keV luminosity during outburst from their net photon flux, by assuming an absorbed blackbody spectrum with a column density of  H = 1.9 × 10 21 cm −2 and a black-body temperature of  = 40 eV, taken from the spectral model of Seq.2 determined by XMM-Newton observations (Henze et al. 2014).For Seq.1, Seq.4,Seq.6 and Seq.7, we examine the source spectrum combining all available ACIS observations.For Seq.1, we further examine the spectrum combining only the three ACIS observations over which the periodic signal was detected.We extract the source and background spectra for each periodic source by utilizing the CIAO tool specextract from the aperture defined in Section 2. For each source, we then generate the combined spectra by CIAO tool combine_spectra, with the corresponding ARFs and RMFs weighted by the effective exposure.The spectra are analyzed using XSPEC v12.12.1, after adaptively binned over 0.5-8 keV such that a minimum of 20 counts and a S/N greater than 2 per bin are achieved.
To fit the spectra, we typically employ an absorbed power-law model, and when necessary, introduce an additional soft component, either diskbb or bbody in XSPEC.The best-fit model parameters are given in Table 3.The unabsorbed 0.5-8 keV luminosity is estimated according to the best-fit model and an assumed distance of 785 kpc for M31 (McConnachie et al. 2005).

CLASSIFYING THE PERIODIC X-RAY SOURCES
Based on the timing and spectral properties of the periodic sources, as well as information from the literature, we attempt to classify them in order below.We note in passing that, while a non-negligible fraction of all X-ray sources detected in the direction of the M31 bulge arise from the CXB (∼0.1 per arcmin −2 ; Voss & Gilfanov 2007b), our recent study based on the Chandra Deep Field South (Bao & Li 2022) suggests that the CXB has a negligible contribution to the periodic signals found here.
Seq.1:This source is the most luminous one among the seven periodic sources, naturally making it an active LMXB.It has been suggested that this source is a BH-XRB candidate (Barnard et al. 2013), based on an argument that its X-ray spectrum at the high hard state is distinctively softer than that of NS-XRBs.We detected its periodic signal, which exhibits an eclipsing or dipping behavior during its high state between 2010-2011.The presence of eclipsing/dipping signatures strongly suggests that the observed signal originates from the binary orbital motion.That this signal is only present in a small fraction of observations is somewhat puzzling, especially in the scope of an orbital modulation.Nevertheless, this signal should be genuine, given the fact that it is detected in both HRC and ASIC observations, disfavoring an instrumental artifact.We speculate that this eclipse/dipping is due to some absorbing structure in the accretion disk, e.g., an accretion disk corona (ADC, White & Holt 1982;Chou 2014).In this scenario, the central source of the XRB system evaporates material from the surface of its accretion disk.If the evaporated material does not escape from the system, it will corotate above and below the disk as a corona and can partly absorb or scatter the Xrays.This may happen even when the viewing angle of the binary system is relatively small.Such an ADC is likely a transient structure depending on the level of the central illuminating source and could be consistent with what is observed in Seq.1.However, we find no significant evidence for an enhanced absorption in the X-ray spectrum of Seq.1 at the epochs when the periodic signal was detected (Figure 2 and Table 3).
Regardless of the nature of the transient dipping, the short period of ∼463 sec, if truly an orbital period, would make Seq.1 a candidate of the shortest-period ultra-compact X-ray binary (UCXB) observed so far.Such an exceptionally tight orbit also indicates that Seq.1 has the potential of becoming an extragalactic gravitational wave source for future detectors.According to He et al. (2023), at least 0.8 merging BH-WD binaries and ∼0.1-0.4 merging NS-WD binaries are expected to be detected in M31 during a 4-yr LISA mission.The gravitational wave of a NS-WD or BH-WD system with  orb = 463.24sec is about 10 3 times below the LISA sensitivity (Babak et al. 2021).Next generation gravitational wave detectors are eagerly awaited to detect extragalactic UCXBs.
Seq.2, Seq.3 and Seq.5:All three sources have been previously confirmed as novae due to their intense outbursts.During the nova outburst, the X-ray luminosity may rise to ≳ 10 36 erg s −1 (Table 3), temporarily offering an X-ray light curve sufficient for probing the periodicity.
Seq.2, also known as M31N 2011-01b, displayed an optical counterpart identified as PNV J00423907+4113258 (Henze et al. 2014).In ObsIDs 13179 and 13180, when the nova burst occurred, a periodic signal with a period of ∼4686 sec is detected.Seq.3 (M31N 2011-11e) with an optical counterpart PNV J00423831+4116313 exhibited a periodic variation of  ∼ 1.3 ± 0.1 h during its outburst (Henze et al. 2014).By utilizing the same HRC observations in 2012 (ObsID 13278,13279,13280), we detect a nearly identical period of 4933.40 sec.Seq.5 was identified as the nova M31N 2004M31N -11b (Henze et al. 2010)).We detect a periodic signal at ∼ 7184 sec during an outburst in 2005 caught by ObsID 5927 and 5928.Judging from the phase-folded light curves (Figure 3), an eclipsing/dipping behavior is seen in at least two sources, Seq.3 and Seq.5, suggesting an orbital modulation.Notably, these periods fall below the so-called orbital period gap (2-3 h) of CVs (Knigge et al. 2011).This might be partly due to a selection effect given the relatively short exposures of the individual Chandra observations.
Seq.4:This source has a period of ∼ 6432.94sec, with an eclipsing behavior strongly suggesting an orbital modulation of a binary system.An orbital period of 6420 sec was previously discovered by XMM-Newton observations in 2000-2002 and Chandra observations in 2001 (Mangano et al. 2004).Here we detect this period  (ObsID 14197, 14198, 13825, 13826, 13827, 13828, 14195, 14196) for Seq.6 (left panels) and Seq.7 (right panels).The top eight rows show the light curves of each individual observation, for which the ObsID and start time of the light light curve are labeled.The botton row overlays the eight light curves, highlighting the eclipsing behavior in both sources.Two eclipses are evident in Seq.7.consistently in all ACIS and HRC observations from 1999 to 2015.Lazzarini et al. (2018) argued that this source is an HMXB, based on the tentative identification of an optical counterpart exhibiting a B-type star multi-band spectral energy distribution.However, given the short orbital period identified here, we argue that this source is more likely an LMXB.In this regard, the previously suggested optical counterpart was probably a mismatch.The X-ray spectrum of this source is best-fitted with a power-law plus a blackbody, consistent with an LMXB.Seq.6:This source exhibits a period of ∼ 22835 sec in both the ACIS and HRC observations.The eclipsing behavior observed in its light curve strongly supports the notion that the period is from the binary orbital modulation.In both the ACIS and HRC phase-folded light curves, there appears a small spike (at phase ∼ 1.0) in the otherwise broad eclipse (Figure 2), which starts around phase ∼ 0.75 and ends at phase ∼ 0.2.However, such a feature is not clearly evident in the unfolded light curves shown in Figure 5, which casts doubt to its reality.The cumulative X-ray spectrum is well described by an absorbed power-law model, with Γ ≈ 1.85 and  X ≈ 6×10 36 erg s −1 .The combination of its periodicity and spectral characteristics is well consistent with an LMXB.
Seq.7:This source exhibits a highly distinctive light curve characterized by two eclipses, a wider one (phase between ∼ 0.4-0.6)followed by a narrow one (centering at phase ∼ 0.7), consistently observed in the ACIS and HRC observations.This dual-eclipse feature is most likely real, because it is directly seen in the unfolded light curves (Figure 5), and the narrower eclipse is almost a full eclipse.The cumulative X-ray spectrum can be well described by a simple power-law, with Γ ≈ 1.72 and  X ≈ 7 × 10 36 erg s −1 .The projected position of this source is closest to the center of M31, raising the possibility that it is an LMXB formed in dynamical encounters in the crowed nuclear region (a projected radius ≲ 1 ′ ; Voss & Gilfanov 2007a).
To our knowledge, a dual-eclipse feature is rarely seen in LMXBs.GRS 1747-312, an LMXB residing in the Galactic globular cluster Terzan 6, is known to exhibit eclipses/dips at different orbital phases (Painter et al. 2024), but these peculiar features are not persistent, unlike the case of Seq.7, in which the dual-eclipse remains highly stable throughout the Chandra observations spanning 16 years.We propose two possible explanations for the behavior of Seq.7.The first explanation involves absorption by an ADC and a nearly edgeon viewing angle of the orbit.Two eclipses will appear in the light curve:a narrow eclipse by the donor star and a wide eclipse by the ADC.A couple of examples of this kind have been discovered in Galactic XRBs (e.g., 4U2129+47, 4U1822-37, White & Holt 1982).The second explanation involves a three-body system, which was also suggested for the case of GRS 1747-312 (Painter et al. 2024).In this scenario, the XRB system consists of three objects: a compact object (NS or BH), a donor star, and a third object (a normal star or a degenerate star), likely sharing the same orbit of the donor but with a quasi-stable phase difference of 1/6 orbit (similar to what is observed in Seq.7).The donor star transfers mass to the compact object, producing the bulk of the X-rays, which is obscured by the donor at a nearly edge-on viewing angle.The presence of the third object can introduce additional obscuration, either by itself or by a mini accretion disk around it.A more quantitative test of this picture is beyond the scope of the present study.

Comparison with previous works
We have detected seven periodic signals from seven X-ray sources in the M31 bulge, among which four are newly discovered.On the other hand, several periodic X-ray sources identified in previous studies, as outlined in Section 1, are not recovered in this work, for various reasons.
XMMU J004252.5+411540 has been identified as the brightest persistent SSS in M31.The XMM-Newton observations revealed a periodic modulation in it with a period of 217.7 sec.Our analysis using the Chandra observations indicates that only one HRC observation displayed similar periodic variation, but with a low GL probability ( GL ∼ 0.5, ObsID=9829).This discrepancy might be attributed to the lower effective area of Chandra compared to XMM-Newton, particularly in the energy range of 0.2-0.5 keV over which the SSS is primarily detected.XMMU J004319.4+411758 is also an SSS with a period of  ∼ 865 sec detected by XMM-Newton.
No point source is detected at the same position by the Chandra observations, probably due to its super-soft spectrum and the transient nature.M31N 2006-04a and M31N 2007-12b are both novae, with their period and outburst detected by Henze et al. (2010); Pietsch et al. (2011).We did not detect their periodicity as they were not in an outburst state during all Chandra observations.XMMU J004314.1+410724 is identified as an LMXB within the globular cluster Bo158, with a period of 2.78 hr detected by Trudolyubov et al. (2002).However, it falls near the edge of the field-of-view of the Chandra observations with a poor PSF, which is not included in our period research.3XMM J004232.1+411314 has been detected with periodic dipping signals of 4.02 hr only during a low-luminosity state ( 0.2−12 keV < 10 38 erg s −1 ) by Marelli et al. (2017).We detected the same period from Chandra data only in one long observation during a low-luminosity state, but with a low GL probability ( GL ∼ 0.72, ObsID=14196).3XMM J004301.4+413017 has been designated as an NS-XRB owing to the presence of a pulsation of 1.2 sec, a timescale shorter than the ACIS frame time.Furthermore, its orbital period of 1.2 days lies beyond our detection limit.
To summarize, in the bulge of M31, there have been eight periodic X-ray signals reported in the literature, which are not recovered by our systematic search with the Chandra data, due primarily to a mismatch of the data capability.These include 2 SSSs, 2 novae and 3 LMXBs (XMMU J004314.1+410724,3XMM J004232.1+411314,3XMM J004301.4+413017).

Implications on the LMXB population of the M31 bulge
As summarized in the previous sections, a total of seven LMXBs with a periodic X-ray signal are found in the M31 bulge.Most, if not all, of these periodic signals can be associated with the binary orbital motion.Together, these sources may provide a useful diagnosis to the underlying close binary populations in the M31 bulge, which in turn allows for a meaningful comparison with the Galactic populations.
In Figure 6, we plot the X-ray (2-10 keV) luminosity versus the putative orbital period of the seven LMXBs of the M31 bulge.The 2-10 keV luminosity of the four Chandra-detected sources are derived from the best-fit spectral model in Table 3, whereas the values of the XMM-Newton-detected sources are taken from the literature (Trudolyubov et al. 2002;Marelli et al. 2017;Zolotukhin et al. 2017).For comparison, Figure 6 also plots Galactic LMXBs with a known orbital period.In particular, we extract from the Galactic LMXB catalog compiled by Avakyan et al. (2023) a total of 64 NS-LMXB candidates and 31 BH-LMXB candidtaes with an orbital period ≳ 80 minutes.This sample is supplied by Armas Padilla et al. (2023), which have classified 21 Galactic UCXBs characterized by their orbital periods shorter than 80 minutes.These sources are shown by different colored symbols in Figure 6.
Certainly, not all intrinsic X-ray periodic signals in the M31 bulge have been identified by the GL algorithm.Following Bao & Li (2020), we evaluate the detection rate using a large set of synthetic time series.Specifically, we adopt a sinusoidal light curve with a 30% variation amplitude as the fiducial input of the periodic variation.Taking a stepwise light curve, which might better mimic an eclipsing behavior but requires an additional free parameter, gives a similar result according to Bao & Li (2020).We try three sets of representative periods (793 s, 7093 s, 17093 s) and six sets of total counts (250,500,1000,2000,5000,10000).A constant light curve with six sets of total counts is also simulated as the no-signal control group.The synthetic light curves, taking into account the Poisson fluctuation and the observing cadence and length same as the true Chandra observations, are then fed to the GL algorithm.For each parameter set, 5000 light curves are generated in order to minimize random fluctuations.A reported period with  GL > 0.9 and within few percents of the input period is considered a valid detection.
The detection rate for a given combination of parameters is taken to be the fraction of the synthetic light curves having a valid detection, as shown in Figure 7.The detection rate of the no-signal control group is nearly zero, showing that false alarm from a constant light curve is negligible, regardless of the source counts.For the periodic light curves, a general trend revealed is that the detection rate increases as the total counts increase, which is naturally understood.Moreover, for the same number of total counts, a higher detection rate is found for a larger period, which is due to the nature of the GL algorithm (Bao & Li 2020).The detection rate for source counts of a few 10 3 (i.e., the case of the M31 sources, Table 2), is close to 100% for all three periods.We point out that for long-period signals, despite the simulation results indicating a detection rate up to 100%, the true signal might be accompanied by numerous false signals due to noise, as illustrated by the simulations presented in Figure 4, a practical effect that would significantly hamper the identification of the true period.
Based on the predicted detection rate and the distribution of source counts among the X-ray sources in the M31 bulge, we can conclude that essentially all sources with a total count greater than 10 3 and an intrinsic periodic variation at a period lower than about 10 4 sec have been identified, whereas a significant fraction of periodic signals from fainter sources might have been hidden in the noise.Combining the estimated detection rate, the observed distribution of source counts, and further assuming 30% of LMXBs having an orbital period below 10 4 sec (i.e., same as the Milky Way population; Figure 6), we expect to have ∼ 45 short-period (< 10 4 sec) LMXBs in the M31 bulge, which can be identified by the GL algorithm.A further factor to consider is the intrinsic fraction of orbital modulation, which 10 4 should occur only when the inclination angle of the orbital plane is high.According to Frank et al. (1987), an inclination angle ranging between 60 • <  < 80 • would lead to eclipsing or dipping signals, which translates to a fraction of ∼ 20% of all LMXBs for a randomly distributed inclination angle.Hence we expect ∼ 9 detectable periodic LMXBs with a period shorter than 10 4 sec, which is several times higher than our actual detections (2 sources).Such a difference implies that the fraction of short-period LMXBs ( ≲ 10 4 sec) in the M31 bulge may be significantly lower than that of the Milky Way LMXB populations, although we caution that selection bias could exist in both the M31 bulge and Galactic LMXB samples.

SUMMARY
We have conducted a systematic search for periodic X-ray signals from several hundreds of X-ray sources in the M31 bulge observed by extensive Chandra observations in a total of 2 Ms spanning a temporal baseline of 16 yrs.Thanks to the GL algorithm that works well for phased-folded X-ray light curves, seven periodic signals are confidently identified, among which four are newly discovered.Three of these sources are novae, the identified periods of which range between 1.3-2.0h and is most likely the orbital period.The other four sources are LMXBs, the identified period of which is also likely orbital due to a clear eclipsing/dipping behavior in the light curve.One of them might be classified as an UCXB due to the very short (463 sec) period.This is the first systematic search for periodic X-ray sources in an external galaxy.Our study demonstrates the potential of using archival X-ray observations to systematically identify periodic X-ray sources in external galaxies, especially LMXBs with short orbital periods, which would provide valuable information about the underlying exotic stellar populations.

Figure 1 .
Figure 1.Counts image of the M31 bulge, combining 64 Chandra HRC observations.The positions of the seven periodic sources are indicated by yellow circles, with source numbering as in Table2.

Figure 3 .
Figure 3. Phased-folded light curves for Seq.2, Seq.3 and Seq.5, whose periodic signals are only detected in a few HRC observations during their outburst.

Figure 5 .
Figure 5. Phase-aligned light curves of 8 long observations(ObsID 14197, 14198, 13825, 13826, 13827, 13828, 14195, 14196)  for Seq.6 (left panels) and Seq.7 (right panels).The top eight rows show the light curves of each individual observation, for which the ObsID and start time of the light light curve are labeled.The botton row overlays the eight light curves, highlighting the eclipsing behavior in both sources.Two eclipses are evident in Seq.7.
Source sequence number.(2) The best-fit model in XSPEC.(3)  2 and degree of freedom of the best-fit model (4) 0.5-8 keV unabsorbed luminosity.*: Spectrum combining only three ACIS observations over which the periodic signal was detected.†: Peak luminosity during the outburst, which is converted from the observed photon flux in the HRC observations.(5) Line-of-sight absorption column density.(6) Photon index of the power-law model.(7) Temperature of diskbb or bbody model.Quoted errors are at the 90% confidence level.

Figure A1 .
Figure A1.Left: the Gregory-Loredo diagram of odds ratio for the seven periodic sources.The -axis is the frequency (/2 ).The -axis is the odds ratio for an optimal number of  bins.Right: the normalized Lomb-Scargle periodogram for the seven periodic sources.Three false alarm levels are denoted by the colored dashed lines.

Table 1 .
Known periodic X-ray sources in M31.Super orbital period.Source located in the globular cluster Bo 158.b Pulsation period of the neutron star.c Pulsation period of the accreting pulsar.Source located in the globular cluster GLC 377.
Padilla et al. (2023)minosity versus orbital period for the seven M31 bulge LMXBs (red circles, with an additional 'x' marker for the four sources discovered in this work).Galactic UCXBs, BH-LMXBs, NS-LMXBs fromAvakyan et al. (2023)and ArmasPadilla et al. (2023)are shown by orange squares, green triangles and blue diamonds, respectively.The vertical dashed line marks the upper bound of the orbital period of UCXBs.
Detection rate as a function of total counts, based on simulated light curves with a sinusoidal variation fed to the GL algorithm.Six values of total counts(250, 500, 1000, 2000, 5000, 10000)and three values of periods (793 sec, 7093 sec, 17093 sec) are simulated.The model with a constant count rate (i.e., no periodic signal) is also plotted.