The High-mass X-Ray Binary Luminosity Functions of Dwarf Galaxies

Drawing from the Chandra archive and using a carefully selected set of nearby dwarf galaxies, we present a calibrated high-mass X-ray binary (HMXB) luminosity function in the low-mass galaxy regime and search for an already hinted at dependence on metallicity. Our study introduces a new sample of local dwarf galaxies (D < 12.5 Mpc and M * < 5 × 109 M ⊙), expanding the specific star formation rates (sSFR) and gas-phase metallicities probed in previous investigations. Our analysis of the observed X-ray luminosity function indicates a shallower power-law slope for the dwarf galaxy HMXB population. In our study, we focus on dwarf galaxies that are more representative in terms of sSFR compared to prior work. In this regime, the HMXB luminosity function exhibits significant stochastic sampling at high luminosities. This likely accounts for the pronounced scatter observed in the galaxy-integrated HMXB population’s L X/SFR versus metallicity for our galaxy sample. Our calibration is necessary to understand the active galactic nuclei content of low-mass galaxies identified in current and future X-ray survey fields and has implications for binary population synthesis models, as well as X-ray-driven cosmic heating in the early Universe.


INTRODUCTION
High-mass X-ray binaries (HMXBs) are binary systems that consist of a compact object, a black hole or neutron star, that is accreting material from a high-mass companion star (M * > 10M ⊙ ).Because of the relatively short lifetimes of the high-mass stars in HMXBs, the total X-ray luminosity of HMXBs serves as a powerful tool for probing star formation rates (SFRs) in galaxies (e.g., Hornschemeier et al. 2000;Grimm et al. 2003;Ranalli et al. 2003;Persic et al. 2004;Lehmer et al. 2010;Mineo et al. 2012;Basu-Zych et al. 2013;Lehmer et al. 2016;Fornasini et al. 2018;Saxena et al. 2021).In addition, there is growing evidence suggesting that HMXBs in dwarf galaxies in the early universe may have contributed significantly to the ionizing radiation during the preheating of the intergalactic medium leading up to the reionization epoch (Warszawski et al. 2009;Madau & Fragos 2017a;Eide et al. 2018).The study of HMXBs also provides an important constraint on binary star and compact object evolution, with important implications for gravitational wave sources (Podsiadlowski et al. 2003;Abbott et al. 2016;Liotine et al. 2023).
For galaxies that are actively forming stars, the HMXB X-ray luminosity function has been observed to predominantly correlate with the overall SFR of the host galaxy as expected (Grimm et al. 2003;Gilfanov 2004;Mineo et al. 2012;Lehmer et al. 2019, L19 hereafter).However, recent observations and theoretical work suggest that the HMXB luminosity function may also depend on factors such as metallicity and star formation history (SFH).Various models have been proposed to describe the dependence of the luminosity function on metallicity (Brorby et al. 2014;Basu-Zych et al. 2016;Ponnada et al. 2020), primarily interpreting the L X /SFR versus metallicity relation.Lehmer et al. (2021) (L21 hereafter) in particular observe that this metallicity dependence causes an excess of sources above 10 38 erg s −1 for low-metallicity galaxies and introduce a framework for modeling the HMXB X-ray luminosity functions as a function of SFR and metallicity.
Thanks to their low metallicities (Tremonti et al. 2004;Mannucci et al. 2010), dwarf galaxies serve as a valuable tool to examine the impact of low-metallicity environments on the HMXB X-ray luminosity function.In dwarf galaxies, another consideration becomes highly relevant.HMXBs are the dominant contributors to the overall X-ray luminosity of any star-forming galaxy without an active galactic nucleus (AGN).In dwarf galaxies, however, the contribution of HMXBs may dominate over that of an active nucleus as well, since dwarf galaxies may harbor very low-mass central black holes (Mezcua 2017;Greene et al. 2020).The X-ray emission from HMXBs becomes a dominant source of con-fusion in the context of intermediate-mass black hole (IMBH) detection in distant dwarf galaxies, where we cannot resolve individual HMXBs with current telescopes (Schramm et al. 2013;Pardo et al. 2016;Mezcua et al. 2018;Halevi et al. 2019).Without a better understanding of the full range of L X per unit star-formation rate from dwarfs, it is not possible to robustly identify AGN candidates using X-rays; while there is a hope to detect and characterize seed black holes in the early universe using next-generation X-ray missions (Natarajan et al. 2017;Barrow et al. 2018;Ricarte & Natarajan 2018;Haiman et al. 2019).A better handle on the role of metallicity in setting the HMXB luminosity function is urgently needed for this purpose.
In this work, we seek to understand this possible metallicity dependence more securely by increasing the sample of dwarf galaxies with measured HMXB luminosity functions.Prior studies suggest that local dwarf galaxies exhibit a metallicity-driven surplus of high-L X sources, consequently affecting the shape of their Xray luminosity functions (Mapelli et al. 2010;Prestwich et al. 2013;Brorby et al. 2014;Douna et al. 2015;Kovlakas et al. 2020).In particular, there appears to be an excess of so-called "ultra-luminous X-ray" (ULX) sources in low-mass and low-metallicity galaxies.ULXs are X-ray point sources with luminosities that exceed the Eddington limit for stellar-mass black holes (10 39 erg s −1 ).Fully characterizing the HMXB X-ray luminosity function for a larger and more representative sample of dwarf galaxies will quantify whether low metallicity leads to this apparent ULX excess.Further, as shown by Fornasini et al. (2020), the metallicity dependence of L X /SFR is not expected to evolve with redshift, allowing us to consider local dwarfs as useful analogs for distant dwarf galaxies.
Past studies have often focused on a limited sample of dwarf galaxies, predominantly favoring those with high specific star-formation rates and the lowest known metallicities (Prestwich et al. 2013).L21 in particular relies on these galaxies to constrain their low-metallicity luminosity functions, but leave a notable gap between 12 + log(O/H) ≈ 7.6 − 8.0 range as a consequence.We fill this gap by broadening our sample to include a wider metallicity range of 12 + log(O/H) = 7.74 − 8.77, thus encompassing typical dwarfs with M * ∼ 10 8 − 10 9 M ⊙ and log sSFR ∼ −10.7 − −8.5.In addition to being targeted by the Chandra X-ray Observatory, these galaxies also benefit from supporting data from GALEX and HST, enabling us to define uniform apertures and obtain reliable measurements of their SFRs.
The structure of this paper is as follows: In Section 2, we discuss our selected sample of galaxies.Section ) and galaxies introduced in this study ("Chandra+HST").Chandra+HST are reprinted by stars while L21 galaxies are marked by square points.The points are colored according to gas-phase metallicity (12 + log(O/H)).Our dwarf galaxy sample supplements the L21 sample in M * , SFR, and metallicity while maintaining a comparable sSFR.The gray circles denote galaxies that satisfy the mass (M * < 5 × 10 9 Msun) and distance (D < 12 Mpc) requirements (in both samples), mainly set by limits in STARBIRDS and LEGUS that our sample is drawn from.The diagonal gray lines show log sSFRs and the vertical line denotes our mass cutoff limit.
3 outlines the data preparation for our analysis.The luminosity function models and the results of our fitting are detailed in Sections 4 and 5, respectively.Finally, Section 6 discusses the L X -SFR-metallicity relation observed in our sample.

GALAXY SAMPLE
Here we present the "Chandra+HST" dwarf sample.We focus on local dwarf galaxies with stellar masses M * < 5 × 10 9 M ⊙ from two Local Volume surveys (D < 12.5 Mpc).The Chandra data are used to identify high-mass X-ray binaries, while the HST data ensure high-fidelity mass and distance measurements for the galaxies.Specifically, we start with a primary sample consisting of two large complementary sets of nearby galaxies, the Legacy ExtraGalactic UV Survey (LEGUS; Cignoni et al. 2018Cignoni et al. , 2019;;Sacchi et al. 2018) and the STARBurst IRregular Dwarf Survey (STARBIRDS; Mc-Quinn et al. 2015).Both LEGUS and STARBIRDS have sensitive multi-band coverage from the Hubble Space Telescope's (HST) WFC3, Galaxy Evolution Explorer Telescope (Martin et al. 2005, GALEX), and the Spitzer Space Telescope's IRAC offering full UV to IR coverage of the galaxy spectral energy distributions.Given the high quality of available data for these sources, they each have extremely accurate (1) direct distances derived from tip of the red giant branch (TRGB), Surface Brightness Fluctuations, Cepheids, and/or SNeII measurements; (2) M * measurements using multi-band photometry; (3) SFRs from HST+GALEX photometry; and (4) gas-phase metallicity measurements using oxygen abundances obtained from spectroscopic follow-up.We use X-ray data from NASA's Chandra X-ray Observatory, which offers the angular resolution and sensitivity needed to quantify the X-ray binary (XRB) populations in nearby galaxies.Among the 52 galaxies in the STARBIRDS and LEGUS datasets with stellar masses M * < 5 × 10 9 M ⊙ , 30 (58%) of these galaxies have publicly available data from Chandra.These observations were taken using the Advanced CCD Imaging Spectrometer (ACIS) from either I (Imaging) or S (Spectroscopy) cameras.Specifically, we label the sample of galaxies that have Chandra ACIS-I/S data as the Chandra+HST sample.
The stellar masses for the Chandra+HST sample were measured by Calzetti et al. (2015) and McQuinn et al. (2010) for the LEGUS and STARBIRDS sub-samples, respectively.Calzetti et al. (2015) obtain stellar masses from extinction-corrected B-band luminosities and color information, using the method described in Bothwell et al. (2009) and based on the mass-to-light ratio models from Bell & de Jong (2001). McQuinn et al. (2010) calculate the total amount of stellar mass from published absolute B-band luminosities of the galaxies, adjusted for extinction from Galactic dust maps published by Schlegel et al. (1998a).We use metallicities compiled by Calzetti et al. (2015) based on direct temperature measurements in the literature.For the Chan-dra+HST sample, we use distances reported by Lee et al. (2009).Galaxies in the Chandra+HST sample span a mass range of log (M * /M ⊙ ) = 6.8 − 9.5 and a distance range of 0.5 − 12.5 Mpc.
Although 30 galaxies in LEGUS+Starbirds have Chandra data, not all of them have sufficient data for our purposes.We remove two galaxies that were only observed with a sub-array, which is unable to cover the full extent of our nearby galaxies.We then remove an additional six galaxies because they do not reach our required depth of log L (erg s −1 ) > 37.5 (see Section 3.4 for details).
The rejected galaxies and the reasons that they were rejected are indicated by flags in the Flag column of Table 1.The total number of Chandra+HST dwarf galaxies in the sample after all cuts is 22. Lehmer et al. 2021 (L21) In addition to our primary sample, after applying our mass (M * < 5 × 10 9 M ⊙ ) and distance (D < 12.5 Mpc) cuts, we obtain 11 local dwarfs from the L21 main sample that are not in the Chandra+HST sample.We identified six additional dwarf galaxies (UGC 05340, NGC 1705, NGC 1569, NGC 5253, NGC 5474, andNGC 7793) that are in both L21 and Chandra+HST.We list the overlapping galaxies under the Chandra+HST category to avoid duplication.For these galaxies we use X-ray photometry and distances from L21 to preserve the X-ray luminosities reported in that work, and any other galaxy properties from measurements made for the Chandra+HST sample (see Section 3).

Additional Dwarfs from
The galaxies in the original L21 main sample span a mass range of log (M * /M ⊙ ) = 7.3 − 10.4 and a distance range of 1.9 − 29.4 Mpc.The metallicities of galaxies in the L21 sample were derived using oxygen abundance measurements either from strong-line calibrations or direct electron-temperature-based theoretical calibration (Lehmer et al. 2021).It is worth noting that for objects in common between the two samples, the metallicities are in good agreement despite being derived from different sources.

ULX Galaxies
Our goal is an unbiased dwarf sample to study X-ray point source distributions.If galaxies were preferentially targeted by Chandra because of known bright X-ray point sources (i.e., ULXs), this could bias the number of detected sources at the bright end of the luminosity function.Many galaxies in the sample were targeted because of the known presence of a ULX.Specifically, we find that four galaxies -NGC 7793, NGC 4490, NGC 4485, and NGC 1313 -were targeted for ULXs.Therefore, we fit the X-ray luminosity functions with and without the ULX-targeted sample to assess the potential bias introduced by their addition.

Summary of Final Sample
We have identified a total of 33 individual galaxies that meet our criteria, 22 from Chandra+HST and 11 from L21.We will refer to this set of galaxies as our final sample.However, we further perform our analysis with and without the four galaxies targeted because they harbor ULXs.The sample containing ULX galaxies will be referred to as the "final+ULX" sample.The galaxies in the final sample span a stellar mass range of log (M * /M ⊙ ) = 6.8 − 9.52, gas-phase metallicity range of 12 + log(O/H) = 7.74 − 8.77, and a distance range of 0.5 − 12.1 Mpc.If the galaxy was targeted by Chandra for ULXs, the galaxy is noted with a black box around its name.The apertures used in this study are plotted as ellipses (magenta).X-ray sources are overplotted as circles and colored according to their X-ray luminosities.For reference, vertical bars of size 1' (blue) and 1 kpc at the galaxy's distance (green) are provided in the lower left and lower right corners of each panel, respectively.

DATA ANALYSIS
In this section, we discuss galaxy apertures (Section 3.1), star formation rate (SFR) measurements (Section 3.2), X-ray photometry (Section 3.3), and recovery (completeness) functions (Section 3.4) used in this study.We compute galaxy sizes, recovery functions, and perform X-ray photometry for all galaxies in the Chan-dra+HST sample (i.e before final sample selection) while utilizing published X-ray catalogs from L21 for galaxies covered in that study.We compute galaxy sizes and associated FUV SFRs for galaxies in the final sample.Table 2 presents X-ray and SFR measurements for galaxies that met all of our selection criteria (i.e final sample + ULX).

Galaxy Apertures Using GALEX
We use galaxy projected footprints as apertures when making photometric measurements and as boundaries when filtering for X-ray point sources.L21 uses apertures defined by 2MASS (Jarrett et al. 2003) 20 mag arcsec −2 isophotal ellipses in the K s -band.However, several dwarf galaxies in this study are fainter in the K s -band than the L21 sample because they are relatively blue, causing the 2MASS K s -band apertures to severely underestimate the sizes of those galaxies, and therefore underestimate the number of HMXBs within them.Taking this into account, we use elliptical apertures from GALEX as a basis to define the projected apertures of galaxies in this study.
In their study, Lee et al. (2009) define the semi-major and semi-minor axes of the outermost elliptical annulus for each galaxy from the GALEX photometry.This boundary is determined based on two criteria: it is either the point where the annular flux error exceeds 0.8 magnitudes or the point where the intensity drops below the level of the sky background.We performed aperture photometry on the GALEX data using full-sized, half-sized, and quarter-sized versions of the apertures reported by Lee et al. (2009).We found that the half-sized apertures minimized cosmic X-ray background (CXB) contamination, closely matched the effective radius of most galaxies in our sample, and that they are similar to the apertures used in L21.We therefore adopt the half-sized apertures as our galaxy footprints.
NGC 5408 and ESO 495-G021 (He2-10) have no available UV-observations with GALEX.We use the 2MASS apertures provided in L21 since these galaxies are sufficiently bright in the K s -band (9.00 and 11.39 mag, respectively, Skrutskie et al. (2006)).The positions and Xray photometry for sources within and surrounding the L21 apertures are presented in Table A1 of Lehmer et al. (2021).To identify sources that fall within our galaxy aperture, we apply a positional filter on the point sources flagged in the L21 table as 1, 3, and 5 which correspond to point sources within the L21 aperture, sources outside the L21 aperture, and sources that are beyond 1.2 times the galaxy boundaries, respectively.We report our adopted apertures in Table 1 and display the apertures of Chandra+HST galaxies in Figure 2 (magenta ellipses).

Star Formation Rates
To ensure consistency and take into account the new apertures discussed in Section 3.1, we recalculate the FUV derived SFRs of each galaxy using GALEX data.
To achieve this, we measure the total flux enclosed within our FUV galaxy apertures and correct for Milky-Way foreground extinction using E(B-V) colors from Schlegel et al. (1998b) (as reported in Lee et al. (2009)).The FUV fluxes are then converted to luminosities using the distance to the galaxy, which in turn is used to measure the SFR.We expect internal FUV extinction in the dwarf galaxies to be small because our combined sample is predominantly comprised of blue dwarf galaxies with negligible dust.We use the scaling relation given by Lee et al. (2009) to convert FUV luminosity to SFR: Where L ν (U V ) is the total FUV luminosity density (erg s −1 Hz −1 ) enclosed.For the two galaxies with no GALEX data, we used the original SFR values from L21.

X-ray Photometry and Source Catalogs
The Chandra ACIS data were retrieved from the Chandra archive and homogeneously reprocessed using chandra repro as part of the CIAO package version 4.14 along with associated calibration files CALDB v4.9.6.As is standard, we applied the latest bad pixel masks and identified afterglow effects to create new masks, as well as applying the very faint processing flag to further clean the particle background for those observations performed in VFAINT mode to create the final reprocessed Level 2 events files and their associated response files.X-ray photometry was further carried out using CIAO.We started by constructing images in the 0.5 − 7.0 keV band from the reprocessed Level 2 events files for the 31 dwarf galaxy candidates in the Chandra+HST sample.For each of the galaxies, we used the CIAO fluximage script to create exposure-corrected images using the reprocessed ACIS graded Level-2 event files in good time intervals (GTI) along with the associated aspect solution files and bad pixel masks.The PSF radius was computed using the mkpsfmap script with the PSF energy set to 1.4967 keV and the encircled counts fraction (ECF) set to 0.9.Using the resulting flux image and PSF map as inputs, we ran wavdetect with wavelet scales of 1.0, 1.414, 2.0, 2.828, 4.0, 5.657, 8.0 pixels and sigthresh set to 10 −6 .This step produced our X-ray source photometry (counts) and catalog (positions).
For each galaxy, we filter out X-ray sources outside of the ACIS chip and the galaxy's aperture.In particular, for ACIS-S observations we filtered out sources that are not in the same chip as the target galaxy, while for front-illuminated ACIS-I observations we kept sources within the four chips (including the gaps).We visually inspected each galaxy to identify the presence of X-ray features that we deemed indicative of an AGN at the center of the flux images.UGC 5139, NGC 5253, NGC 4490, NGC 4605, and NGC 2500 were flagged as having possible AGN in the form of an X-ray point source coincident with their centers.As such, these suspected AGN point sources were filtered out.
We used the PIMMS (Portable, Interactive Multi-Mission Simulator) to calculate the ACIS counts-toflux conversion factors for each galaxy based on its observation cycle.The input energy range was set to Figure 3. Recovery functions within the apertures of the Chandra+HST sample of galaxies.Sub-array mode galaxies are not included in this figure.Galaxies are colored according to relative exposure times (i.e by ranking).The Recovery functions were computed by injecting simulated X-ray sources into Chandra images and measuring the fraction of the synthetic sources recovered during the source detection stage.The minimum LX in the simulations is set by the luminosity, at the galaxy's distance, that correspond to 4 counts.The gray vertical and horizontal lines denote a log luminosity of 37.5 ergs s-1 and 50% recovery respectively.If the ξ(L) profile falls below 50% before 37.5 ergs s-1, we do not include that galaxy in our analysis (Recovery functions of rejected galaxies are plotted as gray lines).We note that the recovery function of NGC 6822 extends to relatively lower luminosities despite its exposure time because of its proximity to the Milky Way.0.5 − 7.0 keV for all galaxies.Assuming a power-law spectral model we set Γ = 1.7 and estimated the galactic neutral hydrogen column density for each galaxy using HEASARC's software.We then calculated the pointsource fluxes by multiplying the counts from wavdetect by the conversion factor we obtained from PIMMS.Furthermore, we calculate the 0.5 − 7.0 keV X-ray pointsource luminosities "L" using the distances to the host galaxy.

Completeness Corrections
The recovery function, denoted by ξ(L), is a statistical measure of the fraction of sources that can be detected at a given X-ray luminosity while addressing source crowd-ing and the galactic backdrop's impact on point source detection.As such, we compute the recovery function of each galaxy in our sample by simulating a mock image and running our detection algorithm on the simulated dataset.We start with the original Chandra dataset as a base and use CIAO's simulate psf function to inject 100 additional synthetic sources into the galaxy aperture within the same point source luminosity bin.simulate psf accounts for streaks and PSF distortions.We run wavdetect on the simulated dataset and measure ξ(L) as the fraction of the injected sources that we are able to recover.We repeat this for each luminosity bin 10 times (a total of 1000 synthetic sources per luminosity bin per galaxy).We defined log L (erg s −1 ) = 39 as our maximum recovery function luminosity, while the lower bound is set to one or two bins below the 4 count limit.We sample in intervals of 0.5 dex for luminosities in the range log L (erg s −1 ) = 37 − 39 and 0.1 dex for log L (erg s −1 ) < 37. We use smaller energy bins for log L (erg s −1 ) < 37 to sufficiently sample the curvature of the recovery function as it begins to fall from unity to zero with decreasing luminosity.The recovery functions of all the galaxies in our sample are displayed in Figure 3.
To ensure the fidelity of our recovery functions, we take the placement of the synthetic sources and the size of each galaxy into consideration.All of the injected point sources are placed such that they are at least one half PSF radius away from each other, even if they are not injected into the same simulated image.This ensures that no two point sources are probing the same area and reduces the chance of source confusion caused by the synthetic sources overlapping with each other.Note that we do allow synthetic sources to overlap with real sources and incorporate such confusion into the recovery function.Since our sample contains dwarf galaxies, their sizes may be too small to allow for a meaningful measurement of ξ(L).To address this, we set the semimajor and semi-minor axes of the galaxy (i.e our sampling area) to a minimum of 3 arcmin.We also restrict the sampling area to be within the ACIS chip projected footprint (except for gaps in ACIS-I images).

MODELING
In this section, we discuss forward modeling of the Xray luminosity functions.Following L19 and L21, we model the observed X-ray luminosity distributions by taking into account the low-end completeness of the observations (Section 4.1), contributions from the cosmic X-ray background (CXB, Section 4.2), and high and low mass binaries (Section 4.3).By modeling each component of the luminosity function separately and combin-ing them into a compound model, we can represent and fit the observed data.The CXB for Chandra has been modeled by Kim et al. (2007) and the completeness of each observation is modeled using simulated recovery functions (see Section 3.4).We find the contributions of LMXBs, which scale with stellar mass, to be negligible for our sample of galaxies and use the L19 LMXB model to represent their magnitude.The last component is the HMXB component, which scales with SFR.We use the Astropy (Astropy Collaboration et al. 2022) models sub-package to implement individual and compound models.We describe each model component in detail below.

Recovery Function Model
We model the fraction of point sources recovered as a function of luminosity L X (i.e completeness) using the recovery functions from L21 and the simulated recovery functions in Section 3.4.We implement the model by using the tabulated recovery fractions as look-up tables, interpolating between data points when necessary.For luminosities that correspond to 4 counts or less, we enforce a recovery fraction of 0. For luminosities above the range of the simulated recovery function, we use a recovery fraction of 1.

Cosmic X-ray Background
The CXB is the combined X-ray flux from all distant bright X-ray sources, such as AGNs, that may be confused with relatively fewer luminous sources within a nearby galaxy.To model the contribution of such sources, we use the extragalactic X-ray point source number counts from Kim et al. (2007).They provide broken power law models, given in Equation 2, for the ChaMP+CDFs (Kim et al. 2004;Rosati et al. 2002) number counts per unit area, with parameters listed in their publication's Table 4 for each galaxy.We scale this model by the area of our apertures and convert the CXB flux values to luminosities using the distance to the galaxy.To cover the same energy range as our observations, we use the 0.5 − −8 keV ChaMP+CDFs models.Given a power-law photon index of 1.4, the difference between the 0.5 − 7 keV and 0.5 − 8 keV energy ranges is negligible (10%).Our adopted model for the flux (S) dependant CXB number counts can be expressed as: where γ 1 and γ 2 are the broken power-law slopes, K is a normalization constant and S ref is a normalization flux set to 10 −15 erg s −1 cm −2 , and S b is the break flux at which the slope changes.The differential number count (dN CXB /dS) as stated in Equation 2 is in units of 10 −15 deg −2 .

LMXB and HMXB Models
We model the LMXB contributions as a function of stellar mass utilizing the L19 broken power-law model (Equation 12in L19) with a cutoff luminosity.We use parameter values obtained by the L19 full sample (Table 4 of L19, column 4).We do not fit any parameters for the LMXB at any point because the contributions from LMXBs are negligible.
For the HMXB, there are two models available.The first model is an SFR-normalized single power-law model from L19, with a cutoff luminosity.The L19 HMXB differential number counts as a function of X-ray luminosity has the following form: (3) where SFR is the host galaxy star formation rate, K HMXB is a normalizing constant, and γ is the powerlaw slope.Following L19, we take L to be in units of 10 38 erg s −1 (denoted by L 38 = L/10 38 ) for the HMXB and LMXB models.The only component that varies across galaxies in this model is the SFR, with B(L) being universal for all galaxies in a given metallicity bin.
The L21 HMXB model introduces a metallicity dependence to the differential number counts (Equations 1, 2, and 3 in L21).The function has the form of a broken power-law with an exponential cutoff.The highluminosity power-law exponent, as well as the cutoff luminosity, have metallicity dependencies.When we include the L21 model for comparison, we do not refit any parameters but use the values reported in Table 2 of L21.

Compound Model
Following L19 and L21, we combine all the elements described above to forward model the underlying HMXB luminosity function (Equation 4) implied by the ensemble measured luminosity distributions, accounting for incompleteness, CXB, and the LMXB population.The underlying HMXB luminosity function is weighted by the SFR of each galaxy contributing to the ensemble.To account for completeness, the contribution of each galaxy is weighted by its recovery function at each luminosity.This provides us with a model that represents the observed total number counts in each L X bin (∆n): (5) Since we are interested in modeling the HMXB contributions, we focus on the SFR independent component B(L) (see Equation 3).We rearrange Equation 5 to separate B(L) from the rest of ∆n as follows: Where: Note again that B(L), which is scaled by the SFR of each individual galaxy, is an underlying model across all galaxies.C 1 (L) and C 2 (L) are components that scale with SFR, projected area, and stellar mass.Since these components do not contain fitting parameters, we can tabulate the values of C 1 (L) and C 2 (L) as a function of luminosity bins and fit for B(L).

LUMINOSITY FUNCTION FITS
In this Section we provide information on how we fit the compound model, discussed in Section 4.4, to the observed luminosity functions of sub-samples of the final sample (see Section 2.3).We summarize our results and compare them to past models in Section 5.5.

Metallicity Dependence of Luminosity Function
The sub-samples of the final sample are differentiated based on two factors: (1) metallicity bins with 12 + log(O/H) values of 7.7 − 8.3 (low), 8.3 − 8.9 (high) and 7.7−8.9(full), and (2) whether the sub-sample contained galaxies that were targeted for ULXs. Figure 4 shows the observed X-ray luminosity functions for the resulting five unique sub-samples1 .

Observed X-ray Luminosity Functions
The observed X-ray luminosity functions are constructed by aggregating the luminosities of all X-ray point sources within all galaxy apertures in a given metallicity bin.X-ray point source luminosities are binned into intervals similar to L21.The luminosity bins range from log L (erg s −1 ) = 35−41.7 and each bin spans 0.1 dex, resulting in 78 bins.When inferring models, the log mid-points of the bins are used.The final sample luminosity functions for each metallicity bin are plotted as black points with error bars in Figure 4.The one-sigma Poisson errors for the number of sources in each luminosity bin are calculated according to Gehrels (1986).

Model Fitting
As discussed above, the observed luminosity function for each sub-sample is constructed by consolidating all X-ray point sources in all galaxies within a metallicity bin.For each metallicity bin, we limit the range of L that we fit to the highest and lowest L bins that contain at least a single source.We initialize the compound model discussed in Section 4.4 using the full sample parameter values in L19 and fix the value of L c to 10 40.7 erg s −1 .The values of C 1 (L) and C 2 (L) are tabulated at the mid point of each luminosity bin.We use a Levenberg-Marquardt algorithm (Astropy's LevMarLSQFitter) to fit K HMXB and γ, which are parameters of B(L), in the compound model (see Equations 3, 4, and 5).The quality of our models' representation of the data was evaluated using the C-statistic, as described in Section 5.4.
We account for the uncertainties in SFR by implementing a Monte Carlo sampling technique, which consists of 10 4 iterations.In each iteration, the luminosity function is refit with the SFR for each galaxy drawn from a normal distribution of SFR with a mean equal to the measured SFR and a standard deviation equal to the corresponding one-sigma error.If the galaxy does not have FUV SFR measurements, we use half the SFR value as the standard deviation to reflect the relatively higher degree of uncertainty.We limit the range of allowed SFRs to 5 standard deviations from the mean and enforce a minimum SFR of 0.0001.The mean of the parameters resulting from all the iterations is taken as our final result.It is worth mentioning that the fluctuations in the SFR do not result in substantial propagated uncertainties in the fitted parameters.The parameter errors of each fit (at each iteration) are estimated by taking the diagonal of the covariance matrix from the Levenberg-Marquardt optimizer for the best SFR.
Consistent with L19, we adopt a cutoff luminosity (L c ) from their comprehensive sample fit.However, it's worth noting that this cutoff luminosity, while empirically motivated, may have significant uncertainty.This choice of fixed L c could bias our fit.Thus, to validate its impact, we conduct a test to estimate the number of sources that our best fit predicts above this cutoff.We integrate the fitted luminosity functions above L c and find that two or fewer sources are expected.This negligible result alleviates concerns about our fit converging on an unphysically flat slope due to artificial truncation.

Goodness of Fit
Following L19, we evaluated the goodness of fit for each of the metallicity bins on a global basis using a modified C-statistic (Cash 1979;Kaastra 2017): where C denotes the C-statistic corresponding to a particular metallicity bin.The sum of the statistic is calculated by iterating through the luminosity bins.n b is the total number of luminosity bins, whereas M i and N i represent the model value and the observed counts for the i th bin, respectively.We mask out bins where the model is equal to zero (C i = 0 if M i = 0).For bins in which the observed number of sources is zero, we set the logarithmic component to zero for that term We adopted the methods outlined by Kaastra (2017) to compute the expected C-statistic (C exp ), along with its variance (C var ), based on Poisson statistics.Subsequently, we evaluated the null-hypothesis probability (P null ) as follows: Where P null is the null-hypothesis probability, C is the C-statistic measured from the data following Equation 9, and erf is the error function.Models with P null < 0.001 can be statistically rejected with greater than 99.9% confidence.Utilizing a modified C-statistic, we assess the goodness of fit across various metallicity bins and this threshold ensures the reliability of the models we test.

Fitting Results and Comparisons
The results of our fits are summarized in Table 3, and the corresponding figures are presented in Figures 4 and  5.The differential number counts of X-ray point sources as a function of luminosity are provided in Figure 4, grouped as described previously, and the best-fit parameters for the luminosity function models are provided in 5.
Our analysis of the full metallicity bin reveals a powerlaw slope of γ = 1.40 ± 0.03.When fitting the same metallicity bin, including galaxies targeted for ULXs, we find a slightly shallower power-law slope of γ = 1.36 ± 0.09.However, this difference is not statistically significant, indicating that the presence of ULX targeted galaxies does not have a meaningful effect on the overall slope of the dwarf HMXB population.
We examine the metallicity dependence of the HMXB luminosity function by fitting and comparing subsamples of galaxies binned by metallicity ranges (Section 5.1).We find that the low metallicity bin is best fit by a shallower slope than the full sample, while the high metallicity bin results in a steeper slope.The shallower slope of the low metallicity bin implies a relative excess of high luminosity sources in that bin.Interestingly, only the higher metallicity bin includes ULX-targeted galaxies.Similar to the full metallicity sample, ULXtargeted galaxies in the high metallicity sample have slightly shallower luminosity functions slopes, but the difference is not statistically significant.An unbiased sample is likely to provide results that lie somewhere between the ULX-included and ULX-excluded sample results.
In Figure 5, we compare the best-fit parameters to the L19 results, which utilizes a broader range of stellar masses, inclusive of galaxies larger than dwarfs.We find that dwarf galaxies exhibit a shallower slope of γ ∼ 1.40 compared with the full L19 sample, 1.65 ± 0.03, equating to an ∼ 8.0σ deviation, suggestive of a mass and/or metallicity dependence on the power-law slope.By contrast, we find no evidence for a dependence on the normalization parameter K HMXB .We also note that the systematic offsets in the normalization parameter due differences in SFR measurements are are found to be minimal.Specifically, the average ratio of our measurements to those in L21 is 1.07, with a standard deviation of 0.5, with a maximum deviation of ∼ 2.4.We test how well these models describe the HMXB population by taking into consideration the null-hypothesis probability derived from their C-stats.We find that all models result in acceptable fits, with the exception of the L19 model in the case of the full metallicity bin with ULX galaxies.

Observed L X /SFR
In addition to looking at the luminosity function for Xray binaries across all the galaxies in the sample, we can also explore how the X-ray luminosity from individual sources varies with the gas-phase metallicity.Since we know that SFR is the primary driver of the L X from high-mass X-ray binaries, we must normalize each L X by the SFR in that galaxy in order to isolate any additional correlation with metallicity.
We consider L X /SFR as a function of metallicity for individual galaxies in Figure 6.For each galaxy in the Chandra+HST sample, the total HMXB X-ray luminosity is computed by summing the observed point source luminosities and subtracting the expected LMXB and CXB contributions.The galaxies with L X below the expected LMXB and CXB contributions have upper-limits that are 3σ above the expected background from these two sources in that galaxy.
With our expanded sample of dwarf galaxies, the vast majority of low-mass systems fall well below the predicted L X /SFR relations.These low L X /SFR values cannot be explained by incompleteness or aperture effects.Many of our deepest exposures have no HMXBs, and even when we double our SFR apertures, we recover the same result.As we will discuss in the next section, this wide range in L X /SFR is actually expected, because at the low SFR of the dwarfs in the Chandra+HST sample (0.01 − 0.1 M ⊙ /yr), galaxies cannot fully populate the HMXB X-ray luminosity function, and instead only a very small number of galaxies are expected to harbor a luminous X-ray binary at any given time due to stochastic sampling (Gilfanov 2004;Lehmer et al. 2021).

THE L X -SFR-METALLICITY RELATION
We have introduced an expanded sample of dwarf galaxies with measured X-ray binary luminosity functions, covering a broader range of sSFR and metallicity than prior work.This additional sample shows evidence for an excess of luminous X-ray sources at low metallicity in the stacked luminosity functions, and also highlights the challenges of systematic studies in dwarf galaxies where the typical star formation rates are low.In order to put our results into context, we first review the theoretical reasons to expect a dependence on metallicity.We then contextualize our measured distribution of L X /SFR with prior work.

Model Predictions
Population synthesis models (Dray 2006;Fragos et al. 2013b,a;Madau & Fragos 2017b) suggest that the observed scatter in the L X -SFR relation can be explained by a secondary dependence of L X on metallicity.
The metallicity dependence is thought to arise because stars with higher metallicity are known to undergo greater angular momentum loss due to stronger winds, which results in orbital expansion that reduces the chance of forming Roche lobe overflow (Fragos et al. 2013b).Moreover, stronger radiatively driven winds in high-metallicity stars cause them to undergo enhanced mass loss before their eventual supernova explosions (Fornasini et al. 2020).This process tends to yield compact objects that are relatively lower in mass and less numerous in high-metallicity galaxies, resulting in diminished X-ray luminosities from HMXBs.This anticorrelation causes us to expect lower-metallicity galaxies to have sources that extend to higher L X .This trend has been seen in observations (Douna et al. 2015;Brorby et al. 2016, L19, L21), however, the same lowmetallicity galaxy sample is in common across nearly all of these studies.These low-metallicity galaxies also have high specific star-formation rates, which correlates with low-metallicity (Mannucci et al. 2010).Thus, our Chandra+HST sample both fills in intermediate metallicities and alleviates the strong bias to the highest starformation rates.

SFR Driven Stochasticity
In agreement with L21, the shallower slopes observed in the stacked luminosity functions (as discussed in Section 5) imply that flatter luminosity functions are expected (i.e., more sources at higher luminosities) with decreasing metallicity.At the same time, we observe a considerable spread in L X /SFR values for individual sources.
We can reconcile the shallower slopes in HMXB luminosity function stacks with the large spread in L X /SFR by considering the low SFRs of the typical dwarfs in our sample.Specifically, the disparity can be largely attributed to stochastic Poissonian sampling effects that  Parameter fits and associated one sigma errors for five sub-samples compared to L19 parameters.The square marker denotes the L19 fitted parameters and onesigma errors for their full sample.The gray dashed lines mark the L19 parameter values to help with comparison.Our sub-samples without ULX galaxies are marked with filled circles while sub-samples with ULX galaxies are marked with empty circles.The metallicity bins of the sub-samples are labeled by name, where All=7.7 − 8.9, Lower=7.7 − 8.3, Upper=8.3 − 8.9.
come into play due to the inherently low SFRs in dwarf galaxies.This is due to the fact that (Grimm et al. 2003;Gilfanov et al. 2004, L21): Where L up and L lo are the upper and lower bounds of the integrated L X .If the power-law γ is less than 2, then the first term dominates, resulting in a highly stochastic L X /SFR distribution through a dependence on the highest X-ray source in the galaxy.In contrast, a γ value greater than 2 results in a stable value.Grimm et al. (2003) and L21 quantitatively explore the impact of stochasticity at low SFR.In Fig. 5 of L21, they present a Monte Carlo simulation of the expected distribution in L X /SFR given their fitted HMXB luminosity function.These simulations show an inherent stochastic scatter in L X /SFR at low SFR and that the distributions become Gaussian at SFR ≥ 2 − 5 M ⊙ /yr, which are much higher than the typical SFR of dwarf galaxies (see Fig. 1).LX/SFR (HMXB) vs metallicity relation for galaxies in the Chandra+HST and all L21 samples.Dwarfs from the Chandra+HST sample are plotted in blue and L21 galaxies are plotted in black.Galaxies from the L21 supplemental sample, which are compact dwarfs, are plotted in gray.The total HMXB Lx values are calculated by summing the observed Lx and subtracting out the LMXB and CXB contributions after correcting for completeness.Upper limits are used for galaxies with few or no sources using the expected LMXB and CXB contributions.
To quantify the spread in L X /SFR for the Chan-dra+HST sample, we conducted Monte Carlo simulations utilizing the L21 HMXB luminosity function, treated as a metallicity-dependent probability density function.For each SFR and 12 + log(O/H) pairing, we initially integrated the HMXB number densities (for log L X > 35) to estimate the expected source count.Subsequently, we performed random sampling of the predicted number of HMXBs from the L21 model.This Monte Carlo approach was iterated 10, 000 times, and the mean L X /SFR for each SFR and 12 + log(O/H) pair was computed to characterize the distribution.We show the resulting probability distributions in Figure 7.The six panels correspond to SFRs of 0.005, 0.015, 0.05, 0.5, 1.0, and 5.0 M ⊙ /yr respectively.We over-plot the galaxies in Figure 6 binned by SFR for reference.
In Figure 7, for the lowest SFR bin (SFR = 0.005), we observe a reduced probability of detecting a galaxy with a high L X /SFR, a consequence of stochastic Poissonian sampling, despite the over-representation of high L X sources in the HMXB luminosity function at low metallicity 12 + log(O/H).This is because in the case of SFR = 0.005, we expect small source number counts (N = 2 sources) which results in near-zero Figure 7.We show the probabilistic distribution of LX/SFR against 12 + log(O/H) for SFRs and gas-phase metallicities similar to the dwarf samples discussed in this work.The Monte Carlo probabilities were computed, for each SFR and 12 + log(O/H) pair, by using the L21 HMXB model as the probability density function.The six panels correspond to SFRs of 0.005, 0.015, 0.05, 0.5, 1.0, and 5.0 M⊙/yr respectively.The greyscale shows the probability of finding a galaxy with a LX/SFR for a given 12 + log(O/H).The red markers show the observed galaxy LX/SFR values and upper-limits, including all galaxies in this work and L21 for reference (not just dwarfs).The red line represents the L21 model prediction for reference.At lower SFRs, the probability of encountering a high LX/SFR galaxy is reduced due to stochastic Poissonian sampling, despite a surplus of high LX sources in the HMXB luminosity function of galaxies with low 12 + log(O/H).With increasing SFR, the probability distribution approaches the L21 prediction (red line).Some observed points may hover slightly above or below the probability evaluated at the SFR displayed, this is due to the scatter introduced by the range of galaxy SFRs.
high-luminosity sources populating and dominating the total L X of the galaxy, thereby inducing substantial stochastic variability in L X /SFR.As the SFR increases (SFR = 0.015 and SFR = 0.05) the probability of encountering a high L X /SFR galaxy is enhanced because of the slightly larger number of sources being sampled from the luminosity functions.If we take a slice at a low 12 + log(O/H) value (vertically) for these SFRs, the L X /SFR probability would be a biomodal distribution due to a single source dominating the L X of the galaxy.At high SFRs (SFR > 1.0), the larger source numbercounts result in the distributions becoming Gaussian, which in turn results in the probabilities approaching the L21 predictions.Thus, at high SFR, the L X /SFR approaches the theoretical value with small scatter.The Prestwich et al. (2013) dwarf sample (L21 supplementary) exhibit higher L X /SFR because they have much higher specific SFR in general.In particular, six galaxies (30% of the supplementary sample) have log L X /SFR above 40 and SFRs that are greater than 0.01.
Returning to one of our main motivations for this work, to understand how and when we may use L X /SFR to search for accreting massive black holes, we have a hopeful message.While it will be necessary to fold in the detailed SFR distributions and (ideally) metallicity distributions for a sample to derive the expected number of detections from HMXBs alone, these distributions mean that the contribution from HMXBs should be low (∼ 3−4%) from dwarfs that lie on the star-forming main sequence and the local mass-metallicity relation.

SUMMARY
In this paper we present the HMXB X-ray luminosity functions of dwarf galaxies within the local volume, and

Figure 1 .
Figure1.Stellar mass versus SFR for the galaxies in the L21 galaxy sample ("Lehmer et al. 21") and galaxies introduced in this study ("Chandra+HST").Chandra+HST are reprinted by stars while L21 galaxies are marked by square points.The points are colored according to gas-phase metallicity (12 + log(O/H)).Our dwarf galaxy sample supplements the L21 sample in M * , SFR, and metallicity while maintaining a comparable sSFR.The gray circles denote galaxies that satisfy the mass (M * < 5 × 10 9 Msun) and distance (D < 12 Mpc) requirements (in both samples), mainly set by limits in STARBIRDS and LEGUS that our sample is drawn from.The diagonal gray lines show log sSFRs and the vertical line denotes our mass cutoff limit.

Figure 2 .
Figure 2.This figure shows the dwarf galaxies in the Chandra+HST final sample.In each panel, we include a DSS image of the galaxy, which is displayed in gray-scale.Galaxy names are displayed at the top of each panel and the corresponding [12 + log(O/H), log(SFR)] values are displayed at the bottom.If the galaxy was targeted by Chandra for ULXs, the galaxy is noted with a black box around its name.The apertures used in this study are plotted as ellipses (magenta).X-ray sources are overplotted as circles and colored according to their X-ray luminosities.For reference, vertical bars of size 1' (blue) and 1 kpc at the galaxy's distance (green) are provided in the lower left and lower right corners of each panel, respectively.

Figure 4 .
Figure4.This figure shows models that were fit to data in three gas-phase metallicity bins.All panels show the mass and distance filtered final sample of galaxies.The first three panels show samples without galaxies that were targeted for ULXs while the last two include them.The metallicity range and the number of galaxies in each panel is given the upper left corner.The first and fourth panels represent the entire dwarf sample within the metallicity range of this study (i.e 12 + log(O/H) = 7.7 − 8.9).Observed distributions of X-ray point source luminosities (∆n = dN/dL) for galaxy sub-samples are plotted with one sigma Poisson error bars (black).Recovery functions have been applied to all models to match the expected completeness of the data.The red and green lines model LMXB and CXB contributions respectively.The gray dashed, solid lines, and blue lines show the outputs of the combined (HMXB + LMXB + CXB) models for L19, L21, our study.The numeric values of the fitted amplitude (K hmxb ) and slope (γ) are given at the top right corner of each panel.
Figure 5.Parameter fits and associated one sigma errors for five sub-samples compared to L19 parameters.The square marker denotes the L19 fitted parameters and onesigma errors for their full sample.The gray dashed lines mark the L19 parameter values to help with comparison.Our sub-samples without ULX galaxies are marked with filled circles while sub-samples with ULX galaxies are marked with empty circles.The metallicity bins of the sub-samples are labeled by name, where All=7.7 − 8.9, Lower=7.7 − 8.3, Upper=8.3 − 8.9.
Figure 6.LX/SFR (HMXB) vs metallicity relation for galaxies in the Chandra+HST and all L21 samples.Dwarfs from the Chandra+HST sample are plotted in blue and L21 galaxies are plotted in black.Galaxies from the L21 supplemental sample, which are compact dwarfs, are plotted in gray.The total HMXB Lx values are calculated by summing the observed Lx and subtracting out the LMXB and CXB contributions after correcting for completeness.Upper limits are used for galaxies with few or no sources using the expected LMXB and CXB contributions.
mode, and RF is recovery function cutoff.(13) Source catalog where LG is LEGUS, SB is STARBIRDS, and L21.

Table 3 .
Fitting results and goodness of fit tests.(1) Gas-phase metallicity.(2-3) The values of the model parameters obtained after fitting.(4) Cash statistics value for fitted parameters.(5-6) Expected value and variance of the Cash statistics for the fitted model.(7) Null hypothesis probability for the fitted model.(8-9) Null hypothesis probability for L19 and L21 models.

Table 4 .
An abbreviated version of the source catalog for the Chandra+HST sample.The full catalog contains 184 sources, and a discritpion is provided in appendix.