A CORRELATION BETWEEN STAR FORMATION RATE AND AVERAGE BLACK HOLE ACCRETION IN STAR-FORMING GALAXIES

We selected Herschel-observed star-forming galaxies in the 9 deg² Boötes field covered by the NOAO Deep Wide-Field Survey (NDWFS; Jannuzi & Dey 1999) to measure the correlation between the SFR and the BHAR in star-forming galaxies. For the redshifts in this study, we primarily used the spectroscopic redshifts in the range 0.25 < z < 0.8 from the AGES Data Release 2 (Kochanek et al. 2012), which covers 7.7 deg² of NDWFS. To maximize the completeness for our sample of IR-selected SF galaxies, we also supplemented the data with photometric redshifts measured using the data from the 8.5 deg² Spitzer ISS (Eisenhardt et al. 2004; Stern et al. 2005) and the 10 deg² SDWFS (Ashby et al. 2009). The photometric redshifts were derived using all four IRAC bands (3.6, 4.5, 5.8, and 8 μm) in ISS and SDWFS, with the algorithm developed by Brodwin et al. (2006, B06 hereafter). We limited our photometric redshifts to the same redshift range as the spectroscopic redshifts, in which the accuracy of the photometric redshifts is σ = 0.06(1 + z) for 95% of galaxies and σ = 0.12(1 + z) for 95% of AGNs (Brodwin et al. 2006).

In addition to the IRAC observations, we also used far-IR and mid-IR data from Herschel and the Multiband Imaging Photometer (MIPS) on board Spitzer, respectively. The far-IR data in this work are based on the publicly available Herschel SPIRE 250 μm observations from the Herschel Multi-tiered Extragalactic Survey (HerMES; Oliver et al. 2012). We re-reduced and mosaiced the Boötes SPIRE observations (S. Alberts et al., in preparation), which include a deep ∼2 deg² inner region in the center of the field and a shallower ∼8.5 deg² outer region. We specifically focused on removing stripping, astrometry offsets, and glitches missed by the standard pipeline reduction. We also convolved the raw maps with a matched filter (see Chapin et al. 2011), which aided in source extraction by lowering the overall noise and de-blending sources. From this, we generated a matched filter catalog with 21,892 point sources above S/N >5. Completeness simulations showed that these catalogs are 95% complete in the inner region and 69% complete in the outer regions above a flux limit of 20 mJy. We also found minimal flux boosting for low signal-to-noise ratio (S/N) sources above these flux cutoffs. In addition, we used the 24 μm flux measurements available from the MIPS for Spitzer GTO observations (IRS GTO team, J. Houck (PI), and M. Rieke) of the Boötes field as a comparison in the source matching and SFR estimation. This catalog covers an area similar to that of the XBoötes and has 52,089 S/N >5 sources with flux >0.15 mJy.

X-ray observations from the XBoötes survey provide the basis of our measurements of AGN luminosity. The X-ray data in this study are drawn from the 9.3 deg² XBoötes survey, which is a mosaic of 126 short (5 ks) Chandra ACIS-I images (Murray et al. 2005; Kenter et al. 2005) covering the entire NDWFS. XBoötes contains 2724 X-ray point sources with at least four counts in the AGES survey region. 362 of the 2724 X-ray point sources are not close to bright stars and matched within 35 to objects with good AGES redshifts at 0.25 < z < 0.8 (Kenter et al. 2005; Hickox et al. 2009). These X-ray point sources have 0.5–7 keV luminosities of 10⁴² < L_X < 10⁴⁵ erg s⁻¹, which are characteristic for moderate- to high-luminosity AGNs and significantly larger than the typical L_X for star-forming galaxies (Ranalli et al. 2003).

To properly estimate the average SMBH accretion rate, it is important to account for obscured AGNs that may have their X-ray flux significantly depressed by photoelectric absorption (e.g., Alexander et al. 2008). We therefore supplemented our AGN sample with AGNs selected with the IRAC color–color selection criteria (Stern et al. 2005) in ISS and SDWFS. These selection criteria has been shown to effectively identify obscured AGNs at moderate redshift (Gorjian et al. 2008; Assef et al. 2010), since mid-IR wavelengths are not as strongly affected by obscuration as the optical or UV (e.g., Lacy et al. 2004; Stern et al. 2005; Donley et al. 2007; Hickox et al. 2007; Goulding & Alexander 2009). In our sample, 1047 AGNs were identified using mid-IR observations (mid-IR AGNs hereafter). Among the 1047 mid-IR AGNs and the 362 X-ray AGNs, 163 of them can be identified as AGNs using both X-ray and mid-IR selection criteria.

Since the large point-spread function (PSF) of SPIRE can lead to spurious matching results, here we discuss our catalog matching procedures. We first matched the entire SPIRE catalog to the 5σ MIPS 24 μm catalog with a matching radius of 5''. We found that ∼80% of the SPIRE sources in the coverage of the MIPS catalog have MIPS counterparts within the 5'' radius. While increasing the radius to 10'' can increase the MIPS detection fraction for the SPIRE sources to ∼92%, the fraction of SPIRE sources with multiple MIPS counterparts would also increase from 2.9% to 14%; therefore, we chose the 5'' matching radius to avoid spurious matching. These matched sources have a minimum 24 μm flux at least 25% larger than the 0.15 mJy flux limit of the 5σ MIPS 24 μm catalog, which ensures that the matched galaxies are star-forming galaxies that are bright in both mid-IR and far-IR, and that our completeness in 250 μm is not strongly affected by the 24 μm flux limit. Then, we matched the coordinates from the MIPS catalog (with 250 μm counterparts) to the B06 photometric redshift catalog with a matching radius of 2''. We also matched the AGES catalog and the AGES-matched optical positions of X-ray AGNs to the photometric redshift catalog with a matching radius of 1''.

To minimize spurious matches, we tested our matching by offsetting the source positions by 1' in a random direction. We found that with radius of 1'', our matching between the AGES and B06 catalog yielded <0.1% spurious matches. The radius of the matching between the 250 μm matched MIPS catalog and the B06 catalog is tested to have less than 4% of spurious matches, which greatly reduced the ∼25% spurious matching rates obtained when directly matching the SPIRE sources to the B06 catalog.

In our sample, there are 1767 galaxies with both SPIRE and MIPS detections. 1112 of these have spectroscopic redshift measurements, while the remainder have photometric redshifts. 121 of the 1767 sources (∼7%) in our sample have been classified as AGNs, in which 34 (∼2%) sources are identified as X-ray AGNs using X-ray selection criteria and 107 (∼6%) are identified as mid-IR AGNs using mid-IR color–color selection criteria. 20 of the 121 AGNs are identified as both X-ray AGNs and mid-IR AGNs. The angular distribution of sources and the coverage of different observations are plotted in Figure 1.

Zoom In Zoom Out Reset image size

In this section, we discuss the methods we used to estimate the SFRs with the SPIRE 250 μm observations.

Extrapolation of total IR luminosities (L_IR, defined as the integrated luminosity in the 8–1000 μm range) from monochromatic fluxes in the near- to mid-IR wavelengths using templates of infrared spectral energy distributions (SEDs) of local star-forming galaxies (e.g., Chary & Elbaz 2001, CE01 hereafter) has been adopted by a number of previous works to estimate the SFR of galaxies. However, it is also shown that AGNs have significant emission at near- to mid-IR wavelengths (e.g., Daddi et al. 2007; Mullaney et al. 2011), which poses challenges in obtaining reliable estimates of the SFR using mid-IR data.

By contrast, far-IR observations have been shown to suffer minimum contamination from AGNs (e.g., Netzer et al. 2007; Hatziminaoglou et al. 2010; Mullaney et al. 2011; Kirkpatrick et al. 2012), which makes far-IR emission a better tracer of the SF-related L_IR in an AGN-hosting galaxy sample. However, recent studies have discovered that estimates using only one SPIRE band and templates based on local star-forming galaxies can substantially overestimate L_IR for z < 1.5 (Elbaz et al. 2010; Nordon et al. 2012) due to the possible lower dust temperature (e.g., Dale & Helou 2002; Pope et al. 2006) that was not accounted for in the local templates. Therefore, in this work, we adopted the composite SED template from Kirkpatrick et al. (2012, hereafter K12). The z ∼ 1 K12 template is derived from the Spitzer and Herschel observations of a sample of star-forming galaxies at 0.4 < z < 1.47, which is comparable to the redshift range of our sample. From the photometric observations available in both K12 and our sample, we have determined that the distributions of the ratio between the observed 250 μm flux (S₂₅₀) and 24 μm flux (S₂₄) in both samples are similar. We found that ∼96% of the star-forming galaxies in our sample have S₂₅₀/S₂₄ that lies within the 2σ range (0.34 dex, adopted from Table 3 in K12) of the K12 template S₂₅₀/S₂₄ distribution derived in the redshift range of our sample. This shows that the K12 template can describe the far-IR to mid-IR color in our sample well, thus we can use this template to estimate the total L_IR in our sample.

In principle, the far-IR part of an SED for star-forming galaxies, which comprises the bulk of the SF-related L_IR, are dominated by the thermal radiation due to cold and warm dust. Thus, the ratio between the monochromatic far-IR flux and the total L_IR of the SED should be very similar for star-forming galaxies with similar dust temperatures. In particular, it has also been shown that the star-forming galaxies with and without strong AGNs have similar cold dust temperature (Kirkpatrick et al. 2012). Recent studies using Herschel observations have also shown that even for AGN host galaxies, the far-IR (⩾100 μm) emissions are still dominated by the cold dust component (e.g., Hatziminaoglou et al. 2010; Kirkpatrick et al. 2012). Thus, we can estimate the SF-related L_IR for our sample by normalizing the total infrared luminosity of this template ($L_{\rm IR}^{\rm T}$) using the 250 μm observations. For each source in our sample, we calculated the ratio between the observed S₂₅₀ and the monochromatic flux of the template at the corresponding observed-frame 250 μm ($S_{250}^{\rm T}$), and derived the total L_IR with the following equation:

$$\begin{equation} L_{\rm IR}=\frac{S_{250}}{S_{250}^{\rm T}}L_{\rm IR}^{\rm T}. \end{equation} \tag{ 1 }$$

For the K12 template we chose for this work, $L_{\rm IR}^{\rm T}=4.2\times 10^{11}\,L_\odot$, which corresponds closely to the median L_IR of our sample.

In this work, the L_IR was derived by assuming that the S₂₅₀/L_IR in our sample is similar to that of the K12 composite. For a population of star-forming galaxies, the observed dispersion in S₂₅₀/L_IR depends on the variations in their SED shapes. This dispersion also depends on the redshift since the observed S₂₅₀ traces the SED at different rest-frame wavelengths. Thus, the uncertainty in our L_IR measurement can be estimated from the variation in the SED shapes between individual galaxies. However, we cannot directly measure the SED shapes for our Boötes sample, because there is only one photometric band of far-IR observations available. Therefore, we took the 39 z ∼ 1 sources from K12, and selected a sub-sample of 25 SF-dominated (with mid-IR AGN fraction less than 10%, see K12 for the details in mid-IR spectral decomposition), far-IR-detected (with at least two bands of SPIRE photometry) sources. This sub-sample spans a redshift range of 0.47 < z < 1.24 and has at least five photometric data points from 24 to 500 μm; thus, it can be used to estimate the dispersion of S₂₅₀/L_IR due to the variation in SEDs for a population of star-forming galaxies.

We first estimated L_IR and the shape of SED for each source in the K12 sub-sample by taking the available photometry and calculated the corresponding rest-frame monochromatic luminosity (L_λ). Combining L_λ for each photometric bands, we calculated a best-fitting spline curve, then integrated along the spline curve in the rest-frame 20–300 μm range. For wavelengths beyond the longest wavelength of the spline curve, we used a linear interpolation by assuming that this part of the SED follows a Rayleigh–Jeans distribution. Since this sub-sample from K12 is selected to be SF-dominated and have far-IR constraints with at least two photometric bands, the best-fitting spline curves can trace the simple shapes of SF-related SEDs at this wavelength range. The integrated 20–300 μm luminosity (L_20–300) probes the thermal radiation from warm and cold dust, and represents the bulk of the total 8–1000 μm L_IR (e.g., for the K12 template, L_20–300 ∼ 0.91L_IR), thus we can use L_20–300 as a good proxy of L_IR. Along the best-fitting spline curves, we calculated the observed-frame 250 μm fluxes at the redshift range of our Boötes sample, 0.25 < z < 0.8, which corresponds to rest-frame wavelengths from 200 to 140 μm. The S₂₅₀/L_20–300 for each source in the K12 sub-sample can therefore be estimated as a function of redshift. For the K12 composite SED, we also calculated the S₂₅₀/L_20–300 in 0.25 < z < 0.8. We compared the S₂₅₀/L_20–300 for each source to that of the K12 template at 0.25 < z < 0.8, and found that the standard deviations in the differences of S₂₅₀/L_20–300 between the K12 sources and the K12 template is ∼0.17 dex at z = 0.25 and ∼0.10 dex at z = 0.8. This shows that the deviations in far-IR spectral shapes across a representative population of star-forming galaxies are reasonably small, and thus we can confidently estimate the SF-related L_IR using Equation (1) from the observed monochromatic 250 μm flux. Even though the dispersion is lower at higher redshift, we conservatively chose 0.17 dex as the uncertainty in our L_IR estimation. This is also consistent with the 0.17 dex uncertainty inferred by K12 for the $L_{\rm IR}^{\rm T}$ of the composite template.

As a check, we compared the L_IR for star-forming galaxies from our method ($L_{\rm IR}^{250}$) with the L_IR obtained by directly fitting the Spitzer MIPS 24 μm flux to the CE01 library ($L_{\rm IR}^{24}$), which has been shown to be able to robustly estimate L_IR for star-forming galaxies out to z ∼ 1 (Magnelli et al. 2009), and found that the average difference between $L_{\rm IR}^{250}$ and $L_{\rm IR}^{24}$ is ∼0.1 dex. However, we note that the difference is much larger for objects with significant AGN contamination (and thus bluer S₂₅₀/S₂₄), highlighting the need for far-IR data to measure SFR in AGN hosts.

The SFRs for our sample were derived from L_IR using the relation from Kennicutt (1998), modified for a Chabrier initial mass function (IMF; Chabrier 2003; Salim et al. 2007):

$$\begin{equation} \frac{{\rm SFR}}{M_{\odot }\;{\rm yr}^{-1}}=1.09\times 10^{-10} \left(\frac{L_{\rm IR}}{L_\odot }\right). \end{equation} \tag{ 2 }$$

Figure 2 shows the distribution in redshift and L_IR of our sample. The comparison between the photometric and spectroscopic samples, and the comparison between the L_IR of AGNs and star-forming galaxies are also shown as the normalized histograms on the side panels. We note that AGNs and star-forming galaxies have similar distributions in L_IR and redshift, suggesting that there is no apparent difference in SF properties between the AGN host galaxies and the star-forming galaxies in our sample.

Zoom In Zoom Out Reset image size

The large observed nuclear luminosities in AGNs are direct manifestations of BH growth through mass accretion (see, e.g., Alexander & Hickox 2012). In this work, we calculated the rest-frame 2–10 keV X-ray luminosity (L_X) for direct comparison with other studies. The k-corrections were calculated based on an X-ray spectral index of 1.7 and a galactic gas column density of N_H ∼ 10²⁰ cm⁻². We chose L_X as a proxy to estimate the BHAR:

$$\begin{equation} \dot{M}_{\rm BH}=0.15\frac{\epsilon }{0.1}\frac{22.4L_{\rm X}}{10^{45}\,{\rm erg\,s}^{-1}}\,M_{\odot }\;{\rm yr}^{-1}. \end{equation} \tag{ 3 }$$

For simplicity, $\dot{M}_{{\rm BH}}$ in Equation (3) is derived from L_X with a constant bolometric correction factor of 22.4 (the mean bolometric correction factor from Vasudevan & Fabian 2007, which is based on a sample of local, L_X = 10^41–46 erg s⁻¹, AGNs.) Here, is the mass–energy conversion efficiency (we used a typical value ∼ 0.1; see Marconi et al. 2004). In the following paragraphs, we describe the methods we used to calculate the average L_X for all of the star-forming galaxies in our sample.

2.6.1. Mid-IR AGNs

2.6.2. X-Ray Stacking of Star-forming Galaxies

2.6.3. Average X-Ray Luminosity

Combining the contributions of X-ray and mid-IR identified AGNs as well as undetected sources, the average L_X for star-forming galaxies can now be calculated in each bin of SFR:

$$\begin{eqnarray} {\left\langle {L_{\rm X}}\right\rangle } &=& \left[\sum _{i=1}^{N_{\rm XAGN}} \left(L_{\rm X}^{\rm XAGN} \right)_i+\sum _{i=1}^{N_{\rm IRAGN}} \left(L_{\rm X}^{\rm IRAGN} \right)_i+N_{\rm SFG} L_{\rm X}^{\rm stacking} \right] \nonumber\\ &&\times [N_{\rm XAGN}+N_{\rm IRAGN}+N_{\rm SFG}]^{-1}. \end{eqnarray} \tag{ 4 }$$

Here, N_XAGN, N_IRAGN, and N_SFG are the total numbers of X-ray identified AGNs, mid-IR identified AGNs without direct X-ray detections, and star-forming galaxies without identified AGNs, respectively. $L_{\rm X}^{\rm XAGN}$ and $L_{\rm X}^{\rm IRAGN}$ are the X-ray luminosities for individual X-ray AGNs and the mid-IR AGNs, and $L_{\rm X}^{\rm stacking}$ is the average X-ray luminosity from the stacking analysis for all star-forming galaxies without direct AGN detections in each bin.

We estimated our errors by propagating the observed uncertainties for $L_{\rm X}^{\rm XAGN}$, and the uncertainties for $L_{\rm X}^{\rm IRAGN}$ and $L_{\rm X}^{\rm stacking}$ estimated in Sections 2.6.1 and 2.6.2 with a bootstrap method. The uncertainties in L_IR were also taken into our bootstrap analysis. In each bootstrapping sub-sample, we first randomly re-sampled our sources with replacements, then replaced the original L_IR for each source using a random normal error with an 1σ value of 0.17 dex. We then re-binned the random sample using the same bins. For each bin, we recalculated the stacked L_X^stacking and the uncertainties for the sources that were not identified as AGNs, then replaced the L_X for every detected AGN in the bin with a new L_X within the normal error of the original AGN L_X. Finally, we recalculated the 〈L_X〉 and the average L_IR for each bin. We repeated the bootstrapping 5000 times, at which the variances in 〈L_X〉 and L_IR converge to finite values. The results are shown in Figures 3 and 4 and are discussed in the next section.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Using Equation (4), we can calculate the average L_X in each bin of SFR. However, it is well known that high-mass and low-mass X-ray binaries (HMXBs and LMXBs) can also generate X-ray luminosity that is correlated with SFR (e.g., Grimm et al. 2003; Ranalli et al. 2003; Gilfanov et al. 2004; Lehmer et al. 2010). To accurately estimate the SMBH accretion rate, we calculated the X-ray luminosities related to star-forming processes in each bin of SFR using the equation $L_{\rm X}^{\rm SF} = \alpha M_\star +\beta {\rm SFR}$, which is the SFR–L_X relation for HMXBs and LMXBs in Lehmer et al. (2010). In this equation, the stellar mass M_⋆ is only weakly correlated with SFR in active star-forming galaxies (e.g., SFR > 5 M_☉ yr⁻¹). Thus, for our sample and the Chabrier IMF we adopted, the equation can be re-written into $L_{\rm X}^{\rm SF} = 10^{26.4} {\rm SFR}^{0.3}\,{+}\,10^{39.3}{\rm SFR}$ (see Equation (3) in Symeonidis et al. 2011, for more details).

In addition, we also tested whether the limited volume of the sample could affect the BHAR–SFR correlation. AGNs with the highest luminosity are rare in this redshift range (and so might not be detected in our survey volume), but may contribute significantly to the 〈L_X〉 of our sample. We estimated the contribution of these rare, extremely luminous AGNs to the 〈L_X〉 using the X-ray luminosity function (XLF) from Aird et al. (2010). We note that the AGN XLF was not designed to represent the X-ray luminosity from the sources without direct X-ray observations, thus we first estimated the effect of limited volume on the sources that were identified as AGNs in our sample, i.e., L_X > 10⁴² erg s⁻¹; then, we calculated the effect on the full population by adjusting the result from detected AGNs based on the AGN detection fraction. In detail, at the average redshift range of each SFR bin of our sample, we first calculated the "intrinsic" average AGN luminosity by directly integrating the XLF at L_X > 10⁴² erg s⁻¹. Then, we estimated the "detected" average AGN luminosity by integrating the XLF with a high-end cutoff luminosity, at which the number of the detected AGNs in the volume of each SFR bin is ⩽1. We found that the difference between the "intrinsic" average AGN luminosity and the "detected" average AGN luminosity is 10%–3% from the first bin to the last bin of SFR in our sample. After the adjustments of the AGN detection fraction in each bin, the corrections on 〈L_X〉 would become 6.4%, 2.5%, 2.2%, and 2.9%, respectively. These corrections are small and do not make a notable difference to our study of the BHAR–SFR correlation in the large volume of the Boötes survey region, but might be important when calculating 〈L_X〉 for a sample with smaller volume. To accurately describe the correlation between the average SMBH growth and SF, we subtracted our 〈L_X〉 with $L_{\rm X}^{\rm SF}$, and also increased our observed 〈L_X〉 values to account for volume effects as described above. The BHAR were then derived using Equation (3).

We have determined that the average BHAR has a correlation to the SFR in our sample. The results are shown in Figure 3. The average X-ray luminosities, 〈L_X〉, as determined in Equation (4) are shown as the red circles. The observed L_X for the AGNs identified through X-ray or IRAC observations are shown as the stars, while the stacked L_X for star-forming galaxies without direct X-ray observations are shown as the downward triangles. For comparison, we present the SFR–L_X relation from Lehmer et al. (2010) in Figure 3. We also show the SFR–BHAR correlation corresponding to the local ratio of M_BH and M_bulge as the green, dashed line on the top of Figure 3. This $\dot{M}_{{\rm BH}}={\rm SFR}/500$ relation is directly derived from the M_BH/M_Bulge ratio found in Marconi et al. (2004). Since the average L_X of detected AGNs is subject to the flux limit in the observations, the fact that our data points for detected AGNs sit on the $\dot{M}_{\rm BH}={\rm SFR}/500$ relation is only a coincidence.

The relation between L_IR and 〈L_X〉 as shown in Figure 3 can be fitted with a linear relation given by

$$\begin{eqnarray} \log (L_{\rm X} [{\rm erg s}^{-1}]) &=& (30.37\pm 3.80) \nonumber\\ &&+(1.05\pm 0.33)\log (L_{\rm IR}/L_{\odot }), \end{eqnarray} \tag{ 5 }$$

which is calculated using the nonlinear least-squares fitting program MPFIT in IDL (Markwardt 2009). The uncertainties in the averages are estimated using bootstrap re-sampling in each bin. The reduced χ² of this fitting is 0.99. An equation correlating SFR and BHAR can immediately be derived using Equations (2), (3), and (5):

$$\begin{eqnarray} &&({\rm BHAR}/M_{\odot }\;{\rm yr}^{-1}) = 10^{(-3.72\pm 0.51)}\left(\frac{{\rm SFR}}{M_\odot \;{\rm yr}^{-1}}\right)^{(1.05\pm 0.33)}. \nonumber\\ \end{eqnarray} \tag{ 6 }$$

For an SFR of 100 M_☉ yr⁻¹, this corresponds to a ratio of BHAR to SFR of log (BHAR/SFR) ∼ −3.6 ± 0.2, where the error is derived from calculating this ratio in each bootstrapping sub-sample when deriving Equation (5).

We also calculated the same correlation using only the sample with spectroscopic redshifts. In this calculation, we adopted the sampling weight w_i to account for the spectroscopic redshift completeness of the AGES sample. The sampling weight is the combination of the sparse sampling weight that accounts for the random target selection incompleteness, the target assignment weight that addresses the fiber-allocation selection, and the redshift weight that accounts for the unsuccessful redshift measurement. The details of the sampling weight can be found in Section 3.1 of Hickox et al. (2009) and Kochanek et al. (2012). Using the same methods, the SFR–BHAR relation for the AGES galaxies with only spectroscopic redshifts can be written as

$$\begin{eqnarray} \log (L_{\rm X} [{\rm erg s}^{-1}]) &=& (29.39\pm 4.72) \nonumber\\ &&+(1.14\pm 0.41)\log (L_{\rm IR}/L_{\odot }). \end{eqnarray} \tag{ 7 }$$

Here, the power-law index is only higher by ∼0.1 compared to Equation (5), which is still consistent with the result from our main sample. Since photometric redshift measurements are not subject to the choices of sampling weights as the spectroscopic sample, we chose the result from Equations (5) and (6) as our primary conclusion in this work.

In Figure 3, we noticed that our stacked L_X is at least ∼0.7 dex higher than the Lehmer et al. (2010) $L_{\rm X}^{\rm SF}$ in the first three bins, suggesting that the L_X contribution from SF in our stacked X-ray luminosity is less than ∼20% in these bins. For the bin with the highest SFR, the stacked X-ray luminosity is still higher than $L_{\rm X}^{\rm SF}$ by a factor of two. This shows that while the galaxies we stacked were not identified as AGNs using the common X-AGN selection criterion of L_X > 10⁴² erg s⁻¹ or the IRAC color–color selection method, a significant fraction of the average L_X for these star-forming galaxies arises from SMBH accretions. We stress that the 〈L_X〉 used in our primary analysis is the average L_X due to SMBH accretion only, which was determined by subtracting the expected SF contribution. We note that in an X-ray stacking study in the same Boötes survey region, Watson et al. (2009) concluded that the spectroscopically selected late-type galaxies have their X-ray luminosities dominated by HMXBs. However, the sample in Watson et al. (2009) are galaxies with lower SFRs at lower redshifts. The stacked L_X in the lowest SFR bin in our work is still in agreement with their sample at comparable SFR. Using our methods, deeper far-IR observations would be required to probe the more moderate SFR galaxies studied by Watson et al. (2009).

We note that there might be observational bias in the BHAR-to-SFR correlation in our sample due to the flux limits. The limited flux would cause galaxies with high SFR to be preferentially found at higher redshift. This bias translates into the different average redshift in each bin (〈z〉 = 0.31, 0.46, 0.67, 0.68 in the four bins of SFR, from low to high.) Studies of redshift evolution of SFR density (e.g., Hopkins & Beacom 2006; Rodighiero et al. 2010) and BHAR density (e.g., Silverman et al. 2009; Aird et al. 2010) suggest that the SFR and BHAR densities are higher at higher redshift (up to z ∼ 2). Thus, when dividing galaxies into bins of SFR in a flux-limited sample, even if the L_X and L_IR are completely uncorrelated in individual galaxies, the redshift evolutions of the XLF and the infrared luminosity function (IRLF) would naturally yield an L_X–L_IR relation.

To account for the selection bias due to this effect and to test whether the redshift coevolution of the SFR density and the BHAR density is the dominant factor driving the L_X–L_IR correlation observed in this work, we examined the sample by directly computing the effect of the redshift evolution of the X-ray luminosity density (XLD; e.g., Aird et al. 2010, A10 hereafter). In principle, if there is no BHAR–SFR correlation in individual galaxies, then the observed difference of the 〈L_X〉 between the bins with the lowest and the highest SFR should be consistent with the pure redshift evolution of the XLD. We found that in our sample, the pure redshift evolution of the A10 XLD in the range of the average redshifts in our bins would translate into a difference in 〈L_X〉 of 0.47 dex.

To address this issue more carefully, we also created a "mock" catalog of galaxies, in which redshift distributions similar to our sample were generated. We generated IR luminosities for the galaxies in the mock catalog based on the IRLF from Rodighiero et al. (2010). Based on the different normalizations of this IRLF and the XLF from A10, we only assigned X-ray luminosities to a fraction of the mock galaxies according to the XLF from A10. For the rest of the galaxies, we assumed their X-ray luminosities to be zero. Since the IR luminosity distribution and X-ray luminosity distribution were derived independently, there is no intrinsic correlation between L_X and L_IR in our mock catalog. To test the effects of the redshift evolution in XLF and the possible Malmquist bias, we took the flux limits in X-ray and far-IR of our sample and applied them to the mock catalog, then repeated the calculations described in Section 3.1 to obtain 〈L_X〉 in bins of SFR. We found a weak correlation between L_IR and 〈L_X〉 in our mock catalog, log L_X = 37.45  ±  1.93 + (0.30  ±  0.17)log L_IR. The difference in 〈L_X〉 between the bins with the highest and the lowest SFR is 0.38 dex, which is similar to the effect of pure XLD evolution we estimated in the previous paragraph. In Equation (5), we found that there is at least 1.31 dex difference in the 〈L_X〉 of the bins with the lowest and the highest SFR in our sample, indicating that most, if not all, of our observed trend is due to the intrinsic correlation between SFR and BHAR.

To examine whether the average SFR–BHAR correlation is subject to the limiting fluxes of the observations, we compared our result with the sample of Herschel-selected star-forming galaxies in the pencil-beam Chandra Deep Field-North (CDF-N; Symeonidis et al. 2011) at redshift z ∼ 1. From Table 2 in Symeonidis et al. (2011), we selected the galaxies with hard (2–10 keV) X-ray detections and L_IR larger than 10¹¹ L_☉, in which the average X-ray luminosity for the X-ray non-detected galaxies have been estimated. For these luminous infrared galaxies (LIRGs) and ultraluminous infrared galaxies (ULIRGs), we calculated the average luminosity and the error in the 2–10 keV X-ray using a bootstrap re-sampling method similar to what we used for the Boötes data. Since there is no stacking signal in the ULIRG bin, we used the lower limit (L_X = 0) for the X-ray non-detected sources in that bin to estimate the error in 〈L_X〉 conservatively. We also estimated the effects of the limited volume in this field using a similar approach described in the second paragraph of Section 3.1, and found that the possible non-detections of the rare, high-luminosity AGNs might decrease the 〈L_X〉 by ∼26% in the LIRG bin, and ∼7% in the ULIRG bin. For the 〈L_X〉 in both our sample and the Symeonidis et al. (2011) sample, we subtracted by $L_{\rm X}^{\rm SF}$ and made adjustments to account for the bias due to limited volume. A comparison of the results are displayed in Figure 4, which shows that in samples of star-forming galaxies with different X-ray flux limits, even though the average L_X for the detected AGNs are different (so are the average L_X for the galaxies without direct X-ray detections), the average L_X to L_IR relations are consistent.