MODELING X-RAY BINARY EVOLUTION IN NORMAL GALAXIES: INSIGHTS FROM SINGS

The main tool we used to perform our PS simulations is StarTrack, a state-of-the-art PS code that has been tested and calibrated using detailed mass-transfer star calculations and observations of binary populations. StarTrack has been applied to numerous interpretation studies of X-ray and radio pulsar binary populations, as well as γ-ray bursts and binary BHs in the context of gravitational-wave sources (see a detailed description in Belczynski et al. 2008, and references therein). In summary, the code incorporates all the important physical processes of binary evolution.

• 1.  
  The evolution of single stars and non-interacting binary components from zero-age main sequence to remnant formation is followed with the use of high-quality analytic formulae (Hurley et al. 2000). Various wind mass loss rates dependent on stellar evolutionary stage are incorporated and their effect on stellar evolution is taken into account.
• 2.  
  Changes in all the orbital properties are tracked. Through numerical integration of four differential equations, the evolution of orbital separation, eccentricity, and component spins is tracked; these depend on tidal interactions as well as angular momentum losses associated with magnetic braking, gravitational radiation, and stellar wind mass losses.
• 3.  
  All types of mass-transfer phases are calculated: stable, driven by nuclear evolution or angular momentum loss, and thermally or dynamically unstable.
• 4.  
  Supernova explosions are treated accounting for mass loss and asymmetries with natal kicks to NSs and BHs at birth; systemic velocities for all binaries are calculated.
• 5.  
  The calculation of mass-transfer rates in binaries with accreting NSs and BHs (driven by stellar winds or RLOF) has been calibrated against detailed mass-transfer sequences, and X-ray luminosities are calculated by incorporating appropriate bandpass corrections and spectral models.

This paper uses results that are based on a recent major revision of the StarTrack code that includes updated stellar wind prescriptions and their re-calibrated dependence on metallicity (Belczynski et al. 2010). Two newer updates have not been taken into account as our simulations were complete before these updates had been implemented. For reference, these are (1) a revised NS and BH mass spectrum, leading to fully consistent supernova simulations (Belczynski et al. 2012; Fryer et al. 2012), and (2) a more physical treatment of donor stars in common envelopes (CEs) via actual λ values (Dominik et al. 2012), where λ is a measure of the donor's central concentration and the envelope binding energy.

Table 2 lists the full set of input parameters used in our PS modeling. We constructed a large PS model grid using a range of values for these parameters as indicated in the table. For a detailed discussion of how each parameter affects the overall XRB population, we refer the reader to F13 (Section 5.1 and Figure 6) and T13 (Section 4.3 and Figure 7).

The input parameters fall into two categories: first, parameters that mostly characterize initial properties of the population, such as initial mass function (IMF), initial binary mass ratio (q), and distribution of initial orbital separations. For these parameters we usually had information from observational surveys of binary stars, constraining the range of values used. The second category comprises parameters that are associated with physical processes that are poorly understood, such as the efficiency, α_CE, of converting orbital into thermal energy that will be used to expel the donor's envelope during a CE phase. In this work, the related parameter varied is λα_CE, which is the product of this efficiency and the donor central concentration, λ.

The complete grid of 288 PS models is the same as that used by F13 and T13. To construct this grid, we vary all parameters known from earlier studies to affect XRB evolution and formation of compact objects (Belczynski et al. 2007, 2010; Fragos et al. 2008, 2009, 2010; Linden et al. 2009). Specifically, we use the following.

• 1.  
  Four λα_CE values (0.1, 0.2, 0.3, 0.5).
• 2.  
  Three stellar wind strengths, η_wind (0.25, 1.0, 2.0). This parameter is used to multiply the stellar wind prescription of Belczynski et al. (2010).
• 3.  
  The distribution of natal kicks for BHs formed through direct collapse (no kicks or 10% of the Hobbs et al. 2005 distribution for NSs).
• 4.  
  A CE-HG flag for systems with a donor in the Hertzsprung gap, either allowing all possible CE events or always imposing a merger (Belczynski et al. 2007).
• 5.  
  Three distributions of binary initial mass ratios, q ≡ M_secondary/M_primary. This distribution specifies the mass of the secondary, whereas the mass of the primary is governed by the IMF. For our full model grid, we use a uniform (flat) distribution, q = 0 → 1, a twins distribution, q = 0.9 → 1, and a mixed distribution, 50% uniform and 50% twins. However, in this paper we do not use the 96 twin q models 97–192, as F13 have already clearly shown that these are very inadequate in reproducing observed XRB populations, because they prevent the production of LMXBs altogether. This limits the models used in this paper to those in the ranges of 1–96 and 193–288, i.e., 192 models in total.

In addition, for each of the 288 models in our full grid, we also use nine metallicities, keeping all other parameters same. In this paper, we only consider solar metallicity models as the best estimates of Moustakas et al. (2010), using two different methods, appear to straddle solar metallicity for all SINGS galaxies in our study. The best-fit SEDs that we convolve with StarTrack models to construct tXLFs for individual galaxies also assume solar metallicities (Noll et al. 2009).

Each model follows the evolution of 5.12 × 10⁶ stars over 14 Gyr. Since we are only using models 1–96, corresponding to a uniform q distribution, and models 193–288, corresponding to a 50%–50% mixed q distribution, we indicate these model ranges in relevant figures for clarity. For reference, the full list of 288 models and associated parameters can be found in Table 4 of F13.

Since our goal is to construct XRB XLFs, we need to identify all binaries in our simulations that become XRBs and register their X-ray luminosities as a function of time. XRBs are mass-transferring binary stellar systems with a compact object accretor, either a BH or an NS. Systems with donors less massive than 3 M_☉ are labeled LMXBs, and vice versa for HMXBs. LMXBs are always RLOF systems, while HMXBs are usually wind-fed, but can also exhibit RLOF behavior. According to whether they undergo thermal disk instability or not, RLOF systems can be either transient or persistent, while wind-fed systems are always persistent.¹³

We follow the methodology described in F08 and Fragos et al. (2009) to identify all model XRB sources and keep track of properties essential for estimating their X-ray luminosity, L_X, as a function of time. These properties include the mass-transfer rate, $\dot{M}$, as well as the mass and radius of the accretor, M_a and R_a, respectively. We also identify BH or NS accretors, transient or persistent sources, evolutionary stages of donors, and donor masses.

3.2.1. RLOF Systems

For persistent systems, L_X is estimated as

$$\begin{equation} L_X = {\rm min} \left(f L_{\rm Edd}, \eta _{\rm bol} \epsilon \frac{G M_a \dot{M} }{ R_a } \right), \end{equation} \tag{ 1 }$$

with f equal to unity. The value of R_a is 10 km for an NS and 3 Schwarzschild radii for a BH, gives a conversion efficiency of gravitational binding energy to radiation associated with accretion onto an NS (surface accretion, = 1.0) or BH (disk accretion, = 0.5), and η_bol is a factor that converts the bolometric luminosity to L_{X, 0.3–10.0 keV}, the X-ray luminosity in the full Chandra energy band, consistent with our observations. We use results from the literature to tabulate the best available estimates of this factor and its 1σ uncertainty for NS and BH accretors in "low–hard" and "high–soft" state systems (see F13 for a detailed discussion and references). Since we are interested in combining our PS models with galaxy SFHs, we define look-back time windows, δt_lb, in which we evaluate the total stellar mass produced in the galaxies (see Section 3.3). In the StarTrack simulations, each model source is associated with its own time-step window, δt_step, which we use to calculate the fractional time, δN, that a source is on in a given galaxy look-back time interval, i.e.,

$$\begin{equation} \delta N \equiv \delta t_{\rm step}/\delta t_{\rm lb}. \end{equation} \tag{ 2 }$$

For transient systems we follow the prescription of F08, calculating the outburst X-ray luminosity depending on the accretor type, either BH or NS, via

$$\begin{equation} L_{\rm X}=\eta _{\rm bol}\epsilon \times \left\lbrace \begin{array}{@{} ll} \! \min {\left(f_{\rm BH} \times L_{\rm Edd}, f_{\rm BH} \times L_{\rm Edd} \left[\frac{P}{10\, \rm h{\rm r}}\right]\right)}, &\rm \! {for\ BH} \\ \ms \ms \! \min {\left(f_{\rm NS} \times L_{\rm Edd}, \frac{GM_a\dot{M}_{\rm crit}^{2}}{R_a\dot{M}_d} \right)}, &\rm \! {for\ NS}. \end{array}\right. \end{equation} \tag{ 3 }$$

Here f_BH and f_NS are factors setting an upper limit for the maximum Eddington luminosity for BHs and NSs, equal to 2 and 1, respectively. $\dot{M}_{\rm crit}$, $\dot{M}_d$, and P are the critical mass-transfer rate for thermal disk instability, the rate at which the donor star is losing mass, and the orbital period, respectively. In addition, transient systems are mostly in a quiescent state and are too faint to be detectable, except when they go into outbursts. The fraction of time they are in outburst defines their duty cycle, DC. Following Fragos et al. (2008), we estimate DC as (Dobrotka et al. 2006)

$$\begin{equation} {\rm DC} = \left(\frac{\dot{M_d} }{ \dot{M}_{\rm crit} } \right) ^ 2 \end{equation} \tag{ 4 }$$

for NSs and 5% for BHs (Tanaka & Shibazaki 1996). For transient systems it then follows that the fractional time a source is on in a galaxy look-back time interval becomes δN ≡ δt_step/δt_lb × DC.

3.2.2. Wind-fed Systems

The StarTrack simulations described thus far link the stellar mass described by a set of physical parameters with an XRB population, but do not contain information to link the XRB population to host galaxy properties. Thus, the stellar mass produced by a given model is in a sense arbitrary. To construct tXLFs corresponding to real galaxies, we need to modify this arbitrary mass by making use of the galaxy's SFH. This in turn modifies the numbers of XRB sources produced by the model. For the galaxies in our sample, SFH estimates exist based on SED fitting by Noll et al. (2009; see Table 3). We use these results to calculate the total stellar mass produced for each SINGS galaxy in our sample in each look-back time window. Noll et al. (2009) used the SED fitting code CIGALE and assumed two exponentially varying SFRs for a young and an old population. Using their symbols (Table 3), the total SFR at a look-back time t_lb is given by

$$\begin{eqnarray} {\rm SFR}_{t_{\rm lb }} &=& \frac{ f_{\rm burst} M_{\rm gal} }{ \tau _{\rm ySP} (e^{ t_{\rm ySP} / \tau _{\rm ySP} } -1)} e^{ t_{\rm lb}/\tau _{\rm ySP} } \nonumber \\ &&+ \frac{(1-f_{\rm burst}) M_{\rm gal} }{ \tau _{\rm oSP} (e^{ t_{\rm oSP} / \tau _{\rm oSP} } -1)} e^{ t_{\rm lb}/\tau _{\rm oSP} }. \end{eqnarray} \tag{ 5 }$$

As in Table 3, M_gal is the total mass of stars and gas that originates from stellar mass loss in the galaxy, f_burst is the mass fraction of the young population, τ_ySP and τ_oSP are the e-folding timescales for the young and old SPs, respectively, and t_ySP and t_oSP are the ages of the young and old SPs, respectively. All of these parameters are SED fit results allowing us to calculate SFR$_{\bm {t_{\rm lb}}}$ via Equation (5), and then stellar masses, used to scale XRB numbers for each galaxy-model pair.

As explained in F13, the bolometric corrections used to convert StarTrack-derived L_X values for model sources to the full Chandra band are empirical and introduce an uncertainty to the estimated X-ray luminosity value. Thus, each L_{X, 0.3–10.0 keV} value can be thought of as the mean of a Gaussian distribution with a standard deviation originating in the uncertainty introduced by the bolometric correction.

Further, as explained earlier, for each model source the fractional time a source is on in a galaxy look-back time window is given by δN. This number can also be considered as the mean of a Poisson distribution that represents the expected number of times that this source will appear in the XLF. For each source, we use this information to draw a random Poisson deviate, giving the number of times that this source will be on. For each case that the source is on, we draw a random Gaussian deviate from the L_X distribution giving an L_X value for that appearance. When this procedure is complete for all sources, we obtain a set of L_X values, which can be converted to an XLF. However, this is a single realization. To obtain a reliable estimate for the mean XLF and its uncertainty, we carry out 500 Monte Carlo realizations of this process for the total XLF, and 100 for subpopulation XLFs (donor/accretor type, LMXB/HMXB). This procedure gives a reliable estimate for the mean number of sources in each L_X bin, as well as ±1σ uncertainties.

Overall, we thus obtain 12 × 192 = 2304 tXLFs (in both cumulative and differential form) combining SFHs for 12 SINGS galaxies and 192 StarTrack models. The tXLFs for the best models (Section 4) are shown as red dotted curves in Figure 3.

Our purpose is to compare our theoretically derived tXLFs for XRB populations to oXLFs of point sources in SINGS galaxies. Whereas all theoretically obtained point sources in the tXLFs are known to be XRBs by construction, this is not necessarily the case with all observed point sources in SINGS galaxies. It is possible that some of the latter are actually background active galactic nuclei (AGNs) rather than galactic XRBs. We correct our pure XRB-based tXLFs by adding a component that takes into account the expected number of background X-ray point sources in the area surveyed for each galaxy.

We correct our differential and cumulative tXLFs as follows. We use the log N–log S results of Kim et al. (2007; Γ = 1.7, broadband B, their Table 3) to estimate the number of background sources expected in each luminosity bin of our tXLFs, taking into account the area observed for each galaxy (either the Chandra S3/I0-I4 detector area or the D₂₅ region, whichever is smaller). We then add these numbers to our "pure XRB" source numbers in each luminosity bin, thus obtaining background corrected tXLFs. In Figure 3, these corrected, final tXLFs are shown as solid red curves for highest-likelihood models (Section 4). A comparison with observational completeness corrections (see the blue curves in Figure 3), which only affect the faint end of the XLF, shows that the background correction affects the XLF over the full range of X-ray luminosities.

For each galaxy, there is a unique, completeness-corrected oXLF coming from our Chandra observations (continuous blue curves in Figure 3). We wish to establish which of the 192 StarTrack models used in this work, after it has been combined with a specific galaxy's SFH and corrected for background AGN, best describes this oXLF. In other words, there are 192 tXLFs for this galaxy that need to be compared with a single oXLF. Thus, we calculate pair likelihoods for 192 o–t XLF pairs.

Given a galaxy and its oXLF, k, and a model tXLF, m, we define the k–mth pair likelihood, ${\cal L}_{{\rm pair,} km}$, as the probability, P, of obtaining an observational set of XLF data points, k, given a theoretical model XLF, m. This is given by

$$\begin{equation} {\cal L}_{{\rm pair,} km} = \prod _i P_{\rm Poisson} (N_{{\rm obs},i}, N_{{\rm th},i}). \end{equation} \tag{ 6 }$$

Here P_Poisson(N_{obs, i}, N_{th, i}) is the Poisson probability of observing N_{obs, i} point sources, in the i^th luminosity bin, treating the theoretically obtained number of point sources, N_{th, i}, in the same bin, as the expectation value for the number of sources at this luminosity. The total pair likelihood is the product of all such probabilities over all luminosity bins. This compares the agreement of the two XLFs in each luminosity bin, and thus its overall shape and normalization.

Note that tXLFs are calculated for 100 equal-sized bins in log L_X, spanning the observational range of luminosity values. The oXLFs are binned to match the tXLF bins. In addition, in order to compare corresponding quantities, N_{obs, i} are the observed numbers before correction for incompleteness, while N_{th, i} are the theoretical numbers after the observationally derived incompleteness information for a given bin has been taken into account.

For a given galaxy, the best model is that for which this procedure produces the highest pair likelihood value. We tabulate highest-likelihood results based on this procedure in Table 4 for all galaxies in our sample. Since the absolute numeric value of ${\cal L}_{{\rm pair,} km}$ obtained via Equation (6) has no particular meaning, each ${\cal L}_{{\rm pair,} km}$ value in Column 10 is normalized to the highest value in the table. This leads to a ranking, shown in Column 9, with the highest ${\cal L}_{{\rm pair,} km}$ value having rank 1. We also plot likelihood results against model number for all o–t pairs (i.e., not just highest-likelihood pairs) in the top panel of Figure 4. For convenience, in this figure, ${\cal L}_{{\rm pair,} km}$ values are normalized so that they range from 0 to 1 (see Section 4.3).

Zoom In Zoom Out Reset image size

The previous procedure determines the best model for the kth oXLF by estimating the probability, ${\cal L}_{{\rm pair,} km}$, of obtaining this XLF with each model. We are also interested in knowing the best model, m, for all oXLFs taken as a whole. To determine this, we calculate the product of all ${\cal L}_{{\rm pair,} km}$ values by defining the global likelihood

$$\begin{equation} {\cal L}_{{\rm global}, m} \equiv \prod _k {\cal L}_{{\rm pair}, km}, \end{equation} \tag{ 7 }$$

where ${\cal L}_{{\rm pair,} km}$ is given by Equation (6) and the index k runs from 1 to 12, corresponding to each of the 12 oXLFs. The results of this procedure are tabulated in Table 5, ranked by the ratio of each global likelihood value to the maximum likelihood value in the table (model 245). In Figure 3, we also show for each galaxy the tXLF that corresponds to the best global model (245) as a dark gray curve.

The specific pair or global numeric likelihood values obtained have no particular meaning, except relative to each other. For plotting purposes (Figures 4 and 5) it is convenient to normalize likelihood values so that the maximum value is 1 and the minimum 0. This is done as follows.

Zoom In Zoom Out Reset image size

Let $\cal L$ stand for either $\cal L_{\rm pair}$ or $\cal L_{\rm global}$. Then, the final normalized likelihood value is given by

$$\begin{equation} \Lambda _{\rm n} \equiv \Lambda _{\rm min} / {\rm max} (\Lambda _{\rm min}), \, \end{equation} \tag{ 8 }$$

where

$$\begin{equation} \Lambda _{\rm min} \equiv \ln [ {\cal L} / {\rm min} ({\cal L}) ]. \end{equation} \tag{ 9 }$$

In what follows we use an additional subscript to specify whether the final normalized likelihood is a pair or a global likelihood (Λ_{n, pair} and Λ_{n, global}, respectively).

A visual inspection of o–t XLF pairs in Figure 3 suggests that in many cases tXLFs (solid red curves) are successfully reproducing oXLFs (solid blue curves) both in shape and normalization. However, there is considerable variation and room for improvement. The best case is NGC 2841 and model 245, which, based on its likelihood estimate, has rank 1 in Table 4. Note that model 245 is also the best global model. The worst case is NGC 3184 and model 277 (rank 12 in Table 4). Since likelihood estimates take into account source numbers in luminosity bins (Equation (6)), the low likelihood estimate in this case is driven by the strong discrepancy in numbers at low luminosities. In many cases, the XLFs for the best pair model (red curve) and best global model (dark gray curve) are very close, but there are also cases where these are in strong disagreement.

The top panel of Figure 4 displays the normalized pair likelihood, Λ_{n, pair}, for all 192 × 12 o–t pairs (all 12 galaxies) versus all models used in this paper. The imposed normalization (see Section 4.3) is such that the maximum normalized likelihood for each o–t pair is 1 and the minimum is 0. Curves plotting normalized likelihood values as a function of model number (1–96 and 193–288) for each galaxy/o–t pair are shown with a different color. It is immediately obvious that for each galaxy there are a number of good models, indicated by maxima in the curves, as well as a number of poor models indicated by minima. Closer inspection of the plot reveals the remarkable fact that, overall, there is strong clustering of local maxima and minima, indicating that good models are good for all galaxies (and vice versa for poor models). This is so despite the fact that for a given galaxy the best models do not always match oXLFs well (Figure 3). This suggests that, to first order, the parameters associated with a good model are useful indicators of XRB physics across all 12 galaxies.

This global trend is corroborated by the lower panel of Figure 4, which shows the normalized global likelihood, Λ_{n, global}. As before, the normalization is such that the maximum is 1 and the minimum 0. Given the similarity of Λ_{n, pair} curves (top panel) with each other, it is not surprising that they are also similar to the global likelihood curve, Λ_{n, global}, as the latter combines, for a given model, all pair likelihoods for all galaxies (Equation (7)).

To understand what highest-likelihood results imply for individual model parameters, we tabulate results in two ways. We first rank highest-likelihood galaxy-model pairs according to pair likelihood value and show the results in Table 4. Clearly, this tabulation favors a λα_CE value of 0.1, a high-end IMF exponent of −2.7, and a mixed initial q distribution. Second, we rank all 192 models according to global likelihood value and show the results for the 15 highest ranked models in Table 5.¹⁴ These correspond to the 15 highest peaks in the lower panel of Figure 4. These results also favor a λα_CE value of 0.1, and to a lesser extent a high-end IMF exponent of −2.7 and a mixed initial q distribution.

By comparing Tables 4 (Column 2) and 5 (Column 1), we note that all five best global models are also among the best pair models. This just underlies the fact that these are the best models overall. Conversely, with the exception of galaxy NGC 5474, all pair models are also among the best global models. It is remarkable that for 11 out of 12 galaxies, one of the best 15 models that represent global averages over all galaxies also describes the individual galaxy XLF.

We note that for galaxy NGC 1291, Luo et al. (2012) use Chandra observations that are deeper than ours to perform a detailed study of the point source population. They show that high-luminosity (>5 × 10³⁸ erg s⁻¹) LMXBs in NGC 1291 are likely associated with a younger SP in the galaxy's ring. Their deeper observations allow them to separate the bulge from the ring population. Even so, we note that our results still support their conclusions, as we find strong LMXB contributions both from our old and young populations.

A key reason for deviations of oXLFs from tXLFs is the presence of two "jumps" at two different luminosities (∼10³⁷ erg s⁻¹ and a few × 10³⁸ erg s⁻¹), which are not expected from observations. There are two reasons for this. First, in our StarTrack models we identify transient and persistent sources by comparing the calculated mass-transfer rate of the binary to a critical mass-transfer rate, below which the thermal instability develops, giving rise to transient behavior. Furthermore, we strictly limit the accretion rate to the Eddington limit. These limiting mass-transfer rates in our modeling are responsible for the jumps in tXLFs. However, in nature there is neither a sharp transition between thermally stable and unstable disks, nor a precise limit to the highest accretion rate possible. Accretion onto a compact object is a nonlinear and much more complex process, which in reality can result in a smoother luminosity distribution with no sharp transitions. Second, the assumption of solar metallicity in practice imposes a maximum BH mass of ∼15 M_☉. If, for instance, part of the SP had a lower metallicity (e.g., 30% solar), then the maximum BH mass would increase to 30 M_☉ or more (Belczynski et al. 2010), which could smooth out the jump at the very luminous end of the XLF.

More generally, it must be stressed that the SFHs we use are very simple. Noll et al. (2009) note that, as SED fitting is computationally intensive, they select a limited set of values for their SED parameter model grid. Thus, the age of the old population is a constant (10 Gyr), and there are only two possible values for the age of the young population, 200 Myr and 50 Myr. The parameter that was allowed to vary the most is "f_burst,"the mass fraction of the young SP at the present time (nine possible values). Further, the detailed analysis of their SED fitting results shows that best-fit e-folding times and the age of the young population are highly uncertain due in turn to photometric errors and uncertainties in SP models. Finally, as already mentioned, their SFHs have uniformly solar metallicities. Although this is based on the best estimates to date, note that these are based on gas metallicities with their own set of significant uncertainties (see Moustakas et al. 2010, where two sets of metallicity results are presented, straddling solar metallicity for all galaxies in this sample). Further progress in matching oXLFs will require more detailed SED fitting.

Table 4 and especially Table 5 suggest likely best values for model parameters in our model grid. As mentioned, both tables suggest λα_CE ≃ 0.1. In Figure 5, we investigate this further by plotting resistant mean¹⁵ values for the normalized global value, Λ_{n, global}, versus the four StarTrack λα_CE parameter values used in our simulations. The calculation of the resistant mean iteratively rejects outliers beyond 3σ. This is useful for highlighting the clustering of likelihood values, which often show a considerable spread.

Figure 5 shows that higher likelihood values are systematically favored for models with λα_CE ≃ 0.1: not only do Λ_{n, global} values converge to their highest value as λα_CE→ 0.1, but they also do so with a decreasing spread, as indicated by the 3σ resistant mean error bars.

Tables 4 and 5 also appear to favor a value of −2.7 for the high-end slope of the IMF. Although a value of −2.7 agrees with some estimates (e.g., Scalo 1986), we note that it is somewhat steep compared to the so-called canonical value of −2.3 established observationally for resolved SPs in the Local Group (Bastian et al. 2010; Kroupa 2012). In spite of the fact that (1) variations with metallicity and SFR density for the integrated galaxy IMF remain hotly debated (Kroupa 2012 and references therein) and (2) taking uncertainties in the canonical value into account, a value of −2.7 is still close to the canonical upper limit, and we do not consider a value of −2.7 to be a robust result. This is because, first, the Noll et al. SED fitting results are assuming the canonical value, and it is these results that have been convolved with our StarTrack models. Second, four out of the fifteen best global models do, in fact, agree with the canonical value.

For the remaining parameters, these tables show that the results are inconclusive. In addition, for most of these parameters we only use two different values in the StarTrack simulations, so we cannot investigate any trends such as those for λα_CE in Figure 5.

Given the caveats for the IMF slope, the results for λα_CE ≃ 0.1, and to some extent for the prevalence of a mixed initial q distribution, provide the most significant constraints for binary star parameters from this work.

As λα_CE is a combination of two parameters, we are unable to set any constraints on either λ or α_CE individually. In their investigation of Galactic merger rates for compact objects, Dominik et al. (2012) set α_CE =1.0 and studied the behavior of λ. The latter is not considered constant throughout the evolution of the donor, but depends on donor parameters such as mass, radius, and evolutionary stage. They found λ values between ∼0.1 and ∼0.2 for NS progenitors, and below ∼0.1 for BH progenitors. With the assumption of α_CE =1.0, our result for λα_CE would thus be largely consistent with the detailed stellar evolution models of Dominik et al. (2012).

F13 and T13 used the same StarTrack models as this paper but combine them with galaxy information from the Millennium II simulation and the semianalytic galaxy catalog of Guo et al. (2011). While F13 investigated the total galaxy-specific X-ray luminosity (L_X/SFR, L_X/M_*) evolution out to z ∼ 20, T13 constructed galaxy XLFs for comparison with observations out to z ∼ 1.4 and made predictions out to z ∼ 20. Thus, our three papers examine XRB formation and evolution in different contexts, and also use different likelihood formulations. To investigate whether there is agreement in highest-likelihood models between the three papers, we also show the ranks from F13 and T13 for our best 15 global likelihood models in Table 5. There is good agreement overall, and in certain cases the agreement is exceptionally good. In particular, our best ranked model (245) is also F13's reference model (their rank 1) and has rank 4 in T13. Models 245, 277, 229, 205, 269, 249, 273, 37, and 85 are all within the 15 highest ranked in all three papers. This is further evidence in support of λα_CE ≃ 0.1 both in the local universe and over cosmic time.

In Figures 6 and 7, we plot the total tXLF for each galaxy and highest-likelihood model (the black continuous curve), together with constituent subpopulation tXLFs.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

We define two sets of subpopulations. The first is shown in Figure 6 where we plot XLFs for LMXBs (red) and HMXBs (blue). In our simulations, XRBs are labeled HMXBs if the donor star has mass M_donor ⩾ 3 M_☉, and LMXBs otherwise. The mixture of populations present in each galaxy is closely tied to the assumed SFH (Table 3). To illustrate this, XRBs originating in the "old" SP are shown by the dashed curve in Figure 6, while those originating in the "young" SP are shown by the solid curve. Note however that the terms old and young can be misleading in this context. Although the two are distinct in terms of age, with the young one appearing several Gyr after the old one, it is possible for the latter to still be actively star-forming, depending on the e-folding timescale τ_oSP. In other words, an "old" population is not necessarily a "red and dead" one.

The trends seen in Figure 6 for each galaxy can be understood qualitatively by referring to (1) the SFH (Table 3), which we combine with StarTrack models to construct tXLFs, and (2) Table 4, which gives the details of each best model for each galaxy. We should also keep in mind the relative evolutionary timescales for LMXBs and HMXBs. As shown in F13, ¹⁶ (their Figure 2), in a single burst population the contribution in X-ray output from HMXBs peaks at about 5 Myr, remaining important up to an age of ∼100–300 Myr, while LMXBs take over at ∼100–200 Myr. In terms of the assumed "young" and "old" SPs, a population's contribution will be more significant as its e-folding timescale is longer and age is younger. A higher young mass fraction will tend to increase the contribution from the young population.

Looking at galaxies NGC 1291 and 2841, we note that they have all best-fit SFH parameters the same, and slightly different M_gal, while the associated StarTrack model is the same. It is evident that the subpopulation XLFs are similar, with the old LMXB SP dominating in both cases.

For galaxies NGC 3184 and 3627, the main SFH difference is a higher τ_oSP for the latter. Based on this, one might expect that in NCG 3627 the old population would be more dominant. However, young SP HMXBs dominate for most of the XLF. This is due to StarTrack model 245 (NGC 3627) having half the stellar wind strength compared to model 277 (NGC 3184). Weaker stellar winds lead to smaller mass loss for primaries that eventually become compact objects, and thus to numerous and more massive BH XRBs. The latter in turn tend to be more luminous than NS XRBs as (1) they can form stable RLOF XRBs with massive companions, and (2) they show higher accretion rates due to the high BH masses. Although a weaker stellar wind also decreases the accretion rate in wind-fed HMXBs, it turns out that this is not the dominant effect (F13).

NGC 3198 has by far the highest τ_oSP among all galaxies and this is reflected in the dominance of the old LMXB population. In contrast, NGC 3521 and 5055 have a much lower τ_oSP (with other parameters except for M_gal identical for all three) and are dominated by young SPs. NGC 4631 is also similar, with a somewhat lesser contribution from the old SP (higher τ_ySP). Compared to NGC 4631, NGC 4826 has lower τ_ySP and f_burst, leading to a relatively stronger contribution from the old SP.

NGC 4736's young SP has the smallest best-fit age (50 Myr) together with a high τ_ySP. This is consistent with the clear dominance of young HMXBs, which are closest to their peak activity timescale (Shtykovskiy & Gilfanov 2007), while LMXBs have not yet had time to contribute significantly to the total X-ray luminosity.

Finally, NGC 5474 has higher τ_oSP but a higher young mass fraction compared to, e.g., NGC 3521, so it is moderately dominated by the young SP.

The second subpopulation set is shown in Figure 7 in terms of donor and accretor stellar type. We classify donors as Main Sequence (MS), evolved, or degenerate. Evolved types have left the MS and can be in any of a number of different stellar evolutionary stages, including the Hertzsprung gap, the red giant branch (GB), core helium burning, early asymptotic GB, thermally pulsating asymptotic GB, helium MS, helium HG, or the helium GB. Degenerate types are helium or carbon/oxygen white dwarfs. On the other hand, accretors are either BHs or NSs.

A number of trends are visible in these plots. Considering the accretor subpopulations, indicated by different line styles in Figure 7, we note that, in general, BH accretors (solid lines) dominate the XLFs at high luminosities. In addition, galaxies with no sources at L_X ≳ 3 × 10³⁸ have no BH accretors. These galaxies are NGC 3184, 3521, and 4631. The reason for this is that, since BH accretors originate in the high-mass end of the IMF, from a statistical point of view, on average there will be fewer such systems. In addition, BH XRBs can have higher luminosities than NS XRBs, but not vice versa. Thus, BH XRBs can only be registered when very luminous sources are present (e.g., see Luo et al. 2013, for NGC 4649), as NS XRBs will always dominate at low luminosities. As a result, for galaxies with very few or no sources above ∼10³⁸ erg s⁻¹, models are unable to constrain the BH XRB population.

Donor subpopulation XLFs, indicated by different colors in Figure 7, are dominated by MS donors at low luminosities (<10³⁷ erg s⁻¹)¹⁷ and by evolved donors at high luminosities. This is a consequence of the fact that evolved donors are mainly giants in RLOF systems, so they all have large accretion disks and long orbital periods. Since $\dot{M}_{\rm crit}$ depends on the size of the accretion disk and orbital period, it follows that in the case of these systems it will always be higher than a threshold value, which corresponds to L_X ∼10³⁷ erg s⁻¹ (King et al. 1996; Dubus et al. 1999). Further, systems with mass-transfer rates less than $\dot{M}_{\rm crit}$ and, thus, with luminosities lower than ∼10³⁷ erg s⁻¹ will be transients, which are mostly quiescent. Thus, in practice, at L_X ≳ 10³⁷ the dominant donor systems will be persistent systems with evolved donors.

It is well established that emission from HMXBs and LMXBs correlates with galaxy-wide SFR and M_*, respectively. This is due to the fact that the former are relatively young (≲ 100 Myr) compared to the latter (≳ 1 Gyr). Lehmer et al. (2010, hereafter L10) parameterize the contribution of LMXBs and HMXBs to the 2–10 keV integrated luminosity in star-forming galaxies as

$$\begin{equation} L_{\rm HX} = L_{\rm HX}^{\rm LMXB} + L_{\rm HX}^{\rm HMXB} \equiv \alpha \ M_{\star } + \beta \ {\rm SFR}, \end{equation} \tag{ 10 }$$

where L_HX ≡ L_{X, 2.0–10.0 keV}. They use a sample of galaxies spanning a large range in star-forming activity observed with Chandra to obtain best-fit values $\alpha = (9.05\pm 0.37) \times 10^{28} {\rm \,erg \, s}^{-1} \,M_{\odot }^{-1}$ and β = (1.62 ± 0.22) × 10³⁹ erg s⁻¹( M_☉ yr⁻¹)⁻¹. We use these results for comparison with our work, noting also that they are consistent with other work both on the L_X–M_* relation in elliptical galaxies (Gilfanov et al. 2004; Zhang et al. 2011; Boroson et al. 2011) and on the L_X–SFR relation in galaxies with high specific SFR, where HMXBs are dominant (Gilfanov 2004; Mineo et al. 2012).

We use our subpopulation XLFs for LMXBs and HMXBs to compare with the L10 results as follows. We obtain galaxy-wide X-ray luminosities by integrating each model XLF over all luminosity bins. Note that for each galaxy we use both the highest pair likelihood model XLF and the XLF that is based on the best global model 245, thus obtaining two sets of galaxy-wide X-ray luminosities. We then use the M_* and SFR values for each galaxy to estimate $L_{\rm HX}^{\rm LMXB}$ and $L_{\rm HX}^{\rm HMXB}$ from Equation (10). Since our data are in the 0.3–10 keV band, we convert L_HX values to L_{X, 0.3–10.0 keV} by using the mean L_{X, 0.3–10.0 keV} to L_HX ratio based on the F13 models. In Figure 8, we plot the integrated X-ray luminosities for LMXBs and HMXBs against M_* and SFR. The gray open circles represent integrated luminosities that use the highest pair likelihood model XLF, while the filled black circles represent model XLFs using the best global model 245. We also show the relations $L_{\rm HX}^{\rm LMXB}\hbox{--}$ M_* and $L_{\rm HX}^{\rm HMXB}\hbox{--}{\rm SFR}$ for the best-fit α and β values of L10, with the associated 1σ scatter.

Zoom In Zoom Out Reset image size

The comparison shows that, when the highest pair likelihood models are used, there is some agreement, especially for LMXBs, but some of these fail to reproduce the L10 relations. In contrast, when the best global model is used, integrated luminosities are in much better agreement with the L10 relations, at least for HMXBs. The residual scatter is likely due to limitations of SFHs and uniformly solar metallicities. The observed agreement suggests that averaging over all 12 galaxies to calculate global likelihoods produces results that, in a statistical sense, are more reliable. On the other hand, individual pair likelihoods inevitably suffer from all of the SFH uncertainties discussed earlier, even though a given tXLF might match a given galaxy oXLF better than the tXLF for model 245. This result is fully consistent with the fact that model 245 is also the top-ranked model of F13. Using an independent likelihood formulation, these authors constrain their best models by comparing with observational work, including the L10 relations.

This result provides strong motivation for future work. One would expect that for a larger data set, such as all 75 SINGS galaxies, the agreement would be further improved.