HMXB luminosity function

# Bright end of the luminosity function of high-mass X-ray binaries: contributions of hard, soft and supersoft sources

S. Sazonov and I. Khabibullin
Space Research Institute, Russian Academy of Sciences, Profsoyuznaya 84/32, 117997 Moscow, Russia
Moscow Institute of Physics and Technology, Institutsky per. 9, 141700 Dolgoprudny, Russia
Max Planck Institute for Astrophysics, Karl-Schwarzschild-Strasse 1, D-85741 Garching, Germany
###### Abstract

Using a spectral analysis of bright Chandra X-ray sources located in 27 nearby galaxies and maps of star-formation rate (SFR) and ISM surface densities for these galaxies, we constructed the intrinsic X-ray luminosity function (XLF) of luminous high-mass X-ray binaries (HMXBs), taking into account absorption effects and the diversity of HMXB spectra. The XLF per unit SFR can be described by a power law  yr from to  erg s, where is the unabsorbed luminosity at 0.25–8 keV. The intrinsic number of luminous HMXBs per unit SFR is a factor of larger than the observed number reported before. The intrinsic XLF is composed of hard, soft and supersoft sources (defined here as those with the 0.25–2 keV to 0.25–8 keV flux ratio of , 0.6–0.95 and , respectively) in 2:1:1 proportion. We also constructed the intrinsic HMXB XLF in the soft X-ray band (0.25–2 keV). Here, the numbers of hard, soft and supersoft sources prove to be nearly equal. The cumulative present-day 0.25–2 keV emissivity of HMXBs with luminosities between and  erg s is  erg s  yr, which may be relevant for studying the X-ray preheating of the early Universe.

###### keywords:
stars: black holes – stars: massive – X-rays: binaries – X-rays: galaxies – galaxies: star formation – early Universe
pagerange: Bright end of the luminosity function of high-mass X-ray binaries: contributions of hard, soft and supersoft sourcesBright end of the luminosity function of high-mass X-ray binaries: contributions of hard, soft and supersoft sources

## 1 Introduction

One of the indicators of recent star formation activity in galaxies is the total number and cumulative luminosity of high-mass X-ray binaries (HMXBs, e.g. Grimm, Gilfanov, & Sunyaev 2003; Lehmer et al. 2010). The X-ray luminosity function (XLF) of HMXBs in nearby galaxies can be described by a power law across more than four decades in luminosity (at –8 keV energies) from to  erg s (Mineo, Gilfanov, & Sunyaev, 2012), with some indication of a cutoff at  erg s. This XLF thus smoothly connects low-luminosity HMXBs with ultraluminous X-ray sources (ULXs, usually defined as off-nuclear point-like objects with  erg s). It is possible that some ULXs, especially the most luminous ones, contain a black hole of intermediate ( ) rather than stellar mass, but the compact object is nevertheless expected to be accreting matter from a massive donor star in a close orbit, hence such systems may still be called HMXBs in a broad sense.

The slope of the HMXB XLF implies that more than 80% of the total X-ray emission from point sources in actively star-forming galaxies is produced by HMXBs with  erg s, with ULXs alone contributing more than half. The X-ray output of the young stellar population is thus strongly dominated by neutron stars and black holes accreting at near- and supercritical accretion rates. Such objects are very interesting but rare, and their population properties can be disclosed by accurately measuring the bright end of the HMXB XLF.

Knowledge of the shape and normalisation (i.e. number of objects per unit star-formation rate, SFR) of the HMXB XLF, as well as its possible dependence on other galactic properties such as metallicity (e.g. Linden et al. 2010; Douna et al. 2015), is also important for studying the ’cosmic dawn’, when X-ray binaries belonging to the first generations of (massive) stars and their remnants might have (together with other mechanisms, such as supernova cosmic-ray heating, Sazonov & Sunyaev 2015) significantly preheated the Universe (e.g. Venkatesan, Giroux, & Shull 2001; Madau et al. 2004; Ricotti & Ostriker 2004; Mirabel et al. 2011) before it was reionised by UV radiation from subsequent generations of stars and quasars.

Previous statistical studies of HMXBs were usually aimed at constructing an observed XLF, and paid litte attention to survey selection effects that may arise due to the diversity of HMXB X-ray spectral shapes and the presence of X-ray absorbing gas along the line of sight. This is a reasonable approach as long as one is dealing with prevalent low- and medium-luminosity HMXBs with  erg s, since most of such objects are X-ray pulsars with characteristic hard X-ray spectra, with the bulk of the luminosity emitted above 2 keV (see Walter et al. 2015 for a recent review). However, the situation is different for luminous ( erg s) HMXBs, as explained below.

First, the vast majority of such high-luminosity HMXBs are powered by Roche lobe overflow rather than by wind, with a significant and increasing with luminosity fraction of systems with a black hole as the compact object. In contrast to neutron-star HMXBs, whose X-ray luminosity is dominated by hard X-ray emission from magnetically and radiatively supported accretion columns on the surface of the neutron star (Basko & Sunyaev, 1976), black-hole X-ray binaries (both low- and high-mass ones) are known to experience dramatic spectral state transitions when the accretion rate is between % and % of the Eddington one and tend to be in a so-called high/soft state (sometimes also in a ’very high’ state) at subcritical accretion rates (see Done, Gierliński, & Kubota 2007 for a review).

In the soft state, most of the X-ray emission originates in a standard thin accretion disk (Shakura & Sunyaev, 1973) and has a multicolour black-body spectrum with characteristic temperature keV, where is the black hole mass and  erg s is the corresponding Eddington luminosity. A photon limited X-ray survey performed at 0.5–8 keV energies might have a better chance to detect such a source than one of the same luminosity but in a hard spectral state (due to the larger number of photons per unit X-ray flux in the former case) in the absence of absorption along the line of sight to the source. However, propagation through the interstellar medium (ISM) in the Galaxy and especially in the host galaxy may significantly suppress the X-ray flux from a source in a soft state and reverse the selection effect in favour of those in hard states.

Secondly, most ULXs are probably stellar-mass black holes accreting matter from a massive donor at a supercritical rate, although a non-negligible fraction of such systems may contain a neutron star rather than a black hole, as witnessed by the well-known pulsars LMC X-4 and SMC X-1 with  erg s and the recently discovered pulsars M82 X-2 (Bachetti et al., 2014) and NGC 7793 P13 (Fuerst et al., 2016; Israel et al., 2016) with  erg s. In this case, accretion onto the black hole proceeds through a thick disk. When viewed through the funnel of the disk with a strong wind outflowing from it, such a system is perceived as a ULX and exhibits a characteristic hard spectrum with a downturn above –10 keV (e.g. Sazonov, Lutovinov, & Krivonos 2014; Walton et al. 2014; Rana et al. 2015; Krivonos & Sazonov 2016). However, if the same system is observed from outside the funnel, through the optically thick wind, only reprocessed emission from the central regions can be visible, which, depending on the viewing direction and accretion rate, may emerge i) in the soft X-ray band or ii) at even longer wavelengths (e.g. Poutanen et al. 2007; Middleton et al. 2015; Urquhart & Soria 2016; Gu et al. 2016; Feng et al. 2016). So-called ultraluminous supersoft X-ray sources (ULSs) – rare objects in nearby galaxies with  erg s and thermal spectra with colour temperature  keV (Di Stefano & Kong, 2003) – may represent case i), while case ii) is probably realised in the Galactic microquasar SS 433 (Fabrika et al. 2015, see also a relevant discussion in Khabibullin & Sazonov 2016). The recent discovery of SS 433-like baryonic jets in a ULS in the M81 galaxy (Liu et al., 2015) provides additional support to this picture.

There may also exist ULSs containing an IMBH accreting in subcritical regime via a standard geometrically thin disk with low temperature  keV (since for a given ).

Measured X-ray fluxes of ULSs are very sensitive to line-of-sight absorption, and X-ray surveys may miss such objects despite their high intrinsic luminosity, due to the presence of interstellar gas in their host galaxies (Di Stefano & Kong, 2003). This raises a question: how many ULSs are there in nature compared to ULXs? The answer could provide strong constraints on models of supercritical accretion.

The purpose of this work is to construct an intrinsic (i.e. unabsorbed) XLF of HMXBs at  erg s and study its composition in hard and soft X-ray sources. To this end, we use a catalogue of X-ray sources detected by the Chandra X-ray observatory in nearby (within 15 Mpc) galaxies. Since its launch in 1999, Chandra has systematically observed many local galaxies and detected thousands of X-ray sources in them. On the other hand, recent UV, infrared, 21 cm and CO surveys have provided high-quality maps of the SFR and atomic and molecular ISM in the same galaxies. We combine this information with the Chandra data to infer the shape and normalisation of the bright end of the intrinsic XLF of HMXBs.

Our study consists of three stages. First, we perform a spectral analysis of bright Chandra X-ray sources to infer their intrinsic luminosities, , in the 0.25–8 keV energy band (and the photoabsorption columns, , affecting the measured spectra). We then use the H, H and SFR maps of the studied galaxies to study how the total SFR ’probed’ by Chandra depends on and source spectral type. Our underlying assumption is that the number of HMXBs follows star formation activity in galaxies. At the final stage, we build an intrinsic XLF of HMXBs by dividing the numbers of Chandra sources of different spectral types in specified bins by the corresponding SFR () functions and adding together the contributions of hard and soft sources.

Making allowance for line-of-sight absorption effects is especially important in the soft X-ray range. We thus also construct an intrinsic HMXB XLF in the 0.25–2 keV band. We use this soft XLF elsewhere (Sazonov & Khabibullin, 2017) for estimating the X-ray preheating of the early Universe.

## 2 Galaxy sample

This study is based on galaxies from The HI Nearby Galaxy Survey (THINGS, Walter et al. 2008), which is a high spatial resolution () 21 cm survey of nearby galaxies performed by the Very Large Array. The original THINGS catalogue comprises 34 galaxies located at distances  Mpc and north of , covering a wide range of Hubble types, star-formation rates, absolute luminosities and metallicities.

Chandra has observed all THINGS galaxies except for NGC 2366. Apart from this object, we also excluded from consideration i) NGC 3077, because of problems with determination of its SFR111A bright star contaminates the GALEX UV image., and ii) five dwarf galaxies, DDO 53, DDO 154, Ho I, M81 DwA and M81 DwB. The total SFR is likely less than 0.1  yr for NGC 3077 and is below 0.01  yr for each of the excluded dwarf galaxies (Leroy et al., 2008; Walter et al., 2008), hence the combined contribution of these galaxies to the HMXB population of the THINGS sample is expected to be negligible222We nevertheless searched, according to our criteria defined below, for bright Chandra X-ray sources in these excluded galaxies and found none..

We have chosen the THINGS galaxy sample for this study not only because of the availability of homogeneous HI data but also because there exist similarly high-quality and sufficiently extended H and SFR maps for these galaxies. The former are provided by the Heterodyne Receiver Array CO Line Extragalactic Survey (HERACLES, Leroy et al. 2009), an atlas of CO emission with angular resolution obtained with the IRAM 30-m telescope, while the latter can be readily constructed from Galaxy Evolution Explorer (GALEX) far-UV (1350–1750 Å) maps and Spitzer/MIPS 24 m maps (mostly from the Spitzer Infrared Nearby Galaxy Survey, SINGS, Kennicutt et al. 2003). The angular resolution (FWHM) of these UV and IR maps is and , respectively (see Bigiel et al. 2008). GALEX and Spitzer data are available for all of our galaxies, whereas HERACLES data exist for all but 7 objects: NGC 1569, NGC 3031, NGC 3621, NGC 4449, NGC 4826, NGC 5236 and NGC 7793. As a result, we use HERACLES H maps for most of our galaxies and construct synthetic H maps from GALEX/Spitzer SFR maps for the 7 galaxies listed above, using a well-established correlation between SFR and H surface density (see Section 4.1 below).

Our sample thus consists of 27 nearby galaxies. Table 1 provides key information on these objects. Distances (), inclinations (), position angles (P.A.) and optical sizes (-band isophotal radii at 25 mag arcsec, ) are adopted from the THINGS description paper (Walter et al., 2008). Galaxy types are taken from the HyperLeda database. Estimates of the total SFR and stellar masses () are adopted from Leroy et al. 2008, except for 7 galaxies absent from their sample. For these, we have added to Table 1 our own estimates of the SFR (obtained, similar to Leroy et al. 2008, using GALEX and Spitzer maps, see Section 4.1) and . The latter were calculated from the -band luminosities given in the Updated Nearby Galaxy Catalogue (Karachentsev, Makarov, & Kaisina, 2013), assuming a mass-to-light ratio , the same as used by Leroy et al. (2008) for the other galaxies. We need total stellar masses to estimate (in Section 6.2) the contribution of low-mass X-ray binaries (LMXBs) to the derived XLFs. Finally, Table 1 includes information on the Galactic HI column density in the direction of the galaxies (Kalberla et al., 2005).

The sample consists mainly of spiral galaxies, although it also includes 5 smaller, irregular galaxies. The total SFRs range from   yr (Ho II) to   yr (NGC 5194, NGC 5236, NGC 5457, NGC 6946 and NGC 7331). The combined SFR of all 27 galaxies is   yr.

## 3 X-ray sources

We used the recently published catalogue of Chandra/ACIS X-ray point sources by Wang et al. (2016). This is the largest publicly available catalogue of Chandra sources, based on the Chandra data archive covering years of observations. Importantly, the sources in this catalogue had been selected based on their detection in the 0.3–8 keV energy band, which reduces bias against soft X-ray sources compared to catalogues based on source detection in harder X-ray bands.

Most of the sources in the Wang et al. (2016) catalogue are weak detections, with just counts in an individual ACIS observation. Since our goal is to estimate unabsorbed X-ray luminosities of sources, which requires measurement of the line-of-sight absorption column from their spectra, we decided to use only sources with at least 100 counts in some observation. A number of galaxies in our sample (Table 1) have been observed more than once by Chandra, and the Wang et al. (2016) catalogue provides information on individual observations for a given source. In what follows, we use the ACIS observation with the highest counts for a given source. However, in Section 6.4 we estimate the influence of this selection procedure on our results using information from all existing Chandra observations for our sources.

We cross-correlated the Wang et al. (2016) catalogue with the positions of the galaxies in our sample and selected X-ray sources located within the 25 mag arcsec isophotes of the galaxies, i.e. at radii in the plane of the galaxy444Note that M51b (NGC 5195) and Holmberg IX, the well-known satellites of NGC 5194 (the Whirlpool Galaxy) and NGC 3031 (M81), respectively, do not fall into their regions and hence are not studied here., but outside . With this choice of the outer boundary we wish to include most of the X-ray sources actually belonging to the galaxy (since usually more than 90% of the total SFR is contained within ) and minimise contamination by background active galactic nuclei (AGN). The condition is introduced to avoid complications related to possible supermassive black hole activity in the nuclei of the studied galaxies and because galactic central regions are often characterized by highest ISM column densities, which is difficult to properly take into account in our analysis. Also, although source confusion plays a minor role in our study (see below) due to the excellent angular resolution of Chandra, it may become problematic in the nuclear regions of the galaxies.

We additionally filtered out sources that are reported by Wang et al. (2016) to be partially confused with another source in a particular Chandra observation. Sometimes, the same source may be reported to be unconfused in another observation, but we nevertheless excluded all such sources from consideration to have better confidence in the results of our X-ray spectral analysis. We also excluded several sources for which the total counts in an X-ray spectrum that we extracted (see below) differed by more than a factor of 2 from the value reported by Wang et al. (2016) for the same Chandra observation. In all such cases, there proves to be another X-ray source in the vicinity, so this discrepancy is yet another manifestation of source confusion. In total, 34 sources (% of all sources) with 100 or more counts have been filtered out.

The remaining sample consists of 330 sources with at least 100 counts, located between and in the 27 studied galaxies.

### 3.1 Spectral modeling

We used all available archival Chandra/ACIS data to extract X-ray spectra for the selected sources. Data search and reduction was performed by means of standard tools from the ciao (v. 4.8) package. Namely, we took advantage of the specextract thread with source and background extraction regions defined as follows: the source region is the -radius circle around the source, while the background region is the -radius circle centred on the source from which all source regions (for the source under consideration and any other sources) are excised. The technique works well sufficiently far from the crowded central regions of the investigated galaxies, which we have already excluded by imposing the condition . We verified (using ciao psfsize_srcs tool) that the -radius extraction region encompasses at least 2/3 of the source photons even for the most off-axis ( arcmin) cases in our sample and more than 90% of the source photons for 97% of a subsample that is used for construction of XLFs below.

Using the extracted spectra along with the corresponding response and aperture correction functions, we performed spectral analysis and measured source fluxes in the 0.25–8 keV energy band, which is only slightly different from the 0.3–8 keV band used by Wang et al. (2016) for source detection. We extended the energy range down to 0.25 keV to minimise loss of information for sources with the softest spectra.

The purpose of our spectral analysis was to i) characterize the spectral hardness/softness of sources, and ii) to determine their intrinsic (i.e. corrected for line-of-sight absorption) X-ray fluxes and luminosities. These tasks are not independent and in fact can only be solved together, since we do not know the positions of our sources with respect to the ISM in their host galaxies along the line of sight. Therefore, we need to determine both the intrinsic spectral shape and the line-of-sight column density, , for a given source from its X-ray spectrum. This is a notoriously difficult task, especially for very soft sources (see a relevant discussion in Devi et al. 2007), even despite the fact that we pre-selected sources with at least 100 counts.

After a number of trials we decided to take the following approach. We fitted the spectrum of each source, using xspec (Arnaud, 1996), by three alternative models: 1) absorbed power law (hereafter, PL), 2) absorbed black-body spectrum (hereafter, BB) and 3) absorbed multicolour disk emission (hereafter, DISKBB) – wabs(powerlaw), wabs(bbody) and wabs(diskbb) in xspec terminology. These models are similar in the sense that each has just 3 parameters: the absorption column density, , a parameter describing the intrinsic spectral shape (slope for PL, temperature for BB and the inner disk temperature for DISKBB) and the normalisation (we actually took advantage of the xspec convolution model cflux which calculates the unabsorbed flux in a given energy band, 0.25–8 keV in our case). We allowed to range between 0 and 10, and from 0.01 to 10 keV, and from to  cm.

As regards the absorption model, we also tried to use tbabs (Wilms, Allen, & McCray, 2000) instead of wabs (Morrison & McCammon, 1983), but found the results to be relatively insensitive to this (using tbabs in conjunction with the abundance set of Wilms, Allen, & McCray 2000 results in % higher estimates), while wabs provides a strong advantage in terms of computational time. We also note that differences in of the same order may arise from variations in the ISM metallicity between different galaxies or within a given galaxy, none of which is taken into account in our analysis.

We used -statistics (cstat in xspec, Cash 1979) for comparing the results of spectral fitting by the different models. For very soft sources, we applied an additional criterion (see below) to finally select the best model.

Our choice of spectral models is not arbitrary. Both DISKBB and PL models are commonly used for spectral analysis of X-ray binaries in their different spectral states, whereas BB often adequately describes the spectra of supersoft sources such as ULSs. Theoretically, DISKBB spectra can be produced by standard accretion disks, PL spectra can naturally arise in coronae of accretion disks and jets, and BB spectra can emerge from supercritical accretion disks with winds. Furthemore, PL should fit well the spectra of AGN, which may contaminate our sample of sources.

However, our single-component modeling is certainly oversimplified for some spectra, especially for brighter sources with counts. We nevertheless consider this approach adequate for the purposes of this study, first because it is practically impossible to reliably deduce more than 3 spectral parameters for the vast majority of our sources, and also because even for the spectra of some brighter sources that seem to require fitting by a multicomponent model (e.g. a combination of PL and BB) it usually turns out (as demonstrated below) that one of these components clearly dominates and its parameters can be inferred with reasonable accuracy by fitting the corresponding single model to the data.

### 3.2 Comparison of different spectral models

As expected, there is strong correlation between the best-fit values of , and for our sources, as shown in Fig. 1. The crucial question is which of the considered spectral models should we prefer for a given source? To address this issue, we show in Fig. 2 the difference in the cstat values, , for the PL, BB and DISKBB models as a function of the best-fitting power-law index.

For % of the softest spectra with , PL provides a worse fit compared to one of the thermal models, with 39 and 11 out of 63 spectra being best fit by BB and DISKBB, respectively. This result is anticipated, since such steep power-law spectra are observed very rarely from astrophysical X-ray sources, with both synchrotron and inverse Compton mechanisms usually giving rise to flatter spectra.

Nevertheless, 13 of the 63 sources with are better described by PL than by BB or DISKBB. However, the improvement in the fit has low statistical significance (), according to standard criteria such as the Akaike information criterion (AIC)555In our case, , where for all three considered spectral models, so a difference between two models implies that one of them is % as probable as the other., except for one source – CXOGSG J203500.1+600908666Hereafter, we use source names from Wang et al. (2016).. In this case, a steep PL model (with ) provides a clear improvement: . We explored the spectrum of this source (see Fig. 3) and found that the addition of a relatively hard () PL component to the supposedly dominant soft thermal (DISKBB) one leads to a fit of similar to the single-component PL model quality. Moreover, if PL is used as a secondary rather than single component, its contribution to the total unabsorbed luminosity proves to be small (%), while the temperature of the dominant thermal component does not change much compared to modeling by DISKBB alone ( vs. 0.31 keV), with the total unabsorbed luminosity increasing by a factor of . In contrast, the PL fit implies unrealistically high values of and .

For the above reasons, we imposed an additional condition on sources: even if PL is the statistically best model, we select the better of BB and DISKBB as our preferred model.

For softer spectra with , there is no clear preference between the models (see Fig. 2): 66 (25%), 48 (18%) and 153 (57%) out of 267 spectra are best fit by PL, BB and DISKBB, respectively. Since all these models are reasonable, we did not introduce any additional criteria for sources and picked their best spectral models based on -statistics alone. However, for relatively soft sources () there is a risk that choosing the wrong model from two or three models with similar values may lead to a large error in the absorption column and hence unabsorbed luminosity of the object. We made a special test to estimate the impact of this systematic uncertainty on our results (see Section 6.3).

### 3.3 Luminous sample

Based on the results of the spectral analysis, we compiled a sample of sources (hereafter, the ’luminous sample’), 219 in total, with estimated unabsorbed luminosities (assuming the sources indeed belong to the corresponding galaxies)  erg s (0.25–8 keV). This sample provides the basis for the subsequent analysis (in Sections 5 and 7). We have not extended our study to lower luminosities partially because the problem of separating HMXBs from other types of sources such as LMXBs and AGN, is greatly exacerbated at  erg s.

## 4 Survey coverage as a function of source luminosity and spectrum

We do not expect our sample of X-ray sources to directly represent the intrinsic distribution of luminosities and proportion of hard and soft sources, because it is based on a photon limited survey. Even in the absence of line-of-sight absorption, detection of an X-ray source with a given luminosity in a given galaxy must depend on the source’s spectral properties and the Chandra/ACIS spectral response. This may bias the resulting sample of sources in favour of soft ones. On the other hand, the presence of interstellar gas in the Milky Way and especially in the host galaxy in the direction of the source may significantly lower its chances of being detected, and this negative bias is expected to mainly affect soft sources. The situation is exacerbated by the fact that HMXBs are usually located close to sites of active star formation, which in turn are correlated with concentrations of cold molecular gas.

As mentioned before, the unique property of the sample of nearby galaxies studied here is that their high-quality Chandra X-ray maps can be superposed on similarly high-quality maps of the SFR and atomic and molecular gas. This makes it possible, under reasonable assumptions about the ISM distribution perpendicular to the plane of the galaxy (see §4.2 below), to evaluate the X-ray source selection effects discussed above and eventually integrate the SFR over those 3D regions of the galaxies where Chandra can detect a source with given intrinsic luminosity and spectral type. In what follows, the so-derived SFR will be sometimes referred to as the SFR probed by Chandra. By dividing the number of actually detected X-ray sources of a given spectral type in a given luminosity bin by the corresponding SFR probed by Chandra, we later construct XLFs of different types of sources per unit SFR (§7).

We assumed that the spatial distribution of HMXBs in a galaxy approximately follows that of its current SFR. As known from HMXB surveys in the Milky Way (Bodaghee et al., 2012; Lutovinov et al., 2013), the two are correlated on a spatial scale of few hundred pc. This is the distance that a HMXB can travel (due to a kick received from the supernova explosion of the progenitor of the relativistic compact object) from its birth site over its lifetime of several million years. By the same argument and as also confirmed by Galactic observations (Lutovinov et al., 2013), HMXBs are expected to be distributed approximately uniformly in the vertical direction within  pc of the central plane of their galaxy. It is worth noting in this connection that the Galactic microquasar SS 433 is located  pc away from the Galactic plane.

The HI component of the host galaxy ISM, again by analogy with the Milky Way (e.g. Nakanishi & Sofue 2016), is expected to have a vertical scale height of a few hundred pc (which may strongly depend on the galactocentric distance). Therefore, the distribution of HMXBs along the normal to the galactic plane should be approximately random with respect to the atomic interstellar gas.

The situation is somewhat different with the colder, H, phase of the ISM. For the Galaxy it is known that a substantial or even dominant fraction of H is concentrated in giant molecular clouds (GMCs), but there is a broad distibution of cloud sizes and masses. However, even GMCs may be regarded as small objects ( pc) compared to the other relevant spatial scales discussed above. The molecular gas is expected to be concentrated near the galactic plane, with a vertical scale  pc (as in the Milky Way, Heyer & Dame 2015). All this suggests that i) one should not expect HMXB positions in a given galaxy to correlate with the positions of individual molecular clouds, although it is natural to expect some correlation between HMXB locations and the large-scale distribution of H in the galaxy (because the latter is correlated with the SFR); ii) to a first approximation, HMXBs are as likely to be located between the H disk of their galaxy and the observer as behind the H disk.

### 4.1 SFR and ISM maps of the galaxies

We used a linear combination of GALEX UV and Spitzer/MIPS 24 m intensity maps to produce maps of the SFR surface density, , for our galaxies. Our analysis closely followed Leroy et al. (2008) and Bigiel et al. (2008), and in fact the maps for 20 of our 27 galaxies have already been presented by Leroy et al. (2008).

Specifically, we calculated according to equation (1) in Bigiel et al. 2008:

 ΣSFR[M⊙ yr−1 kpc−2] = 3.2×10−3I24[MJy ster−1] + 8.1×10−2IFUV[MJy ster−1] ,

where I and I are the 24 m and FUV intensities, respectively, derived directly from Spitzer and GALEX images. We additionally corrected the FUV intensities for the Galactic extinction assuming (Yuan, Liu, & Xiang, 2013) and adopting from Schlafly & Finkbeiner (2011) (using IRSA Galactic Dust Reddening and Extinction Service). Note the significant difference between the ratio adopted here (4.64) and that used by Bigiel et al. (2008) (8.24, as in Wyder et al. 2007); see Yuan, Liu, & Xiang (2013) for a relevant discussion.

For most of our galaxies, we used HERACLES CO line intensity (, K km s) maps (Leroy et al., 2009) to produce maps of the H surface density, . This was done using equation (4) in Leroy et al. (2009), which reduces to

 ΣH2[M⊙ pc−2]=5.5 ICO[K km s−1], (2)

assuming the CO to H conversion factor cm (K km s) and the CO to ratio (see Leroy et al. 2009 and references therein for more details). Note that makes allowance for the contribution of helium.

Seven of our galaxies have not been covered by the HERACLES survey (see Section 2). For these, we constructed synthetic maps from the SFR maps using a well-established correlation between and (Bigiel et al., 2008; Leroy et al., 2008):

 ΣSFR/ΣH2=(5.25±2.5)×10−10 % yr−1. (3)

We used the same method to estimate H surface densities in the outer regions of some galaxies if the HERACLES maps were insufficiently large (usually, is negligibly small in these regions, though). As has been shown by Bigiel et al. (2008), the correlation appears to work universally well in spiral galaxies in a broad range of H surface densities, –50  pc, i.e.  H molecules cm, with the interpretation being that H forms stars at a constant efficiency. However, this correlation may not be applicable both below and above this range. The former case does not present any problem, since usually regions with   pc are dominated by HI (Bigiel et al., 2008; Leroy et al., 2008). The   pc limit pertains to starburst regions, which are usually the central regions of some of our galaxies, and our use of the correlation in this regime is limited because we have excluded the regions from the analysis.

The resulting , and (or ) maps for the studied galaxies, with the positions of the X-ray sources from our clean sample (defined below in Section 6) superposed on them, are presented in Appendix B.

We next constructed radial profiles (with a bin size of ) of , and from the corresponding maps. They are are in good agreement with (Leroy et al., 2008), where most of our galaxies have been analysed before. For the galaxies with available HERACLES data, there is satisfactory agreement between the directly measured H surface densities and the SFR-based estimates (see further discussion in Appendix B), which justifies our usage of synthetic maps as a substitute for missing CO data.

### 4.2 Simulations of X-ray absorption in the ISM

Using the derived SFR and ISM radial profiles for the galaxies, we performed simulations of Chandra observations of X-ray sources of various spectral types with different intrinsic luminosities and locations within the galaxies.

For each galaxy, we chose random celestial positions within each of the five -wide bins (elliptical rings in the sky plane) within and simulated spectra that could be measured in a particular existing Chandra/ACIS observation from X-ray sources located in these positions, depending on the sources’ intrinsic spectral shape and line-of-sight absorption. The 10 different positions per radial bin are needed to take into account possible differences in the Chandra response from one position to another. We used the same three spectral models, wabs(powerlaw), wabs(bb) and wabs(diskbb), as we used before in the spectral analysis of the real Chandra sources in our sample, and ran simulations on a large grid of spectral parameter values: ,  keV,  keV and  cm.

We then assumed that an X-ray source with a given position on the sky can be located quasi-randomly along the line of sight within the ISM of its host galaxy. Namely, motivated by the discussion at the beginning of Section 4, we modeled the HI component of the galaxy as a thick homogeneous slab whereas we assumed the H gas to be concentrated in GMCs located in the central plane of the galaxy (see a schematic in Fig. 4).

In this model, an X-ray source falling into a given elliptical ring of the galaxy, with the ring-averaged HI column density , has a probability to have a column of atomic gas between itself and the observer, with . Here, is the Galactic HI column density in the direction of the source, adopted from Kalberla et al. (2005) (we use the nh tool from heasoft ftools).

The same X-ray source has a 50% chance of being located in front of the galactic plane and thus at least 50% probability to have no GMC between itself and the observer. If the source is located behind the galactic plane, its view can be obscured or not by a GMC, as discussed below.

We assumed that a typical column density through a GMC is  H molecules cm, or, equivalently,   pc (taking into account helium). This is the median mass surface density of Galactic GMCs determined by Heyer et al. (2009) using the Boston University–FCRAO Galactic Ring Survey of CO emission, but a factor of smaller than 170  pc, the value from Solomon et al. (1987) (with a small correction by Heyer et al. 2009) which was for a long time regarded as standard in astrophysics. Heyer et al. (2009) explain a number of reasons (including the use of CO  instead of CO  and improved angular and spectral resolution) that led to this significant revision. However, Heyer et al. (2009) point out that their new reference value for may be underestimated by a factor of .

We implemented the following algorithm. If the radial-bin averaged is smaller than , we posit that a source behind the galactic plane has a chance of having a GMC between itself and the observer. In this case, a fraction of such sources will have an H column toward the observer, while the remaning will have . If , we adopt that and , i.e. use the bin-averaged instead of as the H column density in front of the source. The last assumption may reflect a situation when several GMCs form a larger molecular complex (in a starburst region) and/or there are individual GMCs with column densities higher that our fiducial value of  cm. There is clearly significant uncertainty here, and we performed some tests (see below) to estimate its effects on our results.

Our final assumption is that the total number of X-ray sources of a given kind (i.e. intrinsic luminosity and spectral type) in a given radial bin is proportional to the total SFR () in that bin. This assumption is reasonable, since HMXBs generally follow ongoing star formation activity in galaxies, as mentioned before.

Using the above assumptions and prescriptions, we evaluated for each galaxy in our sample the SFR probed by Chandra as a function of X-ray source intrinsic luminosity and spectral type. Just as for real sources in our sample, a simulated source was considered to be detected if it provided at least 100 counts in any of the available Chandra/ACIS observations for the given position within the galaxy. We integrated the SFR over . We assumed that the innermost region contributes 1/16th of the SFR encompassed within and subtracted this amount from this radial bin. This is unlikely to lead to a significant error.

Figure 5 shows the resulting SFR () dependencies for the PL, BB and DISKBB models, summed over all 27 galaxies.

As could be expected, at very high luminosities  erg s Chandra probes the total SFR of the 27 galaxies of   for all spectral types of sources, except for the softest ones ( keV or  keV). This is because very high ISM column densities are needed to suppress the observed X-ray flux of a luminous source to less than 100 counts, given that all of the considered galaxies are located within 15 Mpc. However, the difference between various spectral types of sources becomes noticeable below  erg s and grows with decreasing luminosity. As also aniticipated, for a given , the integrated SFR first increases going from hard to soft sources (e.g., for the PL model this occurs for ) due to the increasing number of photons per unit flux, but the opposite trend sets in for yet softer spectra because of the strong suppression of soft X-rays in the ISM.

The largest systematic uncertainty in our simulations is associated with the distribution of H in the galaxies. It is clear that our assumption that all molecular gas seen in the CO maps is located in GMCs of fixed column density simplifies the real situation. We therefore tested the sensitivity of our SFR () curves to this assumption.

First, we tried a model where there is no absorbing gas in the direction of the X-ray sources. In this imaginary scenario we, as expected, see (in Fig. 6) a large enhancement of the total SFR probed by Chandra for sources with soft and supersoft spectra (PL with and BB with  keV, respectively) and a noticeable increase (% at  erg s) for hard sources (). As a second test, we retained the atomic gas but removed all molecular gas from the modeled galaxies. In this (also unrealistic) scenario, the impact of absorption on the cumulative SFR is still substantial for soft and supersoft sources.

We next tried using different values of in our ’standard’ model and found that the SFR does not change significantly if is increased from to  cm (see Fig. 6).

Finally, we implemented a more radical change of the model, namely distributed all the molecular gas seen in the CO map of a given elliptical ring homogeneously over the volume of the galaxy subtended by this ring, i.e. essentially allowed H to be mixed with HI. This model leads to a much stronger negative bias for soft sources (see Fig. 6) compared to our baseline model. This has a clear explanation: in our standard model, there is always at least a 50% chance (due to the assumed location of all GMCs in the galactic plane) for sources of any spectral type to have no molecular gas obscuring their view. This probability is in fact often higher, because even if a source is located on the other side of the galactic plane there may still be holes in the molecular layer above it. In the alternative scenario, all sources are immersed in the HI+H ISM and their X-ray fluxes inevitably suffer some attenuation in the molecular gas in addition to absorption in the atomic gas.

Although we find the above ’homogeneous H model’ unrealistic, it indicates that our baseline absorption model is probably too conservative, i.e. it may somewhat underestimate the negative bias associated with soft sources if the molecular gas has a more fragmentary structure than assumed in our baseline model.

## 5 Removal of contaminants

We cross-correlated our luminous sample with standard astronomical databases in an attempt to find out the nature of our sources. Although some published catalogues ascribe detailed classifications, e.g. ’HMXB’ or ’LMXB’, to many of our sources, we considered such classifications unreliable, because they are usually based solely on the location of the source within its presumed host galaxy (e.g. a source within the bulge of a galaxy would be designated a LMXB). Similarly, many of our most luminous ( erg s) sources are marked ’ULX’ in various catalogues, but we do not use this purely empirical classification either.

Specifically, we searched for objects located within 3 arcsec of the Chandra positions of our sources. Since all these sources are bright (more than 100 counts), their X-ray positions are known to better than , except for 3 sources (observed at large offset angles) for which the uncertainties are between and (Wang et al., 2016).

### 5.1 Foreground stars

We looked for the possible presence of foreground stars in the luminous sample and found 3 likely associations (see Appendix A).

### 5.2 Background AGN

We next looked for the presence of AGN in the luminous sample and found 8 likely associations, 7 of which may be considered reliable (see Appendix A). However, the total number of CXB sources contaminating the sample may be higher.

To estimate this number, we adopted the AGN curve in the 0.5–10 keV band from (Georgakakis et al. 2008, their equation (2) and table 2; see also Kolodzig et al. 2013 for a comparison with other published curves) and made a small correction from 0.5–10 keV to our operative energy band of 0.25–8 keV. We assumed all CXB sources to have power-law spectra with , which is approximately equal to the slope of the CXB spectrum and is effectively an average between the spectra of intrinsically unobscured (type 1) AGN and the harder spectra of type-2 objects (see, e.g., Sazonov et al. 2008). We then took into account that the X-ray flux from an AGN located behind any of our studied galaxies will suffer some absorption in its ISM (and in the Galaxy), by performing simulations similar to our modeling of SFR () above.

We found that there are expected to be 19.7 background AGN in our luminous sample. This estimate is quite robust due to the hardness of CXB sources. In fact, we also tried to completely switch off line-of-sight absorption in our model, which led to a negligible enhancement of the expected number: 20.0. Therefore, the luminous sample likely contains background AGN, 8 of which we have already identified. Hence, CXB sources likely remain hidden in the sample.

### 5.3 Measured vs. expected absorption columns, further AGN candidates

It is interesting to compare the columns inferred from our spectral analysis of the Chandra sources with the total column density of gas in their direction: (where and are measured in H atoms cm and H molecules cm, respectively). The latter information can be taken from the same HI and H (or synthetic H) maps that we used before for computing the SFR () dependencies for the galaxies.

Figure 7 shows the result of comparison of with for the luminous sample, excluding the 11 objects associated with either stars or AGN mentioned above, i.e. for 208 sources in total. Disregarding the (significant) statistical uncertainties associated with the measurements, 105 (55%) and 163 (78%) sources lie below the and lines, respectively. Taking the uncertainties into account, 181 (87%) sources are consistent with having .

The above numbers appear to be in fairly good agreement with the expectations of our ISM absorption model (see Section 4.2). Indeed, the total line-of-sight gas column density for most of our sources proves to be dominated by H in the host galaxy: the mean and median values of are 5.5 and  cm, respectively, whereas the corresponding values for are 1.6 and  cm. Hence, our typical source has  cm in its direction, which corresponds to –1 GMC per line of sight (for our adopted  cm). Taking into account that the source may be located either in front or behind the molecular disk of its galaxy, we may expect –50% of our sources to be obscured by molecular gas and hence have .

The satisfactory agreement of the measured X-ray absorption columns with the expectations of our ISM absorption model implies that the absorption evident in the X-ray spectra of most of our sources is indeed due to the ISM in their host galaxies, rather than due to gas intrinsic to the sources. This result fits well into the assumed picture that most sources with  erg s are binary systems where accretion onto the compact object proceeds via Roche-lobe overflow of the (massive) companion, rather than through its wind. If the opposite were true, we would expect to see a lot of additional absorption in the wind. Such absorption, with up to  cm, is indeed observed in many Galactic HMXBs, in particular those discovered with the INTEGRAL observatory (Walter et al., 2015), but all such systems have  erg s.

However, we do see several clear outliers in the vs. diagram. Specifically, there are 8 sources (denoted by squares in Fig. 7 and listed in Appendix A) that are located above the line, taking into account the uncertainties. If they belonged to the explored galaxies (3 objects in NGC 5194 and 5 objects in NGC 5457), the bulk of the absorption observed in their spectra would have to be intrinsic to the sources rather than arise in the ISM (even allowing for metallicity variations). Although this possibility cannot be excluded, we may also consider an alternative explanation that some or all of these objects are background AGN with intrinsic absorption. In fact, the spectra of 6 of them have slopes consistent with if fitted by an absorbed PL model, which is typical for AGN (the remaining two sources have somewhat softer spectra with the 1 lower limits on being 2.3 and 2.5). Since we previously concluded there must remain unidentified AGN in our luminous sample, of which a substantial fraction should have significant X-ray absorption (as is common for AGN), the 8 absorption-diagram outliers are more likely to be these missing AGN compared to other sources in the sample. We therefore excluded these 8 objects from further consideration. This disputable decision has a very limited impact on the results of our study since if affects just % of the entire sample of sources.

## 6 Clean sample

After exclusion of 3 foreground stars, 8 background AGN and another 8 sources suspected to be AGN based on their high X-ray absorption columns, we obtained a ’clean sample’ of 200 sources (Table 2), the vast majority of which presumably belong to the 27 nearby galaxies under consideration.

It is unlikely that all sources in the clean sample are HMXBs (and ULXs). First, it definitely contains a significant number of LMXBs. Unfortunately, it is practically impossible, even in nearby galaxies, to separate HMXBs from LMXBs on a source by source basis. We therefore made an attempt to estimate the contribution of LMXBs in a statistical way, as described below (Section 6.2).

In addition, our clean sample may contain supernovae (SNe) recently exploded in the studied galaxies. We indeed found 4 SNe positionally coincident with our X-ray sources (see Table 2). These are all core-collapse ones, and in fact were the targets of the Chandra observations used in our analysis. Three of these observations (for SN 2002hh, SN 2004dj and SN 2004et) were performed within two months after discovery, whereas the observation of the famous SN 1993J was taken in 2000, i.e. 7 years after the explosion (see Dwarkadas & Gruszko 2012 for the X-ray light curves of these SNe). It is well known that young core-collapse SNe can remain luminous ( erg s) X-ray sources for months, years and sometimes decades after the explosion (e.g. Immler & Lewin 2003; Dwarkadas & Gruszko 2012). For this reason and because the 4 identified SNe constitute just % of our clean sample, we did not exclude them from further analysis despite them being the Chandra targets. It is actually interesting to retain these objects in the sample because young core-collapse SNe are another manifestation of star-formation activity, complementary to HMXBs. As for the X-ray spectra of the 4 identified SNe, two of them (SN 2004dj and SN 1993J) are best fit by an absorbed PL among our simple one-component models, with and , respectively, one (SN 2002hh) by a BB model with  keV and one (SN 2004et) by a DISKBB model with  keV. We also tried to desribe these spectra in terms of optically thin thermal emission, wabs(apec), and obtained worse fits.

Finally, we found 15 possible associations of our sources with supernova remnants (SNRs, see Table 2). All of these SNRs had been identified as such based on their optical (Blair, Winkler, & Long, 2012; Vučetić, Arbutina, & Urošević, 2015) and/or radio (Chomiuk & Wilcots, 2009a, b) emission. However, the SNR is unlikely to be the actual source of the X-ray emission observed by Chandra in most of these cases. First, these SNRs are expected to be at least few hundred years old, while, to our knowledge, there are no well-documented cases of an SNR maintaining an X-ray luminosity in excess of  erg s for such a long time. On the other hand, some core-collapse SNRs are known to be physically associated with a bright HMXB, the best-known example perhaps being the pair of the Galactic microquasar SS 433 (believed to be a ’misaligned’ ULX) and the W50 nebula (SNR G039.7-02.0). A similar situation may well be realised in at least some of our 15 possible SNR–X-ray source associations.

### 6.1 Spectral properties of the sample

According to the results of our spectral analysis (see Section 3.1), the clean sample consists of the following spectral types: there are 49 PL spectra with , 40 BB spectra with  keV (excluding one, exceptionally hard spectrum, the range is  keV), and 111 DISKBB spectra with  keV.

For the following treatment it is convenient to split the sources into three categories: ’hard’, ’soft’ and ’supersoft’. We define these according to the ratio of the unabsorbed flux in the 0.25–2 keV band, , to that in the total (0.25–8 keV) band, :

 Hard: F0.25−2,unabsF0.25−8,unabs≤0.6, (4) Soft: 0.60.95. (6)

This is a purely empirical division, with no direct relation to similar terminology frequently used in the literature to describe different types of X-ray sources. However, our hard class approximately corresponds to X-ray pulsars and black-hole X-ray binaries in their hard state, our soft class to black-hole X-ray binaries in their high/soft and very high states and our supersoft class to ULSs.

In practice, the hard category corresponds to for the PL model,  keV for BB and  keV for DISKBB. For the soft class, the corresponding ranges are ,  keV and  keV; and for the supersoft class they are ,  keV and  keV.

There are in total 119 hard, 50 soft and 31 supersoft sources in the clean sample. Figure 8 (upper panel) shows the distribution of unabsorbed 0.25–8 keV luminosities, , for these classes of sources as well as for the whole sample. We see that the relative fraction of hard sources is highest () in the lowest luminosity bin,  erg s, diminishes to in the next 3 bins and drops to in the highest luminosity bin ( erg s), although there are only 7 sources in total in this last interval. We also see that both soft and supersoft sources are important over the entire studied luminosity range.

As for the supersoft sources, we caution that there is very large uncertainty associated with the two objects of this type present in the  erg s bin. Both have extremely soft spectra (see Fig. 9), which are best fit by a heavily absorbed ( cm) BB model with –0.07 keV. As a result they have huge intrinsic/observed luminosity (0.25–8 keV) ratios of (see Table 2), and these absorption corrections are uncertain by 1–2 orders of magnitude. Both sources are located in the direction of gas rich regions of their presumed host galaxies, with and  cm, so the values inferred from our X-ray spectral analysis are not unexpected, although for one of the sources the lower limit on somewhat exceeds . Moreover, fitting by other models, namely absorbed PL and absorbed DISKBB, leads to even higher . Note that the lower luminosity bins in the histogram are not strongly affected by this kind of systematic uncertainty, since these intervals contain a significant number of supersoft sources with relatively small absorption corrections (see Table 2).

For comparison, the middle panel of Fig. 8 shows the corresponding distribution of observed (i.e. uncorrected for line-of-sight absorption) luminosities, , for the same sources. As could be expected, the relative number of soft and especially supersoft sources is smaller in this case. Finally, the bottom panel of Fig. 8 shows the distribution of soft-band (0.25–2 keV) unabsorbed luminosities, . In this case, we see the opposite and also expected trend: the relative numbers of soft and supersoft sources are higher.

### 6.2 Contribution of LMXBs

Our goal is to obtain the XLF of HMXBs. However, the histograms shown in Fig. 8 must inevitably contain a contribution of LMXBs. The XLF of LMXBs is known fairly well from studies of the bulge of the Galaxy (Revnivtsev et al., 2008), galactic bulges and early-type galaxies (Gilfanov, 2004), i.e. environments nearly free of HMXBs (apart from projection effects). We can use this information to estimate the contibution of LMXBs to our sample of sources.

We adopted the analytic form of the LMXB XLF from (Gilfanov 2004, namely their equations (8) and (9) and the parameters for ’All galaxies’ from table 3). In contrast to HMXBs, the number of LMXBs is proportional to the stellar mass of a galaxy rather than its SFR. We thus multiplied the XLF per stellar mass from Gilfanov (2004) by the stellar masses of our galaxies, given in Table 1. Our clean sample by construction contains objects located within the regions of the galaxies. Because these regions encompass nearly all of the stellar mass of the galaxies, we did not apply any correction to the normalisation, bearing in mind much larger uncertainties associated with the LMXB XLF itself (see relevant discussions in Gilfanov 2004 and Mineo, Gilfanov, & Sunyaev 2012).

The majority of Galactic LMXBs have hard spectra, with –2 if fitted by a power law (e.g. Revnivtsev et al. 2008). However, we are interested here in LMXBs with  erg s, which are very rare in the Galaxy. The majority of such luminous LMXBs are black-hole transients (mostly X-ray novae), which are known to experience transitions between hard and soft spectral states (e.g. Done, Gierliński, & Kubota 2007). Therefore, our clean sample may well contain LMXBs in both hard and soft states. Furthemore, among our lowest luminosity sources ( erg s) there may be present classical supersoft sources, associated with accreting white dwarfs, but such objects are not expected to have  erg s (e.g. Soraisam et al. 2016).

Since it is hardly possible to predict the relative numbers of hard, soft and supersoft sources among the LMXBs contaminating our HMXB sample, and in view of the substantial uncertainties associated with the LMXB XLF, we assumed that for all LMXBs. We further made a flux correction from 0.5–8 keV (the effective energy band in Gilfanov 2004, although the lower boundary of this range is somewhat uncertain in that study) to 0.25–8 keV and simulated the appearence of LMXBs in Chandra observations, as we did before in deriving the SFR ( curves and estimating the AGN contribution to the sample. Specifically, we assumed that LMXBs are distributed uniformly over the face-on image of a given galaxy within (to roughly take into account that LMXBs are more concentrated toward the centre of the galaxy compared to HMXBs, although the results prove to be almost insensitive to the choice of ) and randomly with respect to the ISM in the direction perpendicular to the plane of the galaxy. We thus assumed that LMXBs are subject to the same kind of X-ray absorption as HMXBs. We also assumed that the Gilfanov (2004) LMXB XLF is an intrinsic rather than observed one. This is a reasonable first-order approximation, since most of the Gilfanov (2004) galaxies are elliptical ones, presumably containing very little cold gas, although his sample also contains a few bulges of spiral galaxies (including NGC 3031 and NGC 5457, which are also present in our sample) and the measured fluxes of some of the Gilfanov (2004) X-ray sources may well be affected by absorption in the H disks of these galaxies.

The estimated contribution of LMXBs to the distribution for the clean sample is shown by the dashed line in Fig. 8 (upper panel). It proves to be substantial in the first two bins ( and  erg s) but vanishes at higher luminosities. We also estimated the LMXB contributions to the distributions of observed 0.25–8 keV luminosities (Fig. 8, middle panel) and intrinsic 0.25–2 keV luminosities (Fig. 8, bottom panel).

### 6.3 Robustness of the sample

As mentioned before, selection of a suitable spectral model for a given source in the presence of line-of-sight absorption is an unreliable procedure, which may lead to significant systematic uncertainties in estimating unobsorbed source fluxes.

To evaluate this uncertainty, we constucted an alternative sample of sources as follows. We retained the previous algorithm for very soft sources with (i.e. selected the better of the BB and DISKBB models based on -statistics, see Section 3.2 above), while for the harder sources () we introduced the logarithmic mean between the unabsorbed luminosities for the PL and BB models, , as a replacement for the corresponding to the statistically best spectral model (PL, BB or DISKBB). Note that usually, for a given source, . We thus obtained a sample of sources with  erg s, from which we excluded likely associations with stars and AGN (including dropouts on the vs. diagram), as we did before.

The properties of the alternative sample proved to be not very different from the clean sample. The total number of sources in the new sample is 214 vs. 200. The number of sources with luminosity higher than and  erg s is 46 and 5 vs. 42 and 7. The number of hard, soft and supersoft sources is 128, 55 and 31 vs. 119, 50 and 31.

This test gives us some confidence that the XLF derived below using the clean sample is not strongly affected by the systematic uncertainty in determining both the line-of-sight absorption column and intrinsic source spectral shape from X-ray spectra with limited photon statistics.

### 6.4 Impact of variability

As was explained in Section 3, our spectral analysis is based on the Chandra/ACIS observation providing the most total counts for a given source. Since some of the studied galaxies have been observed by Chandra more than once, our reported measurements of the fluxes of the sources from the same galaxy may represent different observations, and the sum of these fluxes may be biased with respect to the instantaneous total flux of the HMXBs in that galaxy.

We estimated this effect as follows. For each of the 200 sources in the clean sample, we compared the mean of its observed 0.25–8 keV fluxes (derived from fitting by an absorbed PL model), , in all available observations with the observed flux, , measured in the observation that was actually used in our analysis. Figure 10 shows the resulting distribution of ratios. We see that for all of the sources. The mean value of this ratio is 1.26. Of the total sample, 62 sources have only one archival Chandra observation and so have , while for the remaining 138 ones the mean and median values of are 1.41 and 1.18, respectively.

We conclude that our selection procedure (motivated by the goal of achieving the best quality of spectral modelling) led to an average flux increase of % for the studied sources compared to an ideal situation where we could use the same observation for all sources in a given galaxy. The relative smallness of the effect suggests that we often selected the longest of the available observations for a given source rather than one where its flux (i.e. count rate) was highest.

Since the HMXB XLF derived below is fairly flat (has a slope ), the ratio implies that we probably overestimated the XLF normalisation by %. We have thus corrected the resulting XLFs by this factor (i.e. divided them by 1.2).

## 7 X-ray luminosity function

Using the clean sample of sources obtained in Section 6 and the SFR () dependencies for different spectral types computed in Section 4, we constructed the intrinsic XLF of HMXBs at  erg s per unit SFR. Specifically, we calculated the XLF in binned form via the formula

 dNdlogLX,unabs=∑i1SFR (LX,unabs,i,Speci), (7)

where the summation is performed over sources falling into a given bin, and are the source’s luminosity and spectral type (for example a PL with ), respectively, and SFR () is read off from the corresponding SFR () curves (computed on a fine mesh of parameter values (, and ). Using equation (7), we also calculated the XLFs of hard, soft and supersoft sources (as defined before in equation (6)), by performing the summation over sources of the spectral type of interest. Note that the algorithm used here is analogous to the method commonly used in astronomy, with the SFR replacing . Figure 11 (upper panel) shows the resulting intrinsic XLF (0.25–8 keV) and its composition in hard, soft and supersoft sources.

We next subtracted from the derived XLFs the expected contribution of LMXBs. To this end, we assumed that the proportion of hard, soft and supersoft sources among the LMXBs hidden in the clean sample is the same as for the HMXBs in a given luminosity bin. We use this simple approach because of the substantial uncertainty associated with the LMXB XLF and its composition (see Section 6.2). Thorough investigation of this problem is beyond the scope of this paper. As a result, our estimates of the number density of HMXBs in the lowest two luminosity bins of  erg s and  erg s (where the estimated fraction of LMXBs among our sources is –50% and –40%, respectively, see Fig. 8) and the corresponding contributions of hard, soft and supersoft sources should be regarded with caution. At higher luminosities, LMXB contamination is negligible.

The lower panel of Fig. 11 shows the intrinsic XLF, as well as its hard, soft and supersoft components, upon subtraction of the expected LMXB contribution and correction for the variability bias evaluated in Section 6.4. The uncertainties in Fig. 11 combine those associated with application of equation (7) () and the Poisson uncertainties associated with the subtraction of LMXBs.

We derived intrinsic XLF is well fit by a simple power law with a slope of and normalisation given in Table 3). Furthermore, the intrinsic XLFs of hard, soft and supersoft sources are also well fit by power laws with the same (within ) slope, but with larger associated uncertainties. We can thus use the normalising constants given at  erg s in Table 3 to estimate the relative contributions of hard, soft and supersoft sources to the combined XLF, which prove to be %, % and %, respectively. Hence, soft and supersoft sources constitute together about half of all sources over the  erg s luminosity range.