MODELING THE REDSHIFT EVOLUTION OF THE NORMAL GALAXY X-RAY LUMINOSITY FUNCTION

The Millennium Cosmological Simulation is an unprecedented computational effort to simulate the dark matter distribution in the universe (see Springel et al. 2005 for details). In this study, we use the data from the most recent Millennium Run II (hereafter MRII; Boylan-Kolchin et al. 2009). This is an N-body simulation that follows the evolution of 2160³ particles, each of mass 6.9 × 10⁶ h⁻¹ M_☉ within a comoving box with sides each of size 100 h⁻¹ Mpc. The cosmological model used in the simulation is a Λ cold dark matter model with Ω_m = 0.25, Ω_Λ = 0.75, Ω_b = 0.045, and h = H₀/100 km s⁻¹ Mpc⁻¹ = 0.73.

The MRII has 60 snapshots in time that were saved and analysis was done to identify substructures within the dark matter distributions, including dark matter halos. Guo et al. (2011, hereafter G11) use a semi-analytic procedure to track the evolution of the galaxies that exist within these halos. Once subhalos are identified, their merger trees are derived. The evolution of these subhalos provides the base for the galaxy formation model. The models used by G11 build upon the work of De Lucia & Blaizot (2007), making improvements to the treatment of supernova (SN) feedback, reincorporation of ejected gas, galaxy sizes, the distinction between satellite and central galaxies, and the effect of the environment on galaxies. While semi-analytical models do not supply accurate details about individual galaxies, they are very useful for understanding general characteristics of large populations of galaxies. These semi-analytical models, when applied to the MRII simulation, are able to accurately reproduce observed characteristics of galaxy populations, e.g., the abundance and large-scale clustering of low-z galaxies, the Tully–Fisher relation, stellar mass and luminosity functions of low-z galaxies, the halo–galaxy-mass relationship, and the evolution of the cosmic star formation density. However, these models overproduce passive low-mass galaxies and fail to reproduce the observed abundances, clustering, and mass functions of high-redshift (z > 1.0) galaxies. In this paper, we will be comparing with X-ray observations of galaxies out to redshift 1.4, which is in a regime where the G11 model is still fairly accurate.

The result of G11's semi-analytic model is a catalog of the galaxy population at 60 different times between z ∼ 20 (about 13.4 billion years ago) and the present day. These catalogs include properties such as metallicity, stellar mass, bulge mass, the mass of hot and cold gas, rest-frame luminosity magnitudes, etc., for each of the galaxies in the simulation box, as a function of time.

To simulate the XRB populations of the galaxies from MRII, we use StarTrack, a current binary PS code that has been tested and calibrated using detailed binary star calculations and incorporates all the most important physical processes of binary evolution (Belczynski et al. 2002, 2008).

• 1.  
  The evolution of single stars and non-interacting binary components, from their birth, taken as the time of their initial emergence onto the main sequence, to compact remnant formation using evolutionary formulae of Hurley et al. (2000) modified as described in Belczynski et al. (2008). Various wind mass-loss rates and their effect on stellar evolution are also incorporated into the code and have been recently updated (Belczynski et al. 2010).
• 2.  
  The time evolution of orbital properties. StarTrack numerically integrates a set of four differential equations describing the evolution of orbital separation, eccentricity, and component spins, taking into account tidal interactions, magnetic braking, gravitational radiation, and stellar wind mass losses.
• 3.  
  All types of mass-transfer phases. This includes both stable and unstable mass-transfer processes, which are driven by either nuclear evolution or angular momentum loss. Unstable mass transfer is encountered most often as a direct consequence of rapid stellar expansion during nuclear evolution, but angular momentum loss may also lead to instability.
• 4.  
  SN explosions, which are treated by taking into account mass/angular momentum losses as well as SN asymmetries (through natal kicks to neutron stars (NS) and BHs).
• 5.  
  X-ray emission, which is tracked for accreting binaries with compact object primaries (both for wind-fed and Roche-lobe overflowing systems). The resulting X-ray luminosities are calculated from the secular averaged mass accretion rate, but are not calculated for unstable accretion phases because the timescales are very short.

The models in this paper include a recent revision of StarTrack that incorporates updated stellar winds and their recalibrated dependence on metallicity (Belczynski et al. 2010). However, two more recent upgrades have not been incorporated into these results, as the simulations were run long before the changes were made. The first update includes a revised NS and BH mass distribution based on fully consistent SN simulations (Belczynski et al. 2012; Fryer et al. 2012). The second, most recent upgrade improves upon the treatment of donor stars in common envelope (CE) events via usage of the actual value of the λ parameter, the measure of the central concentration of the donor and envelope binding energy, for which usually a constant value is assumed (Dominik et al. 2012).

Table 1 lists the input parameters of our PS models, which can be put into two categories. In the first category there are the parameters that correspond to the initial properties of the binary population. These values are relatively well constrained by the most recent observational surveys. Also in the group are stellar wind prescriptions and natal kick distributions, which can also be constrained by observations. In the second group there are the truly "free" parameters that correspond to poorly understood physical processes, which we are not able to model in detail. One of these truly "free" parameters is the CE efficiency (α_CE), which measures how efficiently orbital energy loss is transformed into thermal energy that will expel the donor's envelope during the CE phase. We note that in our calculations we combine α_CE and λ, the binding energy parameter described above, into one CE parameter. Whenever we mention the CE efficiency α_CE, we refer in practice to the product α_CE × λ, effectively treating α_CE × λ as a free parameter (see Belczynski et al. 2008 for details).

We create a grid of 192 PS models, a subset of those used in F13, run for 9 different metallicities and each simulating 5.12 × 10⁶ stars for 14 Gyr. In this grid, we varied all the parameters known from previous studies to affect the evolution of XRBs and the formation of compact objects in general (Belczynski et al. 2007, 2010; Fragos et al. 2008, 2010; Linden et al. 2009). Specifically, we vary the CE efficiency, initial binary mass ratio distribution, initial mass function (IMF), SN kicks for direct collapse (DC) BHs, and stellar wind strength. We also take into account the possibility of CE inspirals with Hertzsprung gap (HG) donors that could terminate binary evolution barring the subsequent XRB formation (Belczynski et al. 2007).

For all models, we assume a Maxwellian distribution of SN kicks given by Hobbs et al. (2005), with σ = 265 km s⁻¹. For compact objects formed with partial mass fallback, the natal kicks given by the Hobbs et al. (2005) distribution are decreased by a factor of (1 − f_fb), where f_fb is the fraction of the stellar envelope that falls back after the SN explosion. In our standard prescription, DC BHs, BHs formed with f_fb = 1, are given no natal kick. However, due to recent theoretical evidence that even the most massive stellar BHs have probably received small asymmetric kicks (Linden et al. 2010; Valsecchi et al. 2010), in some models used in this work we set a lower limit (0.1) on the amount by which the natal kicks may be decreased due to mass fallback, allowing for small natal kicks to be given to DC BHs.

The mass of the primary star in each binary is determined by the adopted IMF. It is important to note here that, because we sample the IMF with only the primary star, we are only sampling the high-mass end of the IMF because the primary stars that form XRBs must be massive enough to form a BH or NS. The mass of the secondary star is calculated using a distribution function for the binary mass ratio, q = M_secondary/M_primary. We vary the distribution of q between a flat distribution in the range q = 0–1 and a distribution that has 50% of the binaries follow a distribution with q = 0–1 and the other half follow a "twins" distribution, with q = 0.9–1.

The models and their numbers used in this work are the same as those used in F13, except here we exclude the models 97–192, which have all systems following the pure "twins" q-distribution. F13 show that models with the pure "twins" distribution are unphysical, as they fail to reproduce local populations of XRBs. Thus, we exclude them here.

We note that our PS code calculates the bolometric luminosity of each XRB based on the rate of mass transfer. In order to compare our model results with observed data sets we need to estimate the X-ray luminosity in a specific energy band, which in this study is the soft X-ray band of Chandra (0.5–2.0 keV). In order to calculate the bolometric correction, we used two sets of published X-ray spectra from Galactic NS and BH XRBs at different spectral states (Remillard & McClintock 2006; Wu et al. 2010). Following the same procedure outlined in F13, we derive to the bolometric correction factors for different types of XRBs and use these results to estimate the 0.5–2.0 keV X-ray luminosity of our modeled XRB population.

We also note that the PS models used here take into account only binary systems formed in the field and not those formed via dynamical interactions in dense stellar systems. Dynamical formation in globular clusters (GCs), for example, is a significant formation pathway for LMXBs in old, GC-rich elliptical galaxies (e.g., Humphrey & Buote 2008; Zhang et al. 2012). Further, the LMXB populations formed in GCs can be as much as 2–3 times more luminous than the field population in bright elliptical galaxies (Irwin 2005).

For each model and each metallicity value, we calculate the X-ray luminosity as a function of age for single bursts of star formation, where we also take into account the effect of transient XRBs (see Fragos et al. 2008, 2009). Taking into account the assumed IMF and initial mass ratio distribution, we normalize the total X-ray luminosity to a nominal population of 10¹⁰ M_☉. This quantity, L_{X, spec}(t)(erg s⁻¹/10¹⁰ M_☉), is the specific X-ray luminosity of a single-age stellar population as a function of its age. The specific X-ray luminosity coming from our PS models can be convolved with the star formation history (SFH) and metallicity evolution of a galaxy to calculate the total X-ray luminosity of its complex stellar population.

The MRII catalog created by G11 corresponds to 60 snapshots that span a redshift range from z = 0 to ∼20. For each galaxy, we can derive its complete progenitor tree. Each progenitor galaxy has a unique SFR and metallicity, so for every stellar population in a target galaxy we know during what time frame and at what metallicity that population was created. We then convolve the SFHs with L_{X, spec}(t) derived from the PS models for the appropriate metallicity values.

The SFRs given for each galaxy in the G11 catalog are averaged over the entire time step, Δt, between subsequent snapshots so that the total new stellar mass created in a given progenitor galaxy is M_new = SFR_progΔt. In order to account for the possibility of starbursts, we assume that all new stellar mass forms in a 20 Myr burst occurring at a random time between subsequent snapshots, t_i < t_burst < t_{i + 1}, where t_i is the timestamp associated with the snapshot of a given progenitor galaxy. This effect is important only for the HMXBs of young populations, and the effect on LMXBs is minimal since their evolution occurs on timescales much longer than the time steps between snapshots. A 20 Myr duration is reasonable, given that the most cited values for starburst durations are around 3–10 Myr (e.g., Thornley et al. 2000; Harris et al. 2004), though there is evidence for longer burst durations on the order of a few 10⁸ yr (e.g., McQuinn et al. 2010).

By summing the soft-band X-ray luminosities of all the stellar populations in a given galaxy, we derive the integrated X-ray luminosity from XRBs. The end result is a catalog of integrated XRB luminosities of galaxies within the MRII comoving volume from z = 0 to z ∼ 20.

When comparing our results to the observations of T&G08, we want to distinguish between early- and late-type galaxies. For their sample, T&G08 cross-correlate with other surveys to obtain optical counterparts for their X-ray-selected galaxies, which they used for classification.

The classification of a galaxy as early or late type can be based either on its morphology or its spectroscopic properties. In many observational surveys of distant galaxies (e.g., GOODS; Giavalisco et al. 2004), the morphologies of most galaxies that are observed cannot be determined due to inadequate spatial resolution, and therefore colors are used instead.

For color classification of our model galaxies as early or late type, we adopt the method developed by Bell et al. (2004). They showed that it is possible to define the population of early-type galaxies empirically by using the bimodality of the color distribution, which they studied out to z ∼ 1. The MRII database includes absolute rest-frame magnitudes in the Sloan Digital Sky Survey ugr filters, which can easily be transformed to the UBV filters used in Bell et al. (2004). The magnitudes include the effects of dust extinction. Following the Bell et al. (2004) prescription, we define early-type galaxies to be galaxies where 〈U − V〉 ⩾ 1.15 − 0.31z − 0.08*(M_V − 5log₁₀h + 20). Figure 1 shows plots of 〈U − V〉 versus M_V − 5log₁₀h for the MRII galaxies at various redshifts with the cutoff function overlaid on top and the bimodality is clearly present.

Zoom In Zoom Out Reset image size

If we instead define galaxy types based on morphology with late-type galaxies having (M_bulge/M_total) < 0.7 and early-type galaxies (M_bulge/M_total) > 0.7, we find that there is approximately a 1% and 0.5% contamination among the color defined late and early types, respectively. Thus, these methods give nearly identical results, but using colors to define morphology allows us to better simulate observations. The morphology method of classification, since it is independent of color, provides a check on our color classification. The fact that they both yield similar results is encouraging and indicates that the colors provided by the G11 catalog are able to yield accurate morphology classifications.

For our analysis, we only select galaxies with stellar masses greater than 10⁵ solar masses, as galaxies with mass less than this are very unlikely to have X-ray luminosities that are observationally relevant. For instance, the dwarf galaxies in the SINGS sample, with masses ∼10⁷ M_☉, generally have X-ray luminosities below 10³⁷ erg s⁻¹ (0.5–8 keV) with many having no binaries detected at all above 5 × 10³⁶ erg s⁻¹ (0.5–8 keV; L. Jenkins et al. 2013, in preparation). The G11 catalog contains several hundreds of these more massive galaxies at very high redshift (z = 19.9) and over 2 million at z = 0. It should also be noted that we assume that all of the galaxies in our sample are normal galaxies. This is a valid assumption because galaxies with bright active galactic nucleus (AGN) only constitute ∼2%–4% of all galaxies and, therefore, any selection effects on our data would be minimal (Xue et al. 2010; Haggard et al. 2010; Silverman et al. 2009). Lower luminosity AGN have been found in a much higher percentage (∼30%–40%) of LINER galaxies (Ho et al. 1997). However, since LINER galaxies themselves make up only one-fifth to one-third of all galaxies (Ho et al. 1997), this effect is also rather minimal. Additionally, since this is for lower luminosity AGN, it is likely that our more luminous galaxies would still be dominated by hot gas and XRB emission and would be classified as a normal galaxy. For example, Flohic et al. (2006) find that AGN in their sample of LINER galaxies contribute only 60% of the 0.5–10 keV luminosity when only considering the central regions of the galaxies.

To compare our results with observations, we derive the XLF for our galaxies by calculating the number density of galaxies versus their integrated X-ray luminosity. The simulation data are all within a single comoving volume of constant size. For the time slice represented in each snapshot, ϕ(L) is defined as

$$\begin{equation} \phi (L) = \frac{N(L_{{\rm min}},L_{{\rm max}})}{V_{{\rm MRII}}\delta _{L}}. \end{equation} \tag{ 1 }$$

Here, L_min, L_max are the bin limits, V_MRII is the volume of the MRII simulation, which is (100 Mpc h⁻¹)³, and δ_L is the size of the luminosity bin in log space, i.e., δ_L = log₁₀(L_max/L_min). Herein lies a major difference between the theoretical and observational luminosity functions. Observational surveys study a range of redshifts within a light cone. An entire volume of space cannot be observed at a constant redshift so a range of redshifts is explored. Thus, when calculating ϕ(L), observers such as T&G08 use methods like the one found in Page & Carrera (2000) which uses the following definition:

$$\begin{equation} \phi (L,z) = \frac{N}{\int _{L_{{\rm min}}}^{L_{{\rm max}}} \int _{z_{{\rm min}}(L)}^{z_{{\rm max}}(L)} \frac{dV}{dz}\,dz\,dL}. \end{equation} \tag{ 2 }$$

Here, dV/dz represents the rate of change of the survey volume with respect to redshift and z_max(L), z_min(L) are the redshift ranges for a source as a function of luminosity such that it stays within the flux limits of the survey and within the redshift interval. Our simulated galaxies, on the other hand, exist within a comoving volume and our snapshots capture all galaxies that exist at a given redshift. In order to compensate for this difference, we adopt similar redshift intervals as used by T&G08 and calculate the XLF for each redshift individually using Equation (1). This gives us ϕ(L, z). Then, for each luminosity bin centered at L, we take the average value of ϕ(L, z) over all z in the interval, giving us an estimate for ϕ(L) for that bin in that redshift interval.

In addition to XRBs, the hot ISM in a galaxy can have a significant contribution to its overall X-ray luminosity. T&G08 do not distinguish between emission from XRBs and hot gas in their analysis of X-ray bright galaxies. Their analysis is done in the soft X-ray band, so emission from the hot ISM becomes important and needs to be taken into account when calculating the total integrated luminosities of our galaxies.

We use observationally derived relations to estimate the X-ray luminosity of hot gas in early-type galaxies from their K-band luminosity. It has been shown that there is a power-law relationship between the X-ray luminosity of the hot ISM in early-type galaxies and both the K-band luminosity (L_X∝L^1.935_k) of the galaxy and the temperature of its hot gas (T∝L^0.214_X; Boroson et al. 2011). The K-band luminosities of MRII early-type galaxies are easily calculated from mass and age using synthetic stellar population models (Bertelli et al. 2008). With these relations, we estimate both the full-band X-ray luminosity and the temperature of the hot gas in each early-type galaxy. The spectrum of hot diffuse gas is assumed to be that of a collisionally ionized diffuse gas as calculated by the APEC XSPEC model and the ATOMDB code (Foster et al. 2012). The gas temperature estimate from the empirical relations is used as input to the APEC model in order to calculate the luminosity of the hot gas in the soft X-ray band.

For late-type galaxies, we estimate the hot gas X-ray luminosity based on the SFRs given in the G11 catalog and the power-law relationship catalog between the soft-band (0.5–2 keV) X-ray luminosity of the hot ISM and the SFR for late-type galaxies (L_X∝SFR^1.07; Strickland et al. 2004a, 2004b).

In summary, we estimate the total X-ray emission from hot gas in all of the galaxies in the G11 catalog. We add those values to each galaxy's integrated X-ray emission from XRBs, calculated using StarTrack, to obtain the total integrated X-ray luminosity of each galaxy. We find that on average XRB emission contributes to ∼50%–60% of the 0.5–2 keV emission from bright (L_X > 10³⁸) late-type galaxies and ∼40% of the 0.5–2 keV emission from bright early-type galaxies for our best-fitting model (205). Hence, we find that hot gas emission has an appreciable effect on the galaxy XLFs. See Section 4 for more details on our results.

From deep Chandra survey observations, T&G08 present early- and late-type galaxy counts in several luminosity bins and two redshift intervals, 0 ⩽ z < 0.4 and 0.4 ⩽ z < 1.4. They also provide total galaxy counts split into three redshift intervals, but for our analysis we will focus on the early- and late-type counts. Associated with each count is a survey volume that depends on both luminosity and redshift. Let the set of counts be d = {d_{i, j}|i = 1, ...N, j = 1, ...M} and associated volumes be V_{i, j}, where i ranges over the N = 5 luminosity bins and j ranges over the M = 2 redshift bins. T&G08 assume that the d_{i, j} are subject to Poisson statistical errors.

Similarly, for a particular choice of parameters, θ, (for example, see Table 2) our model produces a set of counts n = {n_{i, j}|i = 1, ...N, j = 1, ...M} for galaxies in each luminosity and redshift bin in our (100 Mpc h⁻¹)³ model volume. We assume that the counts n are drawn from a Poisson distribution with (unknown) means λ = {λ_{i, j}|i = 1, ...N, j = 1, ...M} (note that this λ is separate from the one used before to describe the CE efficiency parameter). Because we only observe one particular set of counts, n, we do not measure the rates λ implied by our model directly, but instead must treat λ as a nuisance parameter whose distribution under the observed n must be integrated over.

Zoom In Zoom Out Reset image size

Bayes' rule relates the posterior probability of model parameters θ, p(θ|d), to the likelihood of the data under the model, p(d|θ), the prior probability of the model parameters before the data have been observed, p(θ), and a normalizing constant, p(d), called the evidence, that is independent of θ via

$$\begin{equation} p(\theta | d) = \frac{p(d|\theta) p(\theta)}{p(d)}. \end{equation} \tag{ 3 }$$

Writing the likelihood in terms of the (unknown) true mean λ implied by the model, we have

$$\begin{equation} p(d|\theta) = \prod _{i,j} \int d\lambda _{i,j} \, p(d_{i,j} | \lambda _{i,j}) p(\lambda _{i,j} | \theta ), \end{equation} \tag{ 4 }$$

where

$$\begin{equation} p(d_{i,j} | \lambda _{i,j} ) = \frac{(v_{i,j}\lambda _{i,j})^{d_{i,j}}}{d_{i,j}!} \exp\; (- v_{i,j} \lambda _{i,j} ) \end{equation} \tag{ 5 }$$

is the Poisson probability of drawing d_{i, j} counts in a volume V_{i, j} = v_{i, j}(100 Mpc h⁻¹)³ when the underlying rate is λ_{i, j} per (100 Mpc h⁻¹)³.

The distribution of the underlying rates implied by our model, p(λ_{i, j}|θ), must be estimated from the observed n_{i, j}. Applying Bayes' rule again, we have

$$\begin{equation} p(\lambda _{i,j} | \theta ) = p(\lambda _{i,j} | n_{i,j}(\theta) ) = \frac{p(n_{i,j} | \lambda _{i,j}) p(\lambda _{i,j} )}{p(n_{i,j})}. \end{equation} \tag{ 6 }$$

The counts observed in a model with underlying rate λ_{i, j} are Poisson distributed, so

$$\begin{equation} p(n_{i,j} | \lambda _{i,j} ) = \frac{\lambda _{i,j}^{n_{i,j}}}{n_{i,j}!} \exp\; (- \lambda _{i,j}). \end{equation} \tag{ 7 }$$

We choose a Jeffreys prior¹⁰ on λ_{i, j},

$$\begin{equation} p(\lambda _{i,j} ) = \frac{1}{\sqrt{\lambda _{i,j}}}, \end{equation} \tag{ 8 }$$

whence

$$\begin{equation} p(\lambda _{i,j} | \theta) = \frac{\lambda _{i,j}^{n_{i,j}-\frac{1}{2}}}{\Gamma (n_{i,j} + \frac{1}{2} )} \exp\; (- \lambda _{i,j} ). \end{equation} \tag{ 9 }$$

Combining Equations (9) and (5), we find that

$$\begin{equation} p(d | \theta) = \prod _{i,j} \frac{v_{i,j}^{d_{i,j}}\Gamma \left(\frac{1}{2} + d_{i,j} + n_{i,j}\right)}{(1+v_{i,j})^{\frac{1}{2} + d_{i,j} + n_{i,j}} d_{i,j}! \Gamma \left(\frac{1}{2} + n_{i,j} \right)}. \end{equation} \tag{ 10 }$$

We choose a flat prior on the model parameters, θ, so that the posterior is proportional to the likelihood in Equation (10):

$$\begin{equation} p(\theta | d) \propto p(d | \theta). \end{equation} \tag{ 11 }$$

F13 use the same PS models used in this work (with the inclusion of pure "twins" models) to study the evolution of the overall population of XRBs in the universe. They compare with X-ray observations of local galaxies that give estimates of the specific X-ray luminosity of XRBs in the local universe (Lehmer et al. 2010; Boroson et al. 2011; Mineo et al. 2012). They calculate the likelihood of each model based on these data and find the six highest likelihood models to be, in order of likelihood, 245, 229, 269, 205, 249, and 273 (see Table 2 for model parameters and likelihood values from this study). Figure 3 compares the six highest likelihood models from this work with those of F13. Four out of our top six models are also among the six highest likelihood models from F13, so it is no surprise that the region bounded by the models in this work is very similar to that bounded by the models from F13.

Zoom In Zoom Out Reset image size

Figure 4 shows XLFs for our highest likelihood model (205), with and without hot gas emission, plotted against the data from T&G08. Figure 5 plots XLFs similar to Figure 4, but for the highest likelihood model from F13 (245). These plots show that our models are able to reproduce the redshift evolution of the observed XLFs. Consistent with the analysis of T&G08, the XLF evolution is driven almost entirely by late-type galaxies. Our models also reproduce the shape of the early-type XLF and the normalization of the late-type XLF, though they drastically underproduce bright early-type galaxies and they fail to reproduce the shape of the bright (L_X > 10⁴¹ erg s⁻¹) end of the late-type XLF.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Figures 4 and 5 also show that hot gas can have a large effect on the shape and normalization of the XLF, showing that hot gas emission dominates the integrated X-ray luminosity of the brightest galaxies in our sample. Adding in the hot gas emission suppresses the redshift evolution for the early-type galaxy XLF. For galaxies with L_X > 10⁴⁰ erg s⁻¹ at low (z < 0.8) redshift, emission from XRBs accounts on average for only 1%–5% and ∼15% of early- and late-type galaxy emission, respectively. However, as we will discuss in Section 4.3, the XLFs are still rather sensitive to changes in our model parameters, as seen in the varying likelihood values shown in Figure 2. The important role of hot gas emission on the XLF means that our simplistic prescriptions for hot gas emission add a great deal of uncertainty to our models and could be a major source of the discrepancies between the models and observations, particularly for early-type galaxies. While our method is motivated by observations, it does not take into account the internal characteristics of the gas that contribute to its emission, such as density and metallicity. Further, the relations used for early-type galaxies were derived only from low-redshift sources, which may not be accurate for the high-redshift galaxies studied here.

In addition, the G11 semi-analytic model underproduces massive galaxies at high redshift. This will lead to less large elliptical galaxies, which could explain part of our discrepancy at higher redshift.

Another aspect of our models that can account for the underproduction of bright early-type galaxies is that our PS models only take into account LMXBs formed in the field and not those formed dynamically in GCs. Dynamically formed LMXBs are believed to play a significant role in old, massive, GC-rich elliptical galaxies. (e.g., Humphrey & Buote 2008; Zhang et al. 2012). These LMXB populations can have a significant contribution to the integrated X-ray luminosity of bright early-type galaxies, as they can make up over half of the total number of LMXBs in a galaxy (Irwin 2005). So, including dynamically formed LMXBs in our models could increase the number of bright LMXBs, and therefore the total LMXB luminosity, in early-type galaxies by a factor of ∼3. Changing the q-distribution in model 245 from a 50–50 to a flat distribution is a good proxy for this effect because it increases the LMXB population without changing the distribution of their physical properties (see Section 4.3). F13 show that doing this increases the total luminosity from LMXBs at all redshifts by a factor of two. We find that changing to a flat q-distribution increases the low-luminosity end of the early-type XLF by ∼0.3 dex (see Figure 7). The effect of including dynamically formed LMXBs could have a similar but greater effect, bringing the low-luminosity end of the XLF closer to observations, but having little effect on higher luminosity galaxies. A more detailed calculation will require information on the GC population of each galaxy, which is not included in the G11 catalog.

For younger, star-forming galaxies, our models also have only a very basic formula to simulate starburst activity, which can occur, e.g., due to galaxy mergers. This would have a significant effect on the HMXB populations present in late-type galaxies, and the effect would not necessarily be constant with redshift, as merger rates may evolve with time (e.g., Conselice 2006). Thus, a more detailed SFH is needed to more accurately model HMXB populations of late-type galaxies.

In addition, the higher end of the observed late-type galaxy XLF is more at risk from AGN contamination, even with the efforts of T&G08 to minimize this effect. Since the observations of T&G08, the depth of the X-ray data, combined with better multiwavelength data, has allowed for more accurate classifications of the X-ray sources. Of the 56 1 Ms CDF-S sources used in T&G08, we find 53 counterparts with 4 Ms exposure using a matching radius of 2.5 arcsec. The missing three sources may have been false-positive sources in the 1 Ms data. Of these 53 sources, we find that 25 of them are classified as normal galaxies and 28 of them as AGN according to the 6 criteria highlighted in Section 3.1 of Lehmer et al. (2012). Therefore, it is possible that the T&G08 data points will be lowered by ∼0.3 dex. However, it is difficult to know in detail how this affects the TG08 luminosity functions and recomputing the luminosity functions is beyond the scope of this work.

Figure 6 plots the X-ray luminosity density from normal galaxies as derived from our highest likelihood model, 205. The overall evolution (black line in Figure 6) is very similar to that of the observed SFH of the universe. It is also similar to the evolution of the specific XRB X-ray luminosity of the universe predicted in F13, despite the inclusion here of hot gas emission. This is evidence that XRBs drive the overall evolution of the normal galaxy X-ray luminosities out to at least z = 4. However, our predicted X-ray luminosity density reaches a maximum at z ∼ 2.5, which is lower compared with the XRB models in F13 that reach a maximum at z ∼ 3. This can be attributed to the inclusion of hot gas emission in our models, which has already been shown to have a noticeable effect on the shape and evolution of our XLFs, particularly for early-type galaxies.

Zoom In Zoom Out Reset image size

Splitting the galaxies into three luminosity bins, we find that the evolution of low (10³⁹ < L_X < 10⁴⁰ erg s⁻¹) luminosity galaxy emission is small compared with the evolution for higher luminosity galaxies, which varies by an order of magnitude on the range of z = 0–4. The most luminous galaxies (L_X > 10⁴¹ erg s⁻¹) reach a maximum around z = 3, which also approximately corresponds to the time of maximum SFR density in the universe. Galaxies in the range of 10⁴⁰ < L_X < 10⁴¹ erg s⁻¹ reach a maximum around z = 2. In the local universe, the low-luminosity galaxies dominate the normal galaxy X-ray emission. We do not go to higher redshift here because our hot gas emission prescription relies on galaxy morphology, which becomes harder to classify at higher redshifts (van den Bergh 2002).

Figure 7 shows XLFs for different models compared to model 245. Each model is chosen to encapsulate the effect that each parameter has on the shape of the XLF. Several parameters have significant effects on the shape of the XLFs, while others have only minimal effects. Model 245 was chosen because it has both a high likelihood and is more sensitive to certain changes in parameters, thus better illustrating the different effects on our XLFs.

Zoom In Zoom Out Reset image size

The CE efficiency parameter (α_CE) dictates how efficiently orbital energy is converted to thermal energy that will expel the envelope. A lower efficiency means that it will take more orbital energy to expel the envelope. This parameter mainly affects LMXBs, as most LMXBs formed in the field must go through a CE phase. The CE phase plays an important role in making the orbit close enough to allow for Roche-lobe overflow (RLO), but a lower CE efficiency leads to even more orbital decay and a higher rate of mergers, overall decreasing the rate of LMXB formation. This effect can be seen by comparing models 245 and 248. HMXBs are not as strongly affected by changes in α_CE, as they have other formation channels available that do not involve CE phases (Linden et al. 2010; Valsecchi et al. 2010). This parameter mainly affects the lower luminosity end of the XLF, where a higher CE efficiency increases the number of bright galaxies due to an increased LMXB population.

Wind mass-loss rates affect the evolution of high-mass stars in two major competing ways. Higher wind mass-loss rates will increase the accretion rates of wind-fed HMXBs, increasing their luminosity. On the contrary, lower winds will result in a lower overall mass loss of the primary star and hence increase the formation rate of massive BHs. BH-XRBs tend to be more luminous than XRBs with NS accretors. This is because, on one hand, they can form stable RLO XRB systems with massive companions and, on the other hand, BHs drive higher accretion rates due to their higher mass and, therefore, higher Eddington limits. In this way, weaker stellar winds can increase the luminosities of both LMXB and HMXB populations. Comparing models 245 and 261, we see that weaker stellar winds increase the number of bright early- and late-type galaxies, so the latter effect is dominant. While lower stellar winds will help our highest likelihood model match observations of early-type galaxies, it would overproduce bright late-type galaxies.

Changing the initial binary mass ratio between a flat distribution and a 50% twins and 50% flat distribution affects both the HMXB and LMXB populations. As stated earlier, the binary mass ratio affects the secondary star in the system, which will eventually become the donor star in most cases. All XRBs are accreting mass onto a compact object, which can only come from a high-mass progenitor. Most LMXBs require a high initial mass ratio in order to ensure a high-mass primary star that will evolve into a BH or NS with a lower mass companion. The "50–50" distribution adopted here forces mass ratios close to 1 for half the binary population. This will decrease the LMXB population while increasing the HMXB population, as HMXBs require ratios closer to 1. Changing from "50–50" to a flat distribution (comparing models 245 and 53) increases the lower luminosity end of the early- and late-type galaxy XLFs and slightly decreases the number of very bright late-type galaxies galaxies.

Allowing small natal kicks for DC BHs affects the HMXB and LMXB populations in two competing ways. On one hand, natal kicks can enhance the formation of RLO-HMXBs (Linden et al. 2010). On the other hand, natal kicks inject energy into the binary system and could result in the widening or the complete disruption of the system, thereby decreasing the formation of HMXB and LMXB with BH accretors. Linden et al. (2009) find that imparting small natal kicks to DC BHs is necessary in order to reproduce the lack of observed wide orbit BH-XRBs. Comparing models 245 and 197, not allowing natal kicks increases the high-luminosity end of the late-type galaxy XLF and has little effect on the early-type galaxy population.

Within our grid of PS models, we also have the IMF power-law exponent as a free parameter, allowing it to be either −2.35 or −2.7. It is instructive to note that the IMF referred to in this work represents the integrated galaxy IMF (Weidner & Kroupa 2005; Kroupa & Weidner 2003) and that this work only probes the high-mass region of the IMF, since the primary stars that are created via sampling the IMF must be massive enough to form a BH or NS. The slope of the power law at the high-mass end of the IMF affects the population of XRBs in a way similar to stellar winds, in the sense that a flatter IMF will produce relatively more massive BHs compared to a steeper one. A flatter IMF will result in more bright LMXBs and HMXBs. This effect can be seen by comparing models 245 and 241. As expected, a flatter IMF results in a higher number of bright galaxies, though the very high luminosity end of the early-type galaxy XLF is not strongly affected by this parameter. This is not surprising, as Figure 4 shows that hot gas emission dominates the high-luminosity end of the XLF for early-type galaxies.

Finally, we found that allowing or not for all possible outcomes in CE phases with donor stars in the HG has only a small effect on the galaxy XLFs, slightly increasing the number of bright late-type galaxies (compare models 245 and 269). Thus, although this parameter affects the shape of the XLF of individual XRBs in a galaxy (Luo et al. 2012), it has a negligible effect on the integrated X-ray luminosity of a galaxy.

The models that we find to agree best with observations have DC BH natal kicks, low CE efficiency, steep IMFs, η_wind = 1.0–2.0, and a 50–50 mass ratio distribution. As outlined above, these parameters make for a limited population of both LMXBs and BH-XRBs. The HMXB population is limited by the steeper IMF, BH natal kicks, and higher winds, but also strengthened by the 50–50 mass ratio distribution. However, the latter effect on the XLF is much weaker, as shown in Figure 7.