THE EVOLUTION OF NORMAL GALAXY X-RAY EMISSION THROUGH COSMIC HISTORY: CONSTRAINTS FROM THE 6 MS CHANDRA DEEP FIELD-SOUTH

Chandra studies of local galaxies have yielded remarkable insight into the formation and evolution of populations of X-ray binaries (XRBs; see, e.g., Fabbiano 2006 for a review). Expansive multiwavelength (e.g., from GALEX, Herschel, Hubble, and Spitzer) and Chandra observations of local star-forming and passive galaxy samples have constrained basic scaling relations between the X-ray emission from the high-mass XRB (HMXB) and low-mass XRB (LMXB) populations with star-formation rate (SFR) and stellar mass (${M}_{\star }$), respectively (see, e.g., Grimm et al. 2003; Ranalli et al. 2003; Colbert et al. 2004; Gilfanov 2004; Hornschemeier et al. 2005; Lehmer et al. 2010; Boroson et al. 2011; Mineo et al. 2012a; Zhang et al. 2012); hereafter, we refer to these scaling relations as the ${L}_{{\rm{X}}}$(HMXB)/SFR and ${L}_{{\rm{X}}}$(LMXB)/${M}_{\star }$ relations, respectively. However, the scatters in these relations are factors of ≈2–5 times larger than the expected variations due to measurement errors and statistical noise (e.g., Gilfanov 2004; Hornschemeier et al. 2005; Mineo et al. 2012a), thus indicating that real physical variations in the galaxy population (e.g., stellar ages, metallicities, and star formation histories) likely have a significant influence on XRB formation.

Recently, Fragos et al. (2013a; hereafter, F13a) used measurements of the local scaling relations to constrain theoretical XRB population-synthesis models. The F13a framework was developed using jointly the Millenium II cosmological simulation (Guo et al. 2011) and the Startrack XRB population-synthesis code (e.g., Belczynski et al. 2002, 2008) to track the evolution of the stellar populations in the universe and predict the X-ray emission associated with the underlying XRB populations, respectively. The F13a models follow the evolution of XRB populations and their parent stellar populations throughout the history of the universe, spanning $z\quad \approx$ 20 to the present day, and provide predictions for how the X-ray scaling relations (i.e., ${L}_{{\rm{X}}}$(HMXB)/SFR and ${L}_{{\rm{X}}}$(LMXB)/${M}_{\star }$) evolve with redshift. From this work, F13a identified a "best fit" theoretical model that simultaneously fits well both the observed ${L}_{{\rm{X}}}$(HMXB)/SFR and ${L}_{{\rm{X}}}$(LMXB)/${M}_{\star }$ scaling relations at z = 0. Subsequent observational tests have shown that the F13a best-fit model provides reasonable predictions for (1) the XRB luminosity functions of a sample of nearby galaxies that had simple star-formation history estimates from the literature (Tzanavaris et al. 2013); (2) the metallicity dependence of the ${L}_{{\rm{X}}}$(HMXB)/SFR relation for powerful star-forming galaxies (Basu-Zych et al. 2013a; Brorby et al. 2016); (3) the stellar-age dependence of the ${L}_{{\rm{X}}}$(LMXB)/${M}_{\star }$ relation for early-type galaxies (Lehmer et al. 2014; however, see Boroson et al. 2011 and Zhang et al. 2012); (4) the redshift evolution out to $z\approx 1.5$ of the normal galaxy X-ray luminosity functions in extragalactic Chandra surveys (Tremmel et al. 2013); and (5) the redshift evolution of the total ${L}_{{\rm{X}}}$/SFR relation (i.e., using the summed HMXB plus LMXB emission) for star-forming galaxies out to $z\approx 4$ (e.g., Basu-Zych et al. 2013b).

The F13a theoretical modeling framework, as well as the broad observational testing of its predictions, represent major steps forward in our understanding of how XRBs form and evolve along with their parent stellar populations. Within the F13a framework, the most fundamental predictions include prescriptions for how the ${L}_{{\rm{X}}}$(HMXB)/SFR and ${L}_{{\rm{X}}}$(LMXB)/${M}_{\star }$ scaling relations evolve as a function of redshift (see Figure 5 of F13a). Due to sensitivity and angular-resolution limitations, it is not possible to detect complete and representative populations of cosmologically distant galaxies and separate spatially their HMXB and LMXB contributions. However, with deep (≳1 Ms) Chandra exposures and new multiwavelength databases, several extragalactic surveys now have data sufficient to isolate large populations of galaxies, measure their global physical properties (e.g., SFR and ${M}_{\star }$), and study their population-averaged X-ray emission via stacking techniques (see, e.g., Hornschemeier et al. 2002; Laird et al. 2006; Lehmer et al. 2007, 2008; Cowie et al. 2012; Basu-Zych et al. 2013b). With these advances, we can now conduct the most robust tests to date of the F13a model predictions by examining the XRB emission of galaxies dependence on SFR, ${M}_{\star }$, and redshift.

The Chandra Deep Field-South (CDF-S) survey is the deepest X-ray survey yet conducted. In this paper, we utilize data products based on the first ≈6 Ms of data, which were produced following the procedures outlined for the ≈4 Ms exposure in Xue et al. (2011). An additional ≈1 Ms of data will be added to the CDF-S, eventually bringing the total exposure to ≈7 Ms; these results will be presented in Luo et al. (2016, in preparation). In the ≈6 Ms exposure, 919 highly reliable X-ray sources are detected to an ultimate 0.5–2 keV flux limit of ≈7 × 10⁻¹⁸ erg cm⁻² s⁻¹, including 650 AGN candidates, 257 normal galaxy candidates, and 12 Galactic stars (see Section 3 for classification details). For comparison, the ≈4 Ms CDF-S catalog contained 740 sources down to an ultimate 0.5–2 keV flux limit of ≈10⁻¹⁷ erg cm⁻² s⁻¹, of which 568, 162, and 10 were classified as AGN, normal galaxies, and Galactic stars. In the most sensitive regions of the survey field, the 0.5–2 keV detected normal galaxies rival or exceed the AGN in terms of source density (see, e.g., the Lehmer et al. 2012 analysis of the 4 Ms data).

Source catalogs based on optical/near-IR imaging contain ≈25,000 galaxies within ≈7 arcmin of the CDF-S Chandra aimpoint (e.g., Luo et al. 2011; Xue et al. 2012), indicating that only a small fraction (≲1%) of the known normal galaxy population is currently detected in X-ray emission. In this paper, we utilize X-ray stacking analyses of the galaxy populations within the CDF-S, divided into redshift intervals and subsamples of specific SFR, sSFR ≡  SFR/${M}_{\star }$, which is an indicator of the ratio of HMXB-to-LMXB emission. These measurements provide the first accounting of both HMXBs and LMXBs at $z\gt 0$ to simultaneously constrain the evolution of the ${L}_{{\rm{X}}}$(HMXB)/SFR and ${L}_{{\rm{X}}}$(LMXB)/${M}_{\star }$ scaling relations and provide the most powerful and robust test of the F13a model predictions.

This paper is organized as follows: In Section 2, we define our galaxy sample and describe our methods for measuring SFR and ${M}_{\star }$ values for the galaxies. In Section 3, we discuss the X-ray properties of galaxies and scaling relations of our sample that are individually detected in the ≈6 Ms images. In Section 4, we detail our stacking procedure, and in Section 5 we define galaxy subsamples to be stacked. Results from our stacking analyses, including characterizations of the evolution of X-ray scaling relations, are presented in Section 6. Finally, in Section 7, we interpret our results in the context of XRB population-synthesis models, construct models of the evolution of the X-ray emissivity of the universe, and estimate the cosmic X-ray background contributions from normal galaxies.

Throughout this paper, we make use of the main point-source catalog and data products for the ≈6 Ms CDF-S as will be outlined in Luo et al. (2016, in preparation). The Luo et al. (2016, in preparation) procedure is identical in nature to that presented for the ≈4 Ms CDF-S in Xue et al. (2011). The Galactic column density for the CDF-S is $8.8\times {10}^{19}$ cm⁻² (Stark et al. 1992). All of the X-ray fluxes and luminosities quoted throughout this paper have been corrected for Galactic absorption. Estimates of ${M}_{\star }$ and SFR presented throughout this paper have been derived assuming a Kroupa (2001) initial mass function (IMF); when making comparisons with other studies, we have adjusted all values to correspond to our adopted IMF. Values of H₀ = 70 km s⁻¹ Mpc⁻¹, ${{\rm{\Omega }}}_{{\rm{M}}}$ = 0.3, and ${{\rm{\Omega }}}_{{\rm{\Lambda }}}$ = 0.7 are adopted throughout this paper.

We began with an initial sample of 32,508 galaxies in the Great Observatories Origins Deep Survey South (GOODS-S) footprint as presented in Section 2 of Xue et al. (2012; hereafter X12; see also Luo et al. 2011). The X12 galaxy sample was selected using the HST F850LP photometric data from Dahlen et al. (2010), and contains objects down to a 5σ limiting magnitude of ${z}_{850}\approx 28.1$. We cut our initial sample to the 24,941 objects that were within 7${}^{\prime }$ of the mean ≈6 Ms CDF-S aimpoint, a region where the Chandra point-spread function (PSF) is sharpest and the corresponding X-ray sensitivity is highest. Hereafter, we refer to the resulting sample as our base sample. Of the 24,941 objects in our base sample, 1,124 (4.5%) have secure spectroscopic redshifts. To measure redshifts for the full base sample, X12 used 12–15 band PSF-matched photometric data (including upper limits) covering 0.3–8 μm, and performed spectral energy distribution (SED) fitting using the Zurich Extragalactic Bayesian Redshift Analyzer (ZEBRA; Feldmann et al. 2006). The full redshift range of the sample spans z = 0.01–7. A variety of tests indicated that the resulting photometric redshift estimates are of high quality (normalized median absolute deviation ${\sigma }_{{\rm{NMAD}}}=0.043$) over the range of ${z}_{850}\quad \approx$ 16–26, with a low outlier fraction (fraction of sources with $| {\rm{\Delta }}z| /(1+{z}_{{\rm{spec}}})\gt 0.15$) of ≈7%.

For each of the 24,941 galaxies in the base sample, X12 adopted the best-available redshift and best-fitting SED for that redshift to estimate galaxy absolute magnitudes. In this procedure secure spectroscopic redshifts were used when available and photometric redshifts were used otherwise. For a given galaxy, absolute magnitudes were computed for rest-frame B, V, R, I, J, H, and K bands using the best-fit SED. Stellar masses, ${M}_{\star }$, were calculated following the prescription in Zibetti et al. (2009), using rest-frame B − V colors and K-band luminosities, L_K, along with the following equation:

$$\begin{eqnarray}&&\mathrm{log}{M}_{\star }=\mathrm{log}{L}_{K}+1.176(B-V)-1.39,\end{eqnarray} \tag{ 1 }$$

where ${M}_{\star }$ and L_K are in solar units. The numerical constants in Equation (1) were supplied in Table B1 of Zibetti et al. (2009), which provides mass-to-light ratio estimates for a variety of rest-frame bands and colors. The color term (B − V) in Equation (1) accounts for variations in stellar age (or star-formation history). Therefore, Equation (1) is applicable to the full range of stellar ages and galaxy types. Uncertainties in the above calibration are at the ≈0.15 dex level, and we adopt this uncertainty throughout the rest of this paper. As we discuss below, the stellar masses derived from this method for a large subset of our galaxies are in good agreement with those derived from more detailed SED fitting techniques.

We calculated SFRs for the galaxies in our sample using UV and far-IR tracers. The UV emission is assumed to arise from the young stellar populations, and the portion of the UV light absorbed and re-radiated by dust is assumed to produce the far-IR emission. The majority of the galaxies in our sample have photometric estimates of the continuum emission near rest-frame 2800 Å, thereby allowing accurate estimates of rest-frame 2800 Å monochromatic luminosities, $\nu {l}_{\nu }(2800$ Å) from the best-fit ZEBRA-based SED templates used by Xue et al. (2012) (see discussion above). To determine the broad-band portion of the SED that is associated with only the young UV-emitting population, we constructed an SED that assumes a constant SF history that spanned 100 Myr. For each galaxy, we used this SED to calculate the total observed UV emission for the young population following ${L}_{{\rm{UV,}}\;{\rm{obs}}}^{{\rm{young}}}=C({A}_{V})\nu {l}_{\nu }(2800$ Å), where $C({A}_{V})$ scales $\nu {l}_{\nu }(2800$ Å) to the bolometric luminosity, given an attenuation, A_V, and a Calzetti et al. (2000) extinction curve. In our case, A_V was calculated for all galaxies by X12 in the SED-fitting process. The observed UV luminosity alone provides a poor tracer of the SFRs, since the intrinsic emission from young stars is, in most cases, attenuated by several multiplicative factors due to dust extinction. Therefore, whenever possible, we also measured galaxy IR luminosities (8–1000 μm; ${L}_{{\rm{IR}}}$) to estimate directly the levels of UV power obscured and reprocessed by dust (see, e.g., Kennicutt 1998 & Kennicutt & Evans 2012 for reviews).

To determine ${L}_{{\rm{IR}}}$, we cross-matched our galaxy catalog positions with the far-IR 24–160 μm  GOODS-Herschel catalogs.²¹ These catalogs include data from Spitzer-MIPS 24 μm, as well as the Herschel-PACS 100 μm  and 160 μm  fluxes of the 24 μm  sources; there are no unique sources at 100 μm  and 160 μm  that are not detected by Spitzer  24 μm  imaging (Elbaz et al. 2011). We identified 931 far-IR counterparts (using a constant matching radius of 1'') to the 24,941 galaxies in our base sample (i.e., 3.8%). Using the IR SED template presented by Chary & Elbaz (2001), we estimated ${L}_{{\rm{IR}}}$ for a given IR-detected galaxy by (1) normalizing the template SED to the l_ν value derived from the flux of the galaxy in the reddest far-IR detection band, and (2) integrating the normalized template SED across the 8–1000 μm  band. We tested for systematic differences between ${L}_{{\rm{IR}}}$ values derived using one band versus another, but found no significant offsets or any trends with redshift. The corresponding mean ratios and 1σ scatter between bands were $\mathrm{log}[{L}_{{\rm{IR}}}$(160μm )/${L}_{{\rm{IR}}}$(100 μm)] = 0.08 ± 0.16 and $\mathrm{log}[{L}_{{\rm{IR}}}$(100 μm)/${L}_{{\rm{IR}}}$(24 μm)] = −0.07 ± 0.20, implying the expected error on ${L}_{{\rm{IR}}}$ estimates is on the order of ≈0.2 dex. Of the 931 far-IR detected galaxies in our sample, we derived ${L}_{{\rm{IR}}}$ using 250, 172, and 509 detections at 24 μm, 100 μm, and160 μm, respectively.

For the 931 galaxies with both UV and far-IR detections, we utilized the methods described in Section 3.2 of Bell et al. (2005) and estimated total galaxy SFRs using the following equation:

$$\begin{eqnarray*}&&{\rm{SFR}}({M}_{\odot }\;{{\rm{yr}}}^{-1})=9.8\times {10}^{-11}{L}_{{\rm{bol}}}^{{\rm{young}}}\end{eqnarray*}$$

$$\begin{eqnarray}&&{L}_{{\rm{bol}}}^{{\rm{young}}}\approx {L}_{{\rm{UV,}}\;{\rm{obs}}}^{{\rm{young}}}+{L}_{{\rm{IR}}},\end{eqnarray} \tag{ 2 }$$

where all luminosities are expressed in units of solar bolometric luminosity (${L}_{\odot }=3.9\times {10}^{33}$erg s⁻¹) and ${L}_{{\rm{bol}}}^{{\rm{young}}}$ is the bolometric luminosity associated with the 100 Myr stellar population with constant SFR. For the large fraction of galaxies in our sample (96.2%) that lack far-IR detections, we calculated ${L}_{{\rm{bol}}}^{{\rm{young}}}$ and SFRs using dust-extinction corrected UV luminosities following the Calzetti et al. (2000) extinction law:

$$\begin{eqnarray}&&{L}_{{\rm{bol}}}^{{\rm{young}}}\approx {10}^{0.72{A}_{V}}3.2\nu {l}_{\nu }(2800\;{\rm{\mathring{\rm A} }}),\end{eqnarray} \tag{ 3 }$$

where the factor of 3.2 provides the bolometric correction of the young population SED (see above) to $\nu {l}_{\nu }(2800\;{\rm{\mathring{\rm A} }})$ for an A_V = 0 (i.e., $C({A}_{V}=0)=3.2$). Figure 1 displays the distribution of $({L}_{{\rm{UV,}}\;{\rm{obs}}}^{{\rm{young}}}+{L}_{{\rm{IR}}})/{L}_{{\rm{UV,}}\;{\rm{obs}}}^{{\rm{young}}}$ versus A_V for the 931 galaxies with far-IR detections and overlay the formal extinction-law prediction. It is clear from the distribution of far-IR detected sources that there is substantial statistical scatter in the relation, and a bias is present due to the fact that galaxies with low $({L}_{{\rm{UV,}}\;{\rm{obs}}}^{{\rm{young}}}+{L}_{{\rm{IR}}})/{L}_{{\rm{UV,}}\;{\rm{obs}}}^{{\rm{young}}}$ and low A_V will be excluded, since they would not be far-IR detected. Nonetheless, these data indicate that using A_V and $\nu {l}_{\nu }(2800\;{\rm{\mathring{\rm A} }})$ (without an IR measurement) to estimate ${L}_{{\rm{UV,}}\;{\rm{obs}}}^{{\rm{young}}}+{L}_{{\rm{IR}}}$ and the implied SFR for an object will result in a 1σ uncertainty of ∼0.5 dex, a factor of ∼2 times larger than the uncertainty of a far-IR based measurement of ${L}_{{\rm{IR}}}$. Hereafter, we adopt uncertainties of 0.2 dex and 0.5 dex for SFRs calculated using UV plus far-IR and the UV-only data, respectively.

Zoom In Zoom Out Reset image size

Of the 24,491 galaxies in our base sample, the majority are low mass (${M}_{\star }\lt {10}^{9}$ ${M}_{\odot }$) and will not produce significant X-ray emission per galaxy, nor will they provide substantial contributions to the global X-ray emission of the universe. Going forward, we chose to exclude low-mass galaxies with ${M}_{\star }\lt {10}^{9}$ ${M}_{\odot }$ from further analyses, since (1) they are more than two orders of magnitude below the knee of the stellar mass function of galaxies, which is at ${M}_{\star }\gtrsim {10}^{11}$ ${M}_{\odot }$ out to $z\gt 4$, and therefore not representative of the stellar mass in the universe (e.g., Muzzin et al. 2013); (2) they are optically faint and have large uncertainties in their rest-frame parameters; and (3) they effectively dilute X-ray stacking signals when included with higher-mass galaxies that have, e.g., similar SFR/${M}_{\star }$ values. We note that low-mass galaxies are relatively metal-poor and may indeed have enhanced levels of XRB emission per unit stellar mass or SFR (e.g., Linden et al. 2010; Basu-Zych et al. 2013a, 2016; Prestwich et al. 2013; Brorby et al. 2014, 2016; Douna et al. 2015); however, given their relatively low-mass, these enhancements are not important to the overall average global scaling relations.

After removing the low-mass galaxies, our resulting main sample contains 4,898 ${M}_{\star }\geqslant {10}^{9}$ ${M}_{\odot }$ galaxies within 7${}^{\prime }$ of the CDF-S aimpoint. In Figure 2, we show the distributions of ${M}_{\star }$, SFR, and z for the main sample. Given the low stellar mass limit of our sample ($\approx 2$ orders of magnitude below the knee of the stellar mass function at all redshifts), X-ray scaling relations derived from this sample will be representative of the global stellar populations in the universe. We tested this assumption in the local universe by deriving scaling relations from the Lehmer et al. (2010; hereafter, L10) galaxy sample (see L10 for procedures for deriving scaling constant) that both include and exclude ${M}_{\star }\lt {10}^{9}$ ${M}_{\odot }$ galaxies and found no material differences in derived relations. Specifically, $[{L}_{2-10\mathrm{keV}}({\rm{LMXB}})/{M}_{\star }{]}_{z\quad =\quad 0}^{{\rm{include}}}={2.2}_{-1.1}^{+1.9}\times {10}^{29}$ erg s⁻¹ ${M}_{\odot }$⁻¹ and $[{L}_{2-10\mathrm{keV}}({\rm{LMXB}})/{M}_{\star }{]}_{z\quad =\quad 0}^{{\rm{exclude}}}={2.4}_{-1.2}^{+2.0}\times {10}^{29}$ erg s⁻¹ ${M}_{\odot }$⁻¹, and $[{L}_{2-10}({\rm{HMXB}})/{\rm{SFR}}{]}_{z\quad =\quad 0}^{{\rm{include}}}={1.8}_{-0.4}^{+0.5}\times {10}^{39}$ erg s⁻¹ (M_⊙ yr⁻¹)⁻¹ and $[{L}_{2-10}({\rm{HMXB}})/{\rm{SFR}}{]}_{z\quad =\quad 0}^{{\rm{exclude}}}=(1.7\pm 0.3)\times {10}^{39}$ erg s⁻¹ (M_⊙ yr⁻¹)⁻¹. For consistency, however, when comparing stacking results with the local L10 sample, we make use of scaling relations that exclude ${M}_{\star }\lt {10}^{9}$ ${M}_{\odot }$ galaxies.

Zoom In Zoom Out Reset image size

Within the main sample, 870 galaxies have far-IR based measurements of ${L}_{{\rm{IR}}}$ and SFRs from ${L}_{{\rm{UV,}}\;{\rm{obs}}}^{{\rm{young}}}+{L}_{{\rm{IR}}}$ (Equation (2)), while the remaining 4,028 galaxies have SFRs measured from A_V and $\nu {l}_{\nu }(2800\;{\rm{\mathring{\rm A} }})$ (Equations (2) and (3)). To ensure the properties of our main sample were robust, we compared our values of ${M}_{\star }$ and SFR to those available in the literature. Using the Rainbow catalog galaxy sample builder,²² we constructed a sample of F160W-selected sources from the Cosmic Assembly Near-infrared Deep Extragalactic Legacy Survey (CANDELS; Grogin et al. 2011; Koekemoer et al. 2011) that included stellar mass values from Mobasher et al. (2015) and SFR values following Santini et al. (2015). The 4,898 main sample galaxies were cross-matched to the CANDELS catalog using a matching radius of 10, and 3,284 matches were obtained. We found that the mean ratios and 1σ scatters between our values and those in the CANDELS catalogs were $\mathrm{log}[{M}_{\star }/{M}_{\star }({\rm{CANDELS}})]=-0.04\pm 0.27$ and $\mathrm{log}[{\rm{SFR}}/{\rm{SFR}}({\rm{CANDELS}})]=-0.04\pm 0.60$. This result indicates that there are no systematic differences between methods, despite there being non-negligible scatter. The level of scatter is comparable to our adopted errors on these quantities, which we account for throughout the rest of this paper.

Using a conservative constant 05 matching radius, we cross-matched our main sample of 4,898 galaxies with the 613 X-ray point-sources in the ≈6 Ms CDF-S main catalog that are within 7${}^{\prime }$ of the aimpoint. There are 388 matches to X-ray point sources that are detected in at least one of the Luo et al. (2016, in preparation) standard bands: 0.5–7 keV (FB), 0.5–2 keV (SB), and 2–7 keV (HB). We note that the majority of the 225 X-ray sources without matches within 05 also have optical/near-IR counterparts, but those counterparts are either (1) low-mass galaxies excluded from our sample (see above) or (2) sources with optical offsets >05. We chose to exclude sources with offsets >05 from our analysis to strictly limit false matches. We determined the number of likely false matches by shifting the X-ray source positions by 10'' in four directions and re-matching them to the main sample using the above criterion. This exercise produced 6–14 matches per shift, suggesting a false-match rate of 1.5%–3.6%. We classified ≈6 Ms CDF-S sources with $z\gt 0$ (i.e., not Galactic stars) as either "AGN" or "normal galaxies" using the following five criteria (as outlined in Section 4.4 of Xue et al. 2011):

• 1.  
  A source with an intrinsic X-ray luminosity of ${L}_{0.5-7\mathrm{keV}}\geqslant 3\times {10}^{42}$ erg s⁻¹ was classified as an AGN. We obtained ${L}_{0.5-7\mathrm{keV}}$ using the following procedure. We first estimated the relationship between the 2–7 keV to 0.5–2 keV count-rate ratio and intrinsic column density and redshift for a power-law model SED that includes both Galactic and intrinsic extinction (in xspec ${wabs}\times {zwabs}\times {zpow}$) with a fixed ${{\rm{\Gamma }}}_{{\rm{int}}}=1.8$. For each source, we obtained the intrinsic extinction, ${N}_{{\rm{H,int}}}$, using the count-rate ratio and redshift. For some cases, the count-rate ratios were not well constrained due to lack of detections in both the 0.5–2 keV and 2–7 keV bands; for these sources an effective ${\rm{\Gamma }}=1.4$ was assumed. From our SED model, we could then obtain the absorption-corrected 0.5–7 keV flux, ${f}_{0.5-7\mathrm{keV},{\rm{int}}}$, and luminosity following ${L}_{0.5-7\mathrm{keV}}\;=4\pi {d}_{L}^{2}{f}_{0.5-7\mathrm{keV},\mathrm{int}}{(1+z)}^{{{\rm{\Gamma }}}_{{\rm{int}}}-2}$, where d_L is the luminosity distance.
• 2.  
  A source with an effective photon index of ${\rm{\Gamma }}\leqslant 1.0$ was classified as an obscured AGN.
• 3.  
  A source with X-ray-to-optical (using R-band as the optical reference) flux ratio of $\mathrm{log}({f}_{{\rm{X}}}$/${f}_{{\rm{R}}}\gt -1$ (where ${f}_{{\rm{X}}}={f}_{0.5-7\mathrm{keV}}$, ${f}_{0.5-2\mathrm{keV}}$, or ${f}_{2-7\mathrm{keV}}$) was classified as an AGN.
• 4.  
  A source with excess (i.e., a factor of ≥3) X-ray emission over the level expected from pure star formation was classified as an AGN, i.e., with ${L}_{0.5-7\mathrm{keV}}\gtrsim 3\times (8.9\times {19}^{17}{L}_{\text{1.4 GHz}})$, where ${L}_{\text{1.4 GHz}}$ is the rest-frame 1.4 GHz monochromatic luminosity in units of W Hz⁻¹ and $8.9\times {19}^{17}{L}_{\text{1.4 GHz}}$ is the expected X-ray emission level that originates from star-forming galaxies (see Alexander et al. 2005 for details).
• 5.  
  A source with optical spectroscopic features characteristic of AGN activity (e.g., broad-line emission and/or high-excitation emission lines) was classified as an AGN. Using a matching radius of 05, we cross-matched the ≈6 Ms CDF-S sources with spectroscopic catalogs from Szokoly et al. (2004), Mignoli et al. (2005), and Silverman et al. (2010) to identify these optical spectroscopic AGN.

All sources with $z\gt 0$ that were not classified as AGN from the above five criteria, were classified as normal galaxies. Out of the 388 X-ray detected sources, we classified 141 as normal galaxies, with the remaining 247 sources classified as AGN.

As described in L10), estimates of the ${L}_{{\rm{X}}}$(LMXB)/${M}_{\star }$ and ${L}_{{\rm{X}}}$(HMXB)/SFR scaling relations for a galaxy population can be obtained by empirically constraining the following relation:

$$\begin{eqnarray*}&&{L}_{{\rm{X}}}({\rm{XRB}})={L}_{{\rm{X}}}({\rm{LMXB}})+{L}_{{\rm{X}}}({\rm{HMXB}})=\alpha {M}_{\star }+\beta {\rm{SFR}},\end{eqnarray*}$$

$$\begin{eqnarray}&&{L}_{{\rm{X}}}{({\rm{XRB}})/{\rm{SFR}}=\alpha ({\rm{SFR}}/{M}_{\star })}^{-1}+\beta ,\end{eqnarray} \tag{ 4 }$$

where ${L}_{{\rm{X}}}({\rm{XRB}})$ is the total X-ray luminosity due to the XRB population, and $\alpha \equiv {L}_{{\rm{X}}}({\rm{LMXB}})/{M}_{\star }$ and $\beta \equiv {L}_{{\rm{X}}}({\rm{HMXB}})/{\rm{SFR}}$ are fitting constants. As outlined in Section 4 below, XRBs typically dominate the total galaxy-wide emission at energies above ≈1–2 keV. Therefore, in practice, Equation (4) can be constrained using total galaxy-wide luminosities derived from a hard bandpass (e.g., the 2–10 keV band; see L10 for further details on the $z\approx 0$ normal-galaxy population). However, due to the relatively high background and low Chandra effective area at >2 keV, only 23 of the 141 normal galaxies are detected in the HB,limiting our ability to constrain α and β. Fortunately, at the median redshift of the 116 normal galaxies detected in the FB, ${\langle z\rangle }_{{\rm{median}}}=0.67$, the FB itself samples the rest-frame 0.8–12.8 keV bandpass, which is expected to be dominated by XRBs (see Section 4 below); however, only normal galaxies at $z\gtrsim 2$ are likely to be entirely dominated by XRB emission across the full 0.5–7 keV range. Nonetheless, we can obtain reasonable constraints on α and β for the FB-detected normal-galaxy population.

Figure 3 displays the FB luminosity per unit SFR (${L}_{0.5-7\mathrm{keV}}$/SFR) versus sSFR for the 141 X-ray detected galaxies in our sample (including 116 FB detections). An inverse relation between ${L}_{0.5-7\mathrm{keV}}$/SFR and SFR/${M}_{\star }$ provides a better fit to the data than a simple constant ratio of ${L}_{0.5-7\mathrm{keV}}$/SFR, with an anti-correlation between ${L}_{0.5-7\mathrm{keV}}$/SFR and SFR/${M}_{\star }$ being significant at the >99.99% confidence level (based on a Spearman's ρ rank correlation). Following the functional form presented in Equation (4), we use these data to extract values of α and β for the FB detected sample. The galaxy sample has a median redshift of ${\langle z\rangle }_{{\rm{median}}}=0.67$ with an interquartile range of ${\rm{\Delta }}z$ = 0.41–0.98; the fits produce ${\alpha }_{{\rm{FB}}}(z\sim 0.7)$ = $({7.2}_{-2.0}^{+2.5})\times {10}^{29}$ erg s⁻¹ ${M}_{\odot }^{-1}$ and ${\beta }_{{\rm{FB}}}(z\sim 0.7)={8.5}_{-0.7}^{+0.8}\times {10}^{39}$ erg s⁻¹ (${M}_{\odot }$ yr⁻¹)⁻¹. This value of β is in reasonable agreement with the mean ${L}_{0.5-7\mathrm{keV}}$/SFR $\approx \;7.4\times {10}^{39}$ erg s⁻¹ (${M}_{\odot }$ yr⁻¹)⁻¹ presented by Mineo et al. (2014) from a sample of $z\quad \lesssim$ 1.3 X-ray and radio detected galaxies in the CDF-N and CDF-S that have sSFR $\gtrsim \quad {10}^{-10}$ yr⁻¹ (see Figure 3).

Zoom In Zoom Out Reset image size

If we assume a power-law X-ray SED with photon index ${\rm{\Gamma }}=2$, ${\alpha }_{{\rm{FB}}}$ and ${\beta }_{{\rm{FB}}}$ can be converted to the 2–10 keV bandpass equivalents by dividing by 1.64. When compared with the z = 0 values of α and β, derived for the 2–10 keV band from L10, ${({\alpha }_{z\sim 0.7}/{\alpha }_{z=0})}_{2-10\mathrm{keV}}\approx 4.7$ and ${({\beta }_{z\sim 0.7}/{\beta }_{z=0})}_{2-10\mathrm{keV}}\approx 3.0$.²³ This result suggests that both the ${L}_{{\rm{X}}}$(LMXB)/${M}_{\star }$ and ${L}_{{\rm{X}}}$(HMXB)/SFR scaling relations may increase significantly with redshift, with potentially stronger evolution of the LMXB population; broadly consistent with the F13 population-synthesis predictions.

We caution that the above result is inherently biased in nature, since the galaxy sample includes only X-ray detected galaxies that were identified as normal galaxies using some dependence on X-ray scaling relations appropriate for the local universe (see Section 3.1 of Lehmer et al. 2012). For example, among other criteria, normal galaxies were discriminated from AGN by having ${L}_{0.5-7\mathrm{keV}}/{L}_{\text{1.4 GHz}}$ and X-ray–to–optical flux ratios below some limiting values. These criteria effectively place a ceiling on the maximum ${L}_{0.5-7\mathrm{keV}}$/SFR value for a normal galaxy in the X-ray selected sample. Furthermore, as a result of the X-ray selection, we are more complete to high-SFR/high-sSFR galaxies that are X-ray luminous compared to the low-SFR/low-sSFR galaxy population. This could bias the X-ray selected sample toward galaxies with relatively large values of ${L}_{0.5-7\mathrm{keV}}$/SFR at low-SFR/low-sSFR, which would effectively bias α high.

As discussed in Section 3, the vast majority of the normal galaxies in our sample have X-ray emission below the individual source detection limit of the ≈6 Ms CDF-S. Since scaling relations derived from the X-ray detected normal galaxy population are biased by both the classification of a "normal galaxy" being restricted to sources that are within certain factors of z = 0 scaling relations (see criteria in Section 3), and the fact that the objects include only sources that are X-ray bright and detectable, these relations may not be representative of the broader galaxy population. To mitigate this limitation, we implement X-ray stacking analyses using complete subpopulations of galaxies, selected by physical properties (i.e., SFR and ${M}_{\star }$). In the paragraphs below, we describe our stacking procedure in general terms, as it is applied to a given subpopulation. In the next section (Section 5), we define the galaxy subpopulations for which we apply the stacking procedure. Variations of our stacking procedure have been implemented in a variety of previous studies of distant normal galaxies (e.g., Lehmer et al. 2005, 2007, 2008; Steffen et al. 2007; Basu-Zych et al. 2013b); for completeness, we highlight here the salient features of the procedure adopted in this paper.

Our stacking procedure begins with the extraction of on-source counts, local background counts, and exposure times for all 4,898 galaxies in our main sample using three X-ray bandpasses: 0.5–1 keV, 1–2 keV, and 2–4 keV (see below for justification). Using a circular aperture with a radius of ${R}_{{\rm{ap}}}$, we extracted Chandra source-plus-background counts S_i and exposure times T_i from the X-ray images and exposure maps, respectively. For galaxies with $z\lt 0.7$ we chose to use an aperture with ${R}_{{\rm{ap}}}=2\buildrel{\prime\prime}\over{.} 5$, and for $z\geqslant 0.7$, we used ${R}_{{\rm{ap}}}=1\buildrel{\prime\prime}\over{.} 5$. These apertures correspond to rest-frame physical radii of $\gtrsim$ 11 kpc for z = 0.3–0.7 where ${R}_{{\rm{ap}}}=2\buildrel{\prime\prime}\over{.} 5$, and ≳10 kpc for z = 0.7–4 and ≈8–10 kpc for $z\quad =$ 4–7, where ${R}_{{\rm{ap}}}=1\buildrel{\prime\prime}\over{.} 5$. This choice was motivated by the goal of both encompassing the vast majority of the X-ray emission from XRBs within galaxies while maintaining high signal-to-noise in our stacking (see below). The creation of images and exposure maps is described in Xue et al. (2011). In this process, each aperture was centered on the optical location of each galaxy. Using background maps, we measured the local backgrounds ${B}_{i,{\rm{local}}}$ within a circular aperture with radius ${R}_{{\rm{local}}}\quad =$ 25'' that was centered on each source. In practice, the background maps include X-ray counts from the galaxies that are not detected individually, so our local background circular apertures will include counts from the individually undetected sources of interest. We estimated the expected number of background counts within each on-source extraction aperture B_i by scaling the local background counts to the area of the extraction aperture (i.e., ${B}_{i}={B}_{i,{\rm{local}}}\times {R}_{{\rm{ap}}}^{2}/{R}_{{\rm{local}}}^{2}$).

To avoid contamination from bright unrelated X-ray sources, we excluded galaxies from our stacking analyses that were located within two times the 90% encircled-energy fraction radius of any unassociated X-ray detected sources. However, X-ray detected sources that were classified as normal galaxies by the criteria in Lehmer et al. (2012) were included in stacked subsamples (see Section 3). When stacking a given galaxy subsample, we calculated stacked source-plus-background ($S={\sum }_{i}{S}_{i}$) and background counts ($B={\sum }_{i}{B}_{i}$) and used these quantities to estimate the net counts (S − B) of the subsample. For each stacked subsample, we required that the signal-to-noise ratio (S/N $\equiv \;(S-B)/\sqrt{S}$) be greater than or equal to 3 (i.e., one-sided confidence level of ≈99.9%) for a detection. Due to the fact that our 15 radius stacking aperture encircles only a fraction of the point-source flux²⁴ for sources at relatively large off-axis angles, we calculated off-axis-dependent aperture corrections ${\xi }_{i}$ for each stacked source i. We note that the absolute astrometric match between the optical and X-ray frames (≲01; Luo et al. 2016, in preparation) is much smaller than our stacking aperture size, and therefore no corrections related to errors in the astrometric alignment need to be implemented. Since we are calculating average X-ray counts from the summed emission of sources in differing backgrounds and exposure times, we applied a single, representative exposure-weighted aperture correction, ξ. This factor, which was determined for each stacked subsample, was calculated as follows:

$$\begin{eqnarray}&&\xi \equiv \displaystyle \frac{{\displaystyle \sum }_{i}{\xi }_{i}\times {T}_{i}}{T},\end{eqnarray} \tag{ 5 }$$

where $T={\sum }_{i}{T}_{i}$. We find a range of $\xi \;=$ 1.4–2.0 for all bandpasses and stacking subsamples in this study. Using these corrections, we computed mean X-ray count-rates Φ for each stacked subsample using the following equation:

$$\begin{eqnarray}&&{\rm{\Phi }}=\xi \left(\displaystyle \frac{S-B}{T}\right).\end{eqnarray} \tag{ 6 }$$

For a given stacked subsample, the errors on Φ were calculated using a bootstrap resampling technique. For a given galaxy subsample, we constructed 10,000 "resampled" galaxy lists that were stacked using the above procedure. Each resampled list was constructed by drawing, at random, galaxies from the original list until the resampled list contained the same number of sources as the original list. A given resampled list will typically contain multiple entries of galaxies from the original list, without including all entries in the original list. The 1σ error on the count rate of a given stacked subsample was thus obtained by calculating the scatter in count-rates obtained for the 10,000 resampled lists.

To convert observed count rates to fluxes, we chose to make use of the mean X-ray SED shape of four starburst and normal galaxies from the NuSTAR galaxy program: NGC 253, M83, NGC 3256, and NGC 3310 (Lehmer et al. 2015). These galaxies span sSFR = 0.08–0.8 Gyr⁻¹, SFR = 3–30 M_⊙ yr⁻¹, and ${M}_{\star }=(0.9$–$6)\times {10}^{10}$ ${M}_{\odot }$, values similar to the mean values that we use in our stacked samples from $z\quad \approx$ 0–1.5 (see below). The SEDs for the four galaxies are well constrained over the 0.3–30 keV bandpass via simultaneous Chandra/XMM-Newton and NuSTAR observations (e.g., Wik et al. 2014; Lehmer et al. 2015; Yukita et al. 2016). The galaxy-wide Chandra/XMM-Newton plus NuSTAR spectra for these galaxies were modeled using a combination of two or three thermal plasma (hot gas) components with ${kT}\quad \approx$ 0.3–2 keV (modeled using apec in XSPEC) plus a component associated with XRBs (using power-law or broken power-law models). Both the hot gas and power-law components had varying levels of intrinsic absorption (see references for details); however, here we focus on observed spectra and do not attempt to correct for intrinsic absorption.

Figure 4(a) displays the mean X-ray SED with relative contributions from hot gas and XRBs indicated. On average, hot gas and XRBs dominate below and above $E\quad \approx$ 1.5 keV, respectively. Below ≈6 keV, the mean XRB spectrum is well characterized as an absorbed power-law with ${\rm{\Gamma }}\approx 1.8$–2.0, similar to previous Chandra studies of the mean X-ray spectra of bright XRBs in star-forming galaxies (e.g., Mineo et al. 2012a; Pacucci et al. 2014). Above ≈6 keV, the X-ray spectral shape develops a steeper spectral slope (${\rm{\Gamma }}\approx 2.5$), a feature not widely accounted for in XRB population studies (see, e.g., Kaaret 2014 for a discussion). The spectral turn over is a feature of the brightest XRBs in these galaxies, which include black hole binaries in ultraluminous and intermediate (or steep power-law) accretion states (e.g., Wik et al. 2014). These accretion states have non-negligible components from both an accretion disk and a power-law component from Comptonization (see, e.g., Remillard & McClintock 2006; Done et al. 2007; Gladstone et al. 2009; Bachetti et al. 2013; Walton et al. 2013, 2014; Rana et al. 2015). The Fragos et al. (2013b) population-synthesis models include synthetic spectra that assume both disk and power-law (Comptonization) components. For comparison, we show the Fragos et al. (2013b) synthetic XRB spectra for the global population at z = 0 and z = 7, which brackets the full range of spectral shapes across this redshift range. These synthetic spectra do not differ significantly between each other and are remarkably consistent with the observed XRB spectrum from our local star-forming galaxy SED at $E\gtrsim 1.5\;{\rm{keV}}$. At $E\lesssim 1.5\;{\rm{keV}}$, the synthetic spectra vary somewhat as a function of redshift, primarily due to variations in the assumed intrinsic absorption. Fragos et al. (2013b) assumed that the intrinsic column densities of XRBs varied as a function luminosity following the same distribution as XRBs in the Milky Way. The true intrinsic obscuration of XRBs in high-redshift galaxies is highly uncertain. Nonetheless, the agreement in SED shapes at $E\gtrsim 1.5\;{\rm{keV}}$ for the local star-forming galaxy sample and the synthetic spectra at all redshifts suggests that there is unlikely to be any significant variations in the XRB SED as a function of redshift.

Zoom In Zoom Out Reset image size

The intensities of the hot gas and XRB components have both been observed to broadly scale with SFR, and the hot gas temperature does not appear to vary significantly with SFR (e.g., Mineo et al. 2012a, 2012b). Given these results, and the predicted lack of evolution in the synthetic XRB SED from Fragos et al. (2013b), we do not expect the shapes of the hot gas and XRB components will change significantly with SFR. However, as we reach to higher redshift galaxies, it is expected that the XRB emission components will become more luminous per unit SFR due to declining XRB ages and metallicities (e.g., F13a). These physical changes may have some effect on the hot gas emission per unit SFR, as the supernova rate and mechanical energy from stellar winds are affected by metallicity (e.g., Côté et al. 2015). The exact form of this dependence has yet to be quantified specifically; however, some evidence suggests that the increase in XRB emission with declining metallicity dwarfs the hot gas emission (e.g., see the 0.5–30 keV SED of low-metallicity galaxy NGC 3310 in Figure 6 of Lehmer et al. 2015). We therefore expect that the ratio of hot gas to XRB emission will likely decrease with increasing redshift, leading to a more XRB-dominant SED at high redshift. Unfortunately, our data are not of sufficient quality to constrain such an effect directly.

Since our primary goal is to utilize X-ray stacking of galaxy populations to study the underlying XRB populations, it is imperative that we perform our stacking in a bandpass that will be dominated by XRB emission. Clearly, high-energy bandpasses that sample rest-frame energies $E\quad \gtrsim$ 2 keV would accomplish this goal; however, since the Chandra effective area curve peaks at 1–2 keV and declines rapidly at higher energies, the corresponding S/N of stacked subsamples is prohibitively low for bandpasses above observed-frame $E\quad \approx$ 2 keV. Fortunately, since the rest-frame energy range for a fixed observed-frame bandpass shifts to higher energies with increasing redshift, we can perform stacking in bandpasses near the peak of Chandra's response and still measure directly emission dominated by XRB populations.

Figure 4(b) presents the redshift-dependent fraction of the Chandra counts that are provided by XRBs for observed-frame 0.5–1 keV, 1–2 keV, and 2–4 keV bandpasses, assuming that the SED does not evolve with redshift. The 1–2 keV and 2–4 keV bandpasses probe majority contributions from galaxies dominated by XRB populations beyond $z\approx 0.15$, which corresponds to the lowest redshift galaxies in our sample (see below). In an attempt to measure directly emission from XRBs, while preserving high S/N for our stacked subsamples, we perform stacking in the 1–2 keV bandpass and use these stacking results to measure k-corrected rest-frame 2–10 keV luminosities. Figure 4(a) shows the rest-frame regions of the SED sampled by observed-frame 1–2 keV at z = 0.3, z = 0.6, and z = 1.5, which are, respectively, the lowest, second lowest, and median redshift intervals for stacked subsamples; we define stacking subsamples explicitly in Section 5 below. As a check on our SED assumptions and possible AGN contamination (see Section 6.1 below), we also stacked subsamples in the 0.5–1 keV and 2–4 keV bands. The 0.5–1 keV band stacks sample softer energies that can be used to estimate rest-frame 0.5–2 keV luminosities by applying small k-corrections. Although our focus is on the 2–10 keV emission, we present, for completeness, 0.5–2 keV constraints throughout the rest of the paper.

For a given stacked galaxy subsample, we estimated mean X-ray fluxes using the following equation:

$$\begin{eqnarray}&&{f}_{{E}_{1}-{E}_{2}}={A}_{{E}_{1}-{E}_{2}}{{\rm{\Phi }}}_{{E}_{1}-{E}_{2}},\end{eqnarray} \tag{ 7 }$$

where ${A}_{{E}_{1}-{E}_{2}}$ provides the count-rate to flux conversion within the observed-frame ${E}_{1}-{E}_{2}$ bandpass, based on our adopted SED (see Figure 4(a)) and the Chandra response function. Errors on ${f}_{{E}_{1}-{E}_{2}}$ were calculated by propagating uncertainties on ${A}_{{E}_{1}-{E}_{2}}$ from our starburst and normal galaxy sample and bootstrapping errors on ${{\rm{\Phi }}}_{{E}_{1}-{E}_{2}}$. Finally, X-ray luminosities were computed using X-ray fluxes and k-corrections based on our adopted SED following:

$$\begin{eqnarray}&&{L}_{{E}_{1}^{\prime }-{E}_{2}^{\prime }}={k}_{{E}_{1}-{E}_{2}}^{{E}_{1}^{\prime }-{E}_{2}^{\prime }}4\pi {d}_{L}^{2}{f}_{{E}_{1}-{E}_{2}},\end{eqnarray} \tag{ 8 }$$

where ${k}_{{E}_{1}-{E}_{2}}^{{E}_{1}^{\prime }-{E}_{2}^{\prime }}$ provides the correction between observed-frame ${E}_{1}-{E}_{2}$ band and rest-frame ${E}_{1}^{\prime }-{E}_{2}^{\prime }$ band. Errors on ${L}_{{E}_{1}^{\prime }-{E}_{2}^{\prime }}$ include errors on ${k}_{{E}_{1}-{E}_{2}}^{{E}_{1}^{\prime }-{E}_{2}^{\prime }}$ from our starburst and normal galaxy sample and the errors on ${f}_{{E}_{1}-{E}_{2}}$. Across the redshift range $z\quad =\quad$ 0–4, the k-correction values range from 2.1–2.7 (1.6–2.1) for correcting observed-frame 1–2 keV (0.5–1 keV) to rest-frame 2–10 keV (0.5–2 keV).

When measuring X-ray scaling relations for stacked bins (e.g., ${L}_{{\rm{X}}}$ versus SFR), we made use of mean physical-property values $\langle {A}_{{\rm{phys}}}\rangle$ (where ${A}_{{\rm{phys}}}$ could be the SFR or ${M}_{\star }$). Errors on these mean quantities were calculated following a Monte Carlo approach. For each stacked subsample, we computed 10,000 simulated perturbed mean values of a given physical quantity. For the kth simulation, the perturbed mean was computed following:

$$\begin{eqnarray}&&\langle {A}_{{\rm{phys}},k}^{{\rm{pert}}}\rangle =\displaystyle \frac{1}{{N}_{{\rm{gal}}}}\displaystyle \sum _{i=1}^{{N}_{{\rm{gal}}}}\;{A}_{{\rm{phys}},i}+{N}_{i,k}^{{\rm{pert}}}{\sigma }_{{A}_{{\rm{phys}},i}},\end{eqnarray} \tag{ 9 }$$

where ${N}_{{\rm{gal}}}$ is the number of galaxies in the stack and the value of ${N}_{i,k}^{{\rm{pert}}}$ is a random number drawn from a Gaussian distribution, centered at zero, with a width of unity. The value of ${\sigma }_{{A}_{{\rm{phys}},i}}$ is the 1σ error on the physical property value ${A}_{{\rm{phys}},i}$ measured for the ith galaxy in the stacked subsample. The errors on the SFR and ${M}_{\star }$ were presented in Section 2. Finally, for a stacked subsample, the error on the mean was estimated as the standard deviation of the perturbed mean values:

$$\begin{eqnarray}&&{\sigma }_{\langle {A}_{{\rm{phys}}}\rangle }=\displaystyle \frac{1}{{N}_{{\rm{sim}}}}{\left[\displaystyle \sum _{k=0}^{{N}_{{\rm{sim}}}}{(\langle {A}_{{\rm{phys}},k}^{{\rm{pert}}}\rangle -\langle {A}_{{\rm{phys}}}\rangle )}^{2}\right]}^{1/2},\end{eqnarray} \tag{ 10 }$$

where ${N}_{{\rm{sim}}}=10,000$ is the number of simulations performed for a given stacked subsample. When computing ratios of mean quantities (e.g., $\mathrm{log}{L}_{{\rm{X}}}$/SFR) errors on all quantities were propagated following the methods outlined in Section 1.7.3 of Lyons (1991).

As discussed in Section 1 our primary goal is to measure the redshift evolution of the scaling factors $\alpha \equiv {L}_{{\rm{X}}}$(LMXB)/${M}_{\star }$ and $\beta \equiv {L}_{{\rm{X}}}$(HMXB)/SFR that were introduced in Equation (4). Since the ratio of HMXB-to-LMXB emission is sensitive to sSFR, our strategy for calculating α and β involves stacking galaxy subsamples that are divided by sSFR in a variety of redshift bins. We began by dividing our galaxy sample into 11 redshift intervals that were chosen to both include large numbers (≈200–1000) of galaxies and span broad dynamic ranges in sSFR. Figure 5 presents the SFR versus sSFR for each of the 11 redshift intervals plus the $z\approx 0$ sample studied by L10. In each panel, we have highlighted galaxies with SFRs calculated using UV plus far-IR (black filled circles) and UV-only (gray open circles) data.

Zoom In Zoom Out Reset image size

For each redshift interval, we divided the galaxy samples into bins of sSFR for stacking. Similar to our choice to exclude galaxies with ${M}_{\star }\lt {10}^{9}$ ${M}_{\odot }$, we further placed lower boundaries on the SFR by which galaxies are included in stacking so that low SFR galaxies would not dilute the stacked signals. As displayed in Figure 5, the low-SFR bounds vary with redshift, such that SFR^bound spans ≈0.1–1 M_⊙ yr⁻¹ for the majority of the stacked subsamples ($z\quad =$ 0–1.5), and increases to ≈10–30 M_⊙ yr⁻¹ for the highest redshift intervals ($z\quad \approx$ 1.5–7).

To first order, obtaining average values of the X-ray properties of galaxy populations in bins selected by sSFR alone can be used to provide a good characterization of the average ${L}_{{\rm{X}}}$/SFR as a function of sSFR; however, finer division of the galaxy sample (e.g., by both sSFR and ${M}_{\star }$) would allow investigation of the X-ray emission of galaxy subsamples with varying physical properties (e.g., metallicity and stellar age; e.g., Linden et al. 2010; Basu-Zych et al. 2013a, 2016; Fragos et al. 2013b; Prestwich et al. 2013; Brorby et al. 2014, 2016; Douna et al. 2015). However, the goal of this paper is to establish the average evolution of the scaling relations (i.e., α and β from Equation (4)) with redshift due to changes in global properties. In a future paper, we will investigate how the emission from XRBs varies with further division of the galaxy samples by physical properties.

In each redshift interval, we divided the full sSFR range into sSFR bins that contained nearly equal numbers of galaxies per bin. The number of galaxies per subsample contained between 12 and 60 sources (median of 23 sources); in total, we stacked 63 galaxy subsamples. Figure 5 shows the resulting stacking bins for the subsamples, and Table 1 summarizes the subsample properties. From Figure 5, it is apparent that the sSFR bins cluster around characteristic central sSFR values and are broader for large and small sSFR values. This behavior is due to the densely populated regions along the galaxy "main sequence," where the SFR and ${M}_{\star }$ are correlated (e.g., Elbaz et al. 2007; Noeske et al. 2007; Karim et al. 2011; Whitaker et al. 2014). Despite this effect, we are able to span over 2 orders of magnitude in sSFR for stacked subsamples in each of the redshift intervals for $z\lesssim 2.5$. Therefore, we can place good direct constraints on α and β in individual redshift intervals (see below).

Table 1.  Stacked Subsample Properties

Subsample ID	${z}_{{\rm{lo}}}$–${z}_{{\rm{up}}}$	$\langle z\rangle$	${N}_{{\rm{gal}}}$	${N}_{{\rm{det}}}$	SFR ${}_{{\rm{lo}}}$–SFR ${}_{{\rm{up}}}$	$\langle {\rm{SFR}}\rangle$	$\mathrm{log}{M}_{\star }$	sSFR ${}_{{\rm{lo}}}$–sSFR ${}_{{\rm{up}}}$	$\langle {\rm{sSFR}}\rangle$
					(M⊙ yr−1)	(M⊙ yr−1)	(${M}_{\odot }$)	(Gyr−1)	(Gyr−1)
(1)	(2)	(3)	(4)	(5)	(6)	(7)	(8)	(9)	(10)
1	0.0–0.5	0.39	12	2	0.2–21.7	0.8 ± 0.4	10.7 ± 0.1	0.01–0.02	0.02 ± 0.01
2	⋯	0.35	12	4	0.2–21.7	1.3 ± 0.5	10.5 ± 0.1	0.02–0.06	0.04 ± 0.02
3	⋯	0.37	12	1	0.2–21.7	1.8 ± 0.4	9.9 ± 0.1	0.06–0.31	0.22 ± 0.07
4	⋯	0.32	12	3	0.2–21.7	1.9 ± 0.4	9.7 ± 0.1	0.31–0.49	0.40 ± 0.12
5	⋯	0.33	13	3	0.2–21.7	5.5 ± 1.1	10.0 ± 0.1	0.49–0.69	0.58 ± 0.08
6	⋯	0.33	12	6	0.2–21.7	4.6 ± 1.0	9.7 ± 0.1	0.69–1.03	0.81 ± 0.22
7	⋯	0.36	14	6	0.2–21.7	3.0 ± 0.6	9.4 ± 0.1	1.03–3.00	1.30 ± 0.36
8	0.5–0.7	0.60	21	6	0.2–55.1	0.9 ± 0.4	10.7 ± 0.1	0.01–0.03	0.02 ± 0.01
9	⋯	0.62	21	5	0.2–55.1	1.6 ± 0.6	10.6 ± 0.1	0.03–0.13	0.05 ± 0.02
10	⋯	0.59	23	2	0.2–55.1	3.1 ± 0.5	10.1 ± 0.1	0.13–0.36	0.24 ± 0.07

Note. Col. (1): unique stacked subsample identification number. The subsample ID has been assigned to each subsample and is ordered based on ascending redshift bin, and within each redshift bin, ascending sSFR. Col. (2): lower and upper redshift boundaries of the subsample. Col. (3): mean redshift. Col. (4): number of galaxies stacked in the given bin. Col. (5): number of sources detected individually in each bin. Col. (6): lower and upper SFR boundaries for the subsample in units of M_⊙ yr⁻¹. Col. (7): mean SFR and 1σ error on the mean. Col. (8): logarithm of the mean stellar mass ${M}_{\star }$ in units of ${M}_{\odot }$ and 1σ error on the mean. The stellar mass is bounded by limits on SFR and sSFR, as well as a hard lower bound of ${M}_{\star }^{{\rm{lim}}}={10}^{9}$ ${M}_{\odot }$. Col. (9): lower and upper sSFR boundaries for the subsample in units of Gyr⁻¹. Col. (10): mean sSFR and 1σ error on the mean.

Only a portion of this table is shown here to demonstrate its form and content. A machine-readable version of the full table is available.

Download table as:  DataTypeset image

For the purpose of comparing our CDF-S stacking results with constraints from local galaxies, we calculated average X-ray luminosities for the local galaxy sample culled by L10 using sSFR bins selected following the same basic binning strategy defined above. The L10 local galaxies sample contains data for 66 nearby ($D\quad \lt$ 400 Mpc; median $D\approx 60$ Mpc) galaxies that were all observed by Chandra, and have SFR and ${M}_{\star }$ values calculated following the same procedures outlined in Section 2 (see footnote 23). The full sample of 66 nearby galaxies includes subsamples of normal galaxies from Colbert et al. (2004), as well as LIRGs and ULIRGs from L10 and Iwasawa et al. (2009), and spans a broad range of SFR (SFR = 0.08–200 M_⊙ yr⁻¹) and stellar mass (${M}_{\star }={10}^{8-11}$ ${M}_{\odot }$).

The upper left-hand corner of Figure 5 presents the SFR versus sSFR distribution of galaxies in the L10 sample. A series of seven sSFR ranges were chosen, in which we calculate average properties to be compared with the CDF-S stacking results. As with the CDF-S data, these sSFR ranges were defined to contain nearly the same number of sources per bin; however, mean X-ray luminosities for each bin were calculated by averaging the luminosities obtained for each galaxy. This operation should yield equivalent results to those obtained by our CDF-S stacking procedure, and throughout the rest of this paper we utilize these mean values, in conjunction with our X-ray stacking results, as constraints on the z = 0 population properties. We emphasize, however, that the L10 compilation includes only 2–10 keV galaxy-wide luminosities, which probe directly the XRB emission. In order to compare our CDF-S stacked 0.5–2 keV luminosities, we estimated comparable mean 0.5–2 keV XRB luminosities for each L10 galaxy subsample by applying a bandpass correction to the 2–10 keV XRB emission and adding an estimate of the 0.5–2 keV hot gas emission. The XRB bandpass correction was calculated assuming a power-law spectrum with intrinsic ${N}_{{\rm{H}}}=3\times {10}^{21}$ cm⁻² and ${\rm{\Gamma }}=2.0$ (see Section 3.1 of Mineo et al. 2012a for justification), and the 0.5–2 keV hot gas emission was estimated for each sample using the Mineo et al. (2012b) relation: ${L}_{0.5-2\;\mathrm{keV}}^{\mathrm{hot}\ \mathrm{gas}}$/(erg s⁻¹) $\approx 1.5\times {10}^{39}$ SFR/(M_⊙ yr⁻¹).

6.1. Spectral Properties and AGN Contamination

Applying the stacking procedure described in Section 4, to the CDF-S galaxy subsamples defined in Section 5, we obtained mean X-ray properties for the subsamples, which are summarized in Table 2. Out of the 63 subsamples, we obtained 23, 60, and 11 significant (>3σ) stacked detections in the 0.5–1 keV, 1–2 keV, and 2–4 keV bands, respectively; 3 subsamples were not detected in any of the three bands. Figure 6 displays the (1–2 keV)/(0.5–1 keV) and (2–4 keV)/(1–2 keV) mean count-rate ratios (hereafter, BR1 and BR2, respectively) versus redshift for each of the stacked subsamples and show the expected band-ratios for our adopted SED. For the majority of subsamples, BR1 and BR2 appear to be in agreement with our adopted SED for almost all stacked subsamples, with a few exceptions where 3σ lower limits on BR1 lie above the SED error envelope (at $z\approx 2$) and BR2 measurements are somewhat elevated (at the ≈1σ level) at $z\quad \approx$ 1–2. This broad agreement provides some confidence that our stacked subsamples are not strongly contaminated by obscured AGN and low-luminosity AGN below the detection limit.

Zoom In Zoom Out Reset image size

Table 2.  Galaxy Stacking Results

ID	S/N			Net Counts			BR1	BR2	$\mathrm{log}{L}_{0.5-2\mathrm{keV}}$	$\mathrm{log}{L}_{2-10\mathrm{keV}}$	$\mathrm{log}{L}_{2-10\mathrm{keV}}$/SFR	${f}_{{\rm{AGN}}}^{2-10\;\mathrm{keV}}$
	0.5–1 keV	1–2 keV	2–4 keV	0.5–1 keV	1–2 keV	2–4 keV			(erg s−1)	(erg s−1)	(erg s−1 (M⊙ yr−1)−1)	(%)
(1)	(2)	(3)	(4)	(5)	(6)	(7)	(8)	(9)	(10)	(11)	(12)	(13)
1	4.5	5.0	3.3	55.1 ± 31.1	77.2 ± 53.3	60.8 ± 51.5	1.3 ± 3.2	0.8 ± 2.8	40.2 ± 0.2	40.1 ± 0.2	40.2 ± 0.3	${0.00}_{0.00}^{+0.00}$
2	9.4	17.8	9.1	146.3 ± 71.5	437.4 ± 296.9	199.5 ± 159.3	2.9 ± 6.0	0.5 ± 1.5	40.5 ± 0.2	40.7 ± 0.2	40.6 ± 0.3	${0.00}_{0.00}^{+0.00}$
3	2.4	3.2	0.3	<32.6	46.7 ± 28.3	<50.4	>1.4	<1.1	<39.9	39.8 ± 0.2	39.5 ± 0.2	${0.00}_{0.00}^{+0.00}$
4	6.7	10.1	5.6	90.6 ± 58.0	190.3 ± 156.2	111.0 ± 72.6	2.0 ± 5.9	0.6 ± 3.4	40.2 ± 0.2	40.2 ± 0.3	40.0 ± 0.3	${0.00}_{0.00}^{+0.00}$
5	9.6	13.7	5.9	153.1 ± 115.3	297.1 ± 180.1	122.4 ± 92.0	1.9 ± 7.6	0.4 ± 1.1	40.4 ± 0.2	40.5 ± 0.2	39.7 ± 0.2	${0.00}_{0.00}^{+0.15}$
6	14.1	14.4	6.0	265.6 ± 179.5	315.4 ± 199.3	119.5 ± 73.4	1.2 ± 3.6	0.4 ± 1.1	40.7 ± 0.2	40.5 ± 0.2	39.8 ± 0.2	${0.00}_{0.00}^{+0.00}$
7	5.6	7.9	4.3	76.3 ± 18.4	144.6 ± 22.3	88.0 ± 21.3	1.8 ± 2.4	0.6 ± 0.8	40.1 ± 0.1	40.1 ± 0.1	39.7 ± 0.1	${0.00}_{0.00}^{+0.00}$
8	4.7	9.4	4.0	72.3 ± 19.7	207.0 ± 54.6	96.3 ± 30.9	2.8 ± 3.9	0.5 ± 0.7	40.5 ± 0.1	40.7 ± 0.1	40.7 ± 0.2	${0.00}_{0.00}^{+0.00}$
9	3.7	6.4	2.7	54.4 ± 16.6	128.7 ± 44.8	<70.6	2.3 ± 3.4	<0.6	40.4 ± 0.1	40.5 ± 0.1	40.3 ± 0.2	${0.00}_{0.00}^{+0.01}$
10	2.8	3.7	1.8	<44.2	72.8 ± 29.3	<72.3	>1.6	<1.0	<40.2	40.2 ± 0.2	39.7 ± 0.2	${0.01}_{0.01}^{+0.29}$

Note. Col. (1): subsample identification as defined in Table 1. Col. (2)–(4): signal-to-noise ratio obtained for each stacked subsample in the 0.5–1 keV, 1–2 keV, and 2–4 keV bands, respectively, following the procedure outlined in Section 4. Col. (5)–(7): background-subtracted net counts and 1σ errors obtained for the 0.5–1 keV, 1–2 keV, and 2–4 keV bandpasses, respectively. Col. (8) and (9): band ratios and 1σ errors for BR1 $\equiv$ (1–2 keV)/(0.5–1 keV) and BR2 $\equiv$ (2–4 keV)/(1–2 keV), respectively. Col. (10) and (11): logarithm of the mean stacked 0.5–2 keV and 2–10 keV luminosities, respectively, in units of erg s⁻¹ with 1σ errors on the mean values. Col. (12) Logarithm of the mean 2–10 keV luminosity per mean SFR ratio in units of erg s⁻¹ (M_⊙ yr⁻¹)⁻¹ and 1σ error. Col. (13) Estimate of the fraction of 2–10 keV emission, in percentage terms, due to undetected AGN; error bars are 1σ. These estimates are computed following the methodology outlined in Section 6.1.

Only a portion of this table is shown here to demonstrate its form and content. A machine-readable version of the full table is available.

Download table as:  DataTypeset image

To estimate quantitatively the level by which individually undetected AGN contribute to the stacked subsamples, we employed an approach in which (1) the supermassive black hole accretion distribution (in terms of the Eddington accretion rate) is estimated, using the known AGN population, to very low Eddington fractions, and (2) this distribution is used to estimate the expected total contributions to stacked subsamples from AGN with luminosities that fall below the individual source detection threshold.

First, for each of the 4,898 galaxies in our main sample, we estimated the central black hole mass using the following relation: $\mathrm{log}{M}_{{\rm{BH}}}\approx 8.95+1.40\mathrm{log}({M}_{\star }/{M}_{\odot })$ (Reines & Volonteri 2015), where ${M}_{\star }$ is the stellar mass of the galaxy. For AGN with HB detections, we calculated rest-frame 2–10 keV luminosities in terms of Eddington fraction, $\lambda \equiv {L}_{{\rm{HB}}}/{L}_{{\rm{HB,Edd}}}$, where ${L}_{{\rm{HB,Edd}}}=(1.26\times {10}^{38}$ erg s⁻¹)/${C}_{{\rm{bol}}}$ (${M}_{{\rm{BH}}}/{M}_{\odot }$), and ${C}_{{\rm{bol}}}$ is the luminosity-dependent bolometric correction as presented in Equation (2) of Hopkins et al. (2007). For each galaxy in the main galaxy sample, we derived an Eddington fraction limit, ${\lambda }_{{\rm{lim}}}={L}_{{\rm{HB,lim}}}/{L}_{{\rm{HB,Edd}}}$, below which we would not be able to detect an AGN if it were present. For a given source ${L}_{{\rm{HB,lim}}}=4\pi {d}_{L}^{2}{f}_{{\rm{HB,lim}}}$, where ${f}_{{\rm{HB,lim}}}$ is the HB flux limit, which we extracted from the spatially dependent sensitivity map constructed by HB Luo et al. (2016, in preparation).

Using the above information, Eddington accretion fraction probability distributions, $p(\lambda )$, were calculated by extracting the number of AGN within a bin of λ divided by the number of galaxies by which we could have detected such an AGN (using ${\lambda }_{{\rm{lim}}}$). Past studies (e.g., Rafferty et al. 2011; Aird et al. 2012; Mullaney et al. 2012; Wang et al. 2013; Hickox et al. 2014) have suggested that such a distribution is likely to be SFR dependent (and close to linear) due to the overall correlation between ${M}_{{\rm{BH}}}$ and ${M}_{\star }$. We therefore derived $p(\lambda )$ for two SFR regimes: SFR = 0.1–4 M_⊙ yr⁻¹and SFR > 4 M_⊙ yr⁻¹ to infer the SFR dependence. Figure 7 shows our estimates of $p(\lambda )$ for the two different SFR regimes. We find the distributions to have similar λ dependencies with normalization increasing with SFR. We found that the following parameterization of $p(\lambda ,{\rm{SFR}})$ fits well the observed λ and SFR dependencies:

$$\begin{eqnarray}&&p{(\lambda ,{\rm{SFR}})=\xi {(\lambda /{\lambda }_{{\rm{c}}})}^{-{\gamma }_{E}}\mathrm{exp}(-\lambda /{\lambda }_{{\rm{c}}})({\rm{SFR}}/1{M}_{\odot }{\mathrm{yr}}^{-1})}^{{\gamma }_{s}},\end{eqnarray} \tag{ 11 }$$

where $\xi =0.002$ dex⁻¹, ${\lambda }_{{\rm{c}}}=0.21$, ${\gamma }_{E}=0.31$, and ${\gamma }_{s}=0.57$ are fitting constants. This functional form was motivated by previous estimates of similar distributions presented in the literature (e.g., Aird et al. 2012; Hickox et al. 2014). The curves in Figure 7 show the predicted curves from Equation (11) fixed to the median SFR of each of the two SFR-regimes.

Zoom In Zoom Out Reset image size

To estimate X-ray undetected AGN contributions for a given stacked subsample, we used a Monte Carlo approach. For each galaxy in a stacked subsample that was not detected in the X-ray band, we first drew a value of ${\lambda }_{i}$ probabilistically following the distribution in Equation (11). We then converted ${\lambda }_{i}$ to a HB luminosity following ${L}_{{\rm{HB}},i}({\rm{AGN}})={\lambda }_{i}{L}_{{\rm{HB,Edd}},i}$, where ${L}_{{\rm{HB,Edd}},i}$ is uniquely defined for a galaxy by its black hole mass (see above). Occasionally a random draw will predict a value of ${L}_{{\rm{HB}},i}$(AGN) above the detection limit. In such cases a new estimate of ${L}_{{\rm{HB}},i}$(AGN) is chosen until it falls below the detection limit. This procedure allows us to estimate the total contribution that AGN below the detection limit provide to the stacked emission ${f}_{{\rm{AGN}}}^{{E}_{1}^{\prime }-{E}_{2}^{\prime }}={\sum }_{i}{L}_{{\rm{HB}},i}$(AGN) ${k}_{{\rm{HB}}}^{{E}_{1}^{\prime }-{E}_{2}^{\prime }}$(AGN)/${L}_{{{\rm{E}}}_{1}^{\prime }-{{\rm{E}}}_{2}^{\prime }}$, where ${L}_{{{\rm{E}}}_{1}^{\prime }-{{\rm{E}}}_{2}^{\prime }}$ is the total stacked luminosity (see Equation (8)) and ${k}_{{\rm{HB}}}^{{E}_{1}^{\prime }-{E}_{2}^{\prime }}$(AGN) provides a bandpass correction from the HB to the ${E}_{1}^{\prime }-{E}_{2}^{\prime }$ band for the AGN contribution.

To account for random errors, we computed ${f}_{{\rm{AGN}}}^{{E}_{1}^{\prime }-{E}_{2}^{\prime }}$ for a given stacked subsample using 1000 Monte Carlo trials, allowing us to estimate its most probable value and 1σ range; values of ${f}_{{\rm{AGN}}}^{2-10\;\mathrm{keV}}$ are provided in Column (13) of Table 2. All values of ${f}_{{\rm{AGN}}}^{2-10\;\mathrm{keV}}$ are less than 0.1 (i.e., 10% contribution from AGN), with a median value of 0.007. These values of ${f}_{{\rm{AGN}}}^{2-10\;\mathrm{keV}}$ are much smaller than the errors on the stacked luminosities, and we therefore conclude that AGN below the detection threshold do not significantly impact our results.

6.2. The X-Ray/SFR Correlation

Although our primary goal is to measure α and β as a function of redshift, it is worth exploring first how basic empirical scaling relations that have been studied and widely used for local and distant star-forming galaxies (e.g., ${L}_{{\rm{X}}}\quad \propto$ SFR) compare with those inferred from our stacked subsamples and L10 local comparison.

Since our galaxy subsample selections were based primarily on intervals of sSFR with strict lower limits on ${M}_{\star }$ and SFR, the mean SFRs for our galaxy subsamples span a modest range of SFR (≈2 dex; see Figure 5), which allows direct measurement of the ${L}_{{\rm{X}}}$/SFR correlations for our data. Figure 8 presents ${L}_{0.5-2\mathrm{keV}}$ and ${L}_{2-10\mathrm{keV}}$ versus SFR for the stacked galaxy subsamples (filled black circles) and local L10 comparison subsamples (filled blue triangles), and provides results from additional previous studies for comparison (see discussion below). Consistent with previous investigations, we find strong correlations between ${L}_{0.5-2\mathrm{keV}}$ and ${L}_{2-10\mathrm{keV}}$ with SFR (e.g., Bauer et al. 2002; Grimm et al. 2002; Ranalli et al. 2003, 2012; Persic et al. 2004; Persic & Rephaeli 2007; Lehmer et al. 2008, 2010; Iwasawa et al. 2009; Symeonidis et al. 2011, 2014; Mineo et al. 2012a, 2012b, 2014; Vattakunnel et al. 2012). We performed linear fits to the CDF-S stacked data and local L10 comparison sample to derive best-fit values of the following linear model:

$$\begin{eqnarray}&&\mathrm{log}{L}_{{\rm{X}}}={A}_{1}+\mathrm{log}{\rm{SFR}},\end{eqnarray} \tag{ 12 }$$

where ${L}_{{\rm{X}}}$ is in units of erg s⁻¹ and SFR is in units of M_⊙ yr⁻¹. The best-fit parameters for the 0.5–2 keV and 2–10 keV bands are listed in Table 3 and plotted as solid lines in Figure 8. Despite the strong ${L}_{{\rm{X}}}$/SFR correlation for our data, the linear fits to all stacked bins do not yield statistically acceptable fits; ${\chi }^{2}/\nu$ = 5.78 and 4.23 for the 0.5–2 keV and 2–10 keV bands, respectively, with resulting residual scatter of 0.37 and 0.32 dex. To test whether the poor fits were due to some nonlinearity in the ${L}_{{\rm{X}}}$–SFR relation, we further performed nonlinear fits to our data following the form:

$$\begin{eqnarray}&&\mathrm{log}{L}_{{\rm{X}}}={A}_{2}+{B}_{2}\mathrm{log}{\rm{SFR}}.\end{eqnarray} \tag{ 13 }$$

The best-fit parameter B₂ has some small variation from unity for both the 0.5–2 keV and 2–10 keV bandpasses; however, the fit does not yield a statistically robust characterization of the data (${\chi }^{2}/\nu =2.75$ and 3.38, respectively).

Zoom In Zoom Out Reset image size

Table 3.  Summary of Global Fits to Stacked Data Sets

Model Description	Parameter	0.5–2 keV				2–10 keV
		Param Value	${\chi }^{2}/\nu$	ν	σ(dex)	Param Value	${\chi }^{2}/\nu$	ν	σ(dex)
$\mathrm{log}{L}_{{\rm{X}}}={A}_{1}+\mathrm{log}$ SFR	A1	39.59 ± 0.02	5.78	29	0.37	39.78 ± 0.02	4.23	66	0.32
$\mathrm{log}{L}_{{\rm{X}}}={A}_{2}+{B}_{2}\mathrm{log}$ SFR	A2	40.06 ± 0.05				40.12 ± 0.05
	B2	0.65 ± 0.04	2.75	28	0.29	0.71 ± 0.04	3.38	65	0.29
$\mathrm{log}{L}_{{\rm{X}}}={A}_{3}+{B}_{3}\mathrm{log}{\rm{SFR}}+{C}_{3}\mathrm{log}(1+z)$	A3	39.83 ± 0.07				39.82 ± 0.05
	B3	0.74 ± 0.04				0.63 ± 0.04
	C3	0.97 ± 0.18	1.77	27	0.24	1.31 ± 0.11	1.32	64	0.20
${L}_{{\rm{X}}}={\alpha }_{0}{(1+z)}^{\gamma }{M}_{\star }+{\beta }_{0}{(1+z)}^{\delta }{\rm{SFR}}$ [CDF-S only]	$\mathrm{log}{\alpha }_{0}$	28.87 ± 0.24				29.30 ± 0.28
	$\mathrm{log}{\beta }_{0}$	39.66 ± 0.17				39.40 ± 0.08
	γ	4.59 ± 0.80				2.19 ± 0.99
	δ	−0.10 ± 0.72	0.84	19	0.16	1.02 ± 0.22	0.91	56	0.17
${L}_{{\rm{X}}}={\alpha }_{0}{(1+z)}^{\gamma }{M}_{\star }+{\beta }_{0}{(1+z)}^{\delta }{\rm{SFR}}$ [L10 plus CDF-S]	$\mathrm{log}{\alpha }_{0}$	29.04 ± 0.17				29.37 ± 0.15
	$\mathrm{log}{\beta }_{0}$	39.38 ± 0.03				39.28 ± 0.05
	γ	3.78 ± 0.82				2.03 ± 0.60
	δ	0.99 ± 0.26	0.79	26	0.16	1.31 ± 0.13	1.05	63	0.17
F13a Population Synthesis	Model 245						2.71	67	0.20
F13a Population Synthesis	Model 269						1.56	67	0.19

Note. All global models were fit to a combination of local (z = 0) galaxy subsamples plus stacked high-redshift subsamples from the CDF-S (derived in this work). For the 0.5–2 keV band, we utilized seven local galaxy subsamples compiled from L10 plus 23 stacked subsamples in the CDF-S that were significantly detected in the observed-frame 0.5–1 keV band. For the 2–10 keV band, we utilized seven local galaxy subsamples compiled by L10 plus the 60 stacked subsamples that were significantly detected in the CDF-S in the observed-frame 1–2 keV band. Only detected stacked subsamples (i.e., S/N $\gtrsim \;3\sigma$) were used in the global fits; upper limits were excluded.

Download table as:  ASCIITypeset image

In addition to our local L10 sample, Figure 8 also displays the results derived from multiple local and distant-galaxy samples.²⁵ For local comparisons, we have chosen results from IR-selected galaxies, which span normal star-forming galaxies (${L}_{{\rm{IR}}}\;\approx$ 10⁹–10¹¹ ${L}_{\odot }$) to luminous/ultraluminous IR galaxies (LIRGs/ULIRGs; ${L}_{{\rm{IR}}}\quad \approx$ 10¹¹–10¹² ${L}_{\odot };$ Iwasawa et al. 2011; Symeonidis et al. 2011), normal star-forming galaxy samples selected to have sSFRs $\gtrsim \quad {10}^{-10}$ yr⁻¹ (Mineo et al. 2012a, 2012b), and UV-selected Lyman break analogs (LBAs; Basu-Zych et al. 2013a), which are rare galaxies that have properties similar to z ∼ 2–3 Lyman break galaxies (LBGs, e.g., relatively low metallicity and compact UV morphologies). For distant-galaxy comparisons, we have included results from Basu-Zych et al. (2013b), Mineo et al. (2014), and Symeonidis et al. (2014), which presented results from X-ray stacking of $z\quad \approx$ 1.5–4 LBGs in the CDF-S, correlating X-ray and radio detected CDF-N and CDF-S sources at $z\lesssim 1.3$, and X-ray stacking of Herschel selected star-forming galaxies in the CDF-N and CDF-S at $z\gtrsim 1.5$, respectively.

There is basic agreement in the trends found between these comparison samples and our data. Our data do have higher ${L}_{{\rm{X}}}$/SFR values on average, but the scatter encompasses the majority of the comparison samples. We suspect that the heterogeneity in the comparison sample selections introduces significant scatter in even the local ${L}_{{\rm{X}}}$–SFR relations. For the local galaxies, the scatter may be explained due to sample differences in (1) the relative contributions of LMXBs and HMXBs (e.g., Mineo et al. 2012a, 2012b use only high-sSFR galaxies); (2) metallicity, which results in varying levels in HMXB formation per unit SFR (e.g., Basu-Zych et al. 2013a study low-metallicity galaxies); and (3) X-ray absorption due to galaxy selection (e.g., Iwasawa et al. 2011 and Symeonidis et al. 2011 study IR selected galaxies that may be influenced by absorption; see also Luangtip et al. 2015).

As the typical sSFR, stellar age, and metallicity of the galaxies in the universe evolve with redshift, it is expected that ${L}_{{\rm{X}}}$/SFR will change as the HMXB and LMXB populations respond accordingly (F13a), so redshift-related scatter is expected to be introduced by combining constraints from various redshifts. This behavior is evident in the scatter for the distant-galaxy comparison samples (i.e., the Basu-Zych et al. 2013b; Mineo et al. 2014 and Symeonidis et al. 2014 samples). In fact, for SFR $\gtrsim$ 10 M_⊙ yr⁻¹, it is apparent from Figure 8 that the relatively large scatter in ${L}_{{\rm{X}}}$/SFR for our stacked CDF-S data arises primarily due to a redshift effect (symbol sizes are proportional to redshift). For a given SFR, ${L}_{{\rm{X}}}$/SFR seems to increase with redshift, an indication that the XRB emission per unit SFR is indeed increasing with redshift. For SFR $\lesssim$ 10 M_⊙ yr⁻¹, the large scatter in ${L}_{{\rm{X}}}$/SFR is likely to be due to varying contributions from LMXBs (see below).

6.3. Redshift-Dependent Evolution of XRB Scaling Relations

As discussed above, direct redshift independent scaling relations of X-ray luminosity with SFR are not statistically robust across all galaxy subsamples at all redshifts. From Figure 8, there are qualitative indications for significant redshift evolution in the scaling relations for SFR $\gtrsim$ 10 M_⊙ yr⁻¹. In this section, we investigate how redshift and galaxy physical properties influence ${L}_{{\rm{X}}}$. Hereafter, we focus our discussion on results from the rest-frame 2–10 keV band (probed by observed-frame 1–2 keV), since this bandpass probes directly XRB emission and provides the best constraints on the redshift evolution of X-ray emission due to the relatively large number of significant stacked detections: (60 in the observed-frame 1–2 keV band versus 23 in the observed-frame 0.5–1 keV band). For completeness, however, we provide equivalent measurements and results for the rest-frame 0.5–2 keV emission throughout the rest of this paper.

6.3.1. SFR and Redshift Dependence

Figure 9 displays the ${L}_{2-10\mathrm{keV}}$/SFR ratio of our stacked subsamples versus redshift in four SFR intervals; we include the several comparison samples presented in Figure 8(b). This representation reveals that much of the scatter in the ${L}_{2-10\mathrm{keV}}$/SFR relation and apparent discrepancies between other studies can be reconciled by positive redshift and negative SFR dependence on the ${L}_{2-10\mathrm{keV}}$/SFR ratio. Some previous investigations of X-ray emission from distant galaxy samples (e.g., Lehmer et al. 2008; Ranalli et al. 2012; Vattakunnel et al. 2012; Symeonidis et al. 2014) concluded that ${L}_{{\rm{X}}}$/SFR at $z\approx 1$ is consistent with local scaling relations, and that the relation does not evolve with redshift. Such a conclusion is dependent on both (1) the choice of the local comparison sample, which varies for different galaxy samples (see Figure 8), and (2) the redshift of the sample. For example, comparison of the $z\approx 1$ data from Symeonidis et al. (2014) with the Mineo et al. (2012a) relation at z = 0 (brown dotted–dashed lines in Figure 9) indicates that the two samples have similar ${L}_{2-10\mathrm{keV}}$/SFR. However, the Symeonidis et al. (2014) values are somewhat larger than the z = 0 values from L10 and Iwasawa et al. (2011), which would imply that there is positive evolution of the ${L}_{2-10\mathrm{keV}}$/SFR relation with redshift. In this study, we find that at $z\gtrsim 2$, ${L}_{2-10\mathrm{keV}}$/SFR is higher than both the L10 and Mineo et al. (2012a) relations, indicating that the redshift evolution is robust.

Zoom In Zoom Out Reset image size

Using stacked samples of $z\quad \approx$ 1–4 Lyman break galaxies (LBGs) in the ≈4 Ms CDF-S survey, stacking results from star-forming galaxies in the ≈2 Ms CDF-S/CDF-N surveys from Lehmer et al. (2008), and the local L10 galaxy sample, Basu-Zych et al. (2013b) found that the mean X-ray luminosity of star-forming galaxy samples could be characterized well using the following relation:

$$\begin{eqnarray}&&\mathrm{log}{L}_{{\rm{X}}}={A}_{3}+{B}_{3}\mathrm{log}{\rm{SFR}}+{C}_{3}\mathrm{log}(1+z).\end{eqnarray} \tag{ 14 }$$

Figure 9 indicates the Basu-Zych et al. (2013b) best-fit relation for comparison (dashed curves), and the Mineo et al. (2012a) best-fit local relation. Using the stacking results from this study and the compiled local galaxy samples from L10, we derived fitting constants A₃ = 39.82  ±  0.05, B₃ = 0.63  ±  0.04, and C₃ = 1.31  ±  0.11 (solid curves in Figure 9). These values predict somewhat more rapid redshift evolution than those from Basu-Zych et al. (2012b; ${A}_{3}=39.8$, ${B}_{3}=0.65$, and ${C}_{3}=0.89$), with significant divergence between fits at z ≳ 2.5–3. Differences in the best-fit function are largely driven by 2–3 $z\gtrsim 2$ stacked samples with SFR ≈ 10–35 M_⊙ yr⁻¹ from Basu-Zych et al. (Basu-Zych et al. 2013b; inverted triangles in Figure 9), and our stacked subsamples in the low-redshift ($z\lesssim 2$), low-SFR (SFR ≈ 0.3–10 M_⊙ yr⁻¹) regime, which contain significant scatter. Our model fit to the data produced a best-fit ${\chi }^{2}/\nu =1.32$, which is a substantial improvement over a single ${L}_{2-10\mathrm{keV}}$–SFR scaling relation that does not include redshift evolution (${\chi }^{2}/\nu \approx 4.23$), but is only marginally acceptable. For $\nu =64$ degrees of freedom, there is a 4.4% probability of obtaining ${\chi }^{2}/\nu \geqslant 1.32$.

Despite the statistical limitations, Equation (14) provides a first-order approach for estimating galaxy-wide ${L}_{2-10\mathrm{keV}}$, given values of z and SFR; the resulting best-fit relation has a statistical residual scatter of ≈0.23 dex. We caution, however, that the SFR $\lesssim$ 10 M_⊙ yr⁻¹ galaxies are predicted by XRB population-synthesis models to have ${L}_{2-10\mathrm{keV}}$/SFR that flatten above $z\approx 1.5$, just above the current limits of our survey (see Figure 8 of Basu-Zych et al. 2013b). It is therefore likely that Equation (14) significantly overpredicts ${L}_{2-10\mathrm{keV}}$/SFR for $z\approx 1.5$ galaxies with SFR $\lesssim$ 10 M_⊙ yr⁻¹, and is not appropriate in this regime.

Embeded within the above empirical parameterization is information about the evolution of the underlying physical properties of galaxies (e.g., stellar age and metallicity). For example, Basu-Zych et al. (2013b) made use of the XRB population-synthesis models of F13a, which predicted similar behavior to that of Equation (14), to interpret these trends as being due to the combined evolution of LMXB and HMXB populations. They reported that for a given redshift interval, the low-SFR populations were predicted to have strong contributions from both HMXBs and LMXBs. This conclusion leads to relatively large values of ${L}_{2-10\mathrm{keV}}$/SFR for low-SFR compared with those of high-SFR galaxies, which are expected to be dominated by HMXBs alone. With increasing redshift, declining stellar ages and metallicities produce brighter populations of LMXBs and HMXBs, respectively, thereby leading to a corresponding increase in ${L}_{2-10\mathrm{keV}}$/SFR.

As noted above, there is substantial remaining scatter in the best-fit parameterization from Equation (14). In particular, ${L}_{2-10\mathrm{keV}}$/SFR values for stacked subsamples with SFR $\lesssim$ 10 M_⊙ yr⁻¹, which we argue are expected to have contributions from both HMXBs and LMXBs, have significant scatter for a given redshift. We expect that our selection of subsamples by sSFR broadens the range of LMXB contributions to these stacked subsamples, e.g., compared to a selection by SFR alone.

6.3.2. Cosmic Evolution of LMXB and HMXB Populations

As described in Section 5, our choice to select galaxy subsamples in bins of sSFR was motivated by the expected scaling of LMXB and HMXB emission with ${M}_{\star }$ and SFR, respectively—i.e., to obtain α and β values for a range of redshifts. Following Equation (4), we expect ${L}_{2-10\mathrm{keV}}$/SFR will be inversely proportional to sSFR. Figure 10 presents ${L}_{2-10\mathrm{keV}}$/SFR versus sSFR for the stacked subsamples in each redshift interval. For the $z\lt 2.0$ redshift intervals, it is apparent that ${L}_{2-10\mathrm{keV}}$/SFR declines with increasing sSFR, as expected from Equation (4). Each panel of Figure 10 displays the L10 relation for $z\approx 0$ galaxies (gray curves) with the L10 data provided in the first panel of Figure 10. With increasing redshift, the ${L}_{2-10\mathrm{keV}}$/SFR values become increasingly offset above the L10 relation. In particular, the low-sSFR slopes of the ${L}_{2-10\mathrm{keV}}$/SFR versus sSFR curves become steeper with increasing redshift, indicating that the LMXB luminosity per unit mass (i.e., α) increases with increasing redshift.

Zoom In Zoom Out Reset image size

For the eight redshift intervals that contained $\geqslant 4$ detected subsamples (i.e., all redshift intervals at $z\lt 2.5$), we fit the data using both a constant model, $\mathrm{log}{L}_{2-10\mathrm{keV}}$/SFR $=\;{A}_{1}$, as well as our canonical model from Equation (4), from which we extract best-fit values of α and β. Table 4 provides the best-fit values for A₁ and α and β in each of the eight redshift intervals, and compares the goodness-of-fit for both models. With the exception of the $z\quad =$ 2.0–2.5 bin, our canonical model provides a statistically improved characterization of the stacked data over the constant model at ≈83%–99.99% confidence levels (based on an F-test; see Col. 10 in Table 4). The resulting ${\chi }^{2}$ values for these intervals indicate that the canonical model is generally a good fit to the data (${\chi }^{2}/\nu \quad =$ 0.2–1.2; median ${\chi }^{2}/\nu \;=$ 0.47). Figure 10 shows the resulting best-fit relations for the $z\lt 2.5$ redshift intervals as solid curves, and Figure 11 presents the corresponding best-fit values of α and β versus redshift for the $z\lt 2.5$ redshift intervals as filled black circles with error bars. Figure 11 also indicates the values of α and β measured for the X-ray detected sample (measured in Section 3; orange stars) and L10 (open circles), as well as the value of α derived from LMXBs in local early-type galaxies (open green triangle; Boroson et al. 2011) and β derived from HMXBs in local high-sSFR galaxies (open brown square; Mineo et al. 2012b).

Zoom In Zoom Out Reset image size

Table 4.  2–10 keV Fits to Stacked Data Sets In Each Redshift Bin

${z}_{{\rm{lo}}}$–${z}_{{\rm{up}}}$	${N}_{{\rm{det}}}$	$\mathrm{log}{L}_{2-10\mathrm{keV}}={A}_{1}+\mathrm{log}{\rm{SFR}}$			${L}_{2-10\mathrm{keV}}=\alpha {M}_{\star }+\beta {\rm{SFR}}$
		$\mathrm{log}{A}_{1}$			$\mathrm{log}\alpha$	$\mathrm{log}\beta$
		($\mathrm{log}$ erg s−1 (M⊙ yr−1)−1)	${\chi }^{2}/\nu$	ν	($\mathrm{log}$ erg s−1 (M⊙ yr−1)−1)	($\mathrm{log}$ erg s−1 ${M}_{\odot }^{-1}$)	${\chi }^{2}/\nu$	ν	${F}_{{\rm{prob}}}$
(1)	(2)	(3)	(4)	(5)	(6)	(7)	(8)	(9)	(10)
0.0–0.5	7	${39.83}_{-0.09}^{+0.08}$	1.99	6	${29.62}_{-0.45}^{+0.33}$	${39.66}_{-0.16}^{+0.16}$	1.00	5	0.9536
0.5–0.7	9	${39.71}_{-0.06}^{+0.06}$	5.28	8	${29.83}_{-0.20}^{+0.18}$	${39.42}_{-0.13}^{+0.13}$	0.45	7	1.0000
0.7–0.9	8	${39.73}_{-0.06}^{+0.06}$	2.88	7	${29.90}_{-0.27}^{+0.22}$	${39.55}_{-0.12}^{+0.11}$	0.49	6	0.9990
0.9–1.0	5	${39.89}_{-0.07}^{+0.07}$	0.83	4	${30.38}_{-0.46}^{+0.29}$	${39.74}_{-0.14}^{+0.13}$	0.14	3	0.9797
1.0–1.2	10	${39.85}_{-0.05}^{+0.05}$	1.67	9	${29.74}_{-0.31}^{+0.24}$	${39.74}_{-0.09}^{+0.09}$	0.21	8	1.0000
1.2–1.5	6	${39.99}_{-0.06}^{+0.06}$	0.64	5	${30.13}_{-1.24}^{+0.38}$	${39.92}_{-0.11}^{+0.11}$	0.47	4	0.8293
1.5–2.0	8	${39.90}_{-0.05}^{+0.05}$	1.66	7	${30.54}_{-0.30}^{+0.21}$	${39.76}_{-0.10}^{+0.09}$	1.18	6	0.9012
2.0–2.5	4	${39.97}_{-0.06}^{+0.06}$	2.55	3	${30.78}_{-0.65}^{+0.31}$	${39.86}_{-0.12}^{+0.12}$	3.41	2	0.3311

Note. Model fits to stacked galaxy subsample 2–10 keV luminosities include only stacked detections with observed-frame 1–2 keV S/N $\geqslant \;3$. Only redshift bins with more than five detected galaxy subsamples are included in this table. The redshift range and number of detected galaxy subsamples are provided in Columns (1) and (2), respectively. A constant model ($\mathrm{log}{L}_{2-10\mathrm{keV}}={A}_{1}+\mathrm{log}{\rm{SFR}};$ Col. (3)–(5)) is tested against our canonical model (${L}_{2-10\mathrm{keV}}=\alpha {M}_{\star }+\beta {\rm{SFR}};$ Col. (6)–(9)) for each redshift range. In all cases, except perhaps, for the last bin (i.e., at $z\approx 2.5$), the canonical model provides a better goodness of fit (${\chi }^{2};$ Col. (4) and (8)), and the significance of improvement in terms of the F-test is provided in Column (10).

Download table as:  ASCIITypeset image

From Figure 11 it is clear that α, the LMXB emission per unit stellar mass, evolves rapidly with redshift out to $z\approx 2.5$, while β, the HMXB emission per unit SFR, remains roughly constant over this redshift range, albeit with some evidence for a mild increase with redshift. As discussed in Section 1, and presented in F13a, this result is consistent with the basic expectations from XRB population-synthesis models, which predict that the decline in age and metallicity with increasing redshift would yield changes in respective LMXB and HMXB scaling relations. In the next section, we make direct comparisons of our measurements with the XRB population-synthesis predictions from F13a.

The constraints on α and β discussed above utilize only data with $z\lesssim 2.5$ on a per-redshift interval basis; however, our stacking analyses provide detections for galaxy subsamples out to $z\approx 4.0$. To better characterize the redshift evolution of α and β, we made use of the following global parameterization using the full set of 60 detections out to $z\approx 4$:

$$\begin{eqnarray*}&&{L}_{{\rm{X}}}({\rm{XRB}})(z)={L}_{{\rm{X}}}({\rm{LMXB}})(z)+{L}_{{\rm{X}}}({\rm{HMXB}})(z),\end{eqnarray*}$$

$$\begin{eqnarray}&&{L}_{{\rm{X}}}{({\rm{XRB}}){(z)={\alpha }_{0}(1+z)}^{\gamma }{M}_{\star }+{\beta }_{0}(1+z)}^{\delta }{\rm{SFR}},\end{eqnarray} \tag{ 15 }$$

where ${\alpha }_{0}$, ${\beta }_{0}$, γ, and δ are fitting constants. As discussed above, ${\alpha }_{0}$ and ${\beta }_{0}$ have already been constrained by L10, Boroson et al. (2011), and Mineo et al. (2012a); however, these values are somewhat discrepant due to differences in galaxy sample properties (see discussion in Section 6.1). We therefore chose to obtain values of ${\alpha }_{0}$ and ${\beta }_{0}$ independently by fitting Equation (15) to the CDF-S stacked data alone. We also combine our local L10 galaxy sample (i.e., the average values for seven sSFR bins) with the CDF-S data to improve overall constraints on the global evolution of the 2–10 keV emission from galaxies. By fitting both CDF-S data alone and L10-plus-CDF-S data, we can show how constraints on local galaxies influence the global redshift-dependent solution to Equation (15).

When using the CDF-S data alone, fitting our data to the parameterization provided in Equation (15) provides good overall fits to the average data (${\chi }^{2}/\nu =0.91$, for $\nu =56$ degrees of freedom) and a significant reduction in the resulting spread (0.17 dex; or 48%). Our fits characterize variations in the stacked, population-averaged emission; however, such averaging masks galaxy-to-galaxy variations, which can be significant. The true intrinsic galaxy-to-galaxy spread in ${L}_{2-10\mathrm{keV}}$, within each subsample, is expected to be larger (on the order of ≈0.2–0.4 dex; see, e.g., L10 and Mineo et al. 2012a) and sensitive to variations in metallicity, stellar age, and statistical variations in the XRB populations themselves. We obtain best-fit values of $\mathrm{log}{\alpha }_{0}=29.30\pm 0.28$, $\mathrm{log}{\beta }_{0}=39.40\pm 0.08$, $\gamma =2.19\pm 0.99$, and $\delta =1.02\pm 0.22$. Interestingly, these values of ${\alpha }_{0}$ and ${\beta }_{0}$, which are based solely on the $z\gtrsim 0.3$ CDF-S data, are in excellent agreement with the range of respective values obtained for z = 0: $\mathrm{log}{\alpha }_{0}\;=$ 29.1–29.2 (L10; Boroson et al. 2011) and $\mathrm{log}{\beta }_{0}\;=$ 39.2–39.6 (L10; Mineo et al. 2012a).

When we combine the L10 local comparison values with the CDF-S data and re-fit the ensemble data set to Equation (15), we obtain ${\chi }^{2}/\nu \approx 1.05$ for $\nu =63$ degrees of freedom, and values consistent with the CDF-S only fits: $\mathrm{log}{\alpha }_{0}=29.37\pm 0.15$, $\mathrm{log}{\beta }_{0}=39.28\pm 0.05$, $\gamma =2.03\pm 0.60$, and $\delta =1.31\pm 0.13$. These parameters indicate that the redshift increases of α and β (as noted above from the fits in each of the $z\lesssim 2.5$ redshift intervals) are significant at the ≈3.4σ and ≈10.1σ levels, respectively. However, the significances of evolution in both α and β are dependent on the inclusion of the L10 local galaxy constraints. From the CDF-S data alone, the respective redshift evolution of LMXB and HMXB scaling relations is significant at the ≈2.2σ and ≈4.6σ levels, respectively. The L10, Boroson et al. (2011), and Mineo et al. (2012a) measurements of ${\alpha }_{0}$ and ${\beta }_{0}$ are based on galaxy samples drawn from the archive, which have selection biases that are different from those in our CDF-S galaxy subsamples; therefore, future investigations to measure ${\alpha }_{0}$ and ${\beta }_{0}$ in a local galaxy sample with similar selection to the CDF-S galaxy subsample would help clarify the results presented in this paper.

The above results indicate that the ${L}_{{\rm{X}}}$–SFR correlation is not universal, but rather depends critically on both SFR and ${M}_{\star }$ and evolves with redshift. The top three panels of Figure 12 show the 2–10 keV residuals to the best-fit global relations presented throughout Section 6, which include: (1) a constant "universal" relation $\mathrm{log}{L}_{2-10\mathrm{keV}}$/SFR = A₁ (Equation (12)); (2) a SFR and redshift-dependent relation $\mathrm{log}{L}_{2-10\mathrm{keV}}$/SFR =${A}_{3}+{B}_{3}\mathrm{log}{\rm{SFR}}+{C}_{3}\mathrm{log}(1+z)$ (Equation (14)); and (3) a SFR, ${M}_{\star }$, and redshift-dependent relation ${L}_{2-10\mathrm{keV}}={\alpha }_{0}{(1+z)}^{\gamma }{M}_{\star }+{\beta }_{0}{(1+z)}^{\delta }{\rm{SFR}}$ (Equation (15)). This succession of global models produces significant improvement in goodness of fit (${\chi }^{2}/\nu \approx 4.2$, 1.3, and 1.1, respectively), as well as significant reductions in residual scatter ($\sigma \;=$ 0.32, 0.20, and 0.17 dex, respectively). The key results from this study can be summarized as (1) the LMXB emission per unit stellar mass, $\alpha \equiv {L}_{{\rm{X}}}$(LMXB)/${M}_{\star }$, increases rapidly with increasing redshift $\propto {(1+z)}^{2-3};$ and (2) the HMXB emission per unit SFR, $\beta \equiv {L}_{{\rm{X}}}$(HMXB)/SFR, increases mildly with increasing redshift $\propto (1+z)$. These proportionalities appear to be valid out to at least $z\;\approx$ 2–3; however, as we will discuss below, there are physical reasons to expect that these proportionalities will not extend to even higher redshifts. In this section, we discuss our results in the context of the XRB population-synthesis models from F13a, which have made similar predictions.

Zoom In Zoom Out Reset image size

7.1. What Drives the Evolution of LMXB and HMXB Populations?

The F13a XRB population-synthesis study constructed 288 unique XRB models, which predict the evolution of α and β over the history of the universe (based on the StarTrack population-synthesis code; Belczynski et al. 2002, 2008), accounting for evolution of the star-formation history and metallicity using the Millenium II cosmological simulation with a semi-analytic galaxy evolution prescription (Guo et al. 2011). These models include prescriptions for XRBs that are formed through stellar evolutionary channels only and do not include XRBs that may form due to dynamical interactions in high stellar density environments like globular clusters (see, e.g., Benacquista & Downing 2013 for a review). As noted in F13a, dynamically-formed XRBs are significant to the overall X-ray emission in the most massive elliptical galaxies, but are minority populations globally (see Section 2.1 of F13a for further details).

The six parameters that were varied in the 288 unique scenarios in F13a included common-envelope efficiency (two parameters), wind prescriptions, binary mass ratio distribution, kicks from SNe, and the stellar IMF (see F13a for details). Six out of the 288 models (Models 245, 229, 269, 205, 249, and 273, ordered by decreasing likelihood; see F13a for details on each of these models) provided acceptable predictions of z = 0 observational constraints (from Tzanavaris & Georgantopoulos 2008; L10; Boroson et al. 2011; Mineo et al. 2012a), with the remaining 282 were deemed to be likely unrealistic. The six models with the highest probability had the same prescriptions for common-envelope efficiency and binary mass-ratio distribution, but varied in prescriptions for stellar-wind strength, SNe kick amplitudes, and stellar IMFs. One of the limitations in constraining the F13a models was the lack of strong observational constraints at $z\gt 0$, a limitation that is mitigated significantly by the current study.

We compared each of the 288 F13a model predictions for the redshift evolution of α and β with our combined local L10 and CDF-S stacked sample measurements. Our comparisons were limited to the 2–10 keV constraints, which are expected to directly probe the XRB populations with negligible contamination from hot gas. Each of the 288 models provides a unique redshift, SFR, and ${M}_{\star }$ dependent prediction for every subsample with no free parameters. Although these models are not strictly fit to the data, a goodness-of-fit parameter can be assigned to each model and compared across models.

Table 5 provides a list of the ${\chi }^{2}/\nu$ values for the top ten F13a models, sorted by ascending ${\chi }^{2}$. The highest probability model, Model 269, has a ${\chi }^{2}/\nu =1.56$ for $\nu =67$ degrees of freedom. Given the number of degrees of freedom, the probability of obtaining ${\chi }^{2}/\nu \geqslant 1.56$ is 0.22%, indicating that this model is formally unacceptable; however, given the limited ranges and numbers of parameters that were varied in the F13a models, it is likely that minor tweaks, and the inclusion of additional parameters not varied by F13a (e.g., the fraction of stars in binaries), could yield statistically robust fits to our data. Future generations of XRB population-synthesis models will explore these issues. For further comparison, Table 5 lists the relative rankings of the F13a models from the Tzanavaris et al. (2013) study of XRB luminosity functions in local galaxies and the Tremmel et al. (2013) study of the evolution of normal galaxy X-ray luminosity functions. These investigations, along with the F13a original study, rank Model 269 highly (within the top 12 models). Notably, Model 269 provides a better global prediction to the data than the best-fit constant ${L}_{2-10\mathrm{keV}}$/SFR (see Figure 11 and Table 3) and a comparable fit to the data as the redshift and SFR dependent (Equation (14)) models.

Table 5.  Fragos Model Goodness of Fit

F13a	Rank	Rank	Rank	Rank	${\chi }^{2}/\nu$
Model	(This Study)	(F13a a)	(Tz13b)	(Tr13c)
269	1	3	5	12	1.56
273	2	6	12	6	1.64
77	3	21	19	27	1.82
201	4	10	9	19	1.98
81	5	26	24	20	2.14
53	6	20	11	22	2.16
246	7	25	23	23	2.24
232	8	35	60	28	2.30
230	9	17	35	16	2.33
231	10	37	34	31	2.43

Notes.

^aRank of likelihood based on F13a. ^bRank from Tzanavaris et al. (2013) study of local SINGS galaxies. ^cRank from Tremmel et al. (2013) study of normal-galaxy X-ray luminosity function evolution.

Download table as:  ASCIITypeset image

From a physical point of view, Model 269 is the same in all six parameters as the previous best-fit model from F13a, Model 245, except that Model 269 allows for binary systems to emerge after a common-envelope phase involving the donor stars going through the Hertzsprung gap, while Model 245 assumes that such a common-envelope phase will automatically lead to the merging of the two stars and a termination of any possible subsequent XRB phase. This has the effect of Model 269 having mildly elevated emission from HMXBs (i.e., elevated values of β) over Model 245 across the redshift range studied here. The variation of this parameter is motivated on theoretical grounds and is currently unconstrained by observations. Theoretically, a Hertzsprung gap star does not have a clear entropy jump at the core-envelope transistion (Ivanova & Taam 2004), so once a companion is engulfed within such a star during the common-envelope phase, there is no clear boundary where an inspiral would cease, and consequently a merger is expected (see, e.g., Taam & Sandquist 2000). However, on energetic grounds, it is possible for the system to successfully exit the common-envelope phase without merging. If binaries indeed survive a common-envelope phase where donor stars are in the Hertzsprung gap (i.e., Model 269), then the predicted numbers of gravitational wave sources from double black hole mergers is predicted to be much larger (by up to a factor of ∼500) than if they are destroyed (i.e., the Model 245 case; see F13a and Belczynski et al. 2007 for details). Future studies from Advanced LIGO and Advanced Virgo gravitational wave detectors will likely constrain this effect independently (Belczynski et al. 2015, 2016; Abbott et al. 2016).

Figure 10 displays the Model 269 predicted ${L}_{2-10\mathrm{keV}}$/SFR versus sSFR for all redshift intervals, spanning $z\quad \approx$ 0–7 (dotted red curves). The predictions for Model 269 are remarkably similar to those from our best-fit global model (Equation (15); dashed blue curves in Figure 10) over the redshift range $z\;\approx$ 0–2.5. In general, Model 269 and our best-fit global model are in good agreement; however, differences in the high-sSFR predictions are clear at $z\;\gtrsim$ 1.5–2, with Model 269 generally predicting lower values of ${L}_{{\rm{X}}}$/SFR.

Figure 11 shows the redshift-dependent trajectories of α (Figure 11(a)) and β (Figure 11(b)) for Model 269. As before, Model 269 and the best-fit global model are in general agreement in terms of the predicted evolution of α and β with redshift, with two key exceptions: (1) Model 269 predicts that α flattens above $z\gtrsim 2$, a regime not well constrained by our data; and (2) Model 269 predicts that β increases more slowly with redshift than our data reveal.

The good agreement between the Model 269 and canonical model trends allows for a physical interpretation of the rapid increase in α and mild increase in β observed in our data. The close similarities between Model 269 and the F13a best model, Model 245, imply that the conclusions drawn here are the same basic conclusions drawn by F13a: the rise in α with redshift arises due to the shifting to higher-mass LMXB donor stars, and more luminous LMXBs, as the population age declines with increasing redshift. The predicted turn over in α at $z\gtrsim 2.5$ occurs when the universe was ≈2.5 Gyr old, and the stellar populations were on average ≈1 Gyr in age, which is the peak stellar age for LMXB emission. Therefore, the turn over in α around $z\approx 2.5$ corresponds to the peak stellar age for LMXB formation. Populations with 0.5–1.5 Gyr age ($z\quad =$ 2–5) are expected to produce comparably bright LMXB emission per unit ${M}_{\star }$, due to contributions from massive (≈3 ${M}_{\odot }$) red giant (RG) donors (see, e.g., Kim et al. 2009). ${L}_{{\rm{X}}}$(LMXB)/${M}_{\star }$ then declines rapidly at ages above ≈1.5 Gyr (i.e., $z\lesssim 2.5$), as the brightest RG donor star mass shifts to lower masses ($\lesssim$ 1–3 ${M}_{\odot }$). Some studies of LMXBs in local elliptical galaxies of different ages have provided initial support for this prediction (e.g.,Kim & Fabbiano 2010; Lehmer et al. 2014); however, no statistically conclusive results have been reached based on how LMXBs evolve with stellar-population age using local galaxies (e.g., Boroson et al. 2011; Zhang et al. 2012).

For HMXBs, mild evolution is predicted for β, continuing to high redshifts, due to a decline in metallicity with redshift. As argued by Linden et al. (2010) and F13a, star formation under low-metallicity conditions can yield significantly larger numbers of compact objects (neutron stars and black holes), compared with higher metallicities, due to a reduction in the efficiency of stellar-wind mass loss over the lifetimes of massive stars. Additionally, massive stars at lower metallicity tend to expand to large radii later in their evolution compared to higher metallicity ones. Hence, lower metallicity stars form better defined and more massive cores before entering the common envelope. This in turns allows for the easier ejection of the common-envelope of stars that will produce black holes (Linden et al. 2010; Justham et al. 2015). This increased pool of black holes in low-metallicity environments results in an excess of HMXBs (and thus HMXB emission) per unit SFR compared to higher metallicity environments. Furthermore, the most massive black holes formed in low metallicity environments can be quite large, increasing the potential for the X-ray emission from very luminous X-ray sources to vastly exceed the collective emission from much more numerous low-luminosity XRBs. These sources are thought to provide important contributions to heating the intergalactic medium at $z\gtrsim 10$ (e.g., Mirabel et al. 2011; Fragos et al. 2013b; Kaaret 2014). Studies of the metallicity dependence of ultraluminous X-ray source (ULX) formation and ${L}_{{\rm{X}}}$/SFR have demonstrated this effect empirically, and find dependences similar to that expected from the F13a population-synthesis predictions (e.g., Mapelli et al. 2009, 2010; Basu-Zych et al. 2013a; Brorby et al. 2014; Douna et al. 2015). These findings, along with the results presented in this paper, therefore support the idea that XRB emission was enhanced in the primordial $z\gtrsim 10$ universe, and could provide a non-negligible contribution to the heating of the intergalactic medium.

7.2. Implications for the Cosmic X-Ray Emissivity and X-Ray Background Contribution

Given the potential for normal-galaxy populations to be substantial contibutors to heating of the intergalactic medium at high redshifts, we use our measurements of the evolution of α and β, along with estimates of galaxy stellar mass and SFR density evolution, to provide updated constraints on the evolution of the X-ray emissivity. The stellar mass and SFR density of the universe has been constrained by numerous past studies using radio-to-far-UV emission. The recent review by Madau & Dickinson (2014) provides a coherent census of these constraints and presents analytic formulae describing the redshift evolution of the stellar mass density, ${\rho }_{\star }$, and SFR density, ψ, of the universe that are valid out to $z\approx 8$. Using the Madau & Dickinson (2014) formalism, corrected to our adopted Kroupa (2001) IMF, we estimated the redshift-dependent LMXB and HMXB emissivity as a function of redshift following:

$$\begin{eqnarray*}&&{\rho }_{{\rm{X}}}^{{\rm{LMXB}}}=\alpha (z){\rho }_{\star }(z),\end{eqnarray*}$$

$$\begin{eqnarray}&&{\rho }_{{\rm{X}}}^{{\rm{HMXB}}}=\beta (z)\psi (z).\end{eqnarray} \tag{ 16 }$$

Equation (16) can be applied using α and β values constrained either in individual redshift intervals or via global parameterizations. Figure 13 presents constraints from the individual redshift intervals for the 0.5–2 keV and 2–10 keV bands as filled circles with 1σ error bars. These values were compared with the LMXB and HMXB emissivity predictions from Model 269, extended from z = 0–8 (dashed curves in Figure 13). We must also account for hot gas emission, in particular at 0.5–2 keV, to compare with the data. We computed a first-order estimate of the hot gas emissivity of the universe at all redshifts by assuming a universal hot gas scaling with SFR, of the form:

$$\begin{eqnarray}{\rho }_{0.5-2\;\mathrm{keV}}^{\mathrm{hot}\ \ \mathrm{gas}}/(\mathrm{erg}\;{{\rm{s}}}^{-1}\;{\mathrm{Mpc}}^{-3}) & = & 1.5\times {10}^{39}\psi (z)/\\ & & {\rm{(}}{M}_{\odot }\;{\mathrm{yr}}^{-1}\;{\mathrm{Mpc}}^{-3}),\end{eqnarray} \tag{ 17 }$$

which is based on the Mineo et al. (2012b) relation. We then utilized the hot gas component of our canonical SED (see Figure 4(a)) to scale Equation (17) to the 2–10 keV band; we find,

$$\begin{eqnarray}{\rho }_{2-10\;\mathrm{keV}}^{\mathrm{hot}\ \ \mathrm{gas}}/(\mathrm{erg}\;{{\rm{s}}}^{-1}\;{\mathrm{Mpc}}^{-3}) & = & 3.4\times {10}^{38}\psi (z)/\\ & & {\rm{(}}{M}_{\odot }\;{\mathrm{yr}}^{-1}\;{\mathrm{Mpc}}^{-3}).\end{eqnarray} \tag{ 18 }$$

Figure 13 indicates the breakdown of LMXB, HMXB, and hot gas contributions to the X-ray luminosity density of normal galaxies. All three components provide substantial contributions to the 0.5–2 keV normal galaxy emissivity, with XRBs dominating at all redshifts for both bands. Similar to the results from F13a, LMXBs dominate the X-ray emissivity at $z\;\lesssim$ 1–2, with HMXBs dominating at $z\;\gtrsim$ 1–2, as the stellar mass density of the universe drops precipitously. For comparison, Figure 13 displays the AGN emissivity derived in Aird et al. (2015); these curves are shown to $z\approx 5$, where direct observational constraints are available. Reasonable extrapolations of the AGN curves to $z\gt 5$ indicate that the normal galaxy emissivity should overpower AGN at $z\;\gtrsim$ 6–8, during the epoch of reionization. This conclusion was also reached by Fragos et al. (2013b), based on Model 245, and has important implications for the role of ionizing photons from XRBs in heating the neutral intergalactic medium during the reionization epoch (see, e.g., discussions in Mirabel et al. 2011; McQuinn 2012; Pacucci et al. 2014; Artale et al. 2015).

Zoom In Zoom Out Reset image size

As noted by Dijkstra et al. (2012), the measured unresolved cosmic X-ray background (CXB) provides a firm upper limit on the possible evolution of X-ray scaling relations. From the above estimates of the X-ray emissivity evolution of the universe, it is straightforward to calculate the contributions to the observed CXB from normal galaxy populations throughout the universe. The normal galaxy CXB intensity in bandpass ${E}_{1}-{E}_{2}$, ${{\rm{\Omega }}}_{{E}_{1}-{E}_{2}}^{{\rm{CXB,}}\;{\rm{gal}}}$, can be calculated following:

$$\begin{eqnarray}&&{{\rm{\Omega }}}_{{E}_{1}-{E}_{2}}^{{\rm{CXB,}}\;{\rm{gal}}}={\int }_{0}^{{z}_{{\rm{max}}}}{k}_{{E}_{1}^{\prime }-{E}_{2}^{\prime }}^{{E}_{1}-{E}_{2}}\displaystyle \frac{{\rho }_{{E}_{1}^{\prime }-{E}_{2}^{\prime }}}{4\pi {d}_{L}^{2}}\displaystyle \frac{{dV}}{{dzd}{\rm{\Omega }}}{dz},\end{eqnarray} \tag{ 19 }$$

where ${k}_{{E}_{1}^{\prime }-{E}_{2}^{\prime }}^{{E}_{1}-{E}_{2}}$ provides a k-correction from rest-frame ${E}_{1}^{\prime }-{E}_{2}^{\prime }$ to observed-frame ${E}_{1}-{E}_{2}$ following our canonical SED (defined in Section 4), $\tfrac{{dV}}{{dzd}{\rm{\Omega }}}$ is the cosmology-dependent differential volume element (comoving volume per unit redshift per unit solid angle), and we integrate to ${z}_{{\rm{max}}}=8$.

As reported in Lehmer et al. (2012), the total unresolved X-ray background plus resolved normal-galaxy emission from the ≈4 Ms CDF-S survey has sky intensities of ≈$2.2\times {10}^{-12}$ erg cm⁻² s⁻¹ deg⁻² and ≈$3.2\times {10}^{-12}$ erg cm⁻² s⁻¹ deg⁻² for the observed-frame 0.5–2 keV and 2–8 keV bands, respectively. These values are, respectively, 27% and 18% of the total CXB intensity and constitute the maximum intensities that normal galaxies could contribute to the CXB in these bands. Using Equation (19) with our estimates of the redshift evolution of ${\rho }_{0.5-2\mathrm{keV}}$ and ${\rho }_{2-10\mathrm{keV}}$ (Equations (16)–(18)), we derive normal-galaxy intensities of ${{\rm{\Omega }}}_{0.5-2\;\mathrm{keV}}^{{\rm{CXB,}}\;{\rm{gal}}}\approx 1.0\times {10}^{-12}$ erg cm⁻² s⁻¹ deg⁻² and ${{\rm{\Omega }}}_{2-8\;\mathrm{keV}}^{{\rm{CXB,}}\;{\rm{gal}}}\approx 1.1\times {10}^{-12}$ erg cm⁻² s⁻¹ deg⁻², which are ≈46% and ≈35% of the maximum possible. We therefore find that the redshift evolution of X-ray scaling relations are fully consistent with constraints set by the CXB intensity.

Using the ≈6 Ms depth data in the Chandra Deep Field-South, we have studied the evolution of X-ray emission in normal galaxy populations since $z\approx 7$. Our findings can be summarized as follows:

• 1.  
  Scaling relations involving not only SFR, but also stellar mass (${M}_{\star }$) and redshift (z) provide a significantly better characterization of normal-galaxy X-ray luminosity (${L}_{{\rm{X}}}$) than a single, "universal" ${L}_{{\rm{X}}}$/SFR relation, which has been widely assumed in the literature.
• 2.  
  We deduce that emission from low-mass X-ray binary (LMXB) and high-mass X-ray binary (HMXB) populations drive correlations involving ${M}_{\star }$ and SFR, respectively. We find that out to at least $z\approx 2.5$, the 2–10 keV emission, which probes directly emission from XRB populations, evolves as ${L}_{2-10\mathrm{keV}}{({\rm{LMXB}})/{M}_{\star }\propto (1+z)}^{2-3}$ and ${L}_{2-10\mathrm{keV}}({\rm{HMXB}})/{\rm{SFR}}\propto (1+z)$. This evolution is consistent with basic predictions from XRB population-synthesis models, which indicate that the increases in LMXB and HMXB scaling relations with redshift are primarily due to effects related to declining stellar ages and metallicities, respectively (see Section 7.1 and F13a for details). However, the best-fit XRB population synthesis model is only marginally consistent with our data (at the ≈7% confidence level), calling for minor revisions of the current suite of models.
• 3.  
  When convolving the redshift-dependent scaling relations ${L}_{2-10\mathrm{keV}}({\rm{LMXB}})/{M}_{\star }$ and ${L}_{2-10\mathrm{keV}}({\rm{HMXB}})/{\rm{SFR}}$ with respective redshift-dependent stellar mass and SFR density curves from the literature, we find that LMXBs are likely to dominate the X-ray emissivity of normal galaxies out to $z\quad \approx$ 1–2, with HMXBs dominating at higher redshifts. We find that XRB populations dominate the X-ray emissivity of normal galaxies in both the 0.5–2 keV and 2–10 keV bands; however, hot gas is inferred to provide significant contributions to the 0.5–2 keV emissivity over the majority of cosmic history.
• 4.  
  The total X-ray emissivity from normal galaxies peaks around $z\quad \approx$ 1.5–3, similar to the SFR density of the universe. However, the X-ray emissivity of normal galaxies declines more slowly at $z\gtrsim 3$ than the SFR density due to the continued rise in ${L}_{{\rm{X}}}$/SFR scaling relation with redshift. Extrapolation of our results suggests that normal galaxies will provide an X-ray emissivity that exceeds that of AGN at $z\;\gtrsim$ 6–8, consistent with XRB population-synthesis predictions.

The results presented in this paper could be significantly expanded upon by future investigations with Chandra and planned observatories. In particular, some of the results and interpretations provided in this paper are in need of verification using data for local galaxies. For example, the rise in $\alpha \equiv {L}_{2-10\mathrm{keV}}({\rm{LMXB}})/{M}_{\star }$ and $\beta \equiv {L}_{2-10\mathrm{keV}}({\rm{HMXB}})/{\rm{SFR}}$ scaling relations with cosmic time is thought to be associated with changes in the mean stellar age and metallicity of the universe. However, attempts to prove that α and β vary with stellar age and metallicity in local galaxies, e.g., using Chandra resolved XRB populations, have not yet yielded statistically robust results (see Section 7.1 for discussion of such studies). This situation is due to challenges in (1) measuring accurate stellar ages and metallicities for galaxies, and (2) obtaining Chandra observations for statistically significant samples of galaxies that span broad ranges of these parameters. Future observations with Chandra (and to a lesser extent XMM-Newton) are needed to target significant numbers of galaxies with stellar ages and metallicities that are well constrained by optical/near-IR spectroscopic data and span broad ranges of these parameters.

In addition to local galaxy studies, future stacking investigations of distant galaxy populations separated by all of the most relevant physical parameters (i.e., SFR, ${M}_{\star }$, metallicity, and stellar ages) would allow for direct measurements of how physical properties influence XRB formation, independent of redshift. Statistically significant samples of galaxies could be drawn from the combination of wide and deep Chandra surveys (see, e.g., Brandt & Alexander 2015 for a review).

Future X-ray missions will provide significant gains in our knowledge of XRB populations in extragalactic environments for both the nearby and distant universe. For example, Athena (e.g., Nandra et al. 2013) will provide new spectral and timing constraints for XRBs in a variety of environments within the local galaxy population. Wide-area Athena surveys are planned that will reach depths comparable to the most sensitive regions of a ≈1 Ms Chandra survey, but over ≈1–10 deg² regions (see, e.g., Aird et al. 2013). Such a survey would yield 10,000–100,000 normal galaxy detections with most sources being at $z\gtrsim 0.5$ (see Lehmer et al. 2012). Significant clarification of the detailed evolution of XRB populations would come from future missions with imaging capabilities that supercede those of Chandra, and provide direct detections of X-ray galaxies to high redshifts. The X-ray Surveyor mission concept (e.g., Weisskopf et al. 2015) would reach depths of a few $\times {10}^{-19}$ erg cm⁻² s⁻¹ in the SB for a 4 Ms depth image, and would yield normal galaxy sky densities of ∼500,000 deg⁻¹. Such flux levels are a factor of $\gtrsim 10$ times lower than the mean stacked fluxes for this study, suggesting that most galaxies that are currently studied through stacking would be detected individually by X-ray Surveyor.

We thank the referee for their thorough reading of the manuscript and helpful comments, which have improved the quality of this paper. We thank Myrto Symeonidis and James Aird for graciously sharing data, which have been valuable for comparing with our results. This work has made use of the Rainbow Cosmological Surveys Database, which is operated by the Universidad Complutense de Madrid (UCM), partnered with the University of California Observatories at Santa Cruz (UCO/Lick, UCSC).

BDL gratefully acknowledges financial support from Chandra X-ray Center (CXC) grant GO4-15130B and NASA ADAP grant NNX13AI48G. AEH and ABZ acknowledge funding through CXC program GO4-15130Z and NASA ADAP grant 09-ADP09-0071. WNB and BL  thank CXC grant GO4-15130A and NASA ADP grant NNX10AC99G. TF acknowledges support from the Ambizione Fellowship of the Swiss National Science Foundation (grant PZ00P2_148123). FEB acknowledges support from CONICYT-Chile (Basal-CATA PFB-06/2007, FONDECYT Regular 1141218, "EMBIGGEN" Anillo ACT1101), the Ministry of Economy, Development, and Tourism's Millennium Science Initiative through grant IC120009, awarded to The Millennium Institute of Astrophysics, MAS. YQX acknowledges support of the Thousand Young Talents program (KJ2030220004), the 973 Program (2015CB857004), the USTC startup funding (ZC9850290195), the National Natural Science Foundation of China (NSFC-11473026, 11421303), the Strategic Priority Research Program "The Emergence of Cosmological Structures" of the Chinese Academy of Sciences (XDB09000000), and the Fundamental Research Funds for the Central Universities (WK3440000001).