A Fundamental Plane in X-ray Binary Activity of External Galaxies

We construct a new catalog of extragalactic X-ray binaries (XRBs) by matching the latest Chandra source catalog with local galaxy catalogs. Our XRB catalog contains 4430 XRBs hosted by 237 galaxies within ~130 Mpc. As XRBs dominate the X-ray activity in galaxies, the catalog enables us to study the correlations between the total X-ray luminosity of a galaxy $L_{X,\rm tot}$, star formation rate $\dot{\rho}_\star$, and stellar mass $M_\star$. As previously reported, $L_{X,\rm tot}$ is correlated with $\dot{\rho}_\star$ and $M_\star$. In particular, we find that there is a fundamental plane in those three parameters as $\log L_{X,\rm tot}={38.80^{+0.09}_{-0.12}}+\log(\dot{\rho}_\star + \alpha M_\star)$, where $\alpha = {(3.36\pm1.40)\times10^{-11}}\ {\rm yr^{-1}}$. In order to investigate this relation, we construct a phenomenological binary population synthesis model. We find that the high mass XRB and low mass XRB fraction in formed compact object binary systems is ~9% and ~0.04%, respectively. Utilizing the latest XMM-Newton, and Swift X-ray source catalog data sets, additional XRB candidates are also found resulting in 5757 XRBs hosted by 311 galaxies.


Introduction
Recent discoveries of gravitational waves from merging black holes (BHs) and neutron stars (NSs) have finally opened the multi-messenger astronomy era (Abbott et al. 2016;Abbott et al. 2017). However, the nature of those binary systems is not fully understood yet. In order to elucidate the formation mechanism of those binary systems, various binary population synthesis models are proposed in the literature (e.g., Belczynski et al. 2008;Belczynski et al. 2010;Belczynski et al. 2016a;Belczynski et al. 2016b;Mandel & de Mink 2016;Pavlovskii et al. 2017;Marchant et al. 2017;Kruckow et al. 2018;Mapelli & Giacobbo 2018).
Here, X-ray binaries (XRBs) can provide an independent test of those binary population synthesis models. XRBs are close binary systems found in nearby galaxies and radiate most of their emission in the X-ray band (see, e.g., Done et al. 2007). They are powered by mass accretion from a companion star onto a compact object (i.e., BH and NS) either through Roche lobe overflow or stellarwind fed process. After the bright X-ray phase, XRBs evolve into various compact binary systems such as gravitational wave merger sources and millisecond pulsars.
Among XRBs, the formation of luminous XRBs, socalled ultra-luminous X-ray sources (ULXs; see Kaaret et al. 2017 for recent reviews) would be tightly related to the coalescences of currently observed binary BHs (Abbott et al. 2016;Nitz et al. 2019;Abbott et al. 2020). Most of the observed binary BHs are found to be composed of BHs with > ∼ 20M (Abbott et al. 2016;Nitz et al. 2019), more massive than known Galactic stellar BHs ≈ 10M (Casares & Jonker 2014) 1 . Binary population synthesis models predict that the common envelope scenario (Pavlovskii et al. 2017) and the chemically homogeneous evolution scenario (Marchant et al. 2017) could form those binary BHs through ULX phases. It is also known that the local ULX formation rate is consistent with the measured binary BH merger rate (Inoue et al. 2016;Finke & Razzaque 2017). Therefore, statistical information of XRBs, including ULXs, would be keys to exploring the parameter spaces of binary population synthesis models.
Such statistical studies of XRBs are also useful for the studies of low metal galaxies. In low-metallicity starforming galaxies, nebular He II emission lines are frequently observed. However, normal stellar populations can not reproduce the observed emission (e.g., Shirazi & Brinchmann 2012). This is because the nebular He II emission requires the sources emitting ionizing photons above 54 eV. Schaerer et al. (2019) suggested that X-ray emission from XRBs and ULXs would be the dominant ionizing sources for the nebular He II emission. To understand the contribution of XRBs in He II ionization, quantitative and statistical properties of XRBs in a galaxy and the relation to the host properties such as stellar mass and star formation rate (SFR) are required.
catalogs, Chandra provides the finest angular resolution down to ∼ 0.5 arcsec. In this paper, by cross-matching Chandra detected X-ray sources with nearby galaxy catalogs such as the local volume galaxy (LVG) catalog (Karachentsev et al. 2013) and the catalog of the Infrared Astronomical Satellite (IRAS) survey (Neugebauer et al. 1984), we aim to make a new XRB catalog and revisit the statistical properties of XRBs in the extragalactic galaxies.
This paper is organized as follows. The Chandra Xray source catalog is introduced in §. 2. We describe host galaxy catalogs and their properties in §. 3. The catalog matching procedure and our XRB catalog are presented in §. 4. Statistical properties of XRB host galaxies and XRBs are shown in §. 5 and §. 6, respectively. To understand the statistical properties of XRBs, we describe our binary population synthesis model for comparison with our data in §. 7. Discussion including the Swift and XMM-Newton catalogs and conclusions are given in §. 8 and §. 9, respectively. Throughout this paper, we adopt the standard cosmological parameters of (h, Ω M , Ω Λ ) = (0.7, 0.3, 0.7).

Chandra Source Catalog
In this paper, X-ray source information is extracted from the latest Chandra source catalog. The Chandra X-ray observatory, launched in 1999 July 23, observes the X-ray sky in the energy range between 0.1-10 keV over a field of view of ∼ 60-250 arcmin 2 . Chandra achieves a subarcsecond on-axis point spread function (Weisskopf et al. 2000;Weisskopf et al. 2002). These instrumental capabilities allow to detect sources with low confusion and good astrometry. On 2019 October 24, a new version of the Chandra source catalog (CSC2) has been released (Evans et al. 2010;Evans et al. 2020a) 2 , which contains 317,167 unique X-ray sources in the sky including 1,432,324 per observation data sets. 3 Non-uniform completeness of data sets hamper the construction of unbiased source catalogs. To avoid such biases as much as possible, we select sources from CSC2 as follows. First, multiple observation data would mimic the source properties due to the variability of sources. CSC2 provides both stacked-observation and per-observation detections. The stacked observations compile all the available different-exposure data sets toward sources then stack them, while per-observations provide the information of sources per single observation. Since XRBs are variable sources, we use the perobservation data sets having the longest exposure toward 2 http://cxc.harvard.edu/csc2/ 3 Note for the Referee, we moved the XMM and Swift catalog part to the discussion section. each source.
Second, the sensitivity in a field-of-view is not uniform. Generally speaking, sensitivity decreases toward the edge of the field-of-views. The detection efficiency of Chandra is almost uniform down to the flux threshold of F th = 5 × 10 −7 photons cm −2 s −1 in the 0.5-7 keV band inside of off-axis angle of 10 arcmin for a 125 ks exposure (see Fig. 24 of the first CSC catalog paper; Evans et al. 2010). Thus, we select sources having off-axis angle smaller than 10 arcmin. Then, we also adopt flux threshold to select samples. As CSC2 detection threshold is improved by a factor of 2 from the first CSC 4 , we select sources brighter than F ≥ (125 ks/t obs )F th /2. t obs is the exposure time. The energy band and spectral assumption are corrected accordingly.
Third, source confusion would bias results. To avoid such a source crowding effect, we select sources not flagged as confused sources. Since XRBs should be compact by definition, we further select only compact sources from the catalog 5 .
Lastly, the interstellar medium (ISM) absorption will reduce soft X-ray fluxes (e.g., Strom & Strom 1961;Morrison & McCammon 1983;Balucinska-Church & McCammon 1992;Wilms et al. 2000). Detailed spectral analysis is required to decompose the local and Galactic absorption effect. However, such treatments require a certain spectral quality. In this work, for simplicity, we correct the Galactic absorption column density only (Willingale et al. 2013). This treatment would introduce uncertainty in evaluating intrinsic flux at < ∼ 1 keV. Therefore, we restrict our samples to having the detected photons in the hard band (2-7 keV).

Host Galaxies
We cross-match the observed X-ray sources with two nearby galaxy catalogs: the LVG catalog (Karachentsev et al. 2013) 6 and a catalog of galaxies detected in the IRAS survey (Neugebauer et al. 1984). Karachentsev et al. (2013) compiled 1075 galaxies in the local volume as the LVG catalog 7 . The catalog includes galaxies having radial velocities with respect to the centroid of the Local Group V LG < 600 km s −1 or galaxies within the distance of D host < 11 Mpc. The LVG catalog contains various galaxy characteristics such as angular diameters, axial ratio, apparent magnitudes in various bands, morphological types, and distances. The definition of morphological types of galaxies in LVG is according to de Vaucouleurs et al. (1991). The catalog also provides an estimate of galaxy parameters such as absolute B magnitude and stellar mass estimated via K-band luminosity. We further add infrared flux values at 60 µm and 100 µm measured by IRAS utilizing the NASA/IPAC Extragalactic Database (NED) 8 to estimate the SFR based on infrared (IR) measurements as done in Swartz et al. (2011). 131 LVG galaxies have measurements both at 60 µm and 100 µm.

LVG Catalog
Since we would like to select X-ray sources within a galaxy, the angular diameter, the axial ratio, and the position angle (PA) of a galaxy are needed. The LVG catalog provides the major angular diameter a 26 defined by the Holmberg isophote ∼ 26.5 mag arcsec −2 in B-band and the corresponding apparent axial ratio b/a. Some LVG galaxies lack the information on the axial ratio. For those objects, we assume b/a = 1. PAs are not listed in the original LVG catalog. Therefore, we obtain the PA information from the NED system. The PA information listed in NED is based on the measurement by SDSS-DR6, 2MASS, RC3, and ESO-LV catalogs. When multiple information is available, we adopt the latest available values for each object. As a result, we add PA for 413 LVG galaxies.

IRAS Catalog
The LVG catalog is a complete galaxy catalog, but up to 11 Mpc. To expand the horizon of our catalog, we retrieve IRAS galaxies from the NED database, including information of distance, flux measurements, and size. The IRAS sources in the database are based on IRAS Faint Source Catalog, version 2.0 (Moshir et al. 1990) and IRAS Point Source Catalog, version 2.0 (Joint IRAS Science Working Group 1988). We select sources with 60 and 100 µm flux measurements available in either catalog. The distances of the IRAS sources in the database are taken from various literatures. The distance indicator of ∼ 70% of the IRAS sources with XRB detection used for the analysis in this paper is based on Tully-Fisher relation. The uncertainty of the corresponding redshift is less than 1%. We restrict galaxies within z ≤ 0.03 ∼ 130 Mpc; otherwise, Xray sources would be biased to only luminous XRBs for distant galaxies. The total number of selected IRAS galaxies is 11357. We match these IRAS galaxies with the RC3 catalog (de Vaucouleurs et al. 1991), which provides the information of galaxy morphology. The closest objects are selected within the matching radius of 15 arcsec. We note that some of IRAS galaxies are overlapped with the LVG catalog. In that case, we primarily use the information from the LVG catalog. The major angular diameter, axial ratio, and PA of IRAS galaxies are taken from NED, which are based on the SDSS-DR6, 2MASS, RC3, and ESO-LV catalogs. When multiple information is available, we adopt the latest available values. The major diameters for SDSS-DR6, 2MASS, RC3, and ESO-LV catalogs are defined at r = 25, K s = 20, B = 25, and B = 25 mag arcsec −2 , respectively. In literature, RC3 D 25 , defined by the Holmberg isophote ∼ 25.0 mag arcsec −2 in B-band, is frequently used. In this paper, we also use D 25 . However, available size information depends on galaxies as described above. Fig. 1 shows the ratio between a 26 and various other angular estimation. The median ratio between a 26 and each angular diameter measurement for LVG galaxies is 1.24 (SDSS-DR6), 2.31 (2MASS), 1.07 (RC3), and 1.25 (ESO-LV), respectively. We take these ratios as conversion factors to correct the values to D 25 for the LVG and IRAS galaxies.

Properties of External Galaxies
We estimate the SFR and stellar mass of individual galaxies to compare the number of XRBs against SFRs and stellar masses. In this paper, we adopt the Salpeter initial mass function (IMF; Salpeter 1955).

Star formation Rates
The original LVG catalog provides an estimate of SFRρ in a galaxy using Hα and the Galaxy Evolution Explorer (GALEX) FUV measurements.
Following Kennicutt (1998), Hα SFR is estimated aṡ where d host is the distance to the galaxy in Mpc and F c Hα is the extinction corrected integral Hα line flux in erg cm −2 s −1 . The Hα extinction is A(Hα) = 0.538 is the B-band extinction according to Schlegel et al. (1998) and A i B is the intrinsic B-band extinction of the galaxy according to Verheijen (2001).
Another SFR estimate is based on the measurement of FUV flux of a galaxy. Following Lee et al. (2011), Karachentsev et al. (2013) used the relatioṅ where m c FUV is the extinction corrected FUV magnitude as . As we combine IRAS flux measurements, we also estimate the SFR of a galaxy using IR flux measurements following Kennicutt (1998); Swartz et al. (2011) aṡ where S 60 and S 100 is the IRAS 60 and 100 µm flux measurement in Jy. However, we note that this estimate reflects only the reprocessed emission component by dust.
Lastly, in order to avoid the uncertainty of extinction correction, we also make an estimate of the SFR of a galaxy using FUV and IR measurements aṡ whereρ e ,FUV is the FUV SFR estimate but using the FUV flux measurement before correcting internal attenuation.
For IRAS galaxies, we retrieve GALEX FUV flux measurement through NED and corrected for the Galactic extinction following Schlegel et al. (1998) as in the LVG catalog where we utilzied the dustmap software package (M. Green 2018). We then deriveρ ,tot of the IRAS galaxies in the same manner as for the LVG galaxies.
In this paper, we useρ ,tot as the SFR values (ρ ) of both LVG and IRAS galaxies. If FUV flux information is not available, we adoptρ ,IR . We do not use Hα and FUV SFR estimations as Hα flux and intrinsic attenuation factor are not available for IRAS galaxies. Left panel of Fig. 2 shows the comparison ofρ ,tot withρ ,Hα andρ ,FUV in 129 LVG galaxies , which have both information. As it can be seen, SFR indicators of Hα and FUV may underestimate the total SFR in some galaxies. Right panel of Fig. 2 shows the comparison ofρ ,tot withρ ,IR in 129 LVG galaxies and 4501 IRAS galaxies. SFRs of LVG galaxies are more FUV dominated than those of IRAS galaxies. Even though IRAS galaxies are IR selected,ρ ,FUV may still contribute to the total SFR at some level. However, that effect would be

Fig. 2. Left:
Comparison of SFR estimates of LVG galaxies between Hα and FUV+IR methods (circle). Comparison between FUV and FUV+IR methods are shown in diamond. Right: Same as the Left panel, but comparison between IR and FUV+IR methods for LVG galaxies (circle) and IRAS galaxies (triangle).

Stellar Mass
For the estimation of stellar mass M of a galaxy, we follow the mass-to-light (M /L) relation proposed by Bell & de Jong (2001). We adopt their relation for B-band and K-band measurements. For simplicity, we assume that K and K s band fluxes are the same. The relation is given as where a k = −0.926, b k = 0.205, and c k = 0.966. B and K are absolute magnitudes in the B and K bands. We adopt the formation epoch model with bursts and combine the relations for B − V and V − K in Bell & de Jong (2001). Because the stellar mass estimation by Bell & de Jong (2001) is based on a scaled Salpeter IMF, which is 30% smaller than that based on the Salpeter IMF, we multiply 1.43 to the derived stellar mass for our purpose.

Selection of X-ray sources in external galaxies
We select X-ray sources located in the region of the LVG and IRAS galaxies where the galaxy size is defined by D 25 , b/a, and PA. If PA is not available, we define the galaxy region by the circle with a diameter of the square root mean of the major and minor diameters (D 25 and b/a × D 25 ). We remove largely extended nearby galaxies from our search, LMC (a 26 = 646 ), SMC (380 ), M31 (191 ), and M33 (64.7 ). We also remove NGC 5195 from the LVG catalog, a merging galaxy with NGC 5194, also known as M 51, and D 25 of NGC 5194 includes the coun-terpart galaxy NGC 5195.
By using the distance information of host galaxies, we calculate the intrinsic X-ray flux and luminosity in the 2-10 keV band of XRB candidates. Intrinsic fluxes and luminosities are estimated as follows. The CSC2 provides the observed flux in the 2.0-7.0 keV of each object. We convert those fluxes to the intrinsic luminosities by assuming an absorbed power-law spectrum with photon index Γ = 1.8, which is the typical photon index of XRBs (Swartz et al. 2004;Swartz et al. 2011). Even if we set it as Γ = 1.7 or 1.9, the XLF results do not change significantly. Distances to the objects set as the distance to their associated host galaxies. The absorption due to the Galactic cold interstellar medium is modeled using the TBabs code (Wilms et al. 2000), in which cross sections of dust grains and molecules are taken into account. The absorption hydrogen column density (N H ) of each line-of-sight is fixed to the value estimated by Willingale et al. (2013), in which the contribution of not only neutral hydrogen atoms (N HI ) but also molecular hydrogen (N H2 ) are included.
The X-ray sky is known to be dominated by active galactic nuclei (AGNs; see e.g., Ueda et al. 2014). We expect significant contamination of background AGNs in the X-ray source catalogs. To remove known AGNs from the catalog, we cross-match the X-ray catalogs with available AGN catalogs covering a wide area of the sky. The Wide-field Infrared Survey Explorer (WISE; Wright et al. 2010), observing the IR sky, is one of the ideal missions to identify a huge number of AGNs across the full sky. The WISE AGN catalog contains 4,543,530 AGN candidates with 90% reliability (Assef et al. 2018). AGN candidates are selected by the IR color information. Véron-Cetty & Véron (2010) also provides a catalog of AGNs containing 168,941 AGNs in its 13th Edition. We utilize these two catalogs to reject the contamination of known AGNs. Following the same methods for the X-ray catalog comparison above, we collect all the X-ray sources around each AGN catalog within 60 arcsec. Based on the distributions, we define the tolerance radii, where the chance coincidence becomes 5%. We remove those matched Xray sources, which can make about 7% of the number of X-ray sources. In addition, central objects can be contamination of low-luminosity AGNs or unresolved sources. Therefore, we remove sources whose positional error is within galactic center.
The surface number density of the faintest X-ray sources reaches ≈ 50500 deg −2 , dominated by AGNs, based on the Chandra deep field survey with ∼ 7 Ms observations (Luo et al. 2017), which goes down to 4.2 × 10 −18 and 2.0 × 10 −17 erg cm −2 s −1 at 0.5-2 and 2-7 keV, respectively. The median size of the angular diameter of the selected LVG and the IRAS galaxies are 8.1 and 4.2 , respectively. Thus, there may be contamination of background unknown AGNs even after removing known AGNs. To avoid contamination of such unknown AGNs to the XRB catalog, we remove galaxies whose X-ray source number density is below the AGN surface density fluctuation in all the flux regions. We note that the 1-σ flux fluctuation of the cosmic X-ray background radiation in the 2-10 keV band is 6.49 +0.56 −0.61 % (Kushino et al. 2002). We adopt the AGN 2-10 keV source count distribution of Ueda et al. (2014). Distant galaxies and high-latitude stars would also contaminate. At > 2 keV, however, those populations are insignificant in the source count comparing to AGNs (Lehmer et al. 2012).
The sensitivity of current X-ray observatories allows us to see even extragalactic supernova remnants (SNRs). SNRs may start to become a main X-ray source population at < ∼ 10 35.5 erg s −1 at 0.5-2 keV (e.g., Binder et al. 2012). Although SNRs tend to have soft spectra, we restrict our samples whose X-ray luminosity is ≥ 10 36 erg s −1 to pick up XRBs. We note that some SNRs are known to be as bright as 10 37 erg s −1 in X-rays (Ghavamian et al. 2005;Maggi et al. 2016). Since the highest ULX luminosity is ∼ 10 42 erg s −1 (Farrell et al. 2009), we restrict samples whose X-ray luminosities do not exceed 10 42 erg s −1 .
Finally, after applying the criteria above, our catalog includes 4430 X-ray sources associated with 237 galaxies. Among those XRB candidates, 378 are ULX candidates. Hereinafter, we treat those XRB candidates as XRBs. The catalog is available from the following link http://astro-osaka.jp/inoue/page/exrbcatalog/. The description of the catalog is presented in Table 1. As described in § 8.4, we also include Swift and XMM-Newton sources in this catalog. However, in this paper, we rely on Chandra detected sources, otherwise noted. We also provide a list of extended objects, although they are not used in our study. Figure. 3 shows the cumulative source count distribution of XRBs in the selected LVG and IRAS galaxies in the 2-10 keV band. For comparison, the source count distribution of AGNs is also shown (Ueda et al. 2014).

Properties of XRB Host Galaxies
In this section, we discuss the properties of the 237 XRB host galaxies. Fig. 4 shows the normalized distribution of morphological type T of XRB hosting galaxies. We follow the definition of morphological types in Table.   Vaucouleurs et al. (1991). However, for the plotting purpose, we set that T = 12 corresponds to peculiar galaxies. As clearly seen in the Figure, late-type galaxies (T ≥ 0) dominate the XRB hosting galaxies in our samples. ∼ 85% of the XRB hosting galaxies are late-type. For comparison, we also show the normalized distribution of T of the whole LVG and IRAS galaxies. Since the LVG catalog includes a number of local dwarf galaxies, their distribution becomes different from the others.
We show the relation between stellar mass and SFR of XRB hosting galaxies in Fig. 5. Since we restrict host Late-type (z ∼ 0) Late: 111 galaxies Early: 18 galaxies Fig. 5. Comparison between stellar mass and SFR of galaxies in XRB hosting galaxies. Circle and triangle data points correspond to late-type and early-type galaxies, respectively. Dashed line shows the linear regression fit to the late-type galaxies with 1-σ error region (shaded region). Solid line is the main-sequence relation at z ∼ 0.1 (Salim et al. 2007).
galaxies having both stellar mass and SFR, the number of the total remaining galaxies is reduced to 129 in this Figure. As seen in the plot, heavier XRB host galaxies have higher SFR. This is similar to a well-known galaxy main sequence (e.g., Noeske et al. 2007;Salim et al. 2007), which is defined by a relatively steady star-formation rate in disk galaxies. One would like to quantify the relation between SFR and stellar mass to compare XRB hosting galaxies with other studies. The sequence for the late-type galaxies is best described with a linear fit: logρ = (−0.31 ± 0.05) + (0.55 ± 0.05) × (log M − 10), (6) which is shown in the Figure. This relation is similar to that of galaxies in the nearby Universe (see, e.g., Salim et al. 2007). Slight offset comparing to Salim et al. (2007) for z ∼ 0.1, after the IMF correction, may be due to the redshift evolution. When we include early-type galaxies, the correlation changes as logρ = (−0.46 ± 0.06) + (0.46 ± 0.06) × (log M − 10). (7) The slope index is consistent with that for the late-type galaxies within the uncertainty. This may be due to a small number of XRB hosting early-type galaxies.

X-ray Source Distribution
In this section, we discuss the statistical properties of those XRBs.

Radial Distribution
Radial distribution of XRBs gives an independent critical test for our selection as XRBs are off-nucleus objects. The left panel of Fig. 6 shows normalized cumulative radial distributions of XRB candidates in various X-ray luminosity bins of the 2-10 keV X-ray luminosity, L X . About 0.2% of the sources are located within the scale of 1% of the semi-major axis. As no apparent central clustering exists, most of our samples can be regarded as off-nucleus objects. We perform the Kolmogorov-Smirnov (KS) tests to check whether the distributions among luminosities are statistically different. For all the luminosity bins of 36 ≤ log L X < 39, 39 ≤ log L X < 40, and 40 ≤ log L X < 42, we can not reject the null hypothesis. Further detailed investigations on the spatial distribution of XRBs and ULXs would be required for the understanding of this trend.
It has been discussed that HMXBs are the dominant population in late-type galaxies, while LMXBs are in early-type galaxies (see, e.g., Grimm et al. 2003;Gilfanov 2004). The right panel of Fig. 6 also shows the cumulative radial distributions, but divided into two different galaxy types, late (T ≥ 0) and early (T < 0). The obtained KSvalue is 6.3 × 10 −2 with a p-value of 0.006. This comparison implies statistically different spatial distributions between galaxy types, corresponding to the distributions of HMXBs and LMXBs.

Star Formation Rate Normalized X-ray
Luminosity Function Fig. 7 shows the cumulative source count distribution in luminosity of each member galaxy in our samples. We also show the cumulative distribution of all the XRBs in our catalog. The distribution spreads in wide range, and no clear trends seem to exist. Here, the number of XRBs is known to scale with SFR of the host galaxy (e.g., Grimm et al. 2003;Swartz et al. 2011;Mineo et al. 2012). Fig. 8 shows SFR normalized cumulative luminosity count distributions of each member galaxy in our samples where SFR information is available. The overall shape of cumulative luminosity distribution is similar to each other.
The galaxy having the highest SFR normalized source count is NGC 3379, which is a nearby elliptical galaxy (T = −3) at 11 Mpc away with SFR and stellar mass of 8 × 10 −4 M yr −1 and 8.8 × 10 10 M , respectively. Fig. 9 shows the SFR normalized XLF, d 2 N/d log L X dρ , obtained by taking the median of the combining data of full galaxy samples. Source luminosity is in the range of log L X ≥ 36. We divide our samples into the late-type galaxies, T ≥ 0, and early-type galaxies, T < 0. As most of our galaxy samples are late-type galaxies, the distributions of the late-type and the full samples look comparable.
In order to understand statistical properties of SFR normalized XLFs, we fit the binned SFR normalized XLF data. We adopt a model having a broken power-law with a cut-off for the fit. The XLF model is given in the form of When the broken power-law is not able to fit the binned XLF data, we adopt a single power-law form (i.e., removing the term of γ 1 in the XLF).
Since we can not determine the cut-off luminosity L c from the XLF fit, we fix it as log L c = 40.0. Dependence on L c will be discussed in §. 6.4. We perform Markov chain Monte Carlo (MCMC) fitting in order to constrain parameters by using the emcee package (Foreman-Mackey et al. 2013). We assume flat distributions for priors of parameters. The fitting results are shown in Table. 2 together with likelihood values and Fig. 9. For the error region, we use the highest posterior probability density interval containing 68% of the walker samples.
The XLF of late-type galaxies shows flattening in the low luminosity regime. This flattening might be due to the incompleteness of observatories. In this paper, even though we select sources at off-axis angles smaller than 10 arcmin and sources whose flux is above the flux threshold having uniform completeness , we do not correct the XLF for possible survey incompleteness (i.e., flux dependence of the survey area) of each Chandra observation. Late: 191 galaxies Early: 30 galaxies All: 221 galaxies Fig. 9. Star formation rate normalized X-ray luminosity function of XRBs in nearby galaxies. The circle, triangle, and square points correspond to that of the late-type, early-type, and all galaxies, respectively. Dotted, dashed, and solid line shows the fitted SFR normalized XLF model for the late-type, early-type, and all galaxies, respectively. The shaded region corresponds to its error region. The models for late-type and all galaxies overlay each other.
Later, we discuss a possible cause of this break by comparing it with a phenomenological binary evolution model. We note that, since the total X-ray luminosity is dominated by the high-luminosity end considering the obtained γ 2 , since γ 2 ≤ 1. Therefore, the faint-end slope uncertainty does not affect the estimation of the total luminosity.

Stellar Mass Normalized X-ray Luminosity Function
The number of XRBs is also known to scale with their stellar mass of the host galaxy, especially in early-type galaxies (e.g., Gilfanov 2004;Zhang et al. 2012;Lehmer et al. 2014;Peacock & Zepf 2016). Fig. 10 shows stellar-mass normalized cumulative luminosity count distributions of each member galaxy in our samples. The galaxy sample used for stellar-mass normalized XLFs is partly different from that for SFR normalized XLFs because the samples should have stellar-mass information.
The overall shape of cumulative luminosity distribution is similar to each other. However, there are several outliers. The galaxy having the highest count is the Garland galaxy. Garland is a member galaxy of the M 81 group (Karachentsev et al. 2002). It is an irregular galaxy having a stellar mass of 3.7 × 10 6 M . Their origin and property are still under debate. The lowest count is NGC 0315, an elliptical galaxy, which has a stellar mass of 9.1 × 10 11 M . NGC 0315 is known as a low-ionization nuclear emission-line region AGN having an extended ra-   Same as Fig. 9, Fig. 11 shows the stellar mass normalized XLFs (d 2 N/d log L X dM ). In this plot, the maximum luminosity bin is lower than that in Fig. 9. This is because we utilize the stellar mass available galaxies only. We fit XLFs using the same function as Eq. 8 in the same manner, but replacingρ with M . The fitting results are shown in Table. 2 and Fig. 11.

Fundamental Plane of XRB hosting Galaxies
SFR and stellar mass normalized XLFs are motivated by the fact that SFR traces the young stellar population, i.e., HMXBs, and stellar mass traces the old stellar population, i.e., LMXBs. However, as we see in our Galaxy, XRB populations in a galaxy are expected to be a mixture of HMXBs and LMXBs. Thus, it is naturally expected that bothρ and M are important to characterize the XRB population in a galaxy. In this section, we investigate the Late: 119 galaxies Early: 23 galaxies All: 142 galaxies Fig. 11. Same as Fig. 9, but for stellar mass normalized XLFs. relation among the integrated luminosity of XRBs L X,tot , ρ , and M of each galaxy. Lehmer et al. (2010) pointed out that 2-10 keV X-ray emission of nearby galaxies correlates with bothρ and M using 17 luminous infrared galaxies (LIRGs). Lehmer et al. (2016) further investigated the redshift evolution of that correlation. In this paper, we compile 237 XRB host galaxies. Among them, 129 galaxies have bothρ and M information. Since the previous studies integrated the whole X-ray emission of galaxies, there was about 10% level of flux contamination by galactic hot gas emission (∼ 0.5-1 keV). By extracting individual XRBs, we can directly estimate the total X-ray emission from XRBs. Fig. 12 shows the three-dimensional plot of the relation among L X,tot ,ρ , and M . L X,tot is the summation of X-ray luminosities of XRBs with log L X ≥ 36 in a galaxy. These two panels show the same data but from different viewing angles. As it is clearly seen from the plots, higher Late: 111 galaxies Early: 18 galaxies All: 129 galaxies Fig. 13. Same as Fig. 9, but normalized by µ. ρ and higher M tend to have higher L X,tot . This is because there are three relations: the main-sequence relation betweenρ and M (e.g., Noeske et al. 2007;Salim et al. 2007), the correlation betweenρ and L X,tot (e.g., Grimm et al. 2003;Mineo et al. 2012), and the correlation between M and L X,tot (e.g., Gilfanov 2004). The Spearman correlation coefficients of those relations are ρ S = 0.65, 0.58, and 0.53, respectively, indicating the existence of positive correlations among those parameters.
Considering the relation among these three parameters, here, we propose to introduce a new parameter defined as: The first term of the right-handed side is the star formation rate which represents the young stellar population (i.e., HMXBs), while the second term is linearly related to the stellar mass which represents the old stellar population (i.e., LMXBs). α controls the ratio between these two populations in the unit of yr −1 . The physical meaning of α is discussed in §. 7. Then, from Fig. 12, we have log L X,tot = (38.80 +0.09 −0.12 ) + log µ(ρ , M ), with α = (3.36 ± 1.40) × 10 −11 yr −1 .
Utilizing µ, we can derive the µ normalized XLFs (d 2 N/d log L X dµ). The total number of galaxies is 129 which have bothρ and M . We fit XLFs using the same function as Eq. 8 in the same manner, but replacingρ with µ. The fitting results are shown in Table. 2 and Fig. 13. We fix log L c = 40.0. By introducing µ, the XLFs of the late, early, and all type galaxies show similar shapes indicating a universal trend in the local Universe. The similarity of XLF shapes indicates a single component dominates the XLFs, even though it is a mixture of HMXBs and LMXBs. This suggests that HMXBs and LMXBs may have a similar XLF shape.
Fig. 14 shows the relation between µ and L X,tot , corresponding to the fundamental plane relation. These two parameters show a positive correlation. The Spearman rank correlation coefficient is 0.68, 0.78, and 0.67 for late, early, and all type galaxies, respectively, with p-value of 10 −5 . A linear regression line corresponding to Eq. 10 is also shown as the solid line in the Fig. 14. We note that we can fit the data with higher-order equations. However, as discussed in §.6.2 and §.6.3, the definition of Eqs. 9 and 10 is a physically reasonable choice.
We can reproduce this relation from the µ normalized XLF. We show the integrated luminosity for various µ in Fig. 14 by the dashed line. The minimum and maximum luminosity is set to be 10 36 and 10 45 erg s −1 , respectively. Because of the existence of L c , the choice of the maximum luminosity does not affect the results as far as we set it higher than L c . Different dashed lines correspond to the different L c values. As presented in Fig. 14, log L c = 40.0 almost matches with the regression line. Therefore, we set log L c = 40.0 in this paper.

Phenomenological XRB Population Synthesis Model
The µ normalized XLF reveals the bright-end XLF slope is γ 2 ∼ 0.6 for any galaxies. It is interesting to see whether this γ 2 can be reproduced by considering a binary population synthesis model. Today, various binary population synthesis models are available in the literature (e.g., Belczynski et al. 2008;Belczynski et al. 2010;Belczynski et al. 2016a;Belczynski et al. 2016b;Mandel & de Mink 2016;Pavlovskii et al. 2017;Marchant et al. 2017;Kruckow et al. 2018;Mapelli & Giacobbo 2018). However, each model contains various detailed physical processes, and it is not easy to compare with our current data directly. Therefore, we construct a phenomenological XRB population synthesis model utilizing the latest stellar evolution model. First, we consider the HMXB population. Once we define the star formation rateρ and the IMF dN/dM ZAMS shape at zero-age main sequence (ZAMS) phase in a galaxy, the stellar mass distribution formed in a unit time can be evaluated. M ZAMS is the mass of ZAMS stars. Suppose we have a power-law IMF such as the Salpeter IMF, dN/dM ZAMS ∝ M −αIMF ZAMS , the stellar mass function in a unit time is given as where M ZAMS,max and M ZAMS,min are the maximum and minimum mass of ZAMS stars. We utilize a model by Spera et al. (2015) for the mass spectrum of compact remnants (neutron stars and black holes) after the stellar evolution, including the metallicity dependence (See Fig. 6 and Appendix C in Spera et al. 2015). Here, the metal environment affects the resulting compact object mass distribution (e.g., Spera et al. 2015). It is also known that metallicity of galaxies depends on M andρ (Mannucci et al. 2010;Yabe et al. 2012;Andrews & Martini 2013). In this paper, since we consider the local Universe only, we fix the metallicity to the solar value.
Then, we have the compact object formation rate as where M co is the mass of the compact remnant. Following Spera et al. (2015), we treat a compact remnant as a white dwarf when its final core mass is less than the Chandrasekhar mass (1.4M ). The remnants with masses of 1.4M ≤ M co < 3.0M are treated as neutron stars, while those with M co ≥ 3.0M as black holes. Among produced compact remnants, we require binary systems in order to have accretion disk activity. We describe the binary fraction as f b . Also, we require the fraction of X-ray emitting binaries among those binaries f HMXB (i.e., in the mass transfer phase with an adequate mass accretion rate). Then, the mass function of compact objects having the HMXB activity can be approximated as where t HMXB is the duration of HMXB activity. Here, in our Galaxy, a large fraction of massive stars ( > ∼ 70%, Sana et al. 2012) are members of binary systems since their birth. For simplicity, we set the binary fraction as f b ∼ 0.7. And, although t HMXB is not well constrained, it is expected to be about 0.1 Myr (Mineo et al. 2012).
A study of Galactic accreting stellar-mass black holes indicates that the logarithmic Eddington ratio distribution is given by a normal distribution (Reynolds & Miller 2013;Finke & Razzaque 2017), where Edd is the Eddington ratio. µ is ∼ −2 and σ ∼ 1 for the Galactic XRBs including both HMXBs and LMXBs at all epochs (Reynolds & Miller 2013).
The resulting XLF of HMXBs at an X-ray luminosity of L X can be evaluated as where we set L X = Edd L Edd (M co ). L Edd is the Eddington luminosity ≈ 1.26 × 10 38 (M co /M ) erg s −1 . Next, we consider the LMXB population. We can adopt the same method as in the HMXB population, but LMXBs are expected to follow the stellar mass. Therefore, we rewrite the equations 12 and 14 which are related to star formation as where f LMXB is the LMXB fraction among binaries having compact objects and t LMXB is the duration of LMXB activity, which is expected to be typically 10 Myr. t Gal is the age of the galaxy and we assume a constant star formation activity. By adopting Eq. 17, we can evaluate dN LMXB /dL X as Here, both integral terms in the HMXB and LMXB XLFs are the same, which are described as φ X (L X ) hereinafter. Then, the whole XRB XLF becomes As dN HMXB /dL X ∝ρ and dN LMXB /dL X ∝ M and by adopting Eqs. 17 and20, By comparing to the µ normalized XLF (See Eq. 9), we should have The normalization is determined by f b f HMXB t HMXB . Therefore, α represents the ratio between HMXB and LMXB in a galaxy lifetime. Fig. 15 shows the µ normalized XRB XLF based on our phenomenological binary population synthesis (PBPS) model described above with various σ . The binned XRB XLF is also shown. The model with σ = 1.0 well reproduces the data which is in the consistent range with the Galactic XRB distribution (Reynolds & Miller 2013). For the other parameters except for normalization, we do not change. This result tells that the XLF slopes are determined by the IMF and the Eddington ratio distribution.
In order to match the data, we require the fraction of HMXB in compact binary systems formed in a certain time as This manifests that t HMXB can not be 0.01 Myr otherwise f HMXB > 1. As described above, t HMXB is typically about 0.1 Myr (Mineo et al. 2012). From binary evolution calculations, t HMXB ∼ 0.01 Myr for supergiant systems and t HMXB ∼ 0.1 Myr for Be/X binaries. About 60% of known high-mass X-ray binaries are Be/X systems, while ∼ 30% are supergiant systems (e.g., Liu et al. 2006).
From Eq. 23, the fraction of LMXB in compact binary systems formed in a certain time is Thus, LMXBs are rare objects per star formation activity, but we can observe them a lot because of its long lifetime comparing to HMXBs. The reason for the difference between f HMXB and f LMXB may be due to the required mass ratio. In our own Galaxy, massive stars are known to be members of binary systems whose mass ratio distribution is flat (Kobulnicky & Fryer 2007;Sana et al. 2012;Kobulnicky et al. 2014). Here, LMXBs requires a high mass ratio because the donor star is a low-mass star (e.g., Kalogera & Webbink 1998;Tauris & van den Heuvel 2006). Therefore, it is natural to expect a low f LMXB value. This equation also indicates that younger galaxies will have a lesser LMXB population.

Comparison with Previous Studies of XRB XLFs
Theρ and M normalized XRB XLFs have been studies in literature. In this section, we compare our results with the previous studies on those XLFs.
The left panel of Fig. 16 shows theρ normalized XRB XLFs of our study, Grimm et al. (2003), Swartz et al. (2011), andMineo et al. (2012). We note that Swartz et al. (2011) used ULX data only. As seen from the figure, the slope of the XLFs are mostly consistent in 38 < ∼ log L X < ∼ 40, while the normalization of Mineo et al. (2012) is slightly smaller than the other works. One remarkable difference is seen in the low-luminosity regime. Even though the model of Swartz et al. (2011) is just an extrapolation, the Grimm et al. (2003) and Mineo et al. (2012) covered the luminosity range down to log L X ∼ 35. The reason for the lower estimates of our study is likely due to the incompleteness of the observations. In Grimm et al. (2003) and Mineo et al. (2012), the incompleteness of each observation is corrected, but we did not correct it as explained in § 6.2. However, we note that, as seen in Fig. 15, the PBPS model shows a drop in XLF at lower luminosity due to the Eddington ratio and mass distribution.
The right panel of Fig. 16 shows the M normalized XRB XLF of our study, Gilfanov (2004), and Zhang et al. (2012). As seen from the figure, the shape and normalization of the XLFs are mostly consistent among these three works, but our studies show the existence of luminous XRBs in galaxies. This may be due to our large galaxy sample, roughly 10 times more than those in the previous works. As the main parent galaxy sample class was earlytype galaxy in Gilfanov (2004) and Zhang et al. (2012), their XLFs indicate fewer ULXs in early-type galaxies. However, it is known that elliptical galaxies host numerous ULXs (Swartz et al. 2004), which is qualitatively consistent with our result. Although the observational incompleteness is corrected for Gilfanov (2004) Grimm et al. (2003), Swartz et al. (2011), andMineo et al. (2012), respectively. The number of sampled galaxies is also shown in the panel. Right: The same as in the Left panel, but for M normalized XLFs. Solid, dashed, and dotted line corresponds to our study, Gilfanov (2004), and Zhang et al. (2012Zhang et al. ( ), respectively. et al. (2012, the faint-end slope is similar among these three works. This is because the correction effect of the observational incompleteness is not significant in these studies (see Table. 4 in Zhang et al. 2012).

Galaxy Sample Selection
Some of our galaxy samples do not have information of either ofρ or M . Here, in our sample, the SFR completeness is about 93%, which indicates biases on SFR information is not significant. However, the µ completeness is about 54%, indicating that our µ normalized XLF results can be affected by the selection effect.
The observational status of each galaxy differs one by one. When SFR or stellar mass information is unavailable, it is not easy to get further information. Some of them may not even have measurements. Therefore, it is not straightforward to understand the selection effect. However, we can at least test whether the µ selected samples differ from the full sample using theρ normalized XLFs. Fig. 17 shows the comparison of theρ normalized XLFs of the full galaxy sample and the µ selected samples. As both data overlap, we do not see any clear deviation between these two galaxy samples. Therefore, the selection effect on the µ selected samples is insignificant.

Effect of Background AGNs
AGNs, dominating the X-ray sky, would contaminate our XRB samples. To remove those background AGNs, we subtract known AGNs using available AGN catalogs galaxies whose X-ray source density is higher than AGN surface density. There is a less biased method. In this method, we estimate the number of possible background AGNs at a given flux within an area of a galaxy using AGN flux distribution, then subtract them from observed XRB XLFs assuming all the AGNs at a distance of that galaxy. However, this method requires a correction of detailed position dependence of the detection efficiency. In this paper, which we do not take into account this, although we adopt the flux threshold where almost uniform detection efficiency achieved (Evans et al. 2010).
To investigate the effect of possible background AGN contamination, we recalculate the µ normalized XLFs by subtracting background AGNs using AGN XLF (Ueda et al. 2014)  account the detailed position dependence of the detection efficiency in each observation. Fig. 18 shows the comparison of recalculated XLF with our fiducial XLF. As shown in the figure, although we see a clear deficit at the lowest luminosity band, the two XLFs do not show clear difference at higher luminosity ranges. This may indicate that our AGN subtraction method would not significantly affect our main results. We note that at lower-luminosity ranges, XLF density becomes lower for the AGN XLF method. This is because we do not take into account the detailed detection efficiency which generally drops at fainter fluxes.

Swift and XMM-Newton Source Catalogs
The Neil Gehrels Swift observatory and the XMM-Newton observatory also observe the X-ray. Both observatories recently released their source catalogs the 2SXPS catalog (Evans et al. 2020b) 9 and the 4XMM-DR9 catalog (Webb et al. 2020) 10 , respectively. The 2SXPS and 4XMM-DR9 catalog contains 206,335 and 550,124 unique X-ray sources in the sky, respectively. Although their angular resolution is not as good as that of Chandra, these samples help us to expand our XRB samples. We add these catalogs following the same selection criteria as done for CSC2. However, the 2SXPS provides stacked data only.
There are significant amount of overlapped sources among the catalogs. To avoid such overlaps, we match the catalogs with each other and searched neighboring sources within a tolerance radius r tol (see also Itoh et al. 9 https://www.swift.ac.uk/2SXPS/ 10 http://xmmssc.irap.omp.eu/Catalogue/4XMM-DR9/4XMM_DR9.html 2020). For an X-ray source catalog, we collect all the X-ray sources from the other catalogs within 60 arcsec, which is large enough comparing to the positional uncertainty of the three X-ray satellites. To set the background contamination less than 5%, we set the r tol = 3.3 , 3.5 , and 5.8 for the matching of CSC2 -4XMM-DR9, CSC2 -2SXPS, and 4XMM-DR9 -2SXPS, respectively. We prioritized Xray sources in the order of CSC2, 4XMM-DR9, and 2SXPS, considering their angular resolutions. According to the prioritization, we remove objects matched in other catalogs.
By combining three X-ray catalog and matching with galaxy catalogs, then it includes 5757 X-ray sources associated with 311 galaxies. Among those XRBs, 635 are ULX candidates. We reconstruct the µ normalized XLF using this larger sample data set. Fig 18 also shows the µ normalized XLFs based on the Chandra-Swift-XMM combined catalog and the Chandra only catalog. At log L X < 40.5, both distributions show the similar structure, while the combined catalog extend the high-luminosity end toward log L X = 41.5. This is because of four 2SXPS sources having log L X > 40.5. As described above, the 2SXPS catalog provide the stacked information, which may result in higher fluxes than typical fluxes. Further careful treatment of 2SXPS and XMM-DR9 catalogs including e-ROSITA will help to understand the detailed statistical properties of XRBs.

Conclusion
In this paper, we construct a new catalog of XRBs in external galaxies utilizing the latest Chandra (Evans et al. 2010;Evans et al. 2020a) source catalog. We match these X-ray source catalogs with two nearby galaxy catalogs: the LVG catalog and the IRAS catalog. Our catalog contains 4430 XRBs, including 378 ULXs, associated with 237 galaxies. ∼ 84% of those XRB hosting galaxies are late-type. The XRB host galaxies reproduce the so-called galaxy main-sequence relation seen in the local galaxies (e.g., Noeske et al. 2007;Salim et al. 2007.
> ∼ 99% of our XRB samples are off-nucleus objects. There is no apparent difference in the spatial distributions of XRBs and ULXs, while spatial distributions of XRBs between late and early-type galaxies show difference.
As demonstrated in previous studies, we can reproduce the SFR (ρ ) and stellar mass (M ) normalized XLF relations and both of them are well described by a brokenpower-law XLF model. With our XRB samples, we further investigate a relation among the integrated luminosity of XRBs L X,tot ,ρ , and M in each galaxy. We find that there is a fundamental plane in those three parameters as log L X,tot = (38.80 +0.09 −0.12 ) + log µ, where µ ≡ρ + αM . α is determined as (3.36 ± 1.40) × 10 −11 yr −1 . The µ normalized XLFs of late-type and early-type galaxies are almost equivalent, suggesting that HMXBs and LMXBs may have a similar XLF shape. Furthermore, the µ normalized XLF can reproduce the correlation between µ and L X,tot of our XRB host galaxy samples.
In order to investigate the fundamental plane, we construct a phenomenological binary population synthesis model. We find that the high mass XRB (HMXB) and low mass XRB (LMXB) fraction in formed compact object binary systems is ∼ 9% and 0.04%, respectively. Thus, LMXBs are rare objects per star formation activity, but we can observe many of them because of their long lifetime comparing to HMXBs.
To increase the number of potential XRBs, we further add the 2SXPS catalog (Evans et al. 2020b) and the 4XMM-DR9 catalog (Webb et al. 2020) based on the Swift and XMM-Newton observations respectively. Our XRB catalog also provides these information.
By adding these two X-ray catalogs, the sample size increases to 5757 XRBs associated within 311 galaxies. However, we note that Chandra has better angular resolution than the other two and that the 2SXPS catalog provides the stacked flux information only. Therefore, for the analysis of XLFs, it would be better to use the Chandra sources only.