THE CONTRIBUTION OF X-RAY BINARIES TO THE EVOLUTION OF LATE-TYPE GALAXIES: EVOLUTIONARY POPULATION SYNTHESIS SIMULATIONS

We used the EPS code developed by Hurley et al. (2000, 2002) and updated by Liu & Li (2007, see Appendix A in the paper) and Zuo et al. (2008) to calculate the X-ray luminosity L_X of XRBs, the optical luminosity L_B, and the stellar mass M of the host galaxy, as well as their evolution. The values of the adopted parameters are the same as the default ones in Hurley et al. (2002) if not mentioned otherwise.

Previous works (Shapley et al. 2001; Persic & Rephaeli 2007; Lehmer et al. 2008) have already shown that the galactic X-ray emission is closely related to the SF activity. So we constructed three cases, i.e., constant SF, starburst SF, and cosmic SF cases, to examine their effect. We also examined several key parameters, such as the IMF, CE efficiency parameter, the binary fraction, and metallicity (listed in Table 1 and discussed below), to explore their influence on the X-ray evolution of the galaxies.

1. Constant star formation case

In this case, we adopt a constant SFR =0.25 M_☉ yr⁻¹ for stars more massive than 5 M_☉ derived by Grimm et al. (2003) in our Galaxy and take the star formation duration (SFD) to be 14 Gyr. For each model, we evolve 10⁶ primordial binary systems, with the same grid of initial parameters (i.e., primary and secondary mass, orbital separation) as Hurley et al. (2002). We also evolve 10⁶ primordial single stars, with initial mass logarithmically spaced between 0.1 M_☉ and 80 M_☉. In our basic model (i.e., model M1, listed in Table 1), we assume the binary fraction f to be 0.5 and evolve each binary and single star on the grid. In the following, we describe the assumptions and input parameters in our basic model.

In order to be in parallel with Lehmer et al. (2008), we take the IMF of Kroupa (2001, hereafter KROUPA01) for the mass (M₁) distribution of the primary stars. For the secondary stars (of mass M₂) and binary orbit, we assume a uniform distribution between 0 and 1 for the mass ratio q ≡ M₂/M₁ and a uniform distribution for the logarithm of the orbital separation ln a (Hurley et al. 2002). We fix the metallicity to be solar over the lifetime of the simulated galaxy.

We assume that any system entering Roche-lobe overflow becomes circularized and synchronized by tidal interaction between the binary components (Belczynski et al. 2008). An important parameter in the binary evolution is the CE efficiency parameter α_CE (Paczyński 1976; Iben & Livio 1993), which describes the efficiency of converting orbital energy into the kinetic energy ejecting the envelope (see Section 2.1.1 in Zuo & Li 2010, for details). It can often reduce the orbital separation of the surviving binaries by a factor of ∼100, resulting in different outcomes of binary systems. In our basic model, we adopt α_CE = 0.3, which can best model the luminosity function of the galaxy (Zuo et al. 2008).

We also construct several other models (listed in Table 1) by varying the key input parameters described as follows.

• 1.  
  As stated above, variations of the CE parameter can considerably change the relative numbers of XRBs. However reliable values of α_CE are difficult to estimate due to lack of understanding of the processes involved, although in the literature it is in the range from ∼0.1 to ∼3.0 (e.g., Taam & Bodenheimer 1989; Tutukov & Yungelon 1993; Podsiadlowski et al. 2003). Here we also adopt α_CE = 1.0 (model M2) to examine its effect.
• 2.  
  The IMF determines the percentage of high-mass stars, consequently the number of XRBs produced, and the X-ray luminosity of the galaxy. So, we also make use of the IMF of Kroupa et al. (1993, hereafter KTG, model M5), which is much steeper in the high-mass end than in KROUPA01. For the secondary masses (M₂), we assume the mass ratio q follows a power-law distribution P(q) ∝ q^α, and adopt both the conventional choice of flat mass spectrum, i.e., α = 0 (Mazeh et al. 1992; Goldberg & Mazeh 1994; Shatsky & Tokovinin 2002, our basic model, M1) and α = 1 (model M3), since recent data are more consistent with "twins" being a general feature of the close-binary population (Dalton & Sarazin 1995; Kobulnicky & Fryer 2007).
• 3.  
  Surveys of M dwarfs within 20 pc from the Sun indicate that the binary fraction f may be a function of stellar spectral types (Fischer & Marcy 1992), for example, f>0.5 for G stars and f>0.6 for massive O/B stars in the Cygnus OB2 association (Lada 2006; Kobulnicky & Fryer 2007). So, we also adopt f = 0.8 (model M4) for comparison.
• 4.  
  Observations of the hosts of XRBs revealed that XRBs, especially ultra-luminous X-ray sources, may prefer to occur in galaxies with low metallicities (Mapelli et al. 2009). So, we vary the metallicity to examine its effect on the X-ray luminosity evolution by taking Z = 1.5 Z_☉ (model M6), 0.5 Z_☉ (model M7), 0.1 Z_☉ (model M8), and 0.02 Z_☉ (model M9), in order to compare with our basic model (M1 with  Z_☉). The other parameters in these models are the same as the ones in our basic model.

2. Starburst star formation case

In galaxies like our own Galaxy, continuous SF processes may last for several Gyr, but this is not the case for starburst galaxies. For example, the Antennae galaxies may have experienced SF for the last several hundred Myr with an enhanced SFR of 7.1 M_☉ yr⁻¹ (Barnes 1988; Mihos et al. 1993). To reveal the effect of SFH, we take the following combinations of SFD and metallicities: 100 Myr/ Z_☉ (model M10), 20 Myr/ Z_☉ (model M11), and 100 Myr/0.02 Z_☉ (model M12). We assume that the SF is quenched after the SFD time and set other parameters to be the same as in our basic model.

3. Cosmic star formation case

With improved observations of SF processes, cosmic SFH can be constrained quite tightly within ∼30%–50% up to z∼1 and within a factor of ∼3 up to z ∼ 6 (Hopkins 2004), which makes it possible to investigate the cosmic X-ray evolution of galaxies. Here we adopt the derived expression of the SFH in Hopkins & Beacom (2006),

$$\begin{eqnarray} \dot{\rho }_{\rm SF}(z)&\propto &\left\lbrace \begin{array}{@{}ll@{}}(1+z)^{3.44}&\ \ z\le 0.97 \\ (1+z)^{-0.26}&\ \ 0.97\le z\le 4.48 \\ (1+z)^{-7.8}&\ \ 4.48\le z, \\ \end{array} \right. \end{eqnarray} \tag{ 1 }$$

and scale the SFR at redshift z = 0 to be the same as that of our Galaxy. Moreover, it is known that the cosmic metallicity also evolves strongly with redshift, and galaxies at higher redshift tend to have lower metallicities (Pettini et al. 1999; Prochaska et al. 2003; Rao et al. 2003; Kobulnicky & Kewley 2004; Kewley & Kobulnicky 2007; Kulkarni et al. 2005, 2007; Savaglio et al. 2005, 2009; Savaglio 2006; Wolf & Podsiadlowski 2007; Péroux et al. 2007). So, we adopt an empirical equation Z/Z_☉ ∝ 10^−γz (Langer & Norman 2006) for metallicity evolution with γ = 0.15 (Kewley & Kobulnicky 2007). We also vary the IMF to a steeper one (KTG93) and a shallower one (Baldry & Glazebrook 2003, BG03 for short) to examine its effect in this case. The other parameters are the same as in our basic model.

We adopt the same procedure to calculate the 0.5–8 keV X-ray luminosities of different XRB populations as in Zuo & Li (2010). Mass transfer in XRBs occurs via either Roche-lobe overflow or capture of the wind material from the donor star. We use the classical Bondi & Hoyle (1944) formula to calculate the wind accretion rate of the compact stars. In the case of Roche-lobe overflow, mass is transferred to the accreting star by way of an accretion disk. It is known that accretion disks in LMXBs are subject to thermal instability if the accretion rate is sufficiently low (van Paradijs 1996). We discriminate transient and persistent LMXBs according to the criteria of van Paradijs (1996) for main sequence (MS) and red giant donors and of Ivanova & Kalogera (2006) for white dwarf (WD) donors, respectively. The simulated X-ray luminosity is described as follows,

$$\begin{eqnarray} &&L_{\rm X, 0.5-8\; keV}\nonumber\\ && =\left\lbrace \begin{array}{@{} ll}\eta _{\rm bol}\eta _{\rm out}L_{\rm Edd}&\ \rm transients\ in\ outbursts, \\ \eta _{\rm bol}\min (L_{\rm bol},\eta _{\rm Edd}L_{\rm Edd})&\ \rm persistent\ systems, \end{array} \right. \end{eqnarray} \tag{ 2 }$$

where the bolometric accretion luminosity $L_{\rm bol}\simeq 0.1\dot{M}_{\rm acc}c^2$ (where $\dot{M}_{\rm acc}$ is the accretion rate and c is the velocity of light), the critical Eddington luminosity L_Edd ≃ 4πGm₁m_pc/σ_T = 1.3 ×  10³⁸m₁ erg s⁻¹ (where σ_T is the Thomson cross section, m_p is the proton mass, G is the gravitational constant, and m₁ is the accretor mass in units of the solar mass), and η_Edd is the factor to allow super-Eddington luminosities, taken to be 5 (Ohsuga et al. 2002; Begelman 2002). To transform the bolometric luminosity into the 0.5–8 keV X-ray luminosity, a bolometric correction factor η_bol is introduced (Belczynski & Taam 2004). Generally, its value is ∼0.1–0.5 for different types of XRBs; here we adopt η_bol ≃ 0.1. For transient sources, the X-ray luminosity during outbursts should be larger than the long-term one by a factor η_out. We take η_out = 0.1 and 1 for the short- and long-period systems, and the critical periods are adopted to be 1 day for neutron star (NS) transients and 10 hr for black hole (BH) transients, respectively (Chen et al. 1997; Garcia et al. 2003; Belczynski et al. 2008).

The optical luminosity L_B of a galaxy is mostly from normal stars (both binary and single stars). Assume that the stellar radiation can be reasonably approximated as a blackbody, the B-band luminosity of a star is calculated with $L_{\rm B}=\frac{\int _{\rm B}I_{\lambda } d \lambda }{\int I_{\lambda } d \lambda } \times L$, where L is the total thermal luminosity of the star, and the radiative intensity $I_{\lambda }=\frac{2hc^2}{\lambda ^5}\frac{1}{e^{hc/\lambda kT_{\rm eff}}-1}$, where h is the Planck constant, k is the Boltzmann constant, and T_eff is the effective temperature.

We also examine the contribution of optical luminosity from accretion disks in XRBs, resulting from the reprocessing of X-ray photons. We calculate the optical luminosity L_B from the accretion disk in BH XRBs following Madhusudhan et al. (2008), adopting the same temperature profile (i.e., their Equation (4)) to describe the effective temperature in the disk. For NS XRBs, we find similar results as BH XRBs. Our calculation reveals that optical radiation from accretion disks in XRBs is negligible compared to the overall stellar optical luminosity.

The stellar mass here is the sum of the masses of currently living stars and does not include the contribution from compact stars (WDs, NSs, and BHs), in order to be in parallel with Lehmer et al. (2008), where they used the rest-frame B − V color and K-band luminosity to estimate the masses of the galaxies.

Figure 1 shows the calculated values of L_X, L_B, L_X/M, L_X/L_B, L_X/(M/L_B), and M/L_B against time in the constant SF case. The five panels from left to right correspond to the basic model (M1), models with α_CE = 1.0 (M2), α = 1 (M3), f = 0.8 (M4), and KTG93 IMF (M5), respectively.

Zoom In Zoom Out Reset image size

For X-ray luminosities L_X (Figure 1(a)) we plot the contributions from HMXBs and LMXBs with dotted and dashed lines, respectively. The X-ray luminosity of HMXBs rises rapidly shortly after the first SF and remains nearly constant afterward, because of continuous, constant SF. LMXBs take a much longer time to form, and their number is correlated with the total star mass, resulting in a long-term increasing trend in X-ray luminosity with time. The position of the crossing point of the two lines depends most strongly on the CE parameter α_CE (larger α_CE results in more LMXBs).

The optical luminosity L_B (Figure 1(b)) is contributed by the primary stars (dotted line) and secondary stars (dashed line) in binaries, and by single stars (long-dashed line). When the fraction of binaries is larger (model M4), more XRBs are produced, leading to larger L_X, and hence larger L_X/M, L_X/L_B and L_X/(M/L_B) ratios compared with those in the basic model. The steeper end of IMF (i.e., model M5) implies a smaller number of massive stars (and smaller L_B), resulting in fewer compact objects and smaller L_X. So the L_X/M and L_X/(M/L_B) ratios in this case are smaller than those in models M1–M4.

The L_X/M ratio (Figure 1(c)) has a clear decreasing trend after the age of tens of Myr. Observationally, this phenomenon was regarded as evidence that lower mass galaxies appear to have higher specific SFRs than massive ones (Brinchmann et al. 2004; Bauer et al. 2005; Noeske et al. 2007a; Feulner et al. 2005; Zheng et al. 2007). Our results show that the decrease of L_X/M is a natural consequence of XRB evolution and the stellar mass accumulating process in galaxies—the stellar mass always steadily increases while the X-ray luminosity changes little during most of the evolution.

The L_X/L_B ratio (Figure 1(d)), similar to L_X, rises in the first several Myr, then remains nearly constant afterward. The flattening values are all around 10³⁰ erg s⁻¹L⁻¹_B,☉, but have severalfold changes among different models. They are caused by the diversity in the percentages of both total massive stars and massive stars in binaries, which determine L_B and L_X, respectively. The ratios of the peak values of log(L_X/L_B) are roughly 1:2:0.5:1.5:1 from models M1 to M5.

Figure 2 shows the cumulative X-ray luminosity functions (XLFs) of HMXBs (dotted line) and LMXBs (dashed line) in the top panels. The parameters are the same as in Figure 1. Note that HMXBs and LMXBs dominate at relatively high (>10³⁹ erg s⁻¹) and low (<10³⁹ erg s⁻¹) luminosity ends, respectively. We also show the detailed components of XRBs which contribute the total X-ray luminosity separately in the middle (HMXBs) and bottom (LMXBs) panels of Figure 2. The dotted, dashed, dash-dotted, and dash-dot-dotted lines represent persistent BH XRBs (BHp), transient BH XRBs (BHt), persistent NS XRBs (NSp), and transient NS XRBs (NSt), respectively. It is seen that for HMXBs persistent BH XRBs contribute most to the XLF; for LMXBs, BH-XRBs (both BHp and BHt) dominate the high-luminosity end of XLF, while in the low-luminosity end transient NS-XRBs play a more important role.

Zoom In Zoom Out Reset image size

Figure 3 shows the evolution of L_X, L_B, L_X/M, L_X/L_B, L_X/(M/L_B), and M/L_B with metallicities taken to be 1.5 Z_☉, Z_☉, 0.5 Z_☉, 0.1 Z_☉, and 0.02 Z_☉, corresponding to models M6, M1, M7, M8, and M9 from left to right, respectively. Note that both L_X and L_B have a roughly increasing trend with decreasing metallicity. The values of L_X/M and L_X/L_B are comparable among different models. The values of the peak values are ∼1:1:0.7:1:2 for log(L_X/M) and ∼1:1:0.5:1:1 for log(L_X/L_B). The corresponding cumulative XLF is shown in Figure 4.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

We present the evolution of L_X, L_B, L_X/M, L_X/L_B, L_X/(M/L_B), and M/L_B in the starburst case in Figure 5. Here the metallicity and SFD are assumed to be  Z_☉/100 Myr (left), Z_☉/20 Myr (middle), and 0.02 Z_☉/100 Myr (right). It is seen that the X-ray luminosity L_X (Figure 5(a)) rises to its peak within the SF episodes (peaked at the age of about 100 Myr for models M10 and M12, and at about 20 Myr for model M11), then decreases with time, lasting at least several 10⁸ yr with L_X>10³⁷ erg s⁻¹. Note that the small serration emerging in late evolution is mainly caused by LMXBs (i.e., X-ray transient outbursts due to thermal instability in the accretion disks). The optical luminosity L_B (Figure 5(b)) decreases sharply when the SF process stops, since massive stars contribute significantly during the starburst episode. The L_X/M (Figure 5(c)) (and also L_X/(M/L_B), Figure 5(e)) ratio roughly follows the trend of L_X. However the slope in this case is much steeper than in Figures 1 and 3, because of the rapid decay of the X-ray luminosity. The L_X/L_B (Figure 5(d)) ratio in the three models is all comparable with those in Figures 1 and 3, implying that it is intrinsically not sensitive to the SFH of the galaxies.

Zoom In Zoom Out Reset image size

Figure 6 shows the evolution of L_X, L_B, L_X/M, L_X/L_B, L_X/(M/L_B), and M/L_B with time (left) and redshift z (right), with a cosmic SFH (from Hopkins & Beacom 2006) and a cosmic metallicity evolution history (from Langer & Norman 2006) taken into account. Note that the X-ray luminosity L_X is mainly dominated by HMXBs in this case over the whole cosmic history because of the enhanced SFR with increasing redshift as a whole. This is comparable with the work of Lehmer et al. (2008), where they found that LMXBs on average play a fairly small role in the X-ray emission. The L_X/M ratio has a decreasing trend after the age of tens of Myr (or increases with z), similar to the constant/burst SF cases. This also confirms our previous conclusion that the decrease of L_X/M with time results from XRB evolution and stellar mass growth in the galaxies. The L_X/L_B ratio rises rapidly in the first several Myr, then stays around 10³⁰ erg s⁻¹L⁻¹_B,☉ afterward, similar to the constant/burst SF cases, indicating that the L_X/L_B ratio is not sensitive to the SFH in the galaxies.

Zoom In Zoom Out Reset image size

To compare the theoretical predictions with observations, we replot Figures 6(c) and 6(d) in Figures 7(a) and 7(b), with redshift z ranging from 0 to 2.0. The solid, dotted, and dashed lines represent the modeled results with IMFs of KROUPA01, Kroupa et al. (1993, KTG93), and Baldry & Glazebrook (2003, BG03), respectively. The observational data of log(L_X/M) and log(L_X/L_B) are taken from Shapley et al. (2001, S01) for normal late-type galaxies in the local universe (open symbols), Lehmer et al. (2008, L08) with the stellar masses 10^10.1 <  M/M_☉ < 10^11.2 (squares, Figure 7(a)), B-band luminosities 10^10.5 < L_B/L_B,☉ < 10^11.3 (filled circles, Figure 7(b)) and 10^10.0 < L_B/L_B,☉ < 10^10.5 (filled squares, Figure 7(b)), and Zheng et al. (2007, Z07) with 10^10.0 <  M/M_☉ < 10^10.5 (diamonds, Figure 7(a)). In these samples, the stellar masses are roughly comparable with our simulated ones. To convert the SSFR into L_X/M, we have made use of the local L_X–SFR relation derived by Persic & Rephaeli (2007).

Zoom In Zoom Out Reset image size

It seems that our simulated log(L_X/M) versus z relations match the observations quite well. Our calculations reveal that the modeled stellar masses are similar to each other within 10% when we use different types of IMFs, suggesting that the variation of the L_X/M ratio is mainly caused by the differences in X-ray luminosity L_X. The values of L_X/M can vary by a factor of ∼7 between models (KTG93 versus BG03), as seen in Figure 7(a).

The nearly constant values of L_X/L_B in Figure 7(b) seem to not properly match the observed increase in L_X/L_B with z (Hornschemeier et al. 2002; Lehmer et al. 2008). The discrepancy originates from the fact that in our simulations we have roughly L_X ∝ L_B, giving a flat L_X/L_B − z relation, while observationally it was found that L_X ∝ L^1.5_B (Shapley et al. 2001; Fabbiano & Shapley 2002; Lehmer et al. 2008), leading to increasing L_X/L_B with z. Fabbiano & Shapley (2002) have discussed the X-ray-B-band luminosity correlation and suggested that the nonlinear power-law dependency in disk galaxies is likely to be due to extinction in dusty star-forming regions, which attenuates light from the B band more effectively than it does in the X-ray band. This means that the intrinsic B-band luminosity calculated here is generally larger than the measured one, which suffers local extinction in the galaxies. Supporting evidence for this hypothesis also comes from the strong positive correlations between L_X/L_B and the ultraviolet dust-extinction measure (L_IR + L_UV)/L_UV, and the correlation between L_X/L_B and IR color (Lehmer et al. 2008). Thus, if the increase in L_X/L_B with z in the late-type galaxies is due to an increase in their SF activity (Fabbiano & Shapley 2002; Lehmer et al. 2008), our results are compatible with the observational data at least qualitatively.

Figure 8 shows the corresponding cumulative XLFs. In general, they are similar to those in the case of constant SF.

Zoom In Zoom Out Reset image size

Our simulations are obviously subject to many uncertainties. For current population synthesis investigations, it is difficult to tell confidently which parameter combinations are the best or most realistic by comparison with observations, since there are many uncertainties in the (both explicit and implicit) assumptions and input parameters, and simplifications in the treatment of the detailed evolutionary processes. For example, the simulated X-ray luminosity L_X depends on the adopted values of several parameters, such as the bolometric correction factor η_bol, the CE efficiency parameter, and so on, which may alter its value severalfold. This further affects the values of L_X/M and L_X/L_B since the stellar mass M and the optical luminosity L_B are not sensitive to these parameters. Despite these limitations, it has become possible to investigate the overall evolution of stars in galaxies with population synthesis and to draw useful information by comparing theoretical predictions with observations.