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

In addition to our primary sample, after applying our mass (M_* < 5 × 10⁹ M_⊙) and distance (D < 12.5 Mpc) cuts, we obtained 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, and NGC 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).

The galaxies in the original L21 main sample span a mass range of $\mathrm{log}\ ({M}_{* }/{M}_{\odot })=7.3\mbox{--}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 (L21). 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.

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.

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 $\mathrm{log}\,({M}_{* }/{M}_{\odot })\,=\,6.8\mbox{--}9.52$, gas-phase metallicity range of $12+\mathrm{log}({\rm{O}}/{\rm{H}})=7.74\mbox{--}8.77$, and a distance range of 0.5–12.1 Mpc.

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 the Two Micron All Sky Survey (2MASS; Jarrett et al. 2003) $20\ \mathrm{mag}\ {\mathrm{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 semimajor and semiminor 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 mag 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.

Zoom In Zoom Out Reset image size

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 X-ray photometry for sources within and surrounding the L21 apertures are presented in Table A1 in L21. To identify sources that fall within our galaxy aperture, we apply a positional filter on the point sources flagged in the table in L21 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).

Zoom In Zoom Out Reset image size

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. (1998; 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 the FUV luminosity to an SFR of

$$\begin{eqnarray}&&\mathrm{SFR}=1.4\cdot {10}^{-28}\cdot {L}_{\nu }({\rm{UV}})\end{eqnarray} \tag{ 1 }$$

where L_ν (UV) is the total FUV luminosity density (erg s⁻¹ Hz⁻¹) enclosed. For the two galaxies with no GALEX data, we used the original SFR values from L21.

The Chandra ACIS data were retrieved from the Chandra archive and homogeneously reprocessed using chandra_repro as part of the CIAO package v4.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 and applied 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 along with the associated aspect solution files and bad pixel masks. The point-spread function (PSF) radius was computed using the mkpsfmap script with the PSF energy set to 1.4967 keV and the encircled counts fraction 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, and 8.0 pixels and sigthresh set to 10⁻⁶. 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-to-flux conversion factors for each galaxy based on its observation cycle. The input energy range was set to 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 the HEASARC software. We then calculated the point-source 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 point-source luminosities "L" using the distances to the host galaxy.

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 crowding 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 data set. We start with the original Chandra data set 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 data set 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 $\mathrm{log}L({\rm{erg}}\,{{\rm{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 of $\mathrm{log}L({\rm{erg}}\,{{\rm{s}}}^{-1})=37\mbox{--}39$ and 0.1 dex for $\mathrm{log}L({\rm{erg}}\,{{\rm{s}}}^{-1})\lt 37$. We use smaller energy bins for $\mathrm{log}L({\rm{erg}}\,{{\rm{s}}}^{-1})\lt 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.

Zoom In Zoom Out Reset image size

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 semiminor axes of the galaxy (i.e., our sampling area) to a minimum of 3'. We also restrict the sampling area to be within the ACIS chip projected footprint (except for gaps in ACIS-I images).

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.

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 (Rosati et al. 2002; Kim et al. 2004) 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) dependent CXB number counts can be expressed as

$$\begin{eqnarray}\displaystyle \frac{{{dN}}_{\mathrm{CXB}}}{{dS}}=\left\{\begin{array}{ll}K{\left(\displaystyle \frac{S}{{S}_{\mathrm{ref}}}\right)}^{-{\gamma }_{1}} & S\lt {S}_{b}\\ K{\left(\displaystyle \frac{{S}_{b}}{{S}_{\mathrm{ref}}}\right)}^{({\gamma }_{2}-{\gamma }_{1})}{\left(\displaystyle \frac{S}{{S}_{\mathrm{ref}}}\right)}^{-{\gamma }_{2}} & S\geqslant {S}_{b}\end{array}\right.\end{eqnarray} \tag{ 2 }$$

where γ₁ and γ₂ are the broken power-law slopes, K is a normalization constant S_ref is a normalization flux set to 10⁻¹⁵ erg s⁻¹cm⁻², and S_b is the break flux at which the slope changes. The differential number count (${{dN}}_{\mathrm{CXB}}/{dS}$) as stated in Equation (2) is in units of 10⁻¹⁵ deg⁻².

We model the LMXB contributions as a function of stellar mass utilizing the L19 broken power-law model (Equation (12) in L19) with a cutoff luminosity. We use parameter values obtained by the L19 full sample (Table 4, column (4) of L19). We do not fit any parameters for the LMXB at any point because the contributions from LMXBs are negligible.

There are two models available for the HMXB. 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:

$$\begin{eqnarray}B(L)={K}_{\mathrm{HMXB},38}\begin{array}{cc}{L}_{38}^{-\gamma } & L\lt {L}_{c}\\ 0 & L\geqslant {L}_{c}\end{array}\end{eqnarray} \tag{ 3 }$$

$$\begin{eqnarray}&&\displaystyle \frac{{{dN}}_{\mathrm{HMXB}}}{{dL}}=\mathrm{SFR}\cdot B(L)\end{eqnarray} \tag{ 4 }$$

where SFR is the host galaxy star formation rate, K_HMXB is a normalizing constant, and γ is the power-law slope. Following L19, we take L to be in units of 10³⁸ erg s⁻¹ (denoted by L₃₈ = L/10³⁸) 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)–(3) in L21). The function has the form of a broken power law with an exponential cutoff. The high-luminosity 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.

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):

$$\begin{eqnarray}&&{\rm{\Delta }}n=\displaystyle \sum _{i}^{{N}_{\mathrm{gal}}}{\xi }_{i}(L){\rm{\Delta }}L\left[\displaystyle \frac{{{dN}}_{\mathrm{HMXB},i}}{{dL}}+\displaystyle \frac{{{dN}}_{\mathrm{LMXB},i}}{{dL}}+\displaystyle \frac{{{dN}}_{\mathrm{CXB},i}}{{dL}}\right]\end{eqnarray} \tag{ 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:

$$\begin{eqnarray}&&{\rm{\Delta }}n=B(L)\cdot {C}_{1}(L)+{C}_{2}(L)\end{eqnarray} \tag{ 6 }$$

where

$$\begin{eqnarray}&&{C}_{1}(L)={\rm{\Delta }}L\displaystyle \sum _{i}^{{N}_{\mathrm{gal}}}{\xi }_{i}(L)\cdot {\mathrm{SFR}}_{i}\end{eqnarray} \tag{ 7 }$$

$$\begin{eqnarray}&&{C}_{2}(L)={\rm{\Delta }}L\displaystyle \sum _{i}^{{N}_{\mathrm{gal}}}{\xi }_{i}(L)\left[\displaystyle \frac{{{dN}}_{\mathrm{LMXB},i}}{{dL}}+{\mathrm{CXB}}_{i}(L)\right]\end{eqnarray} \tag{ 8 }$$

Note again that B(L), which is scaled by the SFR of each individual galaxy, is an underlying model across all galaxies. C₁(L) and C₂(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₁(L) and C₂(L) as a function of luminosity bins and fit for B(L).

The subsamples of the final sample are differentiated based on two factors: (1) metallicity bins with $12+\mathrm{log}({\rm{O}}/{\rm{H}})$ values of 7.7–8.3 (low), 8.3–8.9 (high) and 7.7–8.9 (full), and (2) whether the subsample contained galaxies that were targeted for ULXs. Figure 4 shows the observed X-ray luminosity functions for the resulting five unique subsamples. ⁴

Zoom In Zoom Out Reset image size

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 those in L21. The luminosity bins range from $\mathrm{log}L({\rm{erg}}\,{{\rm{s}}}^{-1})=35\mbox{--}41.7$ and each bin spans 0.1 dex, resulting in 78 bins. When inferring models, the log midpoints 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 1σ Poisson errors for the number of sources in each luminosity bin are calculated according to Gehrels (1986).

As discussed above, the observed luminosity function for each subsample 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⁻¹. The values of C₁(L) and C₂(L) are tabulated at the midpoint 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)–(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⁴ 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 1σ 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 is 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.

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):

$$\begin{eqnarray}&&C=2\displaystyle \sum _{i=1}^{{n}_{b}}{M}_{i}-{N}_{i}+{N}_{i}\mathrm{ln}\left({N}_{i}/{M}_{i}\right)\end{eqnarray} \tag{ 9 }$$

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 ith 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 (${N}_{i}\mathrm{ln}\left({N}_{i}/{M}_{i}\right)=0$ if N_i = 0).

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:

$$\begin{eqnarray}&&{P}_{\mathrm{null}}=1-\mathrm{erf}\left(\sqrt{\displaystyle \frac{{\left(C-{C}_{\exp }\right)}^{2}}{{C}_{\mathrm{var}}}}\right)\end{eqnarray} \tag{ 10 }$$

where P_null is the null hypothesis probability, C is the C-statistic measured from the data following Equation (9), and $\mathrm{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.

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 Figure 5.

Zoom In Zoom Out Reset image size

Our analysis of the full metallicity bin reveals a power-law 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, ULX-targeted galaxies in the high metallicity sample have slightly shallower luminosity function 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 utilize 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 of a dependence on the normalization parameter K_HMXB. We also note that the systematic offsets in the normalization parameter due to differences in SFR measurements 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.

In addition to looking at the luminosity function for XRBs 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 HMXBs, 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.

Zoom In Zoom Out Reset image size

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 XRB at any given time due to stochastic sampling (Gilfanov 2004; Lehmer et al. 2021).

Population synthesis models (Dray 2006; Fragos et al. 2013b, 2013a; Madau & Fragos 2017) 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. 2013a). 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 low-metallicity galaxy sample is in common across nearly all of these studies. These low-metallicity galaxies also have high specific SFRs, 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 SFRs.

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 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)

$$\begin{eqnarray}\begin{array}{rcl}{L}_{{\rm{X}}} & = & {\displaystyle \int }_{{L}_{\mathrm{lo}}}^{{L}_{\mathrm{up}}}L\ B(L)\ {dL}={K}_{\mathrm{HMXB}}{\displaystyle \int }_{{L}_{\mathrm{lo}}}^{{L}_{\mathrm{up}}}{L}^{-\gamma +1}{dL}\\ & = & \displaystyle \frac{{K}_{\mathrm{HMXB}}}{2-\gamma }\left[{L}_{\mathrm{up}}^{2-\gamma }-{L}_{\mathrm{lo}}^{2-\gamma }\right]\end{array}\end{eqnarray} \tag{ 11 }$$

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 Figure 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 Figure 1).

To quantify the spread in L_X/SFR for the Chandra+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+\mathrm{log}({\rm{O}}/{\rm{H}})$ pairing, we initially integrated the HMXB number densities (for $\mathrm{log}{L}_{{\rm{X}}}\gt 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+\mathrm{log}({\rm{O}}/{\rm{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 overplot the galaxies in Figure 6 binned by SFR for reference.

Zoom In Zoom Out Reset image size

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+\mathrm{log}({\rm{O}}/{\rm{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 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 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+\mathrm{log}({\rm{O}}/{\rm{H}})$ value (vertically) for these SFRs, the L_X/SFR probability would be a bimodal distribution due to a single source dominating the L_X of the galaxy. At high SFRs (SFR > 1.0), the larger source number counts 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 material) exhibits 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.