Evident black hole-bulge coevolution in the distant universe

TL;DR

This study provides evidence that supermassive black hole growth in bulge-dominated galaxies at high redshift is primarily linked to star formation, supporting a coevolution scenario that explains the tight local black hole-bulge mass correlation.

Contribution

It demonstrates a direct, significant correlation between SMBH accretion rate and SFR in bulge-dominated galaxies at $z=0.5-3$, highlighting coevolution in the distant universe.

Findings

SMBH accretion rate correlates with SFR in bulge-dominated galaxies.

The SMBH-bulge coevolution is evident at high redshift.

The $M_{BH}$-$M_{bulge}$ relation is consistent with lockstep growth.

Abstract

Observations in the local universe show a tight correlation between the masses of supermassive black holes (SMBHs; $M_{\rm BH}$) and host-galaxy bulges ($M_{\rm bulge}$), suggesting a strong connection between SMBH and bulge growth. However, direct evidence for such a connection in the distant universe remains elusive. We have studied sample-averaged SMBH accretion rate ($\overline{\rm BHAR}$) for bulge-dominated galaxies at $z=0.5-3$. While previous observations found $\overline{\rm BHAR}$ is strongly related to host-galaxy stellar mass ($M_\star$) for the overall galaxy population, our analyses show that, for the bulge-dominated population, $\overline{\rm BHAR}$ is mainly related to SFR rather than $M_\star$. This ${\overline{\rm BHAR}}$-SFR relation is highly significant, e.g. $9.0\sigma$ (Pearson statistic) at $z=0.5-1.5$. Such a $\overline{\rm BHAR}$-SFR connection does not exist among our comparison sample of galaxies that are not bulge-dominated, for which $M_\star$ appears to be the main determinant of SMBH accretion. This difference between the bulge-dominated and comparison samples indicates that SMBHs only coevolve with bulges rather than the entire galaxies, explaining the tightness of the local $M_{\rm BH}$-$M_{\rm bulge}$ correlation. Our best-fit $\overline{\rm BHAR}$-SFR relation for the bulge-dominated sample is $\log\overline{\rm BHAR} = \log\mathrm{SFR} - (2.48\pm0.05)$ (solar units). The best-fit $\overline{\rm BHAR}/\mathrm{SFR}$ ratio ($10^{-2.48}$) for bulge-dominated galaxies is similar to the observed $M_{\rm BH}/M_{\rm bulge}$ values in the local universe. Our results reveal that SMBH and bulge growth are in lockstep, and thus non-causal scenarios of merger averaging are unlikely the origin of the $M_{\rm BH}$-$M_{\rm bulge}$ correlation. This lockstep growth also predicts that the $M_{\rm BH}$-$M_{\rm bulge}$ relation should not have strong redshift dependence.

Figures

Tables

Table 1. Table 1: Summary of sample properties

Field	Area	Gal. #	Spec./Photo. #	Galaxy Ref.	B.-D. (X)	Comp. (X)	X. Dep.	X-ray Ref.
(1)	(2)	(3)	(4)	(5)	(6)	(7)	(8)	(9)
GOODS-S	170	1,504	727/777	Santini et al. (2015)	398 (100)	1,106 (241)	7 Ms	Luo et al. (2017)
GOODS-N	170	1,855	391/1,464	Barro et al., in prep.	483 (71)	1,372 (168)	2 Ms	Xue et al. (2016)
EGS	200	2,446	219/2,227	Stefanon et al. (2017)	591 (48)	1,855 (105)	800 ks	Nandra et al. (2015)
UDS	200	2,128	254/1,874	Santini et al. (2015)	549 (42)	1,579 (75)	600 ks	Kocevski et al. (2018)
COSMOS	220	2,369	10∗/2,359	Nayyeri et al. (2017)	603 (22)	1,766 (39)	160 ks	Civano et al. (2016)
Total	960	10,302	1,601/8,701	–	2,624 (283)	7,678 (628)	–	–

Table 2. Table 2: p 𝑝 p -values (significances) of partial-correlation analyses for the bulge-dominated (top) and comparison (bottom) samples.

Bulge-Dominated; $z=0.5\text{–}1.5$
Relation	Pearson	Spearman
$\overline{\mathrm{BHAR}}$-$M_{\star }$	$0.03$ ($2.2\sigma$)	$0.02$ ($2.3\sigma$)
$\overline{\mathrm{BHAR}}$-$\mathrm{SFR}$	${10}^{-18.8}$ ($9.0\sigma$)	${10}^{-13.3}$ ($7.5\sigma$)
Bulge-Dominated; $z=1.5\text{–}3$
Relation	Pearson	Spearman
$\overline{\mathrm{BHAR}}$-$M_{\star }$	$0.26$ ($1.1\sigma$)	$0.44$ ($0.8\sigma$)
$\overline{\mathrm{BHAR}}$-$\mathrm{SFR}$	${10}^{-28.5}$ ($11.2\sigma$)	${10}^{-5.2}$ ($4.5\sigma$)
Comparison; $z=0.5\text{–}1.5$
Relation	Pearson	Spearman
$\overline{\mathrm{BHAR}}$-$M_{\star }$	${10}^{-7.7}$ ($5.6\sigma$)	${10}^{-8.4}$ ($5.9\sigma$)
$\overline{\mathrm{BHAR}}$-$\mathrm{SFR}$	$0.01$ ($2.5\sigma$)	$0.02$ ($2.3\sigma$)
Comparison; $z=1.5\text{–}3$
Relation	Pearson	Spearman
$\overline{\mathrm{BHAR}}$-$M_{\star }$	${10}^{-14.3}$ ($7.8\sigma$)	${10}^{-4.5}$ ($4.2\sigma$)
$\overline{\mathrm{BHAR}}$-$\mathrm{SFR}$	$0.20$ ($1.3\sigma$)	$0.97$ ($0.0\sigma$)

Table 3. Table 3: Properties of each bin in Fig.  5 (left; black points).

Bulge-Dominated; $z=0.5\text{–}1.5$
$\log \mathrm{SFR}$	$\log M_{\star }$	$\log \overline{\mathrm{BHAR}}$	Gal. #	X. #	AGN #
(1)	(2)	(3)	(4)	(5)	(6)
$-1.97$	10.24	$-{4.35}_{-0.17}^{+0.17}$	360	14	0
$-1.24$	10.44	$-{3.67}_{-0.16}^{+0.14}$	360	20	1
$-0.71$	10.48	$-{2.82}_{-0.39}^{+0.20}$	360	29	6
$0.14$	10.25	$-{2.41}_{-0.15}^{+0.12}$	360	46	15
$1.05$	10.19	$-{1.69}_{-0.17}^{+0.12}$	361	69	36
Bulge-Dominated; $z=1.5\text{–}3$
$\log \mathrm{SFR}$	$\log M_{\star }$	$\log \overline{\mathrm{BHAR}}$	Gal. #	X. #	AGN #
(1)	(2)	(3)	(4)	(5)	(6)
$-0.90$	10.62	$-{3.27}_{-0.21}^{+0.15}$	274	7	0
$-0.01$	10.64	$-{2.28}_{-0.18}^{+0.13}$	274	21	4
$1.47$	10.59	$-{1.08}_{-0.11}^{+0.09}$	275	77	41

Table 4. Table 4: p 𝑝 p -values (significances) of partial-correlation analyses for all galaxies.

All; $z=0.5\text{–}1.5$
Relation	Pearson	Spearman
$\overline{\mathrm{BHAR}}$-$M_{\star }$	${10}^{-5.3}$ ($4.6\sigma$)	${10}^{-2.8}$ ($3.1\sigma$)
$\overline{\mathrm{BHAR}}$-$\mathrm{SFR}$	${10}^{-2.5}$ ($2.9\sigma$)	$0.03$ ($2.2\sigma$)
All; $z=1.5\text{–}3.0$
Relation	Pearson	Spearman
$\overline{\mathrm{BHAR}}$-$M_{\star }$	${10}^{-4.4}$ ($4.1\sigma$)	${10}^{-3.0}$ ($3.3\sigma$)
$\overline{\mathrm{BHAR}}$-$\mathrm{SFR}$	${10}^{-2.0}$ ($2.6\sigma$)	$0.04$ ($2.1\sigma$)

Equations

$$\begin{array}{r}\frac{\mathrm{SFR}}{M_{\odot }\mathrm{y}{\mathrm{r}}^{-1}}=1.09\times 10^{-10}\frac{L_{\mathrm{IR}}}{L_{\odot }}.\end{array}$$

$$\begin{array}{r}\overline{L_{\mathrm{bol}}}=\frac{{\Sigma }_{\det }(L_{\mathrm{X}}-L_{\mathrm{X},\mathrm{XRB}})k_{\mathrm{bol}}+(\overline{L_{\mathrm{X},\mathrm{stack}}}-\overline{L_{\mathrm{X},\mathrm{XRB}}})N_{\mathrm{non}}\overline{k_{\mathrm{bol}}}}{N_{\det }+N_{\mathrm{non}}}\end{array}$$

$$\begin{array}{rl}\overline{\mathrm{BHAR}}& =\frac{(1-ϵ)\overline{L_{\mathrm{bol}}}}{ϵc^2}\\ & =\frac{1.58\overline{L_{\mathrm{bol}}}}{10^{46}\mathrm{erg}{\mathrm{s}}^{-1}}{M}_{\text{sun}}{\mathrm{yr}}^{-1},\end{array}$$

$$\log \overline{\mathrm{BHAR}}=\{\begin{array}{ll}(0.88\pm 0.07)\log \mathrm{SFR}-(2.56\pm 0.08),& z=0.5\text{–}1.5\\ (0.89\pm 0.08)\log \mathrm{SFR}-(2.38\pm 0.09),& z=1.5\text{–}3\end{array}$$

$$\begin{array}{r}\log \overline{\mathrm{BHAR}}=(0.92\pm 0.04)\log \mathrm{SFR}-(2.47\pm 0.05),\end{array}$$

$$\begin{array}{r}\log \overline{\mathrm{BHAR}}=\log \mathrm{SFR}-(2.48\pm 0.05).\end{array}$$

Full text

††pubyear: 2018††pagerange: Evident black hole-bulge coevolution in the distant universe–Evident black hole-bulge coevolution in the distant universe

1 Introduction

One of the essential challenges of extragalactic astronomy is to understand the connection between supermassive black holes (SMBHs) and their host galaxies. It is well established that the masses of SMBHs ($M_{\rm BH}$) are tightly correlated with the stellar masses of host-galaxy classical bulges ($M_{\rm bulge}$) in the local universe (the Magorrian relation; e.g. Magorrian et al. 1998; Ferrarese & Merritt 2000; Gebhardt et al. 2000; Häring & Rix 2004; Kormendy & Ho 2013). The intrinsic scatter of the $M_{\rm BH}$-$M_{\rm bulge}$ correlation is only $\approx 0.3$ dex (Kormendy & Ho, 2013). This tight correlation is surprising considering that $M_{\rm BH}$ is only a tiny fraction (a few thousandths) of $M_{\rm bulge}$. Therefore, some fundamental connections between the growth of SMBHs and host-galaxy bulges likely exist over cosmic history. These physical connections are often termed as “SMBH-bulge coevolution”.

Tremendous observational efforts have been made to identify these mysterious connections. It has been found that the cosmic evolution of SMBH accretion rate (BHAR) density and star formation rate (SFR) density are broadly similar, both peaking at $z\sim 2$ (e.g. Aird et al., 2010, 2015; Mullaney et al., 2012; Madau & Dickinson, 2014; Ueda et al., 2014). However, from observations of active galactic nuclei (AGNs), the SFR is a relatively flat function of the observed BHAR at a given redshift (e.g. Harrison et al., 2012; Rosario et al., 2013; Azadi et al., 2015; Stanley et al., 2015; Lanzuisi et al., 2017; Stanley et al., 2017; Dai et al., 2018). This apparent lack of a strong SFR-BHAR connection might be caused by AGN variability. While star formation activity is stable on time scales longer than $\sim 100$ Myr, SMBH accretion could vary strongly on much shorter time scales ($\approx 10^{2}\text{--}10^{7}$ yr; e.g. Martini 2004; Kelly et al. 2010; Novak et al. 2011; Sartori et al. 2018). An intrinsic connection between SFR and long-term average BHAR might be hidden by this strong AGN variability.

To obtain long-term average BHAR, the ideal way is to observe a galaxy for at least millions of years, which is presently infeasible. Practically, it has been proposed to adopt sample-averaged BHAR ($\overline{\rm BHAR}$) as a proxy of long-term average BHAR (e.g. Chen et al., 2013; Hickox et al., 2014). Indeed, a positive $\overline{\rm BHAR}$-SFR connection has been observed (e.g. Chen et al., 2013; Hickox et al., 2014; Lanzuisi et al., 2017; Yang et al., 2017). However, Yang et al. (2017) show, via partial correlation analyses (PCOR), that $\overline{\rm BHAR}$ is actually more strongly related to host-galaxy total stellar mass ($M_{\star}$) than SFR (also see Fornasini et al. 2018 for a similar conclusion). Their results suggest that the apparent $\overline{\rm BHAR}$-SFR relation is only a secondary effect resulting from a primary $\overline{\rm BHAR}$-$M_{\star}$ relation and the star-formation main sequence. Yang et al. (2018b) further show that once $M_{\star}$ is carefully controlled, SMBH accretion is largely independent of the cosmic environment of the host galaxies, consistent with previous AGN clustering studies (e.g. Georgakakis et al. 2014; Leauthaud et al. 2015; Mendez et al. 2016; Powell et al. 2018). Motivated by the important role of $M_{\star}$ in connecting SMBHs and host galaxies, Yang et al. (2018a) quantitatively derived the $\overline{\rm BHAR}$-$M_{\star}$ relation at different redshifts up to $z=4$. Aided by the stellar-mass history from Behroozi et al. (2013), Yang et al. (2018a) predicted the typical $M_{\rm BH}$-$M_{\star}$ relation in the local universe. At the massive end ($M_{\star}\gtrsim 10^{11.2}$ $M_{\odot}$) of their $M_{\rm BH}$-$M_{\star}$ relation, their $M_{\rm BH}/M_{\star}$ ratio is $\approx 1/500$, similar to observed $M_{\rm BH}/M_{\rm bulge}$ values (e.g. Häring & Rix 2004; Kormendy & Ho 2013). This agreement is expected, as the bulge becomes dominant and $M_{\rm bulge}\approx M_{\star}$ for massive galaxies.

Despite the $\overline{\rm BHAR}$-$M_{\star}$ relation being generally supported by observations, it cannot straightforwardly explain the tightness of the $M_{\rm BH}$-$M_{\rm bulge}$ correlation. The key to the origin of the tight $M_{\rm BH}$-$M_{\rm bulge}$ correlation might be related to the morphology of host galaxies, since $M_{\rm BH}$ is only correlated with the masses of classical bulges rather than other galactic components such as pseudo-bulges or disks (e.g. Kormendy & Ho 2013 and references therein). Therefore, SMBH growth might be related to star formation activity of the bulge only. To investigate this potential SMBH-bulge coevolution, one should ideally study the relation between $\overline{\rm BHAR}$ and bulge SFR in the distant universe. However, with current facilities, it is infeasible to separate the bulge SFR from total SFR when disks are present. In this work, we focus on a sample of bulge-dominated galaxies for which bulge SFR $\approx$ total SFR. If SMBHs indeed coevolve with host-galaxy bulges, we expect a strong correlation between $\overline{\rm BHAR}$ and SFR for these bulge-dominated galaxies over cosmic history.

Our bulge-dominated sample is selected from the Cosmic Assembly Near-Infrared Deep Extragalactic Legacy Survey (CANDELS) where deep HST $H$-band observations are available (Grogin et al. 2011; Koekemoer et al. 2011). X-ray emission is a robust tracer of SMBH accretion (e.g. Brandt & Alexander 2015 and references therein). The CANDELS fields also have deep Chandra X-ray observations, allowing us to estimate reliable $\overline{\rm BHAR}$ for any given sample of galaxies (e.g. Yang et al. 2017, 2018b).

This paper is structured as follows. In §2, we describe the data used in this work and define our samples. In §3, we perform data analyses and present the results. We discuss our results in §4. We summarize our work and discuss future prospects in §5.

Throughout this paper, we assume a cosmology with $H_{0}=70$ km s*-1* Mpc*-1*, $\Omega_{M}=0.3$, and $\Omega_{\Lambda}=0.7$. We adopt a Chabrier initial mass function (IMF; Chabrier 2003). Quoted uncertainties are at the $1\sigma$ (68%) confidence level. We express $M_{\rm bulge}$, $M_{\star}$, and $M_{\rm BH}$ in units of $M_{\odot}$, SFR and $\overline{\rm BHAR}$ in units of $M_{\odot}$ yr*-1*. $L_{\rm X}$ indicates AGN X-ray luminosity at rest-frame 2–10 keV and is in units of erg s*-1*.

2 Data and Sample

Our analyses are based on the five CANDELS fields, i.e., GOODS-S, GOODS-N, EGS, UDS, and COSMOS (Grogin et al., 2011; Koekemoer et al., 2011). All these fields have multiwavelength observations from HST, Spitzer, and ground-based telescopes such as Subaru and VLT. These high-quality data sets allow for measurements of galaxy morphology (§2.1), stellar mass ($M_{\star}$; §2.2), and star-formation rate (SFR; §2.2). FIR observations from Herschel are also available in these fields, enabling robust SFR estimation based on cold-dust emission (§2.2). All of the five CANDELS fields have Chandra X-ray observations from which we derive $\overline{\rm BHAR}$ (§2.4). We define our sample in §2.3 and the sample properties are summarized in Tab. 1.

2.1 Morphology

Rest-frame optical/NIR light is essential for morphological measurements (e.g. Conselice, 2014). The HST $H$ band, centered at $\approx 1.6$ $\mu$m, can cover rest-frame optical/NIR wavelengths up to $z\approx 3$. We adopt the $H$-band morphological measurements in Huertas-Company et al. (2015a) that are based on machine learning for CANDELS galaxies with $H<24.5$. The machine-learning technique is chosen to approximate visual morphologies from humans, and is trained with a galaxy sample that has morphological measurements performed by human classifiers (Kartaltepe et al., 2015). This training sample has the same magnitude cut ($H<24.5$) as reliable visual morphological measurements are difficult at fainter magnitudes. For each galaxy, the Huertas-Company et al. (2015a) catalog provides five fractional numbers, i.e., $f_{\rm sph}$, $f_{\rm disk}$, $f_{\rm irr}$, $f_{\rm pt}$, and $f_{\rm unc}$. These fractions represent the probabilities that a hypothetical classifier would have voted for a galaxy having a spheroid, a disk, and some irregularities, being point-like and unclassifiable, respectively. Note that the sum of the fractions might exceed unity, because, for example, a galaxy might have both spheroidal and disky features simultaneously.

A high $f_{\rm unc}$ value indicates that the source might be a spurious detection, e.g. the spikes of a bright star being falsely detected as a source (see, e.g. Fig. 13 of Huertas-Company et al. 2015a). Sources with high $f_{\rm pt}$ value might be stars or broad-line (BL) AGNs. Due to strong light from the AGN central engine, morphology measurements of host galaxies are unreliable for luminous BL AGNs (e.g. §5.3 of Brandt & Alexander 2015). We exclude the $\approx 8\%$ of sources that have $f_{\rm unc}$ or $f_{\rm pt}$ greater than any of $f_{\rm sph}$, $f_{\rm disk}$, and $f_{\rm irr}$. Upon visual inspection, the excluded sources are indeed spurious detections or point-like. Morphological measurements are challenging at high redshift and our work probes up to $z=3$. We discuss some possible redshift-related effects on our results in §3.4.

2.2 Redshift, Stellar Mass, and Star Formation Rate

We obtain redshift measurements from the CANDELS catalogs (see Tab. 1). These measurements are spectroscopic redshifts (spec-$z$) or photometric redshifts (photo-$z$). The photo-$z$ measurements are based on dedicated photometry extracted with careful consideration of PSF sizes and source shapes (e.g. Guo et al. 2013; Galametz et al. 2013). Compared to the available spec-$z$, the photo-$z$ shows high quality, with $\sigma_{\rm NMAD}=0.018$ and an outlier fraction of 2%.111Here, $\sigma_{\rm NMAD}$ is defined as $1.48\times\mathrm{median}(\frac{|\Delta z-\mathrm{median}(\Delta z)|}{1+z_{\mathrm{spec}}})$, and outliers are defined as sources having ${|\Delta z|/(1+z_{\rm spec})>0.15}$. As in §2.1, we discard the 79 spectroscopic BL AGNs reported in the literature (Barger et al. 2003; Silverman et al. 2010; Cooper et al. 2012; Newman et al. 2013; Marchesi et al. 2016; Suh et al., in prep.).

The CANDELS catalogs also provide stellar mass ($M_{\star}$) and star formation rate (SFR) measurements from independent teams. Following Yang et al. (2017), we adopt the median $M_{\star}$ and SFR values from the five available teams (2aτ, 6aτ, 11$a_{\tau}$, 13aτ, and 14a).222For GOODS-N, only three teams are available (2aτ, 6aτ, and 14a) for now. Fig. 1 (top) shows $M_{\star}$ as a function of redshift for $H<24.5$ galaxies that have morphological measurements (§2.1). We limit our analyses to a $M_{\star}$-complete (corresponding to $H<24.5$) sample (§2.3). The limiting $M_{\star}$ ($M_{\rm lim}$) for $H<24.5$ is also displayed in Fig. 1. The $M_{\rm lim}$-redshift curve is derived based on an empirical method (e.g. Ilbert et al. 2013). We first divide our sources into narrow redshift bins with width of $\Delta z=0.2$. For each redshift bin, we calculate $\log M_{\rm lim}^{\rm ind}=\log M_{\star}+0.4\times(H-24.5)$ for individual galaxies in the bin. We then adopt $M_{\rm lim}$ as the 90th percentile of the $M_{\rm lim}^{\rm ind}$ distribution for the redshift bin.

The CANDELS $M_{\star}$ and SFR are based on SED fitting of rest-frame UV-to-NIR photometry using galaxy templates. As demonstrated by previous works (e.g. Luo et al., 2010; Yang et al., 2017; Kocevski et al., 2018), the rest-frame UV-to-NIR light is often predominantly contributed by galaxy component rather than the AGN component for X-ray AGNs in the CANDELS fields. Also, we have removed BL AGNs that might have strong AGN components in their UV-to-NIR SED. Therefore, the AGN SED contribution should not qualitatively affect our results (see §3.4 for other evidence).

The SED-based SFR estimation, which is physically based on obscuration-corrected UV light, is reliable for low-to-moderate levels ($\mathrm{SFR\lesssim 100}\ M_{\odot}\ \mathrm{yr^{-1}}$) of star-formation activity. However, it tends to underestimate SFR in the high-SFR regime, possibly due to strong dust obscuration (e.g. Wuyts et al., 2011; Whitaker et al., 2017; Yang et al., 2017). To alleviate this issue, we adopt SFR from FIR photometry of Herschel when available (Lutz et al. 2011; Oliver et al. 2012; Magnelli et al. 2013). The photometry has been extracted using Spitzer/MIPS 24 ${\mu}$m priors and source-blending issues have been carefully addressed. The FIR-based SFR is more robust than the SED-based SFR, especially in the high-SFR regime (e.g. Chen et al., 2013; Yang et al., 2017). Due to limited sensitivity, Herschel can only detect sources with the highest SFR at a given redshift. There are five Herschel bands available for the CANDELS fields, i.e., 100 $\mu$m, 160 $\mu$m, 250 $\mu$m, 350 $\mu$m, and 500 $\mu$m. We only utilize robust detections with $\mathrm{S/N>3}$. We discard the 100 $\mu$m band at redshifts above $z=1.5$, because the observed 100 $\mu$m corresponds to rest-frame $<40\ \mu$m which might be contaminated by hot-dust emission powered by AGN activity. We adopt the reddest available Herschel band to estimate SFR, since longer wavelengths are “freer” from possible AGN emission. We calculate SFR from FIR flux following the procedure in Chen et al. (2013) and Yang et al. (2017). We first derive galaxy total IR luminosity ($L_{\rm IR}$) from FIR flux based on the star-forming galaxy templates of Kirkpatrick et al. (2012). We adopt the $z\sim 1$ and $z\sim 2$ templates for $z<1.5$ and $z\geq 1.5$ sources, respectively. We then obtain SFR as

[TABLE]

Fig. 1 (middle) shows SFR (based on SED or FIR) as a function of redshift for all $H<24.5$ galaxies.

The comparison sample has a higher fraction of FIR-based SFR measurements than the bulge-dominated sample (32% vs. 9%), because the former generally has stronger star-formation activity than the latter. To investigate whether this difference in SFR measurements could bias our results, we have tested cutting our $z=0.5\text{--}1.5$ ($z=1.5\text{--}3$) sample at $\mathrm{SFR}<10$ $M_{\odot}$ yr*-1* ($\mathrm{SFR}<10^{1.5}$ $M_{\odot}$ yr*-1*), below which the SFR measurements are mostly SED-based (see Fig. 1). Under these cuts, our results do not change qualitatively.333We cannot perform a similar test for high-SFR galaxies with Herschel detections, because the sample size of Herschel-detected sources is too small. We have also tested our results using SED-based SFR only for the entire galaxy sample and our conclusions still hold. It is well known that SFR measurements from SEDs and the FIR do not always agree (e.g. Buat et al., 2010; Rodighiero et al., 2014), and Yang et al. (2017) found that the statistical scatter between these two methods is ${\lesssim 0.5}$ dex. This level of scatter is unlikely to be seriously problematic to our statistical analyses (§3.2), since our SFR bin sizes are typically ${\gtrsim 0.5}$ dex. To verify this point, we perturb our SFR measurements by 0–0.5 dex randomly and our results below in §3 do not change qualitatively after the perturbation. Yang et al. (2017) found the systematic offset between SED-based and FIR-based SFRs is typically small for the general galaxy population (${\approx 0.2}$ dex), and this level of systematic should not change our main conclusions considering our relatively large SFR bin sizes. However, Symeonidis et al. (2016) considered that, for strong AGNs, FIR-based SFR might be systematically overestimated due to the contamination of AGN-heated dust. To assess this potential issue, we test our results by using SED-based SFR only for strong AGNs (${\log L_{\rm X}>43.5}$), and our results below do not change qualitatively. Therefore, our main conclusions should be robust against uncertainties of SFR measurements.

2.3 The Bulge-Dominated and Comparison Samples

Our analyses are based on a bulge-dominated sample and a comparison sample. In this Section, we detail the selections of these two samples. We first select all $H<24.5$ galaxies for which morphology measurements are available (§2.1). As in §2.1 and §2.2, we exclude BL AGNs, stars, and false detections. We then divide these galaxies into two redshift bins, i.e., $z=0.5\text{--}1.5$ and $z=1.5\text{--}3.0$ for our analyses. The relatively broad redshift bins are necessary to guarantee sufficiently large samples for our statistical analyses (§3). We have also tested on narrower redshift bins and found our qualitative results do not change, although the statistical scatter becomes larger due to reduced sample sizes. Therefore, the two wide redshift bins ($z=0.5\text{--}1.5$ and $z=1.5\text{--}3.0$) should not bias our results, and we adopt them throughout this paper.

We select $M_{\star}$-complete samples for the two redshift bins. The limiting $M_{\star}$ at $z=1.5$ and $z=3.0$ are $\log M_{\star}\approx 9.7$ and $\log M_{\star}\approx 10.2$, respectively (see Fig. 1). Therefore, we limit our analyses to $\log M_{\star}>9.7$ and $\log M_{\star}>10.2$ galaxies for the low and high redshift bins, respectively. These $M_{\star}$ thresholds are below the characteristic $M_{\star}$ of the stellar-mass function (SMF), i.e., $\log M_{\star}\approx 10.6$ at $z\approx 0.5\text{--}3$ (e.g. Tomczak et al., 2014; Davidzon et al., 2017). The stellar-mass density above these $M_{\star}$ cuts is $\approx 90\%$ and $\approx 70\%$ of the total for the low and high redshift bins, respectively (calculated with the SMF in Behroozi et al. 2013). After applying the $M_{\star}$ cuts, our sample does not include dwarf galaxies ($\log M_{\star}\lesssim 9.5$). Aside from technical constraints, the exclusion of dwarf galaxies is also motivated by our major science goal, i.e. investigating the origin of the $M_{\rm BH}$-$M_{\rm bulge}$ relation. Since the $M_{\rm BH}$-$M_{\rm bulge}$ relation is mostly established for $\log M_{\rm bulge}\gtrsim 10$ (e.g. Kormendy & Ho, 2013), we should also focus on relatively massive galaxies rather than dwarf galaxies.

The basic properties of the $M_{\star}$-complete sample are summarized in Tab. 1. In the $M_{\star}$-complete sample, we classify a source as bulge-dominated if it satisfies $f_{\rm sph}\geq 2/3$, $f_{\rm disk}<2/3$, and $f_{\rm irr}<1/10$. These empirical criteria are suggested by several previous studies (e.g. Huertas-Company et al. 2015b, 2016; Kartaltepe et al. 2015). As indicated by these criteria, the term “bulge-dominated” refers to galaxies that only display a significant spheroidal component, without obvious disky and/or irregular components. We note that different authors may adopt different terminology for the bulge-dominated galaxies (e.g. “spheroid-like”; e.g. Conselice 2014). If a galaxy does not meet these criteria, we include it in our comparison sample, i.e., galaxies that are not bulge-dominated. The bulge-dominated and comparison samples have $\approx$ 2,600 and 7,700 galaxies, respectively (see Tab. 1). Our analyses in §3 are based on these two samples.

In Fig. 2, we show some random $H$-band cutouts for the bulge-dominated and comparison samples, respectively. The fractions of bulge-dominated galaxies are both $\approx 25\%$ for the low and high redshift bins. Fig. 3 shows the fraction of bulge-dominated galaxies as a function of $M_{\star}$ and SFR, respectively. At the high-$M_{\star}$ end ($\log M_{\star}\gtrsim 11$), the bulge-dominated fraction in the low-redshift bin is much higher than that in the high-redshift bin ($\approx 50\%$ vs. $\approx 20\%$). This is probably due to the fact that galaxy mergers/interactions for massive galaxies are increasingly prevalent toward high redshift, and thus galaxy irregularities are much stronger toward the early universe (e.g. Conselice, 2014; Marsan et al., 2018). The bulge-dominated fraction drops significantly toward high SFR, indicating that bulge-dominated galaxies tend to have low SFR. Similar trends have also been found in previous studies (e.g. Huertas-Company et al. 2015b, 2016). The underlying physical reason might be “morphological quenching”, such that bulges can effectively suppress star formation (e.g. Martig et al., 2009).

Fig. 4 displays the source distributions on the SFR-$M_{\star}$ plane for the bulge-dominated and comparison samples, respectively. The bulge-dominated sample tends to lie below the star-formation main sequence, while the majority of the comparison sample appears to be on the main sequence. However, we note that our morphological classification is essentially different from a star-forming vs. quiescent classification. For example, the quiescent population in our sample is made up of $\approx 55\%$ bulge-dominated galaxies and $\approx 45\%$ comparison galaxies. While the main population of the comparison sample lies on the main sequence, there is a non-negligible fraction (${\approx 20\%}$) of comparison galaxies lying significantly (${\approx 1}$ dex) below the main sequence. We have visually checked the HST cutouts of these low-SFR sources, and found they appear to have significant disk/irregular components. Therefore, the existence of such a low-SFR population among our comparison sample is likely intrinsic, and does not appear to be caused by morphological misclassification.

2.4 Black Hole Accretion Rate

All five CANDELS fields have deep X-ray observations from Chandra. Tab. 1 lists the X-ray depth and number of X-ray detected sources for each CANDELS field. We calculate $\overline{\rm BHAR}$ contributed by both X-ray detected and undetected sources, and thus the resulting $\overline{\rm BHAR}$ should cover essentially all SMBH accretion. This procedure allows us to seamlessly analyze all sources in different CANDELS fields which have different X-ray depths. We have also repeated our analyses but without sources in COSMOS, which has X-ray depth much shallower than other fields (see Tab. 1), and our results do not change qualitatively. For each X-ray detected source, we calculate $L_{\rm X}$ from the X-ray flux from the corresponding X-ray catalog assuming a photon index of $\Gamma=1.7$ (e.g. Yang et al. 2016; Liu et al. 2017). Following Yang et al. (2018b), we choose, in order of priority, hard-band (observed-frame $2\text{--}7$ keV), full-band (observed-frame $0.5\text{--}7$ keV), or soft-band (observed-frame $0.5\text{--}2$ keV) flux to minimize X-ray obscuration effects. Indeed, Yang et al. (2018b) estimated that, under this scheme of band choice, the X-ray flux decrease due to obscuration is typically small ($\approx 20\%$) for bright sources in CDF-S, for which there are enough photons to assess obscuration. We increase the X-ray fluxes of our X-ray sources by 20% to account for the average systematic effect from obscuration.

For X-ray undetected sources, we employ a stacking technique to include their X-ray emission. We perform this process on full-band X-ray images.444For the EGS field, we use the X-ray image from Goulding et al. (2012), since Nandra et al. (2015) did not produce the X-ray image for the entire EGS field. We generally follow the steps in Vito et al. (2016), and we briefly summarize this procedure below. First, we mask X-ray detected sources in the X-ray images. We choose masking radii ($R_{\mathrm{msk}}$) of $2\times R_{\mathrm{90}}$, $2.25\times R_{\mathrm{90}}$, and $2.5\times R_{\mathrm{90}}$ for sources with net counts $<100$, $100\text{--}1000$, and $>1000$, respectively. Here, $R_{90}$ means the radius for a 90% encircled-energy fraction (EEF), and it is a function of off-axis angle (see Appendix A of Vito et al. 2016). After source masking, we derive the net count rate for X-ray undetected sources. To enhance signal-to-noise, we adopt an aperture radius $R_{\mathrm{aper}}=R_{\mathrm{80}}$, $R_{\mathrm{75}}$, $R_{\mathrm{60}}$, and $R_{\mathrm{40}}$ for sources with off-axis angle $<3.5\arcmin$, $3.5\arcmin\text{--}4.25\arcmin$, $4.25\arcmin\text{--}5.5\arcmin$, and $5.5\arcmin\text{--}7.8\arcmin$, respectively. We discard sources whose off-axis angle is $>7.8\arcmin$ and/or whose apertures overlap with masked regions. The counts in the apertures include contributions from both sources and background, and we need to subtract the background counts. We estimate the background counts in an annulus with inner and outer radii of $1.1R_{90}$ and $1.1R_{90}+10\arcsec$, respectively. Due to the limited aperture size, the net counts encircled in the aperture only represent a fraction of the total net counts. We perform an aperture correction for each source, depending on the aperture size adopted. For example, if $R_{\mathrm{aper}}=R_{\mathrm{80}}$ for a source, we divide the aperture net counts by 80% to recover the total net counts. We then obtain the count rate by dividing the total net counts by the exposure time at the position of the source. For sources in each field, we derive fluxes by multiplying the count rates by a constant factor, which is the median flux/count-rate ratio of X-ray detected sources in the field. Finally, for a group of X-ray undetected sources, we can obtain their average X-ray luminosity ($\overline{L_{\mathrm{X,stack}}}$) from the average X-ray flux and redshift, assuming $\Gamma=1.7$. Our derived SMBH accretion power is mostly contributed by the X-ray detected sources, and the stacking procedure typically accounts for less than 20% of the accretion power.

For X-ray detected objects, we have $L_{\rm X}$ for individual sources; for X-ray undetected sources, we have $\overline{L_{\mathrm{X,stack}}}$ for any group of sources. We can then calculate average AGN bolometric luminosity for any sample of sources as

[TABLE]

Here, $N_{\rm det}$ and $N_{\rm non}$ are the numbers of X-ray detected and undetected sources in the sample. $L_{\rm X,XRB}$ is the expected luminosity from X-ray binaries (XRBs) and $\overline{L_{\rm X,XRB}}$ is the average XRB luminosity for the stacked sources. To obtain $L_{\rm X,XRB}$ and $\overline{L_{\rm X,XRB}}$, we adopt model 269 of Fragos et al. (2013) which describes XRB X-ray luminosity as a linear function of $M_{\star}$ and SFR. Model 269 is a theoretical model favored by the observations of galaxies at $z=0\text{--}2$ (Lehmer et al., 2016). The expected X-ray emission from XRBs only accounts for $\approx 15\%$ of the total X-ray power, and thus the uncertainties related to the XRB modelling should not affect our analyses significantly. $k_{\rm bol}$ and $\overline{k_{\rm bol}}$ are the $L_{\rm X}$-dependent bolometric corrections at $(L_{\rm X}-L_{\rm X,XRB})$ and $(\overline{L_{\rm X,stack}}-\overline{L_{\rm X,XRB}})$, respectively. We adopt the bolometric-correction model from Hopkins et al. (2007).555As pointed out in Footnote 4 of Merloni & Heinz (2013), the $k_{\rm bol}$ in Hopkins et al. (2007) appears to be overestimated due to the double counting of IR reprocessed emission. Following Merloni & Heinz (2013), we multiply the $k_{\rm bol}$ in Hopkins et al. (2007) by a factor of 0.7 to address this issue. Assuming a constant radiative efficiency of $\epsilon=0.1$, we can convert $\overline{L_{\rm bol}}$ to $\overline{\rm BHAR}$ as

[TABLE]

where $c$ is the speed of light. The adopted $\epsilon=0.1$ is motivated by observations (see, e.g. §3.4 of Brandt & Alexander 2015). We obtain the $\overline{\rm BHAR}$ uncertainties with a bootstrapping technique (e.g. §2.3 of Yang et al. 2017).

As explained in §1, the $\overline{\rm BHAR}$ quantity is designed to approximate long-term average SMBH accretion rate, and has been widely adopted in the studies of AGN-galaxy relations (e.g. Chen et al. 2013; Hickox et al. 2014; Yang et al. 2017, 2018a, 2018b; Fornasini et al. 2018). Some works proposed to recover the full distribution of BHAR as a function of galaxy properties (e.g. Volonteri et al. 2015; Georgakakis et al. 2017; Aird et al. 2018b, a), and quantities such as $\overline{\rm BHAR}$ and duty cycle can then be derived. However, detailed modelling of the BHAR distribution at given $M_{\star}$, SFR, and morphological type is beyond the scope of this work, and we leave it to future studies.

3 Analyses and Results

In this Section, we study $\overline{\rm BHAR}$ as a function of SFR and $M_{\star}$ (§3.1). We address the question of whether $\overline{\rm BHAR}$ is mainly related to SFR or $M_{\star}$ in §3.2. All these analyses are performed for the bulge-dominated and comparison samples, respectively. In Appendix A, we perform the same analyses for all galaxies. In §3.3, we quantify the $\overline{\rm BHAR}$-SFR relation for the bulge-dominated sample.

3.1 BHAR as a Function of SFR and $M_{\star}$

We plot the $\overline{\rm BHAR}$ as a function of SFR for our bulge-dominated and comparison samples, respectively, in Fig. 5 (black points). In each panel, the bins are chosen to include approximately the same number of sources, and this approach is to reach similar $\overline{\rm BHAR}$ signal-to-noise (S/N) ratios for the bins. Adjusting the bins does not change our conclusions qualitatively, although the statistical scatter of $\overline{\rm BHAR}$ measurements increases.

For the bulge-dominated sample, $\overline{\rm BHAR}$ rises strongly from low to high SFR by a factor of $\approx 400$ ($z=0.5\text{--}1.5$) and $\approx 100$ ($z=1.5\text{--}3.0$). In contrast, for the comparison sample, $\overline{\rm BHAR}$ only increases by a factor of $\approx 10$ ($z=0.5\text{--}1.5$) and $\approx 2$ ($z=1.5\text{--}3.0$) from low to high SFR. We show $\overline{\rm BHAR}$ vs. $M_{\star}$ in Fig. 6 (black points). For the bulge-dominated sample, there is no strong correlation between $\overline{\rm BHAR}$ and $M_{\star}$. For the comparison sample, $\overline{\rm BHAR}$ appears to rise toward high $M_{\star}$ in general. We note that, due to our limited sample size, statistical fluctuations can be strong sometimes. For example, for the black point at $\log M_{\star}\approx 9.8$ in Fig. 6 (right), the $\overline{\rm BHAR}$ is mostly contributed by a single source. These fluctuations inevitably cause some scatter in Figs. 5 and 6.

3.2 Is $\overline{\rm BHAR}$ Mainly Related to SFR or $M_{\star}$?

In this Section, we address the question of whether $\overline{\rm BHAR}$ is mainly related to SFR or $M_{\star}$ for the bulge-dominated and comparison samples, respectively. The analysis methods here are similar to those in Yang et al. (2017). We compare our results with Yang et al. (2017) in Appendix A.

In Fig. 5, we divide each SFR bin into two bins with $M_{\star}$ above and below the median $M_{\star}$ of the SFR bin, respectively. In general, the high-$M_{\star}$ and low-$M_{\star}$ bins have similar $\overline{\rm BHAR}$ for the bulge-dominated sample. However, the high-$M_{\star}$ bins have significantly higher $\overline{\rm BHAR}$ than the corresponding low-$M_{\star}$ bins for the comparison sample. Similarly, in Fig. 6, we also divide each $M_{\star}$ bin into high-SFR and low-SFR bins. The high-SFR bins have much higher $\overline{\rm BHAR}$ than the corresponding low-SFR bins for the bulge-dominated sample. In contrast, the high-SFR and low-SFR bins have similar $\overline{\rm BHAR}$ for the comparison sample.

The results above qualitatively indicate that $\overline{\rm BHAR}$ might primarily depend on SFR rather than $M_{\star}$ for the bulge-dominated sample and that the situation is the opposite for the comparison sample. To further test this point, we perform partial-correlation (PCOR) analyses with pcor.r in the r statistical package (Kim, 2015). We first bin sources based on both SFR and $M_{\star}$ and calculate $\overline{\rm BHAR}$ for each bin, and Fig. 7 shows the results. Following Fig. 5, the bins for the $x$-axis ($y$-axis) include similar numbers of sources. Adjusting the bins does not affect our results qualitatively. We input the median $\log M_{\star}$, median $\log\mathrm{SFR}$, and $\log\overline{\rm BHAR}$ in each bin to pcor.r and calculate the significance levels for the $\overline{\rm BHAR}$-$M_{\star}$ and $\overline{\rm BHAR}$-SFR relations, respectively. The PCOR tests are performed with the Pearson and Spearman statistics, respectively, and the results are summarized in Tab. 2. These results show that, for the bulge-dominated sample, the $\overline{\rm BHAR}$-SFR correlation is significant ($>3\sigma$) while the $\overline{\rm BHAR}$-$M_{\star}$ correlation is not ($<3\sigma$). For the comparison sample, the $\overline{\rm BHAR}$-$M_{\star}$ correlation is significant while the $\overline{\rm BHAR}$-SFR correlation is not. These conclusions are also supported by Figs. 5 and 6. We note that the lack of a significant ${\overline{\rm BHAR}}$-SFR relation for the comparison sample is unlikely to be caused by X-ray obscuration effects, because the effects of obscuration on our $\overline{\rm BHAR}$ measurements are generally small (see §2.4).

Fig. 5 (left) is the key plot in this paper. It displays the strong $\overline{\rm BHAR}$-SFR connection and qualitatively demonstrates that the $\overline{\rm BHAR}$-SFR relation cannot be significantly split by $M_{\star}$. We have also tested dividing each SFR bin by other galaxy properties (instead of $M_{\star}$) such as $f_{\rm disk}$ (§2.1) and rest-frame $U-V$ color, and none of these parameters can significantly split the $\overline{\rm BHAR}$-SFR relation. Therefore, the strong $\overline{\rm BHAR}$-SFR correlation is likely fundamental.

3.3 Quantification of the $\overline{\rm BHAR}$-SFR Relation

In §3.2, we find that $\overline{\rm BHAR}$ is mainly correlated with SFR rather than $M_{\star}$ for the bulge-dominated sample. To quantify this $\overline{\rm BHAR}$-SFR relation, we fit the data points in Fig. 5 (left; black points) with a log-linear model. For convenience, we list the sample properties of each data point in Tab. 3. We adopt a standard least-$\chi^{2}$ fitting method implemented by a python package scipy.optimize.curve_fit. We first fit the data points in the two redshift bins independently, and the results are

[TABLE]

Considering the best-fit parameters are similar for the two redshift bins, we fit all the data points in both redshift bins simultaneously. The best-fit model is

[TABLE]

where the errors are calculated under a 68% confidence level. The reduced $\chi^{2}$ of the fit is 0.8 ($p\text{-value}=53\%$), showing that the fit quality is acceptable. Considering that the slope of the best-fit model is close to unity, we also fit the data with slope fixed to unity. This procedure results in

[TABLE]

The fit quality is also acceptable, with reduced $\chi^{2}$ of 1.2 ($p\text{-value}=32\%$). This best-fit model is displayed in Fig. 5 (left). The best-fit $\overline{\rm BHAR}$/SFR ratio in this model is $10^{-2.48}$.

Our $\overline{\rm BHAR}$ does not include the accretion from BL AGNs (§2). Here, we consider this missed accretion power statistically. We first construct a non-BL AGN sample with $L_{\rm X}$ and redshift matched with the spectroscopic BL AGN sample (§2.2): for each BL AGN, we randomly select a “nearby” non-BL AGN in the $L_{\rm X}$-$z$ plane (within $\log L_{\rm X}\pm 0.15$ and $z\pm 0.2$). We find that $\approx 35\%$ of these non-BL AGNs reside in bulge-dominated galaxies. Assuming, following the unified model, that a similar fraction of BL AGNs have bulge-dominated hosts, we find that the missed accretion power (contributed by BL AGNs) is $\approx 40\%$ of the observed accretion power for bulge-dominated galaxies. Therefore, after including the accretion power of BL AGNs, the $\overline{\rm BHAR}$/SFR ratio might be slightly higher ($\approx 0.15$ dex) than the best-fit value. We note that errors resulting from radiative-efficiency and IMF uncertainties likely exist, and thus the best-fit value of $\overline{\rm BHAR}$/SFR inevitably suffers from systematic uncertainty up to a factor of a few. However, the systematic uncertainties should not affect our main qualitative conclusion, i.e. $\overline{\rm BHAR}$ primarily depends on SFR among bulge-dominated galaxies.

Considering the importance of SMBH-galaxy growth among bulge-dominated galaxies, we also plot AGN fraction as a function of SFR in Fig. 8. Here, we count an X-ray source as an “AGN” if it has $\log L_{\rm X}>42.8$ ($z=0.5\text{--}1.5$) or $\log L_{\rm X}>43.5$ ($z=1.5\text{--}3$). These thresholds are the $L_{\rm X}$ limits at $z=1.5$ and $z=3$, respectively, for a $0.5\text{--}10$ keV flux limit of $8.9\times 10^{-16}$ erg cm*-2* s*-1*. This flux limit is the detection limit of the COSMOS survey (Civano et al., 2016), which is the shallowest CANDELS X-ray survey (Tab. 1). $\overline{\rm BHAR}$ is mainly driven by duty cycle and average accretion rate of AGNs. From Fig. 8,666We cannot derive reliable average AGN X-ray luminosities and thereby accretion rates due to the small AGN sample sizes in most bins (see Tab. 3). AGN fraction increases toward high SFR for both redshift bins. This result indicates that the positive $\overline{\rm BHAR}$-SFR relation is, at least partially, due to the rise of AGN duty cycle toward high SFR. Detailed quantitative analyses of AGN duty cycle and average accretion rate require the full distribution of BHAR (see §2.4), for which we leave to future studies.

3.4 Reliability Checks

Our $M_{\star}$ and SFR measurements are mostly based on SED fitting of rest-frame UV-to-NIR photometry, and we have removed BL AGNs from our sample to avoid strong AGN SED components that might affect our results (§2.2 and §3.3). Also, from Fig. 2, the X-ray detected sources do not appear to have strong central point-like emission, indicating that the presence of AGNs should not significantly affect the $M_{\star}$, SFR, and morphology measurements. On the other hand, if potential AGN SED contamination significantly biased our analyses, we would expect to reach a similar conclusion for bulge-dominated and comparison samples, which is not the case (§3.2). Therefore, we consider that our conclusions are not biased by AGN SED contamination. In §2.2, we have also discussed the SFR uncertainties of SED-based and FIR-based measurements, and found that our main conclusions are unlikely to be affected by those uncertainties.

Our bulge-dominated galaxies at $z=0.5\text{--}3$ are selected utilizing machine-learning morphological measurements based on HST $H$-band imaging (§2.1 and §2.3). Morphological measurements at high redshift are challenging due to effects such as redshifting of photons and image degradation (e.g. Conselice 2014 and references therein). Detailed assessment of these redshift effects on our results requires careful simulations of HST imaging and repetition of the machine-learning measurements on the simulated data. These procedures are beyond the scope of this work. Here, we qualitatively discuss the robustness of our results against such redshift effects.

Due to the redshifting of photons, the same observed-frame wavelength covers different rest-frame wavelengths at different redshifts. The correction for this redshifting effect is named the “morphological $k$-correction”. From multiwavelength observations of local galaxies, Taylor-Mager et al. (2007) found that the morphological $k$-correction is weak in the optical/NIR wavelength range ($\approx 0.36\text{--}0.85$ $\mu$m, where $0.36\ \mu$m corresponds to the Balmer break and $0.85\ \mu$m is the longest wavelength available in their work), especially for elliptical/S0 galaxies. For our work, the observed-frame $H$ band (used for morphological measurements; §2.1) does not reach out to rest-frame UV photons below the Balmer break at $z=0.5\text{--}3$, and it corresponds to the rest-frame NIR ($\approx 1\ \mu$m) and optical ($\approx 0.4\ \mu$m) light at $z=0.5$ and $z=3$, respectively. Also, considering that we only utilize the morphological information for a basic selection of bulge-dominated galaxies rather than, e.g. a quantitative measurement of galaxy size, we conclude that the morphological $k$-correction should not affect our results qualitatively. Another point of support for this conclusion is that, although the $H$ band is sampling different wavelengths for $z=0.5\text{--}1.5$ and $z=1.5\text{--}3$, we have obtained qualitatively the same results for the two redshift bins (§3.2).

At low redshift, the $H$ band samples rest-frame red optical/NIR photons. Since galactic-disk components are generally bluer than bulge components, one might worry that $H$-band imaging could miss disk components. This issue mainly happens for low-$M_{\star}$ faint disky galaxies. Considering that our main focus is relatively massive galaxies (§2.3) and that CANDELS $H$-band data are deep, this issue might not be problematic for our study. However, we still check this issue in the following ways. First, we visually check HST $I$-band cutouts of $\approx 100$ random galaxies in our bulge-dominated samples at $z=0.5\text{--}1$. We do not find significant disk components for these sources. Indeed, the average $I$-band Sérsic index of our low-redshift bulge-dominated sample is $\approx 3.5$ (using the measurements of Scarlata et al. 2007), which is typical for bulge-dominated galaxies (e.g. Buitrago et al., 2013; Conselice, 2014). On the other hand, we visually check the low-redshift disky ($f_{\rm disk}\geq 2/3$; §2.3) galaxies in our sample. We find their disk components do not appear to be significantly weaker in $H$ band than in $I$ band, and we attribute this result to the deep exposure and relatively low extinction of $H$ band. Therefore, we consider that the $H$ band-based morphological classification is robust at low redshift.

The imaging quality generally becomes worse toward higher redshift. Some fine galactic structures (e.g. spiral arms and clumps) might be smoothed out, and thus some disky and irregular galaxies might be classified as the smooth bulge-dominated type (e.g. Mortlock et al. 2013). Therefore, our bulge-dominated sample might be “contaminated”. However, this issue is, at least to some extent, mitigated by the $H$ mag cut ($H<24.5$) applied to our sample (§2.1). This cut guarantees a minimum signal-to-noise ratio (S/N $\approx 80$) of the imaging, with the penalty of a smaller sample size. Also, if our bulge-dominated sample were strongly contaminated, we would observe a similar $\overline{\rm BHAR}$-SFR-$M_{\star}$ relation for both the bulge-dominated and comparison samples. However, the $\overline{\rm BHAR}$ dependences on SFR and $M_{\star}$ are qualitatively different for these two samples (§3.2). We thus consider that image degradation should not be a significant issue for our conclusions.

4 Discussion

4.1 Physical Implications

We emphasize that the $\overline{\rm BHAR}$-SFR correlation only exists for our bulge-dominated sample, while $\overline{\rm BHAR}$ appears to be primarily correlated with $M_{\star}$ for the comparison sample. This difference indicates that SMBHs only coevolve with bulges rather than entire galaxies, consistent with the observations of local systems (e.g. Kormendy & Ho, 2013; Davis et al., 2018). Such SMBH-bulge coevolution might be driven by the amount of cold gas available in the bulge, since both SMBH and bulge growth require cold gas. From the SMBH-bulge coevolution scenario, we expect that $\overline{\rm BHAR}$ is also fundamentally correlated with bulge SFR even when a galactic disk is present.

Earlier studies speculated an intrinsic $\overline{\rm BHAR}$-SFR relation for the overall galaxy population (see §1). However, the scenario of SMBH vs. entire galaxy coevolution leads to $M_{\rm BH}$ being strongly related to $M_{\star}$ rather than $M_{\rm bulge}$ in the local universe, contradicting observations (e.g. Kormendy & Ho, 2013). To reconcile this contradiction, an ad hoc galaxy evolution model was invoked where all stellar mass formed in the distant universe ($z\gtrsim 0.5$) is transformed to bulge mass at $z=0$ (Jahnke et al., 2009; Mullaney et al., 2012). In contrast, the $\overline{\rm BHAR}$-bulge SFR correlation, as revealed by our work, can naturally result in the $M_{\rm BH}$-$M_{\rm bulge}$ relation observed in the local universe, without invoking any unphysical galaxy evolution models. Our findings highlight the critical role of morphological measurements when observationally studying the connections between distant SMBHs and their host galaxies, as the $\overline{\rm BHAR}$-SFR correlation only exists among bulge-dominated galaxies. Without deep HST observations of CANDELS, our discovery would not be possible (see Appendix A).

Some papers attribute the local $M_{\rm BH}$-$M_{\rm bulge}$ relation entirely to a non-causal statistical origin (e.g. Peng 2007; Jahnke & Macciò 2011). If galaxy/SMBH mergers happen frequently enough, the scatter of the $M_{\rm BH}$-$M_{\rm bulge}$ relation could be averaged out. Our results show that there is indeed an intrinsic $\overline{\rm BHAR}$-SFR connection at high redshift that can lead to the $M_{\rm BH}$-$M_{\rm bulge}$ correlation among nearby galaxies (§4.2). Therefore, the non-causal scenarios of merger averaging are not necessary to explain the $M_{\rm BH}$-$M_{\rm bulge}$ relation. Also, recent observations of Yang et al. (2018a) show that frequent mergers will lead to a $M_{\rm BH}/M_{\rm bulge}$ ratio much smaller than the observed values in the local universe.

Kocevski et al. (2017) found that, for compact galaxies, the star-forming population has elevated AGN fraction compared to the quiescent population with matched ${M_{\star}}$ at ${z\approx 2}$. However, for extended galaxies, the star-forming and quiescent populations have similar AGN fractions. Since our bulge-dominated population is morphologically more compact than other galaxy populations in general (see Fig. 2; e.g. Huertas-Company et al. 2015a), our results in Fig. 6 are broadly consistent with the findings of Kocevski et al. (2017). While we consider our results as evidence of SMBH-bulge coevolution, Kocevski et al. (2017) argued that a contraction process might trigger both compact starburst activity and SMBH accretion. In our scenario, bulge SFR is fundamentally correlated with ${\overline{\rm BHAR}}$; in their scenario, compactness is a critical galaxy property linked with SMBH growth. To address the question of which scenario is more physical, one needs to break the degeneracy that bulge-dominated systems are generally compact. We will perform these analyses in a future paper (Ni et al. in prep.).

4.2 Implications for the $M_{\rm BH}$-$M_{\rm bulge}$ Relation

From the best-fit results in §3.3, we have $\overline{\rm BHAR}/\mathrm{SFR}=10^{-2.48}$. This value is similar to the typical observed $M_{\rm BH}/M_{\rm bulge}$ values in the local universe ($\approx 10^{-2.5}\text{--}10^{-2.2}$; Kormendy & Ho 2013). Also, similar to the observed $M_{\rm BH}$-$M_{\rm bulge}$ relation in the local universe, our $\overline{\rm BHAR}$-SFR relation for bulge-dominated galaxies has slope close to unity. These similarities indicate that the observed $M_{\rm BH}$-$M_{\rm bulge}$ relation originates from SMBH-bulge coevolution as revealed by our work, and the $M_{\rm BH}$-$M_{\rm bulge}$ relation is not heavily biased by the possibility that observations tend to select massive SMBHs for $M_{\rm BH}$ measurements (e.g. Shankar et al., 2016).

The strong $\overline{\rm BHAR}$-SFR relation among bulge-dominated galaxies indicates that SMBH and bulge growth are in lockstep. A natural consequence from this lockstep growth is that the $M_{\rm BH}$-$M_{\rm bulge}$ relation should not have strong redshift dependence. Some observations suggest that the $M_{\rm BH}/M_{\rm bulge}$ ratio appears to be higher toward higher redshifts (e.g. Shields et al., 2006; Ho, 2007), contradicting the scenario of lockstep growth. However, this apparent redshift dependence of $M_{\rm BH}/M_{\rm bulge}$ might result from observational biases (e.g. Lauer et al., 2007), because $M_{\rm BH}$ measurements in the distant universe are generally limited to luminous quasars. These luminous quasars are likely the most massive SMBHs accreting at high Eddington ratios, and thus the observed $M_{\rm BH}/M_{\rm bulge}$ should be systematically higher than the typical $M_{\rm BH}/M_{\rm bulge}$ among the entire galaxy population.

4.3 Galaxies that are Not Bulge-Dominated

For our bulge-dominated sample, ${\overline{\rm BHAR}}$ is fundamentally related to SFR. In contrast, for our comparison sample consisting of galaxies that are not bulge-dominated (§2.3), ${\overline{\rm BHAR}}$ is not strongly coupled with SFR (§3.2), likely due to the fact that their total SFR is mostly contributed by non-bulge components. Actually, most (${\approx 80\%}$) of the comparison galaxies are irregular/disk-dominated galaxies with no significant bulge components (${f_{\rm sph}<2/3}$; §2.3). The rest (${\approx 20\%}$) of the population in the comparison sample is bulge-disk systems. For these systems, according to the SMBH-bulge coevolution scenario, ${\overline{\rm BHAR}}$ should be intrinsically correlated with bulge SFR (§4.1).

For our comparison sample, ${\overline{\rm BHAR}}$ is strongly related to ${M_{\star}}$ (§3.2). The implications of this ${\overline{\rm BHAR}}$-${M_{\star}}$ relation are discussed in detail by Yang et al. (2018a),777Although Yang et al. (2018a) focused on the ${\overline{\rm BHAR}}$-${M_{\star}}$ relation for the entire galaxy population, their conclusions should largely hold for our comparison galaxies which are numerically the main population (${\approx 75\%}$; §2.3 and Appendix A). and we only summarized their main points below. Yang et al. (2018a) found that ${\overline{\rm BHAR}}$/${\overline{\rm SFR}}$ rises toward high ${M_{\star}}$, i.e. massive galaxies are more effective in feeding their SMBHs (see their §4.2). This result inevitably leads to a higher ${M_{\rm BH}/M_{\star}}$ ratio for more massive galaxies in the local universe, i.e. the typical local ${M_{\rm BH}}$-${M_{\star}}$ relation should be non-linear (see their §4.3 and §4.4). On the other hand, Yang et al. (2018a) also considered that the local ${M_{\rm BH}}$-${M_{\star}}$ relation might not be tight due to different stellar-mass histories of local galaxies with similar ${M_{\star}}$ (see their §3.4.1). This is because ${\overline{\rm BHAR}}$ is higher toward high redshift, at a given ${M_{\star}}$ (see their Fig. 9). Therefore, for two galaxies with similar ${M_{\star}}$ in the local universe, the one that forms at higher redshift should have a more massive SMBH.

5 Summary and Future Work

We have studied the $\overline{\rm BHAR}$ dependence on SFR and $M_{\star}$ for a bulge-dominated sample and a comparison sample of galaxies, respectively, based on multiwavelength observations of the CANDELS fields. Our main analysis procedures and conclusions are summarized below:

We have compiled redshift, $M_{\star}$, and SFR for galaxies brighter than $H=24.5$ from the CANDELS catalogs (§2.2). The CANDELS $M_{\star}$ and SFR measurements are based on SED fitting. For sources detected by Herschel, we estimate their SFR from FIR photometry. We have applied $M_{\star}$ cuts of $\log M_{\star}>9.7$ ($z=0.5\text{--}1.5$) and $\log M_{\star}>10.2$ ($z=0.5\text{--}1.5$) to our sample to ensure $M_{\star}$ completeness (§2.3). Based on machine-learning morphological measurements (§2.1), we have selected a sample of bulge-dominated galaxies and included the other galaxies in a comparison sample (§2.3). The bulge-dominated galaxies consist of $\approx 25\%$ of the entire galaxy population. 2. 2.

We have measured sample-averaged BHAR for different samples of galaxies based on the deep X-ray observations from Chandra (§2.4). We first measure the $L_{\rm X}$ for each X-ray detected source as well as average X-ray luminosity for undetected sources via a stacking process. From these measurements, we calculate average AGN bolometric luminosity adopting an $L_{\rm X}$-dependent bolometric correction. Finally, we estimate $\overline{\rm BHAR}$ from $\overline{L_{\rm bol}}$ assuming a constant radiation efficiency. 3. 3.

For the bulge-dominated sample, we have shown, with both qualitative and quantitative (PCOR) analyses, that $\overline{\rm BHAR}$ primarily depends on SFR rather than $M_{\star}$ (§3.1 and §3.2). For the comparison sample, the situation is the opposite. The tight $\overline{\rm BHAR}$-SFR connection for bulge-dominated galaxies indicates that SMBHs only coevolve with bulges rather than entire host galaxies (S§4.1). The non-causal scenarios of merger averaging are unlikely the origin of the $M_{\rm BH}$-$M_{\rm bulge}$ relation in the local universe. 4. 4.

Our best-fit $\overline{\rm BHAR}$-SFR relation for the bulge-dominated sample is $\log\overline{\rm BHAR}=\log\mathrm{SFR}-(2.48\pm 0.05)$ where the slope is fixed to unity (§3.3). Our best-fit $\overline{\rm BHAR}$/SFR ratio is similar to the observed $M_{\rm BH}/M_{\rm bulge}$ ratio in the local universe (§4.2). This agreement indicates that our observed $\overline{\rm BHAR}$-SFR relation is indeed responsible for the well-known tight $M_{\rm BH}$-$M_{\rm bulge}$ correlation among local galaxies. On the other hand, our findings support that the observed local $M_{\rm BH}$-$M_{\rm bulge}$ relation is not heavily biased. The strong $\overline{\rm BHAR}$-SFR relation among bulge-dominated galaxies indicate lockstep growth of SMBHs and bulges, predicting that the $M_{\rm BH}$-$M_{\rm bulge}$ relation should not have strong redshift dependence.

This paper probes the redshift range of ${z=0.5\text{--}3}$. Future studies can extend our work down to ${z\approx 0.2}$ using the 2 deg2 COSMOS field, or even to the local universe (e.g. Goulding et al., 2017) based on wide surveys, e.g. XMM-XXL (Pierre et al., 2016), Stripe 82X (LaMassa et al., 2016), XMM-SERVS (Chen et al., 2018), and the Chandra Source Catalog (Evans et al., 2010). Compared to distant systems in deep fields, local sources have the advantages of larger sample sizes and more accurate morphological measurements, and these advantages could reduce the uncertainties of the ${\overline{\rm BHAR}}$-SFR relation significantly. In the near future, we will also investigate whether bulge SFR or galaxy compactness is more tightly linked to SMBH growth (Ni et al. in prep.; §4.1). Future work could furthermore derive the full BHAR distribution as a function of SFR and ${M_{\star}}$ for bulge-dominated and comparison galaxies, respectively, and detailed sample properties such as duty cycle and average accretion rate of AGNs can be further obtained and analyzed (§2.4 and Fig. 8). However, such studies will require a large galaxy sample with reliable morphological measurements, and thus deep HST (or future JWST and WFIRST) imaging over much larger fields than CANDELS is needed.

Since our results indicate that SMBHs grow in lockstep with host-galaxy bulges, we also expect a strong connection between ${\overline{\rm BHAR}}$ and bulge SFR for systems that have both bulge and disk components (§4.1 and §4.3). Future ALMA observations could study the ${\overline{\rm BHAR}}$-bulge SFR connection among these systems. ALMA can cover FIR wavelengths down to observed-frame ${300\ \mu}$m, corresponding to the typical SED-peak wavelength (${\approx 100\ \mu}$m, rest-frame) of cold-dust emission of galaxies at ${z\approx 2}$. Therefore, IR luminosities and thereby FIR-based SFR can be reliably estimated for these systems with ALMA. Since ALMA can reach HST-like resolutions, it should be able to separate reliably bulge SFR from total SFR. The strong ${\overline{\rm BHAR}}$-SFR connection among bulge-dominated galaxies might be physically driven by the amount of cold gas available (see §4.1). To test this idea, one could compare the gas masses of high-SFR vs. low-SFR bulge-dominated galaxies with observations by ALMA. ALMA could measure gas masses by observing the CO lines. ALMA could also observe the Rayleigh–Jeans tail of cold dust emission (${\approx 500\ \mu}$m, rest-frame), which is a reliable tracer of dust masses (e.g. Scoville et al., 2017). Gas masses can then be estimated from dust masses with the assumption of a typical dust-to-gas ratio, although uncertainties inevitably exist in this conversion (e.g. Simpson et al., 2015).

Appendix A Results for All Galaxies

In this appendix, we perform analyses like those in §3 for all galaxies with $H<24.5$, including both the bulge-dominated and comparison samples grouped together. The results are presented in Figs. 9, 10, and 11, and Tab. 4. In Fig. 10, we also compare the $\overline{\rm BHAR}$-$M_{\star}$ relation with that derived in Yang et al. (2018a). The $\overline{\rm BHAR}$-$M_{\star}$ relation in this work agrees with the results of Yang et al. (2018a).

From Tab. 4, $\overline{\rm BHAR}$ is more strongly related to $M_{\star}$ than SFR. This is expected, because the comparison sample is the numerically dominant galaxy population (see §2.3) and $\overline{\rm BHAR}$ is mainly related to $M_{\star}$ for the comparison sample, especially at $z=1.5\text{--}3.0$. This conclusion is also qualitatively consistent with Yang et al. (2017), although their statistical significances of the $\overline{\rm BHAR}$-$M_{\star}$ relation are higher than those in Tab. 4. We attribute this difference to the fact that the dynamic range of $M_{\star}$ probed in Yang et al. (2017) is much wider than that in this work ($\log M_{\star}\approx 8\text{--}11$ vs. $\log M_{\star}\approx 10\text{--}11$), since here we require $H<24.5$ to ensure high-quality morphological information for all galaxies. The narrower dynamic range also results in smaller sample sizes, leading to the relatively large statistical scatter in Figs. 9 and 10 (compared to Figs. 4 and 5 in Yang et al. 2017).

Acknowledgements

We thank the referee for helpful feedback that improved this work. We thank Robin Ciardullo, Mara Salvato, Ian Smail, Yongquan Xue, and Wenfei Yu for helpful discussions. GY, WNB, and QN acknowledge support from Chandra X-ray Center grant AR8-19016X, NASA grant NNX17AF07G, NASA ADP grant 80NSSC18K0878, and the V.M. Willaman Endowment. The Chandra Guaranteed Time D.M.A. acknowledges the Science and Technology Facilities Council through grant ST/L00075X/1. FV acknowledge financial support from CONICYT and CASSACA through the Fourth call for tenders of the CAS-CONICYT Fund. Observations (GTO) for the GOODS-N were selected by the ACIS Instrument Principal Investigator, Gordon P. Garmire, currently of the Huntingdon Institute for X-ray Astronomy, LLC, which is under contract to the Smithsonian Astrophysical Observatory via Contract SV2-82024. This project uses Astropy (a Python package; see Astropy Collaboration et al. 2018).