Monitoring the X-Ray Variability of Bright X-Ray Sources in M33

We present a new five-epoch Chandra X-ray Observatory monitoring survey of the nearby spiral galaxy M33 which probes X-ray variability with time sampling between two weeks and four months. We characterize the X-ray variability of 55 bright point sources outside of the nucleus, many of which are expected to be high-mass X-ray binaries (HMXBs). We detect eight new candidate transients not detected in previous X-ray catalogs of M33 and discuss their possible nature. The final catalog includes 26 known HMXB candidates identified in the literature. We extend the baseline of the X-ray light curves up to 21 yr by including archival X-ray observations of these sources. We compare the detection and nondetection epochs of the sources to suites of simulated source duty cycles and infer that most of our detected sources have duty cycles >30%. We find only four sources whose detection patterns are consistent with having duty cycles below 30%. This large fraction of sources with high duty cycles is unexpected for a population of HMXBs; thus more frequent X-ray monitoring will likely reveal many more low duty cycle HMXBs in M33.


INTRODUCTION
In an X-ray binary (XRB), a compact object accretes material from a stellar companion.Reaching this state requires a complex path of binary, and sometimes triple, star evolution (e.g., Ford et al. 2000;Sana et al. 2012;Duchêne & Kraus 2013;De Marco & Izzard 2017).While our own Galaxy and its satellites contain many XRBs whose physical parameters have been measured in detail (Liu et al. 2007;Antoniou & Zezas 2016;Haberl & Sturm 2016;Tetarenko et al. 2016a), statistical studies of the Milky Way population can be difficult due to our location in the dusty Galactic disk.By extending our studies of these sources to other galaxies where the distance to all XRBs is the same, we expand the sample to other environments and create homogeneous data sets that simplify statistical studies (e.g., Prestwich et al. 2003;Tzanavaris et al. 2013;Haberl & Sturm 2016;Lehmer et al. 2021).
X-ray variability is an important first-order indicator of the accretion mechanism driving HMXB emission.Multiple mass transfer mechanisms are possible depending on the configuration of the system, including stellar winds from the secondary star (e.g., Hirai & Mandel 2021), Roche-lobe overflow (e.g., Chen et al. 2017), passage of the compact primary through a Be-type companion's circumstellar disk (e.g., Coe et al. 2015a), and wind-capture Roche-lobe overflow which may be responsible for powering the most luminous HMXBs (Huarte-Espinosa et al. 2013;El Mellah et al. 2019).Each case results in powerful X-ray emission originating from the compact object.Changes in an HMXB's total X-ray emission are caused by changes in the accretion rate of material onto the compact object from the stellar companion.Several phenomenological X-ray source states for XRBs with a black hole primary describing differences in outburst behavior have been defined (Homan & Belloni 2005;Tomsick 2006; Kalemci et al. 2022).These source states are described by luminosity as well as spectral and variability components, and transitions between them are linked to changes in the accretion mechanism (e.g., Remillard & McClintock 2006).Although our survey is not sensitive enough to obtain spectra of these sources, changes in luminosity alone can indicate state changes.Variability studies over a wide range of timescales are critical to understanding these processes.
One important observational constraint on such variability is the duty cycle (DC), or ratio of time in outburst compared to the object's lifetime.The presence of highly variable sources with some percentage of time in a low X-ray luminosity state make it difficult to empirically determine the total number of HMXBs present in a galaxy (Ducci et al. 2014;Haberl & Sturm 2016;Binder et al. 2017;Lazzarini et al. 2021;Mori et al. 2021).The fraction of a galaxy's HMXBs that always shows high emission (DC = 100%) as a result of persistent accretion is uncertain.Collectively, measuring HMXB DCs is essential to our interpretation of the X-ray luminosity function (XLF) (Belczynski et al. 2004;Paizis & Sidoli 2014;Binder et al. 2017;Lehmer et al. 2019).Within the last decade, the number of HMXB candidates and long-term X-ray observations has sufficiently grown to allow initial characterization of XRB DCs (Binder et al. 2015;Williams et al. 2015;Garofali et al. 2018;Lazzarini et al. 2021).
An additional nearby galaxy with a significant population of HMXBs is M33.At a distance of 817 kpc (Freedman et al. 2001), M33 is the second closest star-forming spiral galaxy.It lies at medium inclination (i = 52 • ; Kam et al. 2015) and has been well-observed across the electromagnetic spectrum (radio: Israel & van der Kruit (1974); Tabatabaei et al. (2022); optical/UV: Williams et al. (2021); γ-ray: Xi et al. (2020)).Previous X-ray surveys have localized point sources with high accuracy using the Chandra X-ray Observatory (Tüllmann et al. 2011, T11 hereafter), as well as with wide area coverage and several epochs observed by XMM-Newton (Misanovic et al. 2006;Williams et al. 2015, M06 and W15 hereafter).M33 has a high star formation rate, making its X-ray binary (XRB) population dominated by HMXBs (e.g.W15) which have secondary stars that are easily detected in HST imaging data (Garofali et al. 2018;Williams et al. 2021, G18 hereafter), and some with NuSTAR detections that allow a tentative classification of the compact object (Yang et al. 2022).These characteristics make M33 an excellent target for studies of XRB variability.
There have been five major surveys of M33 in the X-ray regime, three of which we utilize to study long-term variability of the sources in M33 between 2000-2015.M33 was first observed with the Chandra X-ray Observatory by Grimm et al. (2005), and Grimm et al. (2007) used this survey to both characterize the population's X-ray variability and conduct spectral fitting.M06 conducted a deep XMM-Newton survey from 2000-2003.T11 conducted a deep Chandra survey from late 2005 to late 2006.Williams et al. (2008) analyzed the Plucinsky et al. (2008) first look catalog of the T11 survey to identify seven X-ray transient candidates.W15 revisited M33 with XMM-Newton between 2010-2012.Yang et al. (2022) have also completed a NuSTAR survey of M33 in the hard X-ray band, classifying 28 sources as X-ray pulsars and BH binaries in various source states.
Pairing these surveys with optical data from the Hubble Space Telescope, G18 predicted ∼ 109 HMXBs in M33 assuming a star formation rate of 0.3 M ⊙ yr −1 (Williams et al. 2013) and identified 55 HMXB candidates from archival HST fields that overlapped with the T11 survey.Mass and age estimates for 52 of these candidates have been obtained by G18 by analyzing their host stellar populations.The Panchromatic Hubble Andromeda Treasury: Triangulum Extended Region survey (PHATTER; Williams et al. (2021)) uniformly covers the inner disk of M33, including some regions not covered by the archival HST imaging used by G18.Lazzarini et al. (2023) (L23 hereafter) identified 65 HMXB candidates by comparing the T11 X-ray positions to PHATTER optical counterparts consistent with being massive stars.L23 classified the companion spectral types and fit their spectral energy distributions to measure their physical properties.
In addition to HMXBs, a shallow time-domain X-ray survey of M33 would be sensitive to strong outbursts from low-mass X-ray binaries (LMXBs).Accretion in a large fraction of these systems is attributed to the companion star filling its Roche-lobe (see Bahramian & Degenaar (2023) for a recent overview).This accretion mechanism is highly efficient at converting energy into radiation and peak LMXB X-ray luminosities during outbursts are observed to range from ∼ 10 37 erg s −1 to a few times 10 38 erg s −1 (Yan & Yu 2015) on timescales of days to months.These populations have been studied extensively in the bulge of M31 (e.g., Williams et al. 2006;Peacock et al. 2010;Barnard et al. 2014).While M33's XLF is dominated by HMXBs due to its younger stellar population (Yang et al. 2022), it is likely that some number of LMXBs are also present in the galaxy.Our survey would be sensitive to outbursts of these transient sources.
In this paper we present new Chandra observations of M33 taken over five epochs in 2021, nine years since W15 obtained the last X-ray observations of the galaxy.These new data allow us to characterize the long-term X-ray variability of 55 bright sources, including 26 candidate HMXBs identified by T11, G18, Yang et al. (2022) and L23.In Section 2, we describe the observations and data reduction.In Section 3, we present our catalog.In Section 4, we discuss constraints on the source DCs and how they compare to HMXB populations in other galaxies.

OBSERVATIONS & DATA REDUCTION
We obtained ten 5 ks observations of M33 with Chandra's ACIS-I detector.Five pointings were centered on the Northern half of M33 and five pointings were centered on the Southern half.Table 1 gives the details of each observation.These ten exposures cover ∼66% of the B-band D 25 isophote of M33 with observations taken between November 2020 and June 2021.The observations were designed to probe variability on multiple timescales.The five epochs in each pointing are separated by two weeks, four weeks, two months and four months consecutively.The time sampling is shown in Figure 1.We reprocessed all observations using CIAO v4.11 and CALDB v4.9.4 (Fruscione et al. 2006).Below we describe our technique to measure our point source catalog and compare our catalog to previous surveys.

Source Detection and Point Source Extraction
Source detection was performed with 0.492 arcsec and 0.984 arcsec pixel sizes for completeness before pruning the sources based on the extracted properties.We used the CIAO tool fluximage to create aspect histograms and point spread function files for each observation with both pixel sizes.We conducted source detection on each images with the CIAO tool wavdetect on both pixel sizes and size scales ranging from 1 to 16.We produced separate source location lists for each of the two pixel sizes and merged them to create our candidate source position list.We examined the positions reported for each Observation ID (ObsID) image to ensure that all potential sources were included in the wavdetect lists.We aligned all observations by matching the three brightest source detections in each observation with their literature counterparts in T11 using wcs match.We then used the standard recipes for ACIS-Extract v2021feb4 (AE hereafter; Broos et al. 2010Broos et al. , 2012;;Broos & Townsley 2022) to extract the counts from the locations identified by wavdetect in four energy bands: soft (0.35-1.1 keV), medium (1.1-2.6 keV), hard (2.6-8.0 keV) and full (0.35-8.0 keV), chosen to match those in the T11 catalog.AE extracts source properties from each ObsID and merges them into a single multi-ObsID extraction with multiple improved position estimates.We visually inspected all position estimates produced for the data in SAOImageDS9 (DS9 hereafter; Joye & Mandel 2003) and chose the most accurate estimate to be the new source position.We ran AE a second time with the updated positions using the 0.492 arcsec pixel size for the final point source extractions.
We made final catalog selections based on the PROB NO SOURCE statistic (PNS) calculated by AE, which describes the probability of the extracted region not being a source.Figure 2 shows that there were many sources in the initial list with high PNS in the multi-ObsID extraction, with a gap at about 10 −10 .We visually inspected the event lists of sources with PNS > 10 −5 in DS9 and compared the source lists based on cuts between PNS > 10 −8 and PNS < 10 −3 .Our final catalog cut required PNS < 10 −4 (corresponding to ≳ 2 × 10 −4 ct s −1 ) in at least one energy band in either the multi-or single-ObsID extraction, to include sources with high variability.The final source list is insensitive to change down to PNS > 5 × 10 −5 and includes all position matches to T11 (described in Section 4).For observations of some sources, AE reports zero or negative extracted counts, indicating that there was no detection during those particular observations.
To estimate the sensitivity of our survey we created sensitivity maps for each ObsID (see Figure 3) using the exposure corrected background images produced by wavdetect with a cutoff signal to noise ratio (S/N) > 2. The average of the median sensitivities in all ten observations is 8 × 10 −4 ct s −1 .Of the final catalog sources, 76.8% have count rates larger than the average median sensitivity in the full 0.35-8.0keV band, and all source count rates exceed this sensitivity when considering their upper errors.The S/N cutoff > 2 aligns with sources that pass the multi-ObsID PNS cut as seen in Figure 4, and no sources with PNS > 10 −4 have S/N > 2. 41 sources have a S/N > 2 while 56 sources have PNS < 10 −4 due to using a flat S/N-based sensitivity which does not account for lower sensitivity at high off-axis angles.
The nucleus of M33 contains the known ultraluminous X-ray source (ULX) M33 X-8 (West et al. 2018;Krivonos et al. 2018), which is known not be a transient source.We include this source in our tables for completeness but in this work we focus on the non-ULX HMXB population.We have excluded the nuclear source from all figures.Our observations may be of interest for future dedicated studies of M33 X-8.

Matching to Previous X-ray Surveys
The rich observational history of M33 allows us to compare our observations to multiple archival X-ray surveys of the galaxy to study the variability of sources over a longer baseline, as well as to calibrate our astrometry.We cross matched our catalog to source positions in M06, T11 and W15 using the catalog matching tool NWAY (Salvato et al. 2018).We used a 10 arcsec search radius for completeness of matches to faint off-axis sources in our catalog with high positional uncertainties.38 sources in the catalog have matches to M06, with a median offset of 1.3 arcsec and maximum offset of 6.6 arcsec.46 sources in the catalog have matches to T11 with a median offset of 0.5 arcsec and maximum offset of 3.6 arcsec.45 sources in our catalog have matches to W15 with a median offset of 1.0 arcsec and maximum offset of 5.3 arcsec.Sources with archival counterparts are indicated in Table 3 and the catalog source identifiers are given in the full table online.We discuss eight sources which do not have matches to any archival sources as new transients in Section 3.2.
The soft-band sensitivity of the ACIS instrument on Chandra has been found to decrease over the observatory's lifetime (Plucinsky et al. 2020).AE reports flux estimates that take the instrument's response files into account (which are available in the full catalog table online), but for this shallow monitoring study we are primarily interested in whether or not a source is detected.To simplify comparison to archival observations, we found standard conversion factors between archival fluxes and modern Chandra count rates.We converted the archival fluxes of T11 to modern-day Chandra sensitivity using the web interface of the Portable, Interactive Multi-Mission Simulator (WebPIMMS; Mukai (1993)) by assuming a power law spectrum with photon index 1.7 and Galactic N H of 6.0×10 20 cm −2 , which are the same assumptions used by M06 and W15 for conversions between count rates and fluxes.We found a conversion factor of 0.608 between count rates observed in Chandra Cycle 7 when T11 observations were taken and Cycle 22 when our data was obtained.The fluxes measured by the two XMM-Newton surveys must also be converted to their equivalent Chandra count rate due to differences in the sensitivity of the two observatories.We again used WebPIMMS and found a conversion factor of 1 equivalent Chandra count s −1 = (6.1×10 10 )×F XM M for the two XMM-Newton surveys, assuming the same spectral model.
We calculated the median residual count rate between the archival surveys and our data for sources with PNS < 10 −10 to estimate the uncertainty in our calibration.The median residual compared to T11 is -21%.The median power law index for the 38 T11 sources with best fit power laws is Γ = 1.65 ± 0.17, corresponding to a range of conversion factors between 0.589 to 0.646.About 20% of the cross matched sources to T11 have soft power law spectral fits with Γ > 1.83, which accounts for some of the overall difference in calibration between our survey and T11.The comparison between our measurements and T11 is similar to typical comparisons between X-ray data sets.For example, the median difference from W15 is 10.3%, and for M06 it is 16%.With these data, we provide a final catalog of 55 sources plus the nucleus, discuss eight new transients detected in this survey, and present X-ray light curves with our detections and archival observations.The catalog positions and survey coverage area are shown in Figure 5.
We report a summary of the catalog in Table 3 with source positions and archival counterpart identifiers.Table 4 gives source likelihood statistics (PNS, S/N) and count rates for the merged and Right Ascension (J2000.0) Right Ascension (J2000.0)Declination(J2000.0) 1h34m00.0s30.0s 33m00.0s32m30.0s30.0s 35m00.0s30.0s 1h34m00.0s30.0s 33m00.0s32m30.0s30.0s 35m00.0s30.0sTable 2. Sensitivities of archival surveys compared to in this paper.Limiting unabsorbed luminosity in each respective band converted to flux assuming a distance of 817 kpc to M33 and to modern Chandra count rate using the same assumptions described in Section 2.2.The average single-ObsID sensitivity of our survey is given in the last row.We describe the calculation of this sensitivity in Section 3.

Survey
Limiting L (erg s −1 ) Energy band (keV) Limiting F (erg cm −2 s −1 ) Chandra count rate (ct s individual ObsID extractions in all four energy bands.Errors on the extracted net counts are given at the 1σ confidence level using the Gehrels approximation for Poisson statistics (Gehrels 1986) as implemented in AE.The full table online includes all extraction results reported by AE.Our final catalog includes 55 point sources, eight of which do not appear in previous X-ray surveys and are discussed in Section 3.2.The limiting luminosity of our survey is determined by the lowest positive count rate in a single observation among all point sources in the catalog.The limiting count rate is 4 × 10 −4 ct s −1 .Assuming a power law spectrum with photon index 1.7, and Galactic N H of 6.0 × 10 20 cm −2 , we used WebPIMMS to find a conversion of 1 ct s −1 = 2.11 × 10 −11 erg cm −2 s −1 , yielding a limiting unabsorbed flux in the full energy band of 8.32 × 10 −15 erg cm −2 s −1 .Assuming a distance of 817 kpc to M33, the limiting luminosity of our survey is 6.6 × 10 35 erg s −1 .For sources with soft X-Table 3. Preview of the final source catalog.The full catalog online contains 56 sources.Theta describes the average off-axis angle of the source on the detector.T11 Source Numbers refer to the T11 Final Catalog designations, not to be confused with the Plucinsky et al. (2008) first look catalog source numbers.Positional uncertainties reported by AE as less than 0.5 arcsec were set to 0.5 arcsec to match the minimum astrometric uncertainty of the T11 survey which we used for alignment.ray spectra the limiting luminosity is higher due to Chandra's decrease in soft-band sensitivity, as discussed above.

X-ray Light Curves
We show two examples of X-ray light curves for our catalog sources in Figures 6 and 7. Light curves for all of the sources are included as supplemental data to the article.Long-baseline light curves incorporating the archival detections and upper limits from M06, T11, and W15 and our merged count rate are shown in the left panels.The duration of our monitoring campaign is indicated by the shaded region.Our survey's sky coverage lies within the exposure area of each of the archival catalogs, so in cases where our sources do not have archival counterparts we take the limiting sensitivity of each survey as the upper limit on the source's flux at that time.The archival upper limits used are given in Table 2.
Short-term light curves using our own observations alone are given in the right panels.Observations from Northern and Southern pointings are marked by different colors.Nondetections when the source position was covered by the observation are represented in the light curve figures as downward arrows to zero ct s −1 .We find the upper limits for these points by adding the upper error to the extracted count rate.The extracted count rates are unchanged in the reported photometry online.

New Candidate Transients
Galactic transient XRBs in quiescent states have X-ray luminosities less than 10 33 erg s −1 and reach values from 10 36 to 10 39 erg s −1 in outburst (Degenaar et al. 2012), most of which peak around 5 × 10 36 erg s −1 (Binder et al. 2017).Given our detection limit of 6.6 × 10 35 erg s −1 , we expect to detect XRB sources in outburst and high-luminosity active accretion states only.
Within our catalog, we find eight sources to be new X-ray detections with no counterparts in archival X-ray catalogs of M33.Each new candidate transient has a strong detection in at least one ObsID as summarized in Table 5, and a few have positive extracted NET CNTS in a second ObsID that are close to our detection limit.We calculated minimum flux ratios between the detected count rate with the highest 1σ lower limit and the T11 survey sensitivity, which was the deepest archival survey and covered all source positions in our catalog.Previous studies define candidate Xray transients as sources with a minimum flux ratio between five and ten (e.g.Williams et al. 2008;Laycock et al. 2017).Seven of our eight previously undetected sources have minimum flux ratios greater than ten, and one source is just below this threshold.The minimum flux ratios for these new candidate X-ray transients are given in Table 5.We stress that these are candidate transients because the large uncertainties on these source fluxes near our relatively shallow detection limit can lead to an overestimate of the number of transient sources in the catalog.Brassington et al. (2012) address this issue in their study of LMXB transients by utilizing a Bayesian approach to determine the minimum flux ratios with uncertainty estimates on the number of transient sources.If we take a more conservative approach to calculate the minimum flux ratio and divide the lower limit on the maximum detected count rate by the upper limit on the nondetection survey sensitivity, then three of the eight sources have minimum flux ratios greater than five (013314.01+303839.0,013324.37+303323.0, and 013406.22+304042.9).Nevertheless, the low background of the ACIS instrument and uncrowded field still make these detections possible transient sources, particularly the three sources which are detected primarily in the soft 0.35-8.0keV band despite its degraded sensitivity.The variability characteristics and DC constraints for these sources are discussed in Sections 4.1 and 4.2.
We checked the positions of each source in SIMBAD as well as by eye in optical HST imaging (Williams et al. 2021), and found no bright blue stars (F475W magnitude < 23) or background galaxies in seven of these positions (see Fig. 8 for an example).However, 1.16 arcsec away from 013420.89+304947.9 is located a blue supergiant IFM-B 1809 reported in Ivanov et al. (1993).This star is within 3σ of our source position.The proximity and classification of the star make it a strong HMXB candidate if it is the true counterpart to the candidate transient X-ray source we detected.Spectroscopic follow-up of this source is under preparation by another team to classify the companion and confirm its association with the candidate X-ray transient (M.Lazzarini, private communication).
The lack of optical counterparts for the remaining seven candidate transients may suggest outbursts arising from LMXB systems, where the companion is a low-mass star that would be too faint to detect in current optical surveys.The low count rates of each candidate transient detection are consistent with this idea, as the brightest new candidate transient has luminosity of order 10 36 erg s −1 .If these detections are true transient sources, they could belong to the faint and very faint Xray transient classes which are both phenomenological classes of likely LMXB systems (Bahramian et al. 2021).Many previous studies of extragalactic LMXBs focus on sources in globular clusters (e.g.Bildsten & Deloye 2004;Hunt et al. 2023).These sources are not, however, expected to be observed as transient X-ray sources since they are mostly in compact orbits that enable persistent accretion (Armas Padilla et al. 2023).We checked the candidate transient source positions against the PHATTER Star Cluster Catalog (Johnson et al. 2022) to see if any sources might be associated with a globular cluster position, and did not find any associations within 5σ of the X-ray source positions.These sources could therefore be field LMXBs, which is indeed where most transient XRBs are expected to be found if they form according to the "standard" binary formation model (see Brassington et al. (2012) and references therein).A handful of transient X-ray sources in M33 have been reported by previous surveys (see Section 4.2), most of which are identified as likely HMXBs (Williams et al. 2008).Validating any of the new candidate transients presented here as LMXBs with deeper X-ray monitoring observations would constitute a significant increase to the fraction of LMXB transients in M33.

DISCUSSION
We now discuss how we use our catalog of bright X-ray sources to characterize the amplitudes of the changes in detected X-ray emission within our survey and across archival observations, simulate DCs which are consistent with the observed light curves, and discuss the variability results and hardness ratios of sources which have counterparts already identified as candidate and confirmed HMXBs in the literature.

X-ray Variability Amplitudes
We look for variability on two timescales: short-term (shorter than the roughly seven month baseline of our survey) and long-term (∼1 year or longer).Thus we consider only our survey data when investigating short-term variability, and compare to previous catalogs to investigate long-term variability.
We first consider source variability on short timescales.To quantify variability, we define the variability amplitude Dev max for each source as the maximum difference from the mean count rate divided by the error bar of the measurement most different from the mean count rate in the full 0.35-8.0keV energy band.We calculate the maximum difference from the mean count rate as ∆F max = max(|F − F mean |) where F is the 84% upper or lower limit on each count rate measurement Table 5. Summary of measured counts for new candidate transients detected in this survey.The ObsID of the strongest detection and the 0.35-8.0keV full energy band PNS and net counts are listed for each source.The narrow energy band with the strongest detection is given with the extracted counts.The minimum flux ratio F max /F lim between the 1σ lower limit on the strongest detection count rate and the T11 survey sensitivity is given in the last column.The three sources with F max /F lim > 20 also exceed a minimum flux ratio of five when calculating this value with the lower limit on the maximum detected count rate.(corresponding to a Gaussian 1σ limit), depending on whether the measurement is above or below the mean count rate.The mean count rate F mean is weighted by the inverse of the variance.We calculate the variability amplitude as Dev max = ∆F max /σ, where σ is the upper error bar on the measurement that determines ∆F max when the maximum measured flux is higher than the mean count rate (such as in Fig. 6), or the lower error bar for sources with a "maximum" measured flux lower than the mean count rate (such as in Fig. 9).We also consider the variability index definition of η = (F max − F min )/ (σ max low ) 2 + (σ min upp ) 2 used by T11, where σ is the upper or lower error bar on the minimum and maximum measured count rates.The resulting distribution of variability statistics based on our survey observations alone is given in the top panels of Figure 10.
We investigated the long-term variability of these sources by considering the archival count rates for each source.The maximum baseline for a source covered by all available data increases from our survey's 7 month duration to about 21 years.We calculated the same Dev max and η statistics described above from the merged-ObsID count rates in our survey combined with the modern Chandra equivalent count rates for archival X-ray flux data in M06, T11 and W15 (converted as described for the long-term light curves above in Section 2.2).The distribution of variability statistics for the extended baseline is given in the lower panel of Figure 10.
The top panels of Figure 10 shows that most sources in our catalog exhibit only low amplitude variability within the seven month duration of our observations, with peaks at Dev max = 2 and η = 2.The same is true for the extended baseline variability statistics shown in the lower panels of Figure 10, although several sources vary with much higher significance than in the short-term calculations.The distribution of Dev max for both short-term and long-term data falls off at Dev max = 3, thus we chose to define sources with Dev max > 3 as significantly variable.Within our observations only, eight of our sources pass this threshold, while 23 of our sources pass the threshold when including archival data, showing that significant variability is more common at longer timescales than shorter ones in our sample.The sources displaying variability amplitudes greater than Dev max = 3 on longor short-timescales are given in Table 6.
The distribution of η for short-term variability also peaks just below η = 2 and has a tail that drops off considerably by η = 3, which we again use as the threshold for flagging sources as significantly variable according to this test.Seven sources pass this cut for the short-term calculation of η.For the long-term calculation, the η distribution has a much longer tail.Nearly half of the sources in our catalog (26 sources) pass the threshold for variability η > 3 on the extended baseline, 18 of which also satisfy η > 5 which is the threshold that T11 uses to flag sources likely to be significantly variable.A total of 27 sources are indicated as variable on at least one timescale according to the thresholds for η.
Our definition of the variability index Dev max is sensitive to a different variability pattern than the index η used by T11, but still identifies similar numbers of variable sources among the bright X-ray source population in the direction of M33.For long-term variability, 32 sources are indicated to be variable by at least one of the two statistics, with 17 of these sources flagged by both Dev max and η.For the short-term variability within our survey's duration, nine sources are flagged as highly likely to be variable with at least one index, with agreement from both on six sources.One notable source is 013311.76+303843.1.In the top left panel of Figure 10, this source appears at the very start of the tail with a value of Dev max = 3.29.The variability increases significantly to Dev max = 58 when considering literature values.The error bars for both literature measurements and in our observations have little to no overlap with each other except in four of our observations for which the observed count rate is essentially constant.Thus, the calculated variability amplitudes provide a robust description of the true variability of the source.For 013311.76+303843.1 it is likely that the fluctuations in count rate in our observations are part of a more lengthy variability cycle.
Sources 013343.28+304630.9 and 013342.02+303848.6 are known variable foreground stars in the direction of M33 (Hatzidimitriou et al. 2006;Tüllmann et al. 2011), which were also the only sources identified as foreground stars in T11 detected in our survey.Our variability statistics both find that these two sources are significantly variable on the extended baseline (see Table 6).Without a known nonvariable foreground source to compare to, it is difficult to define a threshold for variability intrinsic to sources in M33.Even so, a more conservative threshold of η > 5 yields 18 sources with significant variability on at least one timescale.
The numerous nonvariable sources in our catalog serve as an additional validation tool for the data reduction techniques used for this survey.If the images from each ObsID are well aligned, then sources that are persistent across timescales should show low variability within our observations and when including the archival observations.There are 15 sources that show persistent flux within our observations and long-term, and 23 sources that do not vary over the duration of our observations but do vary long-term according to at least one variability index threshold.Since over one quarter of the sources in our catalog are persistent, it is likely that our images are well aligned and the 23 sources varying only long-term are sources with variability cycles extending beyond the duration covered by our survey.
The fraction of persistent sources in our catalog is consistent with trends observed in the Milky Way and Magellanic Clouds.Two thirds of Galactic X-ray binaries are persistent sources (Tanaka & Shibazaki 1996).The recurrence timescale of X-ray transients is on the order of years (Degenaar et al. 2012), which is longer than the duration of this survey but encompassed by the sources showing longterm variability with the extended baseline.Supergiant fast X-ray transients, a subclass of HMXBs, flare on timescales of a few ks repeatedly over timescales of days to months (Sidoli et al. 2019).Although these sources reach relatively low peak X-ray luminosities (rarely exceeding L X ∼ 10 36 erg s −1 ), the strongest of these types of sources could be detectable in our survey.

Duty Cycle Constraints
A source's DC describes the fraction of time it is in a high emission state.Sidoli & Paizis (2018) found that each HMXB subclass typically follows a different range of DCs, with soft X-ray transients displaying low DCs (< 5%) and wind-fed HMXBs with supergiant companions displaying DCs above 10%.To constrain the X-ray DCs of sources in our catalog, we simulated light curves for a grid of DCs and timescales and determined the frequency with which the observed light curves match the simulated light curves.This duty cycle analysis can only be performed for sources that are considered "off" (by our criteria described below) in at least one observation, which is a sub-sample of the total population detected and analyzed in this work.
We attempted to constrain the DCs of sources in our catalog which were not detected in all observations that covered their positions.If the X-ray sources in the direction of M33 all had 100% DCs, we would expect to detect all 94 T11 sources covered by our survey region that were previously above our 10 −4 ct s −1 detection limit and did not have very soft Γ > 3 spectral index fits (which our survey would likely not detect due to the deterioration of the ACIS instrument's soft-band sensitivity; see Section 2.2).Instead, we detected 48 of these sources.Thus, about half of these sources, many of which are expected to be HMXBs, are likely variable on long timescales.The archival survey luminosity limits in Table 2 show that previous X-ray catalogs would have detected any sources in our catalog if they had been in outburst during those epochs.
To simulate the short-term variability within the duration of our survey (∼ seven months), we produced 100 simulated light curves for each combination of DCs and periods of variability (which we refer to as "timescale").The DC grid ranged from 5-100% at 5% intervals and the timescale grid included 10-200 day timescales at 10 day intervals.Each simulated light curve is a series of quasiperiodic cycles in which the source is assigned as "on" (in a high-emission state) for the number of days equal to the cycle length multiplied by the DC.All days the source is on in each cycle are sequential, and the day in the cycle that the source turns on is random.The time between our earliest observation in the survey presented here to the end of our observations is 228 days.We determined the number of cycles needed to construct each curve by dividing the total observation time by the timescale and rounding up, then adding one extra cycle, allowing us to randomize the cycle start time.

% Match
Figure 11.Representative duty cycle simulation result for a source (013314.01+303839.0)marked as on in all our observations and off in at least one of M06, T11, and W15.The plotted quantity is the fraction of simulated light curves that match the observed light curve for each combination of DC and variability timescale.The DCs for these sources are constrained to be greater than 30% with good agreement across all timescales we could probe within our extended baseline.
We then mapped our observed light curves to a modified light curve to directly compare to the simulations, marking the source as "on" or "off" (in a low-emission state) for each ObsID.We define a source to be off if the upper error is lower than 1% of the lower error on the faintest detection.The 1% buffer region is included to mark a single observation of source 013350.49+303821.2 as on so that its DC may be constrained, despite the upper error being slightly below the faintest detection.Next, we compared the observed light curves to the simulated curves by randomly choosing a day in the first cycle of the simulated curve to set as the first day.We considered a simulated light curve to match an observed curve if the source activity was consistent on all days of our observations.For each combination of timescale and DC, we recorded the fraction of simulated curves that matched each source.The results are represented as 2D histograms as in Figure 11.
The distribution of the number of off epochs by peak count rate is shown in Figure 12.Within our survey, only 11% of our sources have any off epochs.When including literature results, still only 34% of our sources have any off epochs.To further break down the variability of the observed sources, we compare the number of off epochs observed to the peak count rate recorded.In our observations, 28% of sources with peak count rates less than or equal to 10 −3 ct s −1 (corresponding to L X ≲ 1.7 × 10 36 erg s −1 using the same assumptions to derive our limiting luminosity in Section 3) have at least one off epoch and 22% have at least four off epochs.For peak count rates greater than 10 −3 ct s −1 and less than or equal to 10 −2.5 ct s −1 (1.7 × 10 36 erg s −1 < L X < 5.3 × 10 36 erg s −1 ), 4% of sources have at least one off epoch and none have more than three.There are also no sources with peak count rates at values greater than 10 −2.5 ct s −1 and any off epochs.The decreasing number of sources with off epochs as peak count rate increases suggests that faint sources are more likely to turn off than bright sources, although fluctuations near our survey sensitivity limit may exaggerate this trend.When including past survey measurements, 93% of sources with a peak count rate less than or equal to 10 −3 ct s −1 have at least one off epoch, and 29% have at least four.For peak count rates greater than 10 −3 ct s −1 and less than or equal to 10 −2.5 ct s −1 , 21% of sources have at least one off epoch and none have greater than three.There are no sources with any off epochs at peak count rate values larger than 10 −2.5 ct s −1 .The long-term data also supports the trend of the number of off epochs decreasing as peak count rate increases.The larger percentage of sources with off epochs including long-term observations is also consistent with the finding from the variability statistics in Section 4.1 that sources in our data set tend to vary more long-term than short-term.
We next simulated additional light curves with an extended baseline to match to the long-term X-ray light curves including our observations and those from M06, T11 and W15.The time from the earliest survey observation to our last observation is just over 7646 days, so we used 50-1000 day timescales with 50 day steps, and 5-100% DCs in increments of 5%. to create the modified on-off light curves of our observations and archival observations, we considered the source on for all individual observation days in the archival surveys where the source was detected.We checked for matches to transient sources in the archival surveys for which the assumption of being on at all survey times would not be valid.Only one source in the match to M06 (their Source Number 251, our source 013431.97+303454.2) was flagged as having both dropped below their detection threshold and being significantly variable.Our catalog includes one additional match to a transient source identified from the T11 data set by Williams et al. (2008), designated as XRT-2 by that paper (T11 Source Number 210, our source 013332.24+303955.4).W15 detects XRT-2 but does not flag it as transient since it was a faint source that did not have to vary by a factor > 10 to explain the past nondetections.We discuss the potential impact on our DC results for these two known transient sources below.
We matched our observed long-term light curves to the new simulated curves, considering a match to be where source activity was consistent for all days in our observations and at least one day in each literature survey.Most sources had the highest fraction of simulated light curves in agreement with the observed light curves for shorter timescales, from 50 to 200 days.Thus, we created a new set of simulated curves with timescales every 10 days from 10 to 200 days and compared to the data again.
We performed this DC analysis for the 19 sources in our catalog which could be classified as off in at least one observation when including archival data.Of these, strong constraints can be made on six sources that were definitely not detected in at least one observation within our survey.The other 13 sources were on in all of our relevant observations and only off in at least one of M06, T11, W15 (six of which are new candidate transients presented in Section 3.2).These 13 sources are consistent with DCs > 30% without any strong constraint on timescale (see Fig. 11), so we expect these sources to have DCs that operate on timescales longer than the seven month duration of our survey.This is consistent with common intervals between XRB transient outbursts of multiple years (Degenaar et al. 2012).These 13 sources with DCs > 30% that operate on timescales longer than those probed in this monitoring survey are 013314.01+303839.0,013324.37+303323.0,013336.04+303333.1, 013337.94+303837.4,013356.09+303024.8,013357.09+304621.7,013358.03+303201.0,013358.09+303438.0,013401.12+303242.2,013402.87+304151.4,013406.22+304042.9, 013416.32+304319.2, and 013420.89+304947.9.
For the six sources which have at least one off observation in our survey and an archival survey, the fraction of simulations that match the observed light curves are shown in the density plots in Figure 13.We report the timescale and DC most consistent with the observed light curves as well as the strength of the peak for these sources in Table 7.The upper and lower errors are defined by the highest and lowest DC and timescale that match the observed light curves in at least 1% of the simulations.Of these six sources, two have DCs above 50% (013350.49+303821.2 and 013426.47+304446.4).Source 013426.43+304446.4 varies significantly over our observations, while both sources with DCs higher than 50% meet the criteria for being significantly variable when considering the survey observations each was present in.Source 013350.49+303821.2 was on the detector for all ten of our observations and 013426.43+304446.4 for eight.Meanwhile, the remaining four sources were only on the detector for five observations each.Two of these sources are new candidate transients, 013347.76+303300.0 (see Fig. 7) and 013353.05+304710.4,and are measured as being on in only one observation, indicating very short outburst lengths.Although these sources have lower DCs, these small outbursts indicate they may be transients with shorter cycle durations.
Our survey detects the known transient M06 Source Number 251 as on in all four observations where it was on the detector.This source is marked as on in each of the archival surveys as well, so no DC information could be obtained when using the assumption that the source was on during all days in the archival surveys.We ran an additional set of duty cycle simulations for this source which included the M06 nondetection with an upper limit below the observed flux of the source as an off epoch and found that the DC must be higher than 30% with agreement across timescales, similar to the 13 sources discussed above.We also detect the known transient XRT-2, which is listed as an HMXB candidate by T11 and G18 (T11 Source 210; see Section 4.3 below).We ran an additional suite of duty cycle simulations for this source including the off epoch within the T11 ObsIDs and obtained the same result of DC > 30%.
We   (Mori et al. 2021).We can identify high DC sources with the present data set, but more frequent X-ray monitoring would likely enable constraints on many more low-duty DCs.

Confirmed and Candidate HMXBs
Our catalog includes observations of 26 confirmed and candidate HMXBs previously identified in the literature by T11, G18, Yang et al. (2022) and L23.These are listed in Table 8 and included in the full version of Table 3 available online, including the nucleus source M33 X-8 (013350.88+303936.5) and the eclipsing HMXB M33 X-7 (013334.14+303211.1;see Fig. 9).Every HMXB candidate is associated with a T11 source by the study that identifies them.We find that 17 of these known HMXB candidates show significant variability on at least one timescale.Six of these sources are persistent within our survey and only off in an archival survey, consistent with the DCs constrained to > 30% discussed above in Section 4.2.The remaining eight sources show persistent emission.
G18 identified 55 HMXB candidates from a combination of the T11 catalog and archival HST and Spitzer imaging.Fourteen of the G18 candidates were observed in our survey.There is one source in the G18 list that we do not detect at all in our survey (J013410.69+304224.0) that falls within our survey region and was above our detection limit in the past, which was classified as a supernova remnant by Long et al. (2010).
Twelve of the 28 sources from the Yang et al. (2022) candidate list observed with NuSTAR are detected in our survey.The positions of these sources on a hardness-intensity diagram in that paper (Figure 3) indicate the nature of the compact object in each XRB candidate.Most of our detected sources with matches to the Yang et al. (2022) catalog are identified as likely containing a BH, with the exceptions of 013328.70+302724.4 (M33 X-6; Y22 Source 2), 013342.55+304253.5 (Y22 Source 18) and possibly 013436.36+304713.9 (Y22 Source 13) which are more consistent with where the pulsar population lies.
Our survey detects 17 of the 65 HMXB candidates reported in L23, three of which are classified as HMXB candidates for the first time.Multiple optical counterparts from PHATTER (Williams et al. 2021) are listed for nine of these sources.Among the top ranked optical counterparts identified by L23 for the 17 sources we detect in this survey, about half are main sequence B stars.This is a lower fraction than in the full L23 catalog, which reports 75% of top-ranked optical counterparts to T11 sources to be main sequence B stars.We also detect two of the four sources with main sequence O star companions, all three of the sources with giant O star companions, and the single source with a giant B star companion.This is consistent with our systematic limitation of detecting only the most energetic systems with our shallow survey, and we also expect some intrinsic variability due to a few sources being at different phases in their activity cycles during the T11 observations and during our survey.Seven sources of this list are "off" in at least one observation in our survey which suggests that these sources are significantly variable and therefore stronger HMXB candidates.All sources from that list are covered by our survey's detection area and four of the undetected sources were above our detection threshold in the past (T11 IDs: 274, 277, 391 and 452).T11 ID 233 was also above our detection limit in the past, but was found in T11 to have a soft power law spectrum with Γ = 2.09 so this source would be below our detection limit due to the deterioration of Chandra's soft-band sensitivity.The other four sources would have been detectable by our survey at their luminosities observed by T11 which suggests that they were in a low X-ray emission state at the times of all of our observations and are also significantly variable sources.T11 277 has a supergiant B star optical companion, and T11 391 has a main sequence B star.T11 274 has two possible optical counterparts which are a main sequence B star and a main sequence O star.T11 452 has three possible counterparts, two of which are main sequence B stars.
Eight of the sources in our catalog with cross matches to T11 are identified there as XRB or related classes.Only one of these does not make it into the three other lists discussed, and was identified as a quasar by Neugent & Massey (2011).We include this source (013445.07+304924.5) in Table 8 for completeness.We detect all sources marked as XRB or related classes that were above our detection limit in the T11 catalog and covered by our survey region.
HMXBs with active accretion are hard X-ray sources (Di Salvo et al. 2004;Yukita et al. 2016;Yang et al. 2022).We attempted to characterize the hardness of all sources in our catalog using the Bayesian Estimation of Hardness Ratios program (BEHR; Park et al. (2006)).Due to the poor soft-band response of the ACIS instrument, hardness ratios (HRs) using counts in the 0.35-1.1 keV energy band are severely biased.We therefore only consider a HR between the medium 1.1-2.6 keV energy band (M) and the hard 2.6-8.0 keV energy band (H).Seven sources resulted in unconstrained HRs due to low counts.We show a plot of HR versus count rate for the 48 sources in our catalog that BEHR found HRs with constraining 1σ error bars for in Figure 14.We find that the HRs of 17 out of the 21 HMXBs with constrained HRs are as hard or harder than the well-studied HMXB X-7, indicated in green in Figure 14.The remaining four HMXB candidates (Sources 013337.94+303837.4,013340.04+304323.0,013350.49+303821.2 and 013407.80+303553.8)have HR upper limits that are softer than X-7, but still may not necessarily be soft sources since we cannot calculate a HR with the soft band counts.

CONCLUSIONS
We have presented a new five-epoch survey of 55 bright X-ray sources in M33 probing multiple timescales, including eight previously undetected candidate transients and 26 known HMXB candidates outside of the nucleus.We characterized the short-and long-term variability amplitudes of 32 significantly variable sources, identified 15 nonvariable sources, and found constraints on the DCs of 21 sources in our catalog.We found a higher fraction of sources with DCs > 30% than previous HMXB literature suggests due to the infrequent observation history of M33.There are many observations that would be useful to classify the variable X-ray sources in this work, including further X-ray monitoring and inspection of the optical counterparts to these sources.While a majority of the bright X-ray sources in M33 are predicted to be HMXBs, there are likely also background AGN present in the sample as contaminants.A careful multi-wavelength selection is needed to get the cleanest intrinsic XRB sample for targeted timing studies.Once the population of HMXBs in M33 is well defined, deeper X-ray observations across the entire population are necessary to improve our understanding of these energetic systems as a whole.

Figure 1 .
Figure 1.Time sampling of ObsIDs in our survey.The lower points (ObsIDs 23603 to 23607) correspond to pointings centered on the Southern half of M33, and the upper points (ObsIDs 23608 to 23612) correspond to pointings centered on the Northern half of M33.

Figure 2 .
Figure 2. Distribution of PNS values in the aligned multi-ObsID source extraction for each energy band.The vertical dashed lines indicate the maximum value of PNS = 10 −4 used to select sources for the final catalog.Any source with PNS = 0 was omitted in the plot(s) for the bands in which that measurement occurred.

Figure 3 .
Figure 3. Sensitivity map for ObsID 23603 in the 0.35-8.0keV band.The colorbar is in units of ct s −1 .The feature in the chip gap is at the position of the bright nucleus source M33 X-8.

Figure 4 .
Figure 4. Log of the signal to noise ratio (SRC SIGNIF reported by AE) plotted against the log of the PNS for all sources.The full band (0.35-8.0 keV) data is used if the source passes the PNS < 10 −4 cut in the full band; otherwise the data is for the band where the PNS is lowest.Orange points are from the merged multi-ObsID extraction.Sources which only pass the PNS cut in an individual observation are represented by green points.Blue points show sources which do not pass the PNS cut in any observation or energy band, so the lowest PNS value between all merged and individual values is used.The horizontal dashed line indicates a signal to noise ratio of two.The vertical dashed line marks the PNS = 10 −4 threshold used to select our final catalog as used in Figure 2.

Figure 5 .
Figure 5. Left panel: Outline of the total sky coverage of the ten pointings in this survey (solid black line) and positions of our catalog sources (red crosses) overlaid on a smoothed GALEX UV image of M33 (Lee et al. 2011).Right panel: Stacked exposure maps on a square root scale in the 0.35-8.0keV band.
Figure 6.X-ray light curve of Source 013314.01+303839.0, a new candidate transient detected by our survey with the highest net counts.The left panel shows the long-term light curve with archival upper limits from the M06, T11 and W15 surveys along with our merged-ObsID count rate.The blue shaded region indicates the duration of time covered by our survey compared to the full baseline.The right panel shows the short-term light curve with count rates from individual ObsIDs within our survey only.Observations from pointings centered on the Northern half of M33 are marked in blue, and observations from pointings centered on the Southern half of M33 are marked in orange.Upper limits are shown with downward arrows where the source position was on the detector but no counts above the background were extracted.

Figure 7 .
Figure 7. X-ray light curve of Source 013347.76+303300.0, a new candidate transient detected by our survey with a duty cycle constrained by our simulations between 10-60% on an outburst recurrence timescale between 10-200 days.The dashed black line on the right panel shows an example iteration of the simulated DC and timescale combination with the peak match to the observed light curve (DC = 50% on 110 day timescale).The off times are represented at zero ct s −1 and the on times are represented at 1 × 10 −3 ct s −1 .

Figure 8 .
Figure 8. Stacked 10 ′′ × 10 ′′ (∼ 40 × 40 pc at the distance of M33) HST image with the red F160W, green F814W, and blue F475W bands used to check for optical counterparts to the new candidate transient 013347.76+303300.0.The 1σ and 3σ error circles (0.4 ′′ and 1.2 ′′ respectively) for the X-ray detection are shown, and no bright blue star or background galaxy is seen.The same visual inspection was performed for the other seven new candidate transients as discussed in Section 3.2.

Figure 10 .
Figure 10.Distribution of sources by variability indices Dev max (left) and η (right).The top panels show the distributions based on short-term variability (i.e.within our ten ObsIDs covering ∼seven months) and the bottom panels show the distributions based on long-term variability over both our observations and archival surveys which cover a maximum baseline of 21 years.

Figure 12 .
Figure 12.Number of epochs sources are labelled as off binned by the log of the peak count rate in ct s −1 .The top panel shows the number of off epochs for short-term light curves.The lower panel shows the number of off epochs for long-term light curves including our observations and archival data.

Figure 13 .
Figure 13.Simulated duty cycle match grids for the sources with at least one "off" epoch within our survey.
identify 17 sources with DC > 30% out of the 21 sources with at least one off epoch.Observed HMXB populations in both the Milky Way and Magellanic clouds generally show low fractions of high DC sources.Of the sample of 56 sources observed in the Milky Way by INTEGRAL, Sidoli & Paizis (2018) reports only 23% as having DCs above 25%.Similarly, Tetarenko et al. (2016b) report

Figure 14 .
Figure 14.Hardness ratio (HR) between the hard 2.6-8.0 keV energy band the medium 1.1-2.6 keV energy band versus merged-ObsID 0.35-8.0keV count rate.HMXB candidates are plotted in orange, and the well-studied HMXB M33 X-7 is plotted in green.Other catalog sources which we could constrain the HRs for are plotted in blue.Most HMXB candidates are as hard as or harder than X-7.

Table 1 .
Summary of Chandra observations taken for this survey.The Pointing column indicates whether the ObsID was centered on the Northern (N) or Southern (S) halves of M33.The last column gives the effective exposure time after pipeline filtering.Pointing ObsID Ob.Start R.A. (J2000.0)Dec. (J2000.0)Roll Angle ( • ) Exp.Time (ks) Source detected in the 0.5" pixel size wavdetect run only (see Section 2.1).Ten additional single-pixel size detections are indicated in the full table online.

Table 4 .
Preview of the final catalog source likelhood and count rates in every energy band considered for each individual ObsID and the merged extraction.The first source is an example of a source primarily located in Southern pointing observations, and the second source is an example of one primarily located in Northern pointing observations.Blank entries appear where a source was not within the field of an ObsID.1σ upper limits are given when a source's position was on the detector for an ObsID but was not detected.The full version of this table can be found online, as well as an additional table with one source per row including all extracted source properties returned by AE including source significance and background levels.

Table 6 .
Sources in our catalog with variability indices η > 3 and Dev max > 3 for longterm variability across archival and current observations and short-term variability within our observations.

Table 7 .
Results from the most successful duty cycle simulation constraints.Source 013350.49+303821.2 had two equally likely simulated eruption timescales (70 d and 80 d) with the same duty cycle, indicating that its best match was on the edge of the pixel at 75 d.The upper and lower errors are defined by the highest and lowest timescale/duty cycle values where more than 1% of the simulated light curves matched the observed light curve.the majority of observed transient black hole XRBs in both the Milky Way and Magellanic Clouds to have DCs less than 10% (most of these are LMXBs).About 75% of BH-LMXBs in the Milky Way turn on once every 50 years

Table 8 .
HMXB candidates previously identified in the literature with matches to sources in our catalog.