Multiwavelength Characterization of the High Mass X-ray Binary Population of M33

We present multi-wavelength characterization of 65 high mass X-ray binary (HMXB) candidates in M33. We use the Chandra ACIS survey of M33 (ChASeM33) catalog to select hard X-ray point sources that are spatially coincident with UV-bright point source optical counterparts in the Panchromatic Hubble Andromeda Treasury: Triangulum Extended Region (PHATTER) catalog, which covers the inner disk of M33 at near infrared, optical, and near ultraviolet wavelengths. We perform spectral energy distribution (SED) fitting on multi-band photometry for each point source optical counterpart to measure its physical properties including mass, temperature, luminosity, and radius. We find that the majority of the HMXB companion star candidates are likely B-type main sequence stars, suggesting that the HMXB population of M33 is dominated by Be-XRBs, as is seen in other Local Group galaxies. We use spatially-resolved recent star formation history (SFH) maps of M33 to measure the age distribution of the HMXB candidate sample and the HMXB production rate for M33. We find a bimodal distribution for the HMXB production rate over the last 80 Myr, with a peak at $\sim$10 Myr and $\sim$40 Myr, which match theoretical formation timescales for the most massive HMXBs and Be X-ray binaries (Be-XRBs), respectively. We measure an HMXB production rate of 107$-$136 HMXBs/(M$_{\odot}$ yr$^{-1}$) over the last 50 Myr and 150$-$199 HMXBs/(M$_{\odot}$ yr$^{-1}$) over the last 80 Myr. For sources with compact object classifications from overlapping NuSTAR observations, we find a preference for giant/supergiant companion stars in BH-HMXBs and main sequence companion stars in neutron star HMXBs (NS-HMXBs).


INTRODUCTION
HMXBs, allowing us to measure population demographics that can constrain theoretical models of massive binary stellar evolution.
The properties of HMXB populations scale with galaxy-wide properties such as star formation rate (SFR; Ranalli et al. 2003;Gilfanov et al. 2004;Antoniou et al. 2010;Mineo et al. 2012;Antoniou & Zezas 2016;Lehmer et al. 2019;Antoniou et al. 2019), metallicity (Basu-Zych et al. 2013Brorby et al. 2016), stellar mass (Lehmer et al. 2010;Antoniou et al. 2019), and the number of O-and B-type stars (Antoniou et al. 2019). In nearby galaxies, we can place HMXBs in context of both their host galaxy and their more local parent stellar population.
Massive stars (≳8 M ⊙ ) play a critical role in many astrophysical processes. During their lives they inject ionizing radiation into the interstellar medium (ISM) and at the end of their lives impact the evolution of galaxies via supernova feedback (e.g., Dalgarno & Mc-Cray 1972;Oey 1999;Telford et al. 2023). Most massive stars form in binary or higher order multiple systems and a majority of these systems have separations small enough for the stars to interact at some point during their lives (e.g., Sana et al. 2012;Moe & Di Stefano 2017). These interactions, including mass transfer and common envelopes, impact the evolution of the stars in the binary. Some massive stellar binaries go on to form binary compact objects which can merge, depending on their separation, within a Hubble time generating detectable gravitational waves (e.g., Tauris et al. 2017). HMXBs are also thought to be a major contributor of ionizing radiation in the early universe, likely playing a key role in cosmic reionization (e.g., Justham & Schawinski 2012;Mesinger et al. 2014;Madau & Fragos 2017;Greig & Mesinger 2018).
Given their importance in many fields of astrophysics, HMXBs have been studied extensively in the Milky Way and other galaxies. Populations of HMXBs can be studied in nearby galaxies with relatively few telescope pointings, and the proximity of these galaxies allows us to resolve their individual stars, which is necessary for identifying the companion star to the accreting compact object in the HMXB system. The populations of HMXBs in the Large and Small Magellanic Clouds have been studied extensively at both X-ray and optical wavelengths. These studies have identified the companion star spectral type and compact object type for over 200 HMXBs in both galaxies combined, providing an insight to HMXB formation at low metallicity and high SFR (e.g., Liu et al. 2006;Antoniou et al. 2009Antoniou et al. , 2010Antoniou & Zezas 2016;Yang et al. 2017a;Haberl & Sturm 2016;Antoniou et al. 2019;Haberl et al. 2022). In M31, overlapping Hubble Space Telescope (HST) and Chandra observations have been used to measure HMXB population demographics in a large spiral galaxy of similar mass and metallicity to the Milky Way Lazzarini et al. 2018Lazzarini et al. , 2021.
M33 is an obvious site to study HMXBs in the Local Group. It is the third largest spiral galaxy in the Local Group and has a measured metallicity gradient; the center of the galaxy is slightly super-solar and the metallicity drops to roughly that of the LMC at the edges of the disk (Carrera et al. 2008;Cioni 2009;Magrini et al. 2009Magrini et al. , 2010Beasley et al. 2015;Toribio San Cipriano et al. 2016;Lin et al. 2017). Although it is about ten times less massive than M31 or the Milky Way (e.g., Quirk et al. 2022), it has a higher SFR intensity (SFR per area; Verley et al. 2009;Lewis et al. 2015), making it an excellent site to study young stellar systems. Because of its high SFR, M33 also hosts a large population of supernova remnants (SNRs), which have been studied at X-ray and optical wavelengths (e.g., Long et al. 2010;Garofali et al. 2017;Koplitz et al. 2023).
M33 has been observed extensively at X-ray wavelengths with the Einstein Observatory (Long et al. 1981;Markert & Rallis 1983;Trinchieri et al. 1988), ROSAT (Schulman & Bregman 1995;Long et al. 1996;Haberl & Pietsch 2001), XMM-Newton (Pietsch et al. 2004;Misanovic et al. 2006;Williams et al. 2015), Chandra (Grimm et al. 2005;Tüllmann et al. 2011), and NuS-TAR (Yang et al. 2022). The Chandra ACIS Survey of M33 (ChASeM33) observed 70% of the area within the D 25 isophote of M33, providing a deep catalog of X-ray point sources with high angular resolution, which is critical for identifying the optical counterparts to accreting compact objects in HMXBs. A uniform survey of M33 at optical wavelengths was not available until recently with the Panchromatic Hubble Andromeda Treasury: Triangulum Extended Region (PHATTER), an HST survey of the inner disk of M33 in six photometric bands spanning near-infrared through near-ultraviolet wavelengths . The survey's photometry catalog contains six-band photometry for over 20 million individual stars and has been used to study the galaxy's recent star formation history (SFH; Lazzarini et al. 2022), structure (A. Smercina et al. 2023, in prep.), and population of star clusters Johnson et al. 2022).
Previous studies of the HMXB population of M33 have led to both detailed characterization of a few individual sources and properties of the population. Garofali et al. (2018) performed a population study using archival HST observations and found a peak in the M33 HMXB age distribution at <5 Myr and another peak at ∼40 Myr, with a valley between these two peaks with little to no HMXB production. Previous studies of HMXBs in M33 focused on individual systems. The nucleus of M33 hosts the nearest ultraluminous X-ray source (ULX) (Long et al. 2002), X-8. The X-ray emission is powered by super-Eddington accretion, although whether the compact object is a black hole or neutron star has not yet been firmly determined (Krivonos et al. 2018;West et al. 2018). There are also two eclipsing X-ray binaries in M33 (Pietsch et al. 2004(Pietsch et al. , 2009. One of the two, M33 X-7, comprises a massive stellar-mass black hole in orbit with an O-type stellar companion (Orosz et al. 2007;Ramachandran et al. 2022). The exact orbital separation and masses of the black hole and companion star in M33 X-7 have been debated, but the system has provided interesting constraints for massive binary stellar evolution models. The other eclipsing binary in M33 (Pietsch et al. 2004) has been less well studied and is not covered in the footprint of our analysis in this paper.
In this paper, we leverage the deep, uniform coverage of the inner disk of M33 at near-infrared/optical/nearultraviolet wavelengths to systematically identify HMXB companion star candidates. We perform spectral energy distribution (SED) fitting on the HMXB companion star candidates to infer their physical properties including temperature, radius, luminosity, and mass. We use the spatially-resolved recent SFH of M33, measured using the optical PHATTER photometric catalog, to place constraints on the age distribution of the HMXB population of M33 and its HMXB production rate over time.
We describe the X-ray and near-infrared/optical/nearultraviolet observations used in this analysis in Section 2. In Section 3 we describe how we identify and evaluate the quality of the HMXB candidate sample and our methodology for SED-fitting of the companion star candidates and measurement of the HMXB candidate ages using the SFH maps. We include a discussion of our results in Section 4 and summarize our findings in Section 5. We assume a distance to M33 of 859 kpc, or a distance modulus of 24.67, throughout our analysis (de Grijs & Bono 2014).

DATA
We use X-ray and optical/near IR/near UV imaging and photometric catalogs of the inner disk of M33 to identify our HMXB candidate sample. We use the final source catalog from the ChASeM33 survey (Tüllmann et al. 2011) and the optical/near IR/near UV catalog from the PHATTER survey . In Figure 1 we show the outline of each survey. We provide more details on each survey below.

X-ray Data
The ChASeM33 ACIS survey covered ∼70% of the D 25 area of M33 in the 0.35−8.0 keV band (Tüllmann et al. 2011). The resulting catalog includes 662 sources, 202 of which lie within the PHATTER survey footprint in the inner disk of M33. The catalog has a limiting unabsorbed luminosity of ∼ 2.4 × 10 34 erg s −1 in the full 0.35−8.0 keV band and includes positions, fluxes, and net counts for each source in multiple energy bands.
We use the X-ray source positions from this catalog, along with their associated errors, to identify the optical counterpart candidates for each X-ray source, thus creating our sample of HMXB candidates. We also use count rates and fluxes, and spectral fits for sources where these are available, from the ChASeM33 catalog when evaluating the quality of our HMXB sample. We provide more details on how we identify optical counterparts to X-ray sources in Section 3.1 and how we evaluate the quality of our HMXB candidates in Section 3.4.
In Table 1 we include a list of HMXB candidates. We include each source's identification number, X-ray source position, and 0.35−8.0 keV luminosity from (Tüllmann et al. 2011). We calculated the luminosity using the count rates published by Tüllmann et al. (2011). For this calculation we assume a power law index of 1.7 and a foreground N H of 5×10 20 cm −2 (Lebouteiller et al. 2006). We used the Chandra proposal planning toolkit PIMMS tool to derive a conversion factor of 1.087×10 −11 from count rate (counts s −1 ) to flux , which we then convert to luminosity.
We also cross-matched each X-ray source with other X-ray surveys of M33 including Grimm et al. (2005, Chandra), Williams et al. (2015, XMM-Newton), and Yang et al. (2022, NuSTAR). We include the counterpart of each source in the aforementioned catalogs in Table 1 and include their classifications (if provided) in each catalog. We describe our cross-matching methodology in Section 3.5.

Optical Data
We use the photometric catalog and imaging from the PHATTER survey in our data analysis ). The PHATTER survey covered an area of 14 kpc 2 in the inner disk of M33, providing six-band photometry for over 22 million stars. In this study, we use the mosaic imaging from the PHATTER survey and the photometric catalog  to identify the optical counterpart for Tüllmann et al. (2011) Xray sources within the PHATTER survey footprint. We also use the spatially resolved recent SFH maps derived from the PHATTER optical photometry by Lazzarini

Astrometric Alignment
We use the X-ray source properties and positions from the ChASeM33 survey (Tüllmann et al. 2011) and optical imaging and photometric catalogs from PHAT-TER ) in our analysis. The ChASeM33 catalog is aligned to 2MASS (Cutri et al. 2003;Tüllmann et al. 2011;Garofali et al. 2018) and the PHATTER catalog is aligned to GAIA. To ensure that our X-ray and optical data were aligned, we selected sources in the PHATTER catalog with a F160W magnitude <15.5, to match the 99.9% completeness limit of the 2MASS survey in the H band, which closely overlaps with the F160W filter on HST. We then matched the 2MASS catalog to the PHATTER sources. We performed this task in an iterative fashion so that false matches could be culled from our list of matches.
To confirm our matches, we plotted the PHATTER photometry and the 2MASS photometry for each source to create a spectral energy distribution (SED). We visually inspected the SED for each match to confirm that the SED shape looked smooth at the transition from the HST to 2MASS filters. We also plotted the positions of all stars within 5 arcsec that fit our PHATTER criteria (F160W<15.5) and all 2MASS stars that lie within the same area to confirm that the sources look like a good match. Due to our magnitude cuts, the source density of PHATTER sources is greatly reduced, reducing the chances of an erroneous match.
We find an offset between the two catalogs of ∆R.A. 0.04 arcsec and ∆Dec. 0.004 arcsec, with a RMS of 0.1 in both R.A. and Dec. The offset in both R.A. and Dec. is of the same order of magnitude as the mean 2MASS positional errors for the sources we use as matches in our alignment. Because the offset between the two catalogs is less than the mean 2MASS positional error for all sources, we did not update the positions of the sources in the ChASeM33 catalog (aligned the 2MASS) or the PHATTER catalog.

ANALYSIS
We use near-UV/optical and X-ray observations to identify 65 HMXB candidates in M33. We perform SED fitting to determine the likely physical properties of candidate HMXB companion stars. Using maps of the spatially resolved recent SFH of M33, we measure the age distribution and production rate for the HMXB candidate sample. We use multi-wavelength information including Chandra hardness ratios, UV/optical SED shapes, local star formation rates, and measurements of extinction to evaluate our HMXB candidate sample and identify the highest quality HMXB candidates. Lastly, we cross-match our HMXB candidate sample with previous catalogs of X-ray sources and HMXB candidates in M33.

Identification of HMXB Candidate Sample
We identify HMXB candidates using HST optical/near-UV imaging and photometry from the PHATTER survey  and Chandradetected X-ray source positions and fluxes/count rates from the ChASeM33 survey (Tüllmann et al. 2011). As described in more detail in this section, we identify 64 HMXB candidates using the methodology described here. We include M33 X-8, the nuclear ULX, in our sample when analyzing the HMXB ages and production rates, which brings our total HMXB candidate sample of 65 sources.
To identify our sample of HMXB candidates, we started with all X-ray point sources from Tüllmann et al. (2011) that were within the PHATTER survey footprint that were not classified as AGN based on optical spectroscopy (Tüllmann et al. 2011), supernova remnants based on their X-ray, radio, and optical properties (Long et al. 2010), or colliding wind binaries (Garofali et al. 2019). We also removed X-ray sources from our sample if the optical counterpart identified in the PHATTER imaging was a resolved galaxy. We include the nuclear HMXB (X-8; ChASeM33 318) in our HMXB candidate sample for age and production rate measurements, although we cannot resolve its optical counterpart within the nuclear star cluster of M33.
We selected optical counterpart candidates from within 1.5σ of the X-ray source position, using the errors and X-ray source positions presented in the ChASeM33 survey catalog. This radius was chosen to maximize the recovery of the HMXB population in M33, while minimizing chance superpositions between X-ray sources and optical counterpart candidates. To settle on this distance, we started with a 1σ error circle and expanded the radius iteratively by 0.1σ. We expanded out to 1.5σ, which allowed us to recover the HMXB candidates identified by Garofali et al. (2018) that fall within the PHAT-TER footprint but that were not otherwise disqualified (i.e. likely SNR, resolved galaxy in PHATTER imaging, colliding-wind binary).
When only selecting companion star candidates from within 1σ, we identified 40 unique HMXB candidates, and when we select from within 1.5σ we identify 64 unique HMXB candidates. As described in more detail in Section 3.1.1, with a match radius of 1.5σ we expect 3−4 false HMXBs due to a chance superposition of a star in the PHATTER catalog that meets our selection criteria and an X-ray source in the ChASeM33 source catalog. With a match radius of 1σ we expect 1−2 false HMXBs due to a chance superposition. When we expand our match radius from 1σ to 1.5σ we gain 24 HMXB candidates, suggesting that the majority of the gained HMXB candidates are not false matches. We expect the X-ray source to be located within the 1σ error on its positions 68% of the time, and when we expand to 1.5σ we expect the X-ray source to be located within this 1.5σ positional error 87% of the time. We find that the trade off of a slight increase in the number of chance superpositions is outweighed by the increase in our sample of HXMB candidates.
We selected optical counterpart candidates that had colors and magnitudes in the PHATTER photometric catalog consistent with being massive (M≳8 M ⊙ ) main sequence, giant, or supergiant stars. We chose this mass cut off because it is the known lower limit for the masses of Be stars in Be-XRBs (Shao & Li 2014;Coe & Kirk 2015). We used Padova stellar models (Marigo et al. 2008) to define the ranges of colors and magnitudes of massive main sequence and giant/supergiant stars. To define the ranges in color and magnitude, we started with Padova stellar tracks and plotted them onto optical and UV color-magnitude diagrams (CMDs). We then applied extinction to each photometric band for our model stars using the coefficients for reddening in the HST bandpasses used in the PHATTER survey (Schlafly & Finkbeiner 2011) and applied 1 and 2 magnitudes of dust extinction (A V ). To ensure that we selected a complete sample, we included stars within 0.5 mag of the lowest expected magnitude for an 8 M ⊙ star. This margin makes us potentially sensitive to unreddened counterpart candidates down to ∼6 M ⊙ .
We used four sets of criteria to select likely massive stars associated with X-ray point sources. We list the exact color and magnitude cuts we used to select HMXB companion star candidates in Table 2. Counterparts were selected if they were within 1.5σ of an X-ray point source position and had colors consistent with being a massive (M≳8 M ⊙ ) main sequence or giant/supergiant star with up to two magnitudes of dust extinction, A V .
Using these criteria, we identify 64 HMXB candidates. The 65th HMXB we include in our sample is the ULX in the nuclear cluster of M33, X-8 (ChASeM33 318), which does not have an identified companion star candidate due to high stellar density in the central region of the galaxy. In our identified HMXB candidate sample, 23 X-ray sources have more than one star that meets one of the aforementioned criteria. We list the HST photometry for all potential HMXB companion star candidates in Table 3 and list the PHATTER position and which  ChASeM33 source the star is associated with in Table  1.
In Figure 2 we present color-magnitude diagrams showing the location of the optical counterparts we identified using the criteria described above. While the location of the HMXB companion star candidates spans about 2 magnitudes in F475W−F814W color (left CMD) we see that all of these stars lie on or close to the main sequence in the optical-UV (upper right) and UV-only (lower right) CMDs, suggesting these stars are likely to be massive stars.

Chance Superposition
We calculated the number of spurious HMXB companion star candidates we expect based on a chance superposition of an OB star within the 1.5σ error circle of our X-ray sources. We applied the selection criteria outlined above to the full PHATTER photometric catalog to find the total of stars in the PHATTER catalog that met our criteria, and turned this into a density of stars (arcsec −2 ) when divided by the total area of the PHATTER survey. We then summed the area of the 1.5σ error circles for all ChASeM33 X-ray sources (regardless of optical counterpart type) that fall within the PHATTER survey footprint. We multiplied this area by the density of OB stars that meet our selection criteria in the PHATTER catalog to find the total number of chance superpositions we expect. We find that there should be 3−4 chance superpositions, so we expect that <10% of our HMXB candidates could be spurious.

AGN Contamination
In the process of identifying HMXB companion star candidates, we visually inspected the PHATTER imaging at the location of each ChASeM33 X-ray source within the PHATTER survey footprint. We identified 29 resolved background galaxies that were spatially associated with ChASeM33 X-ray sources.
There are 202 ChASeM33 X-ray sources that lie within the footprint of the PHATTER survey. Using the logN−logS relation presented in Tüllmann et al. (2011) for AGN, we expect ∼100 AGN within the area of M33 that is covered by both the ChASeM33 and PHATTER surveys. There are 28 X-ray sources identified as SNRs by Long et al. (2010) using narrow band imaging. We identify 65 HMXB candidates using the methods described in Section 3.1.
We note that while we expect there to be ∼100 background AGN in the area of the PHATTER survey in M33, we only identify 29 through visual inspection of the PHATTER HST images. This difference is likely due to the limited depth of the PHATTER imaging and difficulty detecting and resolving the disk of a galaxy seen through the disk of M33, particularly in regions with high stellar density and/or dust extinction. We also expect our quality cuts to remove background AGN with point source optical counterparts from our sample, as the optical/UV SEDs for these sources should differ in shape from massive stars in the disk of M33. We describe this series of quality cuts that we perform on our sample in more detail in Section 3.4.

SED Fitting
We used the Bayesian Extinction and Stellar Tool (BEAST; Gordon et al. 2016) software to fit SEDs for the identified point source optical counterparts. The BEAST fits observed SEDs with theoretical SEDs from the Padova/PARSEC single-star stellar evolution models (Marigo et al. 2008;Bressan et al. 2012;Marigo et al. 2017) using Bayesian methods (Gordon et al. 2016).
We use the four optical to near-ultraviolet bands from the PHATTER survey in our fitting: F275W, F336W, F475W, F814W bands, with central wavelengths 2750, 3375, 4750, and 8353Å, respectively. When performing our SED fits, we assume that all sources are in the disk of M33, at a fixed distance modulus of 24.67, or 859 kpc (de Grijs & Bono 2014), and are moving with a constant radial velocity of -179.2 km s −1 (McConnachie 2012).
The BEAST fits for parameters including age, mass, metallicity, distance, dust column density (A V ), average grain size (R V ), and f A -a parameter describing the distribution of different types of dust observed in the Local Group. From these parameters, the BEAST derives the luminosity, effective temperature, radius, and surface gravity for each HMXB companion candidate.
The BEAST imposes a Kroupa initial mass function (Kroupa 2001) as a prior on stellar mass, a uniform prior for A V , R V , f A , stellar age, and stellar metallicity (Gor-don et al. 2016). The BEAST then maps the initial mass and log(t) onto log(T ef f ) vs. log(L) and log(T ef f ) vs. log(g) diagrams to produce priors on the other stellar physical parameters-including luminosity, effective temperature, radius, and surface gravity-which reflect the expected distribution of known stars (Gordon et al. 2016). The BEAST returns the parameters for the bestfit stellar and dust model and probability distributions for each parameter. We use these probability distributions to report the 16th and 84th percentile errors for each parameter.
We allow the log(t) to range from 6 to 8 (10 6 -10 8 yr) in steps of 0.1. The metallicity grid was constrained to Z initial = [0.03, 0.019, 0.012, 0.008, 0.004, 0.001,0.0004], which covers the slightly supersolar to solar range as seen in M33 (Cioni 2009). We allow A V values ranging from 0.0 to 5.0 mag in steps of 0.25 mag. We chose this range to be consistent with the values of A V measured by Lazzarini et al. (2022) when recovering the spatially resolved recent SFH of M33, but allowing for higher A V to account for the fact that HMXBs may be in areas of their local stellar environment with higher extinction than average.
We determine the most likely spectral type for each source based on its best-fit effective temperature, luminosity and stellar radius. In Figure 3, we plot the Hertzsprung-Russell diagram for all HMXB companion star candidates using the best-fit luminosity and effective temperature. For HMXBs with only one companion star candidate, we plot the star with a diamond-shaped point. For HMXBs with multiple companion star candidates, we plot each star with a square shaped point. The error bars represent the 16th and 84th percentile values corresponding to the ±1σ. If an error bar only extends in one direction, for example if a point does not have a lower error in log(T ef f ), the best-fit value for the given parameter lies outside the 16th−84th percentile range. The size of each point scales with its best fit stellar radius. We also include background patches to show the expected ranges in effective temperature and luminosity for O-and B-type main sequence, giant, and supergiant stars (Lamers & Levesque 2017). We include isochrones from the Padova stellar models at different ages (Marigo et al. 2008). We list the ranges of effective temperature, radius, and luminosity for each type of star in Table 4.

Age Determination with Spatially Resolved SFHs
We use the spatially resolved recent SFH of M33 derived from the PHATTER optical-only photometry catalog (Lazzarini et al. 2022) to determine the age distribution and production rate for HMXBs in M33 within the PHATTER survey footprint. Lazzarini et al. (2022) measured the spatially resolved recent (back to ∼630 Myr) SFH of M33 with colormagnitude diagram (CMD) fitting. They performed CMD fitting with the software package MATCH (Dolphin 2002) using the optical-only (F475W and F814W bands with central wavelengths of 4750 and 8353Å, respectively) PHATTER catalog as input. The optical-only PHATTER catalog covers an area of ∼38 kpc 2 . To derive the spatially resolved SFH, the area of the opticalonly PHATTER catalog was divided into 2005 roughly 100 pc by 100 pc (24 ′′ on a side) regions for which the SFH was measured independently.
The SFH was measured from log(t yr −1 ) of 6.6 to 10.15 with a step size of 0.1. Their CMD fits focused on main sequence stars, optimizing for recovery of the recent SFH. They excluded older, redder stars from their fits, including the vast majority of red giant branch and red clump stars. The recent SFH is reliable back to ∼630 Myr, which more than includes the expected mainsequence lifetimes of stars with masses greater than or equal to 8 M ⊙ , which are the focus of our analysis in this paper.
The SFR measurements from the Lazzarini et al. (2022) include errors that were measured using hybrid Monte Carlo (MC) fitting and represent the range of all SFHs that could recreate the observed CMD in that spatial region of M33. As described in Section 3.2 in Lazzarini et al. (2021), this fitting can result in measurements of 0 SFR in certain time bins with upper errors, creating asymmetrical errors on SFR measurements, resulting in asymmetrical errors on our age distribution.
We determine the age distribution for our HMXB candidate sample using the spatially resolved recent SFH within the last 80 Myr in the region of the SFH maps containing each of our HMXB candidates. For the region containing each HMXB candidate, we calculate the stellar mass formed in each age bin back to 80 Myr by multiplying that bin's measured SFR by its width in years. We then divide this stellar mass formed in each time bin by the total stellar mass formed within the last 80 Myr its local region of the SFH maps to get a probability distribution for its age. We sum these probability distributions for all HMXB candidates to get an age distribution for the population, which we show for various sub-samples of our HMXB candidate sample in the plots in the left column of Figure 4.
We perform the same analysis for an equal number of regions from the Lazzarini et al. (2022)   . Hertzsprung-Russell diagram with the SED-fit derived effective temperatures and luminosities for our HMXB companion star candidates. Each HMXB companion star candidate is plotted as an individual point with errors. The shape of the point denotes the number of companion star candidates, diamond-shaped points represent HMXBs with only one companion star candidate and square-shaped points represent HMXBs with multiple companion star candidates. The size of the point scales with the best-fit radius of the source. The opacity of the point scales with its "quality" according to the evaluation criteria described in Section 3.4. The background shaded regions indicate the expected ranges in log(T ef f ) and log(Luminosity) for O and B type main sequence (MS), giant (G), and supergiant (SG) stars from Lamers & Levesque (2017), also listed in Table 4. Isochrones from the Padova stellar models (Marigo et al. 2008) are over-plotted.
for regions across the galaxy. We perform this random selection of regions 10,000 times and plot the median age distribution as a green line in Figure 4. We list the identification number of the region in the Lazzarini et al. (2022) SFH maps that contains each HMXB candidate in Table 1 so that the SFH can be retrieved for any individual HMXB. We also include the median age of the stars in the region of the SFH maps containing each HMXB candidate that formed in the last 80 Myr in Table 1. The median age is the age at which 50% of the cumulative stellar mass in that region was formed, with errors showing the ages at which 16% and 84% of the cumulative stellar mass was formed.

Determining Best HMXB Candidate Sample
One major goal in our analysis is to identify a clean sample of HMXB candidates in M33 for which both HST and Chandra data are available. We want to remove any non-HMXB sources including background galaxies, su-pernova remnants, and foreground Milky Way stars. We assess the quality of our HMXB candidates using each source's X-ray and ultraviolet/optical/infrared properties, which we discuss in more detail in the following subsections.
We use a system of flags to assess the quality of each HMXB candidate. In this system, no single flag completely removes a source from our HMXB candidate sample, but we note that it may be more likely to be a non-HMXB source. We present the values used to evaluate each source in Table 5 and list each source and the flags that it raised in Table 6.

Flag: Soft Chandra Hardness Ratios
HMXBs are known to have hard X-ray spectra, so we can use a source's hardness ratio from the Tüllmann et al. (2011) Chandra catalog to flag sources with soft hardness ratios, suggesting that they may not be HMXBs. We use the two hardness ratios defined in Plucinsky et al. (2008): To meet the criteria for a soft Chandra hardness ratios flag, a source must have HR 1 < −0.4 and HR 2 < 0.1. These cutoff values were determined using the hardness ratio diagrams in Plucinsky et al. (2008) and are intended to flag potential foreground stars and supernova remnants. We find that two of the HMXB candidates in our sample raise the soft hardness ratio flag (ChASeM33 IDs 363, 542) suggesting that they are more likely to be foreground stars or supernova remnants.

Flag: Foreground IR Colors
Infrared colors can be used to flag potential foreground stars in our sample of HMXB companion star candidates. Foreground stars create a nearly vertical sequence on the F160W versus F110W-F160W CMD with F110W-F160W colors between 0.4 and 0.8 (see Figure  19 in Williams et al. 2014).
We flag HMXB companion star candidates in our sample that have IR colors in this range, which could suggest that the optical counterparts may be foreground stars, rather than HMXB companion stars in the disk of M33. We emphasize that while a source may have an IR color in this range, that does not definitively indicate that it is a Milky Way foreground star, which is why we use this series of checks to evaluate the quality of our HMXB companion star candidates. We find that 15 stars in our HMXB companion star candidate sample raise the IR color flag. We list the F110W-F160W color for each HMXB companion star candidate in Table 5 and denote whether or not each HMXB companion star candidate raised the "IR color flag" in Table 6.

Flag: Flat SED Shape
Background AGN are known to have fairly flat spectra in the range of optical wavelengths that are covered by our photometric measurements from PHATTER (Brown et al. 2019) while we expect a star to have a more curved, black-body type spectrum. When we run the BEAST, it creates SED plots for each source showing the flux at the central wavelength of each filter included in the input photometry. Two co-authors visually inspected these SEDs and flagged sources for which the SED shape appears flat. We note that this flattening of the spectrum could be due to large amounts of dust extinction on a bright blue star and thus re-emphasize the fact that we do not rule out a HMXB companion star candidate if it raises just one flag. We find that the sources for which the SED shape appears flat tend to have higher best-fit A V values, between A V =3.6 and 4.6.
We find nine HMXB companion star candidates that have flat SED shapes, which have a value of 1 listed in the "FLAG: Flat SED Shape" column in Table 6.

Flag: No Recent SF
Due to their young ages, we expect HMXBs to be in regions of the galaxy that have had active star formation within the time scale of the lifetimes of massive stars, within the last ∼ 80 Myr. For each source, we identify the region of the M33 spatially resolved recent SFH maps (Lazzarini et al. 2022) that include the X-ray source position and flag sources where the surrounding 100 pc by 100 pc region has a best-fit SFH consistent with zero star formation in the last 80 Myr. This could indicate that our source is not actually a HMXB, which we would expect to be spatially correlated with recent SF, but could instead be a background galaxy or foreground star, neither of which should have positions correlated with SF regions in M33.
One of the HMXB candidates in our sample is located in regions of M33 with 0 measured star formation within the last 80 Myr (ChASeM33 ID 205). We list the mean SFR over the last 80 Myr for the region containing each source in Table 5.

Flag: Mis-matched NH and AV
For HMXBs in the disk of M33, we expect the absorption measured in the X-ray spectrum, N H , to correlate with the absorption we measure for the companion star at optical wavelengths, A V , because the X-ray source and optical source are located within the same binary stellar system. We can use this assumption to flag sources that may not be HMXBs in the disk of M33: i.e. sources where the N H measured for the X-ray source does not correlate with the measured A V for the optical counterpart. Not all X-ray sources we identify as HMXB candidates have measured N H values, thus we only apply this flag to sources for which this measurement is available.
We adopt N H values for sources from Tüllmann et al. (2011). They performed spectral fits for sources with a S/N of at least 2.0 in each bin and a minimum of 8 spectral bins. The best-fit model for each source is listed in Table 9 of Tüllmann et al. (2011). The N H value listed for each source in Tüllmann et al. (2011) is the intrinsic N H in M33, with a fixed Galactic N H of 0.06 × 10 22 cm −2 . We list the N H value for each of our HMXB candidates for which a fit was performed in Table 5.
We compare these N H values with measurements of A V from BEAST SED-fits of the companion star candidates and local measurements of A V that were derived by Lazzarini et al. (2022) while measuring the spatially resolved recent SFH. See Section 3.2 for a detailed description of the SED fitting. Lazzarini et al. (2022) measured the foreground extinction, A V , and differential extinction, dA V by fitting optical color-magnitude diagrams over a grid of A V and dA V values, using maximum likelihood estimation.
To compare the X-ray and optical extinction, we use published relationships between N H and A V (Predehl & Schmitt 1995;Güver &Özel 2009) to calculate the A V that would be expected for a source given the N H measured in its X-ray spectrum. We then compare this expected A V against the A V we measure with SED-fitting and from the SFH maps.
We flag sources for which the A V expected for the source based on its measured N H differs significantly from the measured A V from SED-fitting and from the SFH maps by more than 3 times the 90% confidence errors reported on the N H measurements. For sources without reported errors on the measured N H , we assumed 10% errors. We perform this calculation with both the Predehl & Schmitt (1995) and Güver &Özel (2009) correlations and find that the outliers are the same for both sets of calculations.
We only find one source that is a clear outlier for which the A V expected based on the X-ray source's N H exceeds the A V measured with SED-fitting of the companion star candidate and the A V measured for the area surrounding the source from the SFH maps. This is ChASeM33 ID 207 and its associated optical counterpart candidate. One possible explanation for this mismatch may be that the X-ray source we are observing is not an accreting black hole or neutron star in the disk of M33, but is instead a background AGN viewed through M33's disk. When we fit for the A V with SED-fitting, we assume sources are stars in M33, which could account for the mis-match. We expect that the N H measured for a background galaxy would not match the A V from the SFH maps near the position of the background galaxy on the sky.
We do not find any sources for which the A V expected based on the X-ray source's measured N H is significantly less than the A V we measure for the source with SEDfitting and the A V in the local area of M33 from the SFH maps. A mis-match of this type could suggest that the X-ray source we are looking at is a foreground Milky Way star rather than an HMXB in the disk of M33. As we described for the previous set of mis-matches, the fact that our SED-fits report a higher A V could be due to the fact that we fix the distance of the companion star candidate to M33 when we perform the fitting. Foreground stars should be very bright and red, which the BEAST could interpret as being a massive star in M33 with significant dust extinction.

Flag: Spectral Type not consistent with massive star
In Section 3.2 we describe how we inferred the most likely spectral type for our HMXB companion star candidates based on their best-fit effective temperature, luminosity, and radius. In Figure 3 we plot the values we fit for each star along with the expected ranges in effective temperature and luminosity for massive main sequence, giant, and supergiant stars, which we also list in Table 4. We flag HMXB companion star candidates for which the best-fit effective temperature, luminosity, and radius do not fit into the expected range for O or B type main sequence, giant, or supergiant stars. We find five HMXB companion star candidates that raise this flag. This could indicate that the point source for which we are performing the SED fitting is not a star in M33, which could mean that it is a foreground Milky Way star or a background galaxy viewed through the disk of M33.

Cross-Match with Previous Catalogs
We cross-match our sample with previously published X-ray catalogs of M33 (Grimm et al. 2005;Williams et al. 2015;Yang et al. 2022) and previous catalogs of HMXB candidates in M33 (Garofali et al. 2018).
We identify counterparts to our HMXB candidates in the Yang et al. (2022) NuSTAR survey of M33. NuS-TAR detected sources in that paper are listed with their counterpart IDs in the ChASeM33 survey (Tüllmann et al. 2011), which we use to identify counterparts between the two catalogs. To identify counterparts in the Williams et al. (2015) XMM-Newton catalog and the Grimm et al. (2005) Chandra catalog, we positionally cross-match sources within 5 ′′ . We list the ID number of the counterparts in the Yang et al. (2022), Williams et al. (2015), and Grimm et al. (2005) catalogs in Table  1.

Sample
We performed a detailed cross-check against the HMXB candidate sample presented in Garofali et al. (2018), which also used the ChASeM33 catalog to identify HMXB candidates using overlapping optical data from archival HST fields. Because of the non-uniform coverage of the HST fields used in their analysis, not all HMXB candidates identified by Garofali et al. (2018) fall within the PHATTER survey footprint, and thus are not included in this analysis. Garofali et al. (2018) presented 55 HMXB candidates, 40 of which lie in the footprint of the PHATTER survey footprint and 15 of which lie outside the footprint of the PHATTER survey. We identify 15 HMXB candidates in this study that fall within the footprint of the HST observations analyzed by Garofali et al. (2018) but that are not included in their catalog. Many of the HST observations used by Garofali et al. (2018) were shallower than the PHATTER data used here, which may account for some of this discrepancy. Additionally, many of the HST fields used in their analysis did not cover the full photometric range of the data we use here, and many did not have coverage at UV wavelengths, which we used to identify HMXB companion star candidates.
There are 8 HMXB candidates identified by Garofali et al. (2018) that fall within the PHATTER footprint but that are not included in our HMXB candidate sample. Two of these X-ray sources (ChASeM33 sources 210, 337) do not fit our selection criteria and did not have the same photometric coverage in the data analyzed by Garofali et al. (2018). There are 4 HMXB candidates (328,357,362,427) that were included by Garofali et al. (2018) for completeness, but that we do not include in our sample because they were classified as supernova remnants by Long et al. (2010). Source 416 is identified as an HMXB candidate by Garofali et al. (2018) but we do not include it in our sample because there is a resolved background galaxy at the location of the X-ray source in the PHATTER imaging. We do not include source 535 in our sample, although Garofali et al. (2018) included it, as it was later classified as a colliding-wind binary by Garofali et al. (2019).
The remaining 32 HMXB candidates identified by Garofali et al. (2018) that fall within the footprint of the PHATTER survey were independently identified as HMXB candidates in this analysis. In Table 1 we include a column that indicates which sources are a match between the two catalogs and list the ChASeM33 name used to identify each source in the Garofali et al. (2018) catalog for matches.

RESULTS & DISCUSSION
In this section we discuss the SED-inferred physical properties of the HMXB companion star candidates in our sample, the age distribution and HMXB production rate we measure, compare the characteristics of the M33 HMXB population with those measured for the population of HMXBs in M31, which was studied with an analogous approach in Lazzarini et al. (2021), and discuss the hard X-ray properties of HMXB candidates in our sample that were also detected in a recent NuSTAR survey of M33 (Yang et al. 2022).

SED-inferred properties of HMXB Companion
Star Candidates in M33 As described in Section 3.2, we classify HMXB companion star candidates as O-or B-type main sequence, giant, or supergiant stars based on their best-fit effective temperature, luminosity, and radius. There are several known classes of HMXBs, based on the nature of the companion star and the accretion mechanism. Supergiant HMXBs (sgHMXBs) are systems where the black hole or neutron star accretes material from a supergiant companion star, either via stellar winds or Roche lobe overflow. In Be X-ray binaries (Be-XRBs), the companion star is an Oe/Be star, which is rapidly rotating and forms an equatorial decretion disk that exhibits characteristic Hydrogen emission lines in the star's spectrum. Accretion onto the compact object occurs when the compact object, typically a neutron star, passes through or very close to this decretion disk. The Oe/Be companion stars in Be-XRBs have been observed to have spectral types ranging from late O-type stars (O9) to early Btype stars (see Reig 2011, for a review).
Classifying the companion star type in HMXB systems is only possible for nearby galaxies and the Milky Way, where optical spectroscopy is readily available. In the Milky Way sgHMXBs make up roughly 30% of the observed HMXB population, Be-XRBs make up about 50%, and the remaining 20% is a combination of giant, main sequence, and unknown companion star types. In contrast, in the Magellanic clouds, out of over 120 HMXB systems, all but two are known to contain Be companion stars (Liu et al. 2006;Maravelias et al. 2014;Walter et al. 2015;Haberl & Sturm 2016;Yang et al. 2017b).
We can attempt to compare the breakdown of likely spectral types of the HMXB companion stars in our observed sample with the populations in the Milky Way and Magellanic Clouds. However, we note that this is challenging because 23 of our HMXB candidates have more than one candidate companion star.
In systems with only one potential HMXB companion star candidate, we find that 31 (76%) are likely B-type main sequence stars and 2 (5%) are likely O-type main sequence stars. Of the remaining systems, 1 is a likely B-type giant, 2 are likely O-type giants, and 2 are likely B-type supergiants.
In systems with more than one candidate companion star, we find that 14 (61%) only have B-type main sequence stars as potential companions (i.e. all candidate companion stars are likely B-type main sequence stars). The remaining systems with multiple candidate companion stars have a mix of O and B-type main sequence, giant, and supergiant stars.
Because there are multiple candidate companion stars for many of our HMXB candidates, we cannot state a firm fraction of the sample that host different spectral type companion stars. However, because most of the systems contain a likely O-or, more frequently, B-type main sequence companion star, our findings suggest that most of the HMXBs in M33 are Be-XRBs. Future optical spectroscopy is required to confirm the true spectral type of the companion stars in these systems, but these initial findings suggest that the fraction of Be-XRBs for M33 lies somewhere between the fraction in the Milky Way and the fraction in the Magellanic Clouds. Metallicity has been suggested to affect the fraction of giant/supergiant HMXBs to Be-XRBs in a population of HMXBs, due to lower rates of mass loss via line driven winds at lower metallicities (Linden et al. 2010).
We note that our SED fits use single stellar models and these models do not include a component to model the potential decretion disk around Be stars. Be stars exhibit infrared excesses compared to normal O/B type stars. This effect is orientation dependent, and is most significant when the system is viewed pole-on. This excess has a small (<0.5 mag) effect at optical wavelengths but can dominate at infrared wavelengths (e.g. Rivinius et al. 2013). We include the two optical and two ultraviolet bands from the PHATTER photometric catalog in our SED-fits, so we do not expect infrared contamination from the decretion disk to affect our measurements. There is a chance that we get contamination in the ultraviolet bands from the compact object's accretion disk, which may bias our SED fits towards higher masses. This is a challenging system to disentangle with purely photometric data, and we plan to tackle this with forthcoming spectroscopic measurements (Lazzarini et al., in preparation).

Age Distribution of HMXBs in M33
As described in Section 3.3, for each HMXB candidate, we used the SFH within the last 80 Myr in the 100 pc by 100 pc region surrounding the source to calculate a probability distribution for its age. We performed this calculation for all HMXB candidates in our sample and then sum their probability distributions to create an age distribution for the population. The subsamples represented in each row of Figure 4 reflect the HMXB candidates that raise any number of flags, at most one flag, and zero flags in the top, middle, and bottom rows, respectively. We include the known nuclear HMXB (X-8) in our best HMXB candidate sample although we cannot identify its optical counterpart due to the high density in the nuclear cluster of M33, which brings the total number of HMXBs in our full sample to 65.
We present this age distribution, for different subsamples, with the black histogram in the plots in the left column of Figure 4. The left plot in the top row of Figure 4 shows the age distribution for our full sample of HMXB candidates (65 sources), the middle row shows the sub-sample that raised one or fewer flags (63 sources), and the bottom row shows our highest quality sub-sample that raised 0 flags (53 sources), using the evaluation criteria described in Section 3.4. We include error bars on the age distribution for each sub-sample. We calculate the error bars by generating a PDF for each source from a random sampling within the errors on its SFR measurements. We performed this random selection for all sources in the sub-sample 10,000 times and present the 16th and 84th percentile of the age distribution derived from this sampling as the lower and upper errors on the age distribution, respectively.
On each age distribution plot, we also include a green line which shows the results when we perform the same analysis on an equal number of randomly selected regions from the SFH maps that have a measured SFR>0 M ⊙ yr −1 over the last 80 Myr. We run this random selection 10,000 times and the green line represents the median age distribution, which is a proxy for the overall SFH in M33.
As shown in the bottom row, left panel of Figure 4, we find more HMXB candidates associated with star formation between 0 and 10 Myr ago and between 40 and 50 Myr ago than would be expected based on the green random line, which represents sampling from the SFH maps agnostic of the positions of HMXB candidates. This was also seen in a previous study of the HMXB population of M33 by Garofali et al. (2018) where the age distribution of HMXBs in M33 was measured using SFHs derived with CMD-fitting of archival HST observations. Garofali et al. (2018) suggested that these peaks repre-sent two different HMXB formation channels: a prompt channel (0-10 Myr) and a more delayed formation channel (40-60 Myr).
Different types of HXMBs have been shown to form via different channels, with associated timescales (e.g.; Linden et al. 2009Linden et al. , 2010. The 40−50 Myr peak in the HMXB production rate that we observe is the expected formation timescale for Be-XRBs. The 0−10 Myr peak in the HMXB production rate we observe is characteristic of more massive binaries, where the primary star forms a compact object and starts accreting within a few Myr. Previous studies in other galaxies including the SMC (Antoniou et al. 2009), NGC 300, and NGC 2403 (Williams et al. 2013), found that the preferred age of HMXBs in those galaxies is ∼40−55 Myr. In M31 the preferred ages for HMXBs is ∼10−50 Myr Lazzarini et al. 2018Lazzarini et al. , 2021, a broader range that overlaps with the findings in the SMC, NGC 300, and NGC 2403. In the LMC, by contrast HMXBs are associated with younger star formation episodes. The SFH is peaked at ∼6.3 Myr in the regions surrounding sgHMXBs and black hole HMXBs, ∼12.6 Myr in the regions surrounding Be-XRBs and X-ray pulsars, and ∼25 Myr in the region surrounding a candidate white dwarf-Be X-ray binary (Antoniou & Zezas 2016). Our findings suggest that we are observing an HMXB population in M33 that is dominated by two different formation channels with timescales characteristic of Be-XRBs (40−50 Myr) and more massive HMXBs (0−10 Myr).

Masses and Ages
To better understand formation channels for our observed HMXB population, we explored the relationship between the mass of the HMXB companion star candidates that we infer from SED-fitting as described in Section 3.2 and the ages of the systems measured using the SFH maps from Lazzarini et al. (2022).
In Figure 5 we plot the HMXB companion star candidate mass against two measures of likely age for the system. In both panels points are color-coded according to the likely spectral type of the companion star. Diamonds indicate HMXBs for which there is only one companion star candidate, while squares represent HMXBs with multiple companion star candidates. The opacity of each point scales with its quality as defined in Section 3.4. The gray line in the background shows the maximum stellar mass expected at each age in single stellar evolution, using Padova isochrones (Marigo et al. 2008). We also include outlines around points for HMXB candidates associated with a compact object with a classification based on its NuSTAR colors and luminosity, as described in Yang et al. (2022) and listed in Table   1. Circles indicate sources that are likely black holes (in any accretion state) and stars indicate sources that are likely neutron stars (either pulsars or lower magnetic field neutron stars).
In the right panel we plot the median age of stars that formed in the last 80 Myr in the 100 pc by 100 pc region of the SFH maps that contains each HMXB candidate. The median age is defined as the age at which 50% of the total stellar mass in the region has formed, and we include error bars that extend to the ages at which 16% and 84% of the stellar mass was formed, respectively. The median age is not generally a good representation of the likely age of the HMXB candidate because of contamination from older stars that are not necessarily associated with the star formation event that formed the HMXB candidate. As Figure 4 shows, some regions may have both an early and a later peak in their SFH, which would cause the median age to be somewhere in between these two peaks, when the SFR was actually lowest and the systems are less likely to have formed.
We looked at the age of the most recent significant star formation to try and remove contamination from older star formation episodes in the surrounding region that are not associated with the star formation episode that produced the HMXB. We defined significant star formation as a SFR that would produce 1000 M ⊙ during the 10 Myr time bin, which assuming a Kroupa (2001) initial mass function, should produce ∼10 stars with M>8 M ⊙ . This threshold SFR is greater than the noise in the SFR measurements in the SFH maps, which is ∼1.0×10 −5 M ⊙ yr −1 after re-binning the spatially resolved recent SFH maps from Lazzarini et al. (2022) into 10 Myr time bins.
As is shown in Figure 5, the age of the most recent significant star formation correlates with the mass of the companion star derived via SED-fitting, although these two measurements are completely independent. We see that systems containing our most massive HMXB companion star candidates are found in regions of M33 with very recent (between 0 and 20 Myr ago) significant SF. We find that for HMXB companion star candidates reside in regions with less recent episodes of significant star formation, the masses of the companion stars are lower.
There are two outliers from this relationship with ages of most recent significant star formation older than expected given the SED-derived mass of the HMXB companion star candidate. These systems could indicate an extension of the stellar lifetime due to binary interactions (e.g., Zapartas et al. 2017) or a system with an abnormally high velocity, but additional characterization of these HMXB systems is needed. We have an ongoing program to obtain optical spectroscopic observations of the M33 HMXB companion star candidates presented in this paper, which is needed to confirm the spectral types, and thus masses, of the HMXB companion star candidates (M. Lazzarini, et al., in preparation).

HMXB Production Rate
The production rate of HMXBs in M33, or the number of HMXBs produced per unit stellar mass or unit SFR, can be compared with measurements from other galaxies to explore the relationship between HMXB production and other galaxy-wide properties such as metallicity, stellar mass, and star formation rate. We present two measures of HMXB production rate in the following subsections. First we present the HMXB production rate measured in the number of HMXBs produced per unit SFR (HMXBs/(M ⊙ yr −1 )). We also measure a time-resolved HMXB production rate, or the number of HMXBs produced per stellar mass, over the last 80 Myr. We provide more details on each calculation in the following subsections and discuss the HMXB production rate in M33 in comparison with other galaxies for which this measurement has been done.
Our HMXB production rate measurements are a reflection of the number of HMXBs that were active at the time of the ChASeM33 survey Chandra observations. This number is impacted by the duty cycle for HMXBs, or the fraction of time that they spend in the active HMXB phase. Empirical studies of mostly Galactic HMXBs have revealed duty cycles ranging from around 10% (for transient Be-XRBs) (Sidoli & Paizis 2018) to greater than 10% (for SG HMXBs), which means that our HMXB production rates are a lower limit on the true HMXB production rate because a significant fraction of HMXBs in M33 were likely not active at the time of observation.
We calculate the HMXB production rate per unit SFR over 50 and 80 Myr timescales. 50 Myr is the expected maximum age for massive stars from single stellar evolution, however this lifetime may be extended for stars in binary systems due to rejuvenation from mass transfer via binary interactions (e.g., Zapartas et al. 2017).

Calculation of HMXBs per Unit SFR
We can use calculations of the HMXB production rate per unit SFR to compare with measurements in other galaxies where time-resolved SFH measurements are not available. We calculate the number of HMXBs formed per unit SFR (HMXBs/M ⊙ yr −1 ) over the last 50 and 80 Myr, to account for the lifetimes of single massive stars and potential extension of that lifetime due to binary stellar interactions.
To calculate the number of HMXBs formed per unit SFR, we integrate the age distribution of HMXBs shown in the left column in Figure 4 over 50 and 80 Myr. We then divide this number of HMXBs by the mean SFR over the full PHATTER footprint over the same timescale. The whole PHATTER footprint is covered by the ChASeM33 survey, which means that we are sensitive to detecting HMXBs over its full area. We obtain a lower and upper limit on this measurement from the error bars on the SFR measurements. We present the HMXB production rates per unit SFR for various subsamples in Table 7.
We aim to compare this HMXB production rate with rates measured in other galaxies. However, this direct comparison is difficult due to differences in the measurement of SFR in these galaxies and the way that HMXB samples were generated including the sensitivity of Xray observations. Due to their proximity, M33 and M31 both has spatially resolved SFH maps. The LMC and SMC also have spatially resolved SFH maps (Harris & Zaritsky 2004, 2009a) that have been used to associate their HMXBs with star formation episodes, and thus ages and production rates ( Harris & Zaritsky (NGC 55;2009b). They used the number of known HMXB candidates in each galaxy and the mean SFR of each galaxy to calculate their global HMXB production rates. These HMXB production rates show a general trend in increasing HMXB production rate with decreasing metallicity. The lowest metallicity galaxies in their sample, NGC 55 and the SMC, have HMXB production rates of 181−240 and 240−880 HMXBs/(M ⊙ yr −1 ). They measure an HMXB production rate of 96−256 and 52−86 HMXBs/(M ⊙ yr −1 ) in the LMC and Milky Way, respectively. For our full HMXB candidate sample, we measure an HMXB production rate of 181−240 HMXBs/(M ⊙ yr −1 ) over the last 80 Myr, which places the HMXB production rate of M33 higher than that of the Milky Way, similar to the HMXB production rate of the LMC, and lower than the measured rate in the SMC and NGC 55.

Rate
We can also calculate a time-resolved measure of the HMXB production rate in M33 using the measured HMXB age distribution and the SFH of M33 in the PHATTER footprint. The age distribution indicates the number of HMXB candidates that likely formed in each 10 Myr time bin. To derive the number of HMXBs produced per stellar mass, we must first calculate the total stellar mass formed in each 10 Myr time bin in the area covered by both the PHATTER and ChASeM33 surveys in M33. We find the total stellar mass formed in each time bin by multiplying the SFR in each time bin by its width, 10 Myr. We can then divide the HMXB age distribution by the stellar mass formed in each time bin to produce a measurement of the number of HMXBs formed per stellar mass. We present this time-resolved HMXB production rate in the right column of Figure  4, with the full sample of HMXBs on the top row, the sample of HMXBs which raise one or fewer flags in the middle row, and our highest quality sample of HMXBs in the bottom row. The green line in the right column of Figure 4 indicates the same analysis performed with an equal number of randomly selected regions from the Lazzarini et al. (2022) SFH maps with measured SFR over the time range analyzed (mean SFR over the last 80 Myr >0).

Comparison with M31 HMXB Population
We can compare the HMXB populations of M33 with the HMXB population of M31 to understand how differences between their histories of star formation are reflected in the HMXB populations we observe today. Previous studies of the HMXB population of M31 used the PHAT dataset and overlapping Chandra-PHAT survey ) to identify HMXB candidates, characterize the companion stars with SED-fitting, measure their age distribution and HMXB production rates with spatially resolved recent SFH maps Williams et al. 2018;Lazzarini et al. 2021). The PHAT and PHATTER surveys were designed to be very similar, making the population demographics of their HMXB populations inferred from these surveys directly comparable. Both the Lazzarini et al. (2021) study of M31 and this one use Chandra observations with similar depths.
We find that the HMXB populations in M31 and M33 are remarkably similar. Despite known differences in the area of the surveys and variation in the SFH of each galaxy, we find a roughly equal number of HMXB candidates in each study. We identify 65 HMXB candidates in our full sample in M33 and Lazzarini et al. (2021) identified 57 HMXB candidates in M31 using very similar methodology. While at first glance, this may seem surprising, upon closer inspection this similarity is supported by known correlations between SFR, metallicity and the X-ray luminosity function.
While M33 is known to exhibit a decreasing metallicity gradient with radius, at the radii probed by the PHATTER survey in M33 and the PHAT survey in M31, the gas-phase metallicities of these galaxies are comparable. The gas-phase metallicity of M33 within the PHAT-TER footprint ranges from 12 + log(O/H) ∼ 8.4 − 8.7 while within the M31 PHAT survey footprint, the gasphase metallicity ranges between 8.6 and 9.0 , and references therein). While the luminosity of HMXBs is known to scale with metallicity, this effect is dominated by high luminosity X-ray sources (L 0.5−8.0keV > 10 38 erg s −1 ) and below this limit, the Xray luminosity function has been shown to be a nearly universal power-law distribution (Lehmer et al. 2021). There is only one identified HMXB in M31 and M33 with a measured L 0.5−8.0keV > 10 38 erg s −1 -M33 X-8 -which suggests that the difference we see in the HMXB populations cannot be explained by metallicity alone.
The total SFR measured in the areas of M31 and M33 covered by PHAT/Chandra-PHAT and PHAT-TER/ChASeM33, respectively, are remarkably similar, despite covering significantly different de-projected areas. The area covered by overlapping PHAT and Chandra-PHAT observations in M31 is ∼328 kpc 2 , which is almost 9 times larger than the ∼38 kpc 2 area covered by the overlapping PHATTER and ChASeM33 surveys in M33. We summarize the mean SFR in the areas of M31 and M33 covered by both HST and Chandra, used in this analysis and Lazzarini et al. (2021) in Table 8 over the last 50, 80, and 100 Myr. M33 has a slightly higher SFR on all timescales.
It is also interesting to contrast the time resolved HMXB production rate measurements. In Figure 6 we plot the time-resolved HMXB production rates for the full sample of HMXB candidates in M31 and M33 on the same axes. See Section 4.3 for a detailed description of how the time resolved HMXB production rates were calculated in both galaxies. Briefly, to calculate the time resolved HMXB production rate the number of HMXBs expected to have formed in each time bin, from the population age distribution, is divided by the total SFR of the galaxy within the analysis region during the same time span. Because the total stellar mass formed in the region of each galaxy covered by HST and Chandra surveys is used in the calculation of the timeresolved HMXB production rate, and the total stellar mass is derived from the SFH, differences in the SFH of M31 and M33 should not affect the time resolved HMXB production rates. The HMXB production rates in both galaxies agree within errors in most time bins, but M33   . Age distributions and time-resolved production rates for our HMXB sample. Each row represents a different subsample based on the evaluation criteria described in Section 3.4, the top row is our full sample, the middle row represents HMXB candidates that raise one or fewer of our quality flags, and the bottom row shows our highest quality HMXB candidate sample. The black histograms in the left column represent the age distribution for each sub-sample and the green line represents a random sample. For more detail, see Section 4.2. In the right column, the histograms represent the time-resolved HMXB production rate for each sample and the green line represents a random sample for comparison. See Section 4.3 for more details.  Table 1.
On the y-axis we plot the mass of the HMXB companion star candidate derived via SED-fitting. On the x-axis we plot the age bin (i.e. we plot a point at 5 Myr for an age bin of 0−10 Myr) of the most recent significant star formation in the SFH map region including the HMXB candidate. We define significant star formation in more detail in Section 4.2.1. Points are color-coded according to the likely spectral type of the companion star. Diamonds indicate HMXBs for which there is only one companion star candidate, while squares represent HMXBs with multiple companion star candidates. The opacity of each point scales with its quality as defined in Section 3.4. The gray line in the background shows the maximum stellar mass expected at each age in single stellar evolution. We outline HMXB candidates with compact object classifications from NuSTAR observations, described in Section 4.5. Right: We plot the HMXB companion star candidate masses on the y-axis and the median age of stars formed in the last 80 Myr in the SFH map region including the HMXB companion star candidate. Point colors and shapes follow the same conventions as the left panel. has a slightly higher HMXB production rate at ages beyond ∼50 Myr.

Hard X-ray Properties of M33 HMXB Candidates
Comparing the optical and hard X-ray properties of HMXB candidates can be useful for understanding the binary evolution that led to the systems we observe today. Most of the PHATTER survey footprint has also been covered by a hard X-ray survey of M33 with NuSTAR (Yang et al. 2022). In their catalog, Yang et al. (2022) identified the counterpart for each NuSTAR source in the ChASeM33 survey catalog. Yang et al. (2022) detected 28 hard X-ray sources in M33, 7 of which are HMXB candidates identified in this paper. Due to the large point spread function (PSF) of NuSTAR compared to Chandra, there are several NuSTAR-detected sources with multiple associated ChASeM33 sources. Two of our HMXB candidates (ChASeM33 sources 321 and 398) are associated with NuSTAR sources for which there are multiple ChASeM33 counterparts, thus the NuSTAR source may be a blend of multiple ChASeM33 sources.
Hard X-ray observations with NuSTAR can be used to constrain the type and/or accretion state of an accreting compact object. Due to spectral differences between accreting black holes and neutron stars at hard X-ray wavelengths (E>10 keV), the hardness ratios and luminosities of accreting compact objects of unknown type can be used to classify them by comparing to Galactic compact objects for which the compact object type/accretion state is known (e.g., Wik et al. 2014;Yukita et al. 2016;Vulic et al. 2018;Lazzarini et al. 2018Lazzarini et al. , 2019Yang et al. 2022). We list the source ID and compact object classification for HMXB candidates in our sample that were detected in the Yang et al. (2022) NuSTAR survey of M33 in Table 1.
While the number of HMXB candidates for which NuSTAR detections and classifications are available is small, we observe some trends when we compare the likely compact object type inferred with NuSTAR and the likely spectral type of the companion star candidate based on its optical/UV photometry. We find that of the 7 matches, the 3 of the 4 that are classified as black holes, soft or intermediate accretion state, have a likely giant or supergiant companion star candidate. The fourth has a best fit SED-derived luminosity, effective, temperature, and radius that is consistent with being an B supergiant within errors (source 299). One of the matched sources is classified as a ULX and the companion star candidate for this source is a likely Btype main sequence star. The remaining two matches are both classified as neutron stars, one as a pulsar and one as Z-type, and both have likely main sequence type companion stars.
It is interesting that we find supergiant/giant companion star candidates associated with compact objects classified as black holes based on their NuSTAR colors and luminosities but not neutron stars. Supergiant/giant companion stars have been observed in HMXB systems in other galaxies with both black hole and neutron star accretors. Black hole HMXBs (BH-HMXBs) such as Cyg X-1 (Orosz et al. 2011), LMC X-1 (Orosz et al. 2009), and M33 X-7 (Orosz et al. 2007;Ramachandran et al. 2022, which is included in our sample as Source ID 225), all comprise black holes accreting from supergiant companion stars with masses ≳30 M ⊙ , suggesting they are very young systems because 30 M ⊙ stars have maximum lifetimes ≤ 10 Myr. Neutron star HMXBs (NS-HMXBs) with supergiant companion stars are also thought to be very young systems. The companion stars in these systems have been observed with masses ≳15 M ⊙ and companion star ages of <12 Myr (Reig et al. 2016, and references therein). It is important to emphasize that the classification of the compact object type based on its NuSTAR colors and luminosity is not conclusive because we do not have dynamical mass measurements for the accreting compact objects in these HMXB systems. More detailed spectroscopic observations are needed to 1. confirm the true spectral types and masses of the companion stars and 2. constrain the compact object masses dynamically. We also do not wish to over interpret this result due to the potential impact of small number statistics in our sample with 4 likely black hole supergiant HMXBs.

CONCLUSIONS
We present a study of the HMXB population of M33 using a combination of high resolution X-ray observations from Chandra (ChASeM33 Tüllmann et al. 2011) and near-IR/optical/near-UV observations from HST (PHATTER Williams et al. 2021). We performed SEDfitting on the HMXB companion star candidates to infer their likely spectral types. We used the spatially resolved recent SFH of M33 (Lazzarini et al. 2022) to measure the HMXB age distribution and production rate for M33. We discuss our findings in context of previous studies of HMXB populations in nearby galaxies and previous work in M33. We list our main findings below: • We identify 65 HMXB candidates in M33, including X-8 in the nuclear star cluster, using multiband photometry from the PHATTER survey of the inner disk of M33.
• Using SED-fitting, we infer that the majority of our HMXB companion star candidates are likely B-type main sequence stars, suggesting that a majority of HMXBs in M33 are Be-XRBs, although spectroscopic observations are required to confirm these classifications.
• We find that the production rate for HMXBs in M33 is peaked between 0 and 10 Myr ago and between 40 and 50 Myr ago. The 0−10 Myr formation timescale is associated with BH-HMXBs with massive companion stars and the ∼40 Myr formation timescale is associated with Be-XRBs, we detect signals of these two formation channels in the time-resolved HMXB production rate of M33. • We calculated the time-resolved HMXB production rate (the number of HMXBs formed per stellar mass) in M33 and find that it ranges from ∼ 3 × 10 −6 M −1 ⊙ to ∼ 1 × 10 −6 M −1 ⊙ over the last 80 Myr.
• We also calculate the number of HMXBs formed per unit SFR over the last 50 and 80 Myr. For our full HMXB candidate sample, we measure an HMXB production rate of 120−153 HMXBs/(M ⊙ yr −1 ) over the last 50 Myr and 181−240 HMXBs/(M ⊙ yr −1 ) over the last 80 Myr.
• We find that the HMXB production rates (HMXBs M −1 ⊙ ) of M33 and M31 agree within errors over the last 50 Myr and only differ slightly between 50 and 80 Myr ago.
• We find that for the 7 HMXB candidates in our sample with compact objects classified as black holes or neutron stars based on NuSTAR observations (Yang et al. 2022

B-MS
Note-We present X-ray and optical information for each HMXB candidate in this study. 23 of the X-ray sources from Tüllmann et al. (2011) have multiple companion star candidates. For X-ray sources with multiple candidate companion stars, the X-ray source properties are listed once, and the star properties are each listed in a separate row. Column descriptions follow, in order: T11 ID, T11 R.A., T11 Dec., and T11 ePos are the ID number, R.A., Dec. and error on the X-ray source position presented in Tüllmann et al. (2011). The PHATTER Name is the position of the HMXB companion star candidate in the PHATTER catalog . The CP Sep. column denotes the separation between the PHATTER star's position and the ChASeM33 X-ray sources position, in arcseconds. We list each source's 0.35−8.0 keV luminosity, derived from its count rate in the ChASeM33 catalog (Tüllmann et al. 2011). The log(T ef f ), log(L), R, A V , M, and SED Age columns come from the SED fits performed for each HMXB companion star candidate, described in Section 3.2 and list the star's best-fit effective temperature, luminosity, radius, dust extinction, mass, and age. We include the median age of stars formed in the last 80 Myr in the region from the Lazzarini et al. (2022) SFH maps of M33 within which each HMXB resides. We also include the number of the region in the Lazzarini et al. (2022) SFH maps that contains each HMXB candidate for ease in retrieving the full SFHs.
We include the most likely spectral type for each companion star candidate, based on the methodology described in Section 3.2. Lastly, we present the source's ID and classifications from previous X-ray surveys of M33. We present the classification for each source listed in the Tüllmann et al. (2011, T11). For sources that also appear in the Garofali et al. (2018, G18) HMXB catalog, we list the source name in their catalog. We list the source ID and classification from the Yang et al. (2022, Y22) NuSTAR survey of M33. Yang et al. (2022) classified NuSTAR-detected sources based on their hard X-ray colors and luminosities. The abbreviations are as follows: ULX=ultra-luminous X-ray source, S=soft accretion state black hole, I=intermediate accretion state black hole, H=hard accretion state black hole, P=pulsar, Z=Z-type neutron star. For sources in Yang et al. (Y22; that were consistent with multiple compact object types, the most likely type is listed first, with other possible types listed in parentheses. We present the source ID for sources in the XMM-Newton survey of M33 by Williams et al. (W15; and in the Grimm et al. (G05; Table 3 continued Note-Color and magnitude cuts used to select massive main sequence (MS) and giant/supergiant (G/SG) stars within the 1.5σ positional errors of ChASeM33 X-ray source positions. We plot the stars selected using these criteria on a CMD in Figure 2.  Table 3 continued Note-Multi-band photometry for each identified HMXB companion star candidate. We list the name of each star (derived from its position) in the PHATTER catalog and the ChASeM33 ID number of the X-ray source each HMXB companion star candidate is associated with for ease in cross-referencing Table 1. Photometry comes from the full PHATTER photometric catalog presented in Williams et al. (2021). Note-List of the ranges of effective temperatures, luminosities, and radii used to classify companion stars. We plot these ranges as colored patches on the Hertzsprung-Russell diagram in Figure 3. Values come from Appendix B of Lamers & Levesque (2017).  Table 5 continued   Table 5 continued

0.26
Note-We present the multi-wavelength properties used to evaluate the quality of each HMXB candidate. For a full discussion of how we determined the highest quality HMXB candidate sample, see Section 3.4. The Chandra hardness ratios listed here were calculated using count rates from Tüllmann et al. (2011). The F110W-F160W color is calculated using the HST photometry for each HMXB companion star candidate, listed in Table 3. We include the mean SFR over the last 80 Myr, calculated from the recent SFH maps of M33 generated with the PHATTER optical photometry, presented in Lazzarini et al. (2022). Spectral fits were performed on X-ray sources in the ChASeM33 catalog with sufficient counts, and we present the resulting best-fit N H for these sources from Tüllmann et al. (2011). The SED-derived A V value for each HMXB companion star candidate comes from our SED fits, and is also listed in      Table 6 continued  Note-For each identified HMXB companion star candidate, we list whether it raised each of our flags used to evaluate the quality of HMXB candidates. The values used to determine the flags are listed in Table 5 and a discussion of criterion is in Section 3.4. Note-HMXB production rates calculated per unit SFR for our HMXB candidate sample. The rows correspond to subsamples based on our evaluation system described in more detail in Section 3.4. The samples used for the calculations in each row match the samples in Figure 4. Note-Mean SFR over various timescales in M31 and M33 in the areas of each galaxy covered by overlapping HST and Chandra observations. The SFR values come from spatially resolved recent SFHs of each galaxy (Lewis et al. 2015;Lazzarini et al. 2022).