X-rays across the galaxy population - II. The distribution of AGN accretion rates as a function of stellar mass and redshift

TL;DR

This study uses deep X-ray imaging to analyze the distribution of black hole accretion rates across different galaxy types, masses, and redshifts, revealing stochastic AGN fueling and self-regulation at high redshift.

Contribution

It provides the first detailed measurement of the distribution of specific black hole accretion rates across a wide range of galaxy properties and cosmic time, using a Bayesian approach on deep X-ray data.

Findings

AGN activity depends on galaxy stellar mass and redshift.

Probability of hosting an AGN varies between star-forming and quiescent galaxies.

High-redshift quiescent galaxies have an AGN duty cycle of about 20%.

Abstract

We use deep Chandra X-ray imaging to measure the distribution of specific black hole accretion rates ($L_X$ relative to the stellar mass of the galaxy) and thus trace AGN activity within star-forming and quiescent galaxies, as a function of stellar mass (from $10^{8.5}-10^{11.5} M_\odot$) and redshift (to $z \sim 4$). We adopt near-infrared selected samples of galaxies from the CANDELS and UltraVISTA surveys, extract X-ray data for every galaxy, and use a flexible Bayesian method to combine these data and to measure the probability distribution function of specific black hole accretion rates, $\lambda_{sBHAR}$. We identify a broad distribution of $\lambda_{sBHAR}$ in both star-forming and quiescent galaxies---likely reflecting the stochastic nature of AGN fuelling---with a roughly power-law shape that rises toward lower $\lambda_{sBHAR}$, a steep cutoff at $\lambda_{sBHAR} \gtrsim 0.1-1$ (in Eddington equivalent units), and a turnover or flattening at $\lambda_{sBHAR} \lesssim 10^{-3}-10^{-2}$. We find that the probability of a star-forming galaxy hosting a moderate $\lambda_{sBHAR}$ AGN depends on stellar mass and evolves with redshift, shifting toward higher $\lambda_{sBHAR}$ at higher redshifts. This evolution is truncated at a point corresponding to the Eddington limit, indicating black holes may self-regulate their growth at high redshifts when copious gas is available. The probability of a quiescent galaxy hosting an AGN is generally lower than that of a star-forming galaxy, shows signs of suppression at the highest stellar masses, and evolves strongly with redshift. The AGN duty cycle in high-redshift ($z\gtrsim2$) quiescent galaxies thus reaches $\sim$20 per cent, comparable to the duty cycle in star-forming galaxies of equivalent stellar mass and redshift.

Equations

$$\log \mathrm{SF}{\mathrm{R}}_{\mathrm{cut}}[{\mathcal{M}}_{\odot }\mathrm{y}{\mathrm{r}}^{-1}]=-8.9+0.76\log \frac{{\mathcal{M}}_{\ast }}{{\mathcal{M}}_{\odot }}+2.95\log (1+z).$$

$${\lambda }_{\mathrm{sBHAR}}=\frac{k_{\mathrm{bol}}{L}_{\mathrm{X}}}{1.3\times 10^{38}\mathrm{erg}{\mathrm{s}}^{-1}\times 0.002\frac{{\mathcal{M}}_{\ast }}{{\mathcal{M}}_{\odot }}}$$

$$f({\lambda }_{\mathrm{sBHAR}}>0.01)={\int }_{-2}^{\infty }p(\log {\lambda }_{\mathrm{sBHAR}}|{\mathcal{M}}_{\ast },z)d\log {\lambda }_{\mathrm{sBHAR}}$$

$$⟨{\lambda }_{\mathrm{sBHAR}}⟩=\frac{{\int }_{-2}^{\infty }p(\log {\lambda }_{\mathrm{sBHAR}}|{\mathcal{M}}_{\ast },z){\lambda }_{\mathrm{sBHAR}}d\log {\lambda }_{\mathrm{sBHAR}}}{{\int }_{-2}^{\infty }p(\log {\lambda }_{\mathrm{sBHAR}}|{\mathcal{M}}_{\ast },z)d\log {\lambda }_{\mathrm{sBHAR}}}.$$

$$f_{\mathrm{total}}(\lambda )=A_1{f}_{\mathrm{gal}}(\lambda )+A_2{f}_{\mathrm{AGN}}(\lambda )$$

$$k_i=\kappa \times \eta (z_i)\times t_i\times \frac{{\mathcal{M}}_i}{{\mathcal{M}}_{\odot }}$$

$$\mathcal{L}(N_i|{\lambda }_{\mathrm{sBHAR}},L_{\mathrm{G}},b_i,t_i,z_i)=$$

$$\displaystyle\sum_{C_{\mathrm{A}}=0}^{N_{i}}\sum_{C_{\mathrm{G}}}^{N_{i}-C_{\mathrm{A}}}\Bigg{(}$$

$$\frac{(l_i{L}_{\mathrm{G}}{)}^{C_{\mathrm{G}}}}{C_{\mathrm{G}}!}{e}^{-{\text{l}}_i{L}_{\mathrm{G}}}$$

$$\displaystyle\frac{b_{i}^{C_{\mathrm{B}}}}{C_{\mathrm{B}}!}e^{-b_{i}}\Bigg{)}$$

$$C_{\mathrm{A}}+C_{\mathrm{G}}+C_{\mathrm{B}}=N_i$$

$$\pi (L_{\mathrm{G}}|\widetilde{{\mathcal{M}}_i},\widetilde{{\mathrm{SFR}}_i})=\frac{1}{\Gamma (\alpha ){\theta }^{\alpha }}{L}_{\mathrm{G}}^{\alpha -1}{e}^{-L_{\mathrm{G}}/\theta }$$

$$\mathcal{L}({\mathbf{D}}_{\mathrm{bin}})=$$

$$\displaystyle\prod_{i=1}^{n_{\mathrm{source}}}\int_{0}^{\infty}\int_{0}^{\infty}\bigg{[}$$

$$\displaystyle\pi(\lambda\;|\;\mathcal{M}_{\mathrm{bin}},z_{\mathrm{bin}})\bigg{]}\;dL_{\mathrm{G}}\;d\lambda$$

$$\pi (\lambda |{\mathcal{M}}_{\mathrm{bin}},z_{\mathrm{bin}})d\lambda \equiv p(\log {\lambda }_{\mathrm{sBHAR}}|{\mathcal{M}}_{\ast },z)d\log {\lambda }_{\mathrm{sBHAR}}$$

$$\mathcal{L}({\mathbf{D}}_{\mathrm{bin}})=\prod _{i=1}^{n_{\mathrm{source}}}\left[\sum _{j=1}^K{A}_j{\mathcal{W}}_{ij}\right]$$

Full text

††pagerange: X-rays across the galaxy population - II. The distribution of AGN accretion rates as a function of stellar mass and redshift–C††pubyear: 2017

1 Introduction

The properties of galaxies and the supermassive black holes that lie at their centres appear to be fundamentally connected, but the physical origin of this relationship remains unclear (see e.g. Kormendy & Ho, 2013). Supermassive black holes grow primarily through intense periods of mass accretion, during which they are observed as Active Galactic Nuclei (AGNs). The total star formation rate density and the total AGN accretion density appear to follow similar patterns with cosmic time (e.g. Boyle et al., 1998; Delvecchio et al., 2014; Aird et al., 2015), indicating that the overall build up of galaxies and their central black holes proceeds in a coherent manner. However, prior studies have found little evidence for a direct correlation between the current levels of black hole growth (traced by the AGN luminosity) and the wider-scale properties of a galaxy (such as the total stellar mass or the current rate of star formation, see e.g. Azadi et al. 2015 and references therein). To determine the physical mechanisms that promote AGN activity—and determine the impact of AGN feedback on the evolution of galaxies—requires a detailed mapping of the distribution of AGN activity within the evolving galaxy population.

X-ray surveys provide an efficient method of identifying AGNs over a wide range of redshifts and down to relatively low luminosities where the host galaxy dominates the observed light at other wavelengths (see Brandt & Alexander, 2015, for a recent review). A number of studies have thus taken X-ray selected samples of AGNs and examined the properties of the galaxies they lie in. Initial studies found that X-ray selected AGN samples predominantly lie in galaxies with red optical colours, indicating a possible connection between the presence of an AGN and the quenching of star formation throughout a galaxy (e.g. Nandra et al., 2007). Subsequent studies, however, demonstrated that the majority of X-ray AGNs are found in galaxies with moderate-to-high stellar masses ($\mathcal{M}_{*}\gtrsim 3\times 10^{10}\mathcal{M}_{\odot}$), although lower-mass hosts are also identified (e.g. Brusa et al., 2009; Xue et al., 2010). When compared to samples matched in stellar mass, X-ray AGNs either show no preference in terms of host colour or a mild preference to be found in star-forming galaxies (e.g. Silverman et al., 2009; Rosario et al., 2013).

Aird et al. (2012, hereafter A12) took a different approach, starting with a sample of galaxies out to $z\sim 1$ (from the Prism Multi-object Survey: Coil et al., 2011; Cool et al., 2013) and determining the probability of finding an X-ray AGN within such galaxies as a function of the galaxy properties: stellar mass, redshift, and colour. A12 showed that the probability of hosting an AGN can be described by a power-law distribution of specific black hole accretion rates (hereafter, $\lambda_{\mathrm{sBHAR}}$, the rate of black hole accretion normalised relative to the host stellar mass). The power-law distribution was found to be consistent over a broad range of stellar masses ($9.5\lesssim\log\mathcal{M}_{*}/\mathcal{M}_{\odot}\lesssim 12$) but with a normalization that evolves strongly with redshift, indicating that the probability of hosting an AGN with a given $\lambda_{\mathrm{sBHAR}}$ is independent of stellar mass but drops rapidly between $z\sim 1$ and the present day. Thus, the predominance of higher mass host galaxies for samples of X-ray AGN can be ascribed to a selection effect: it is simply harder to identify AGNs with an equivalent $\lambda_{\mathrm{sBHAR}}$ in lower mass galaxies than in higher mass galaxies (see also Aird et al., 2013). Bongiorno et al. (2012) confirmed these findings and extended the analysis to higher redshifts ($z\sim 2.5$). A12 also showed that the probability of hosting an AGN (at a given $\lambda_{\mathrm{sBHAR}}$) is mildly enhanced (by a factor $\sim 2-3$) in galaxies with bluer optical colours. This relationship was explored further by Azadi et al. (2015), who demonstrated an increase in the normalization of the distribution of $\lambda_{\mathrm{sBHAR}}$ with increasing star formation rate (measured relative to the “main sequence of star formation”: Noeske et al. 2007; Elbaz et al. 2011). Georgakakis et al. (2014) also measured the distribution of $\lambda_{\mathrm{sBHAR}}$ for stellar-mass-limited samples of X-ray AGNs, split according to their rest-frame colours, confirming a broad distribution of $\lambda_{\mathrm{sBHAR}}$ and the dominance of AGNs in star-forming galaxies.

These studies have had an important impact on our understanding of the relationship between galaxies and the growth of their central black holes. The broad distribution of accretion rates for samples of galaxies with similar physical properties indicates that the instantaneous level of AGN activity can vary over short timescales relative to the wider-scale properties of the galaxy (Mullaney et al., 2012; Aird et al., 2013; Hickox et al., 2014; Schawinski et al., 2015). Thus, the processes that bring gas into the centres of galaxies and trigger AGN activity appear to be stochastic in nature (e.g. Hopkins & Hernquist, 2006; Novak et al., 2011). This variability of AGN activity can blur out any correlation between the overall levels of accretion and the galaxy properties, such as the star formation rate (e.g. Chen et al., 2013; Stanley et al., 2015). Thus, to determine the underlying relationship between AGN activity and galaxy evolution requires large, statistical samples of galaxies, enabling accurate measurements of the distribution of accretion rates within different types of galaxies throughout cosmic history.

Despite this important progress, a number of open questions remain. A universal, power-law distribution of $\lambda_{\mathrm{sBHAR}}$ is difficult to reconcile with the overall shape and evolution of the total luminosity function of AGNs. A turnover at high $\lambda_{\mathrm{sBHAR}}$ (at a point roughly corresponding to the Eddington limit) and a steep power-law tail is required to recover the break in the X-ray luminosity function (at luminosity $L_{*}$) and the steep bright-end slope (e.g. Aird et al., 2013; Caplar et al., 2015). Additional evolution in the $\lambda_{\mathrm{sBHAR}}$ distribution (e.g. a shift of any break to higher $\lambda_{\mathrm{sBHAR}}$ with increasing redshift) may also be required to recover the overall evolution of $L_{*}$ of the X-ray luminosity function (e.g. Aird et al., 2010, 2015). Furthermore, some level of stellar-mass-dependent evolution of the $\lambda_{\mathrm{sBHAR}}$ distribution may be required to account for changes in the overall shape of the luminosity function toward higher redshifts (e.g. Aird et al., 2013; Bongiorno et al., 2016). Ultimately, accurate measurements of the probability distribution of $\lambda_{\mathrm{sBHAR}}$ in different galaxy populations as a function of stellar mass and redshift are required to build a coherent picture of galaxy–AGN co-evolution.

In this paper, we build on the approach pioneered by A12 to measure the distribution of accretion rates within the galaxy population. We construct a large, stellar-mass-limited sample of $\sim 120,000$ galaxies spanning out to $z\sim 4$ and with a wide range of stellar masses ($8.5\lesssim\log\mathcal{M}_{*}/\mathcal{M}_{\odot}\lesssim 11.5$). We extract X-ray data from deep Chandra imaging for every galaxy in our sample and use a sophisticated Bayesian methodology to probe below the nominal detection limits of the Chandra imaging and recover the intrinsic underlying distribution of $\lambda_{\mathrm{sBHAR}}$. Our method is non-parametric and allows us to divide our galaxy sample according to redshift, stellar mass, and galaxy type (star-forming and quiescent) and track how the incidence of AGNs and the distribution of their accretion rates changes across the evolving galaxy population. Our approach, starting with a large galaxy sample, contrasts with other recent studies that start with an X-ray selected sample of AGNs, determine the properties of their host galaxies, and use either parametric (Bongiorno et al., 2016) or non-parametric (Georgakakis et al., 2017) methods to reconcile the observed X-ray AGN population with independent measurements of the galaxy stellar mass function and thus recover the underlying distribution of $\lambda_{\mathrm{sBHAR}}$ (see also Aird et al., 2013).

This paper is the second in a series, building on the analysis of Aird et al. (2017, hereafter Paper I). In Paper I, we measured the intrinsic distribution of X-ray luminosities for samples of star-forming galaxies, binned by stellar mass and redshift. In each bin we identified two different origins for the X-ray emission:

1. the total emission from X-ray binaries and hot gas, tracing the average star formation rate of our galaxy sub-samples; and
2. the emission from active galactic nuclei (AGN), tracing supermassive black hole accretion activity. In this paper (hereafter Paper II) we focus on the AGN emission. We present measurements of the intrinsic distributions of $\lambda_{\mathrm{sBHAR}}$ within all galaxies, as well as within the star-forming and quiescent galaxy populations, as a function of stellar mass and redshift. We compare the shapes of these distributions and analyse their evolution with redshift and their stellar mass dependence, providing crucial insights into the extent and distribution of black hole growth across the galaxy population.

Paper III (Aird et al. in preparation) will expand this analysis further, presenting measurements of the distribution of $\lambda_{\mathrm{sBHAR}}$ within galaxies as a function of star formation rate (relative to the main sequence of star formation) and thus exploring the connection between the growth of galaxies and their central black holes in greater detail.

Section 2 briefly describes our data sets and the sample selection for this paper. In Section 3 we present our measurements of the probability distribution functions of $\lambda_{\mathrm{sBHAR}}$ as a function of stellar mass and redshift for all galaxies (Section 3.2) as well as within the star-forming and quiescent populations separately (Section 3.3). We also use these distributions to estimate the AGN duty cycle (the fraction of galaxies with black holes growing above a certain $\lambda_{\mathrm{sBHAR}}$) and the average accretion rate of these active black holes, providing further insights into the dependence of AGN activity on stellar mass and redshift (Section 3.4). Section 4 discusses our findings and we summarize our overall conclusions in Section 5. We adopt a flat cosmology with $\Omega_{\Lambda}=0.7$ and $H_{0}=70$ km s*-1* Mpc*-1* throughout this paper and assume a Chabrier (2003) stellar initial mass function (IMF).

2 Data and sample selection

The results in this paper are based on the same datasets used in Paper I, described fully therein and briefly summarized here.

Our work is based on near-infrared (NIR) selected samples of galaxies taken from four of the CANDELS survey fields (GOODS-S, GOODS-N, AEGIS and COSMOS: Grogin et al., 2011; Koekemoer et al., 2011) and the larger UltraVISTA survey of $\sim$1.6 deg2 of the COSMOS field (McCracken et al., 2012). We adopt photometric catalogues from Skelton et al. (2014) and Muzzin et al. (2013) for CANDELS and UltraVISTA respectively.

We match the NIR-selected catalogues with X-ray detected sources in the Chandra imaging of our fields (Alexander et al., 2003; Xue et al., 2011; Nandra et al., 2015; Civano et al., 2016), where all the Chandra data have been reprocessed using our own procedures (see Laird et al., 2009; Georgakakis et al., 2014; Nandra et al., 2015; Aird et al., 2015). We also extract X-ray data (total counts, background counts, and effective exposures) at the positions of the remaining NIR-selected sources that are not matched to a significant X-ray detection, We extract counts within a circular aperture that corresponds to a 70 per cent enclosed energy fraction for the exposure-weighted Chandra point spread function at each source position. These extracted X-ray data contain information on the sensitivity of the Chandra imaging as well as providing information from sources that fall below the nominal X-ray detection threshold that is used in our analysis (see Section 3 and Appendix B below).

In this paper, we use X-ray information from the hard (2–7 keV) band rather than the 0.5–2 keV band adopted in Paper I. While Chandra is less sensitive in the 2–7 keV band, using a harder band reduces biases due to intrinsic absorption and allows us to reliably estimate the X-ray luminosity of AGNs, assuming absorption columns of $N_{\mathrm{H}}\lesssim 10^{23}$ cm*-2*. We estimate intrinsic rest-frame 2–10 keV X-ray luminosities based on the observed count rates in the 2–7 keV band. Following Aird et al. (2015), to convert the observed count rates to an X-ray luminosity we assume an X-ray spectrum with a fixed photon index of $\Gamma=1.9$, a reflection component (using the pexrav model) with fixed strength of $R=1$, and Galactic absorption only. We note that the impact of intrinsic absorption on the observed 2–7 keV band depends on redshift and is thus reduced at $z\gtrsim 1$, although we will still severely underestimate the luminosities of Compton-thick sources ($N_{\mathrm{H}}\gtrsim 10^{24}$ cm*-2*). A full and accurate characterisation of intrinsic absorption—including the resulting selection biases—requires advanced techniques (e.g. Buchner et al., 2015; Aird et al., 2015) and is beyond the scope of this paper.

We compile high-resolution spectroscopic redshifts (see Skelton et al., 2014; Muzzin et al., 2013; Aird et al., 2015, 2017; Marchesi et al., 2016, and references therein for full details), low-resolution spectroscopic redshifts (i.e. from PRIMUS, Coil et al. 2011, or 3DHST, Momcheva et al. 2016) and photometric redshifts (including high quality AGN photo-$z$, where available, from Hsu et al., 2014; Nandra et al., 2015; Marchesi et al., 2016) for sources in our NIR-selected catalogues. Rest-frame colours are calculated using EaZY (Brammer et al., 2011). We use the FAST code (Kriek et al., 2009) to fit the UV-to-NIR spectral energy distributions (SEDs) of our sources, as described in the appendix of Paper I. Appendix A of this paper (see below) describes our modifications of FAST to account for an AGN component in addition to a galaxy component when fitting the SEDs. We only allow for an AGN component for objects with significant X-ray detections. We also estimate star-formation rates (SFRs) based on either the sum of the UV and IR emission (for sources with 24µm detections) or the dust-corrected SFR estimate from the SED fit (see Paper I for details). These SFRs are corrected for any AGN contribution as described in Appendix A below.

To define our sample of galaxies, we apply the magnitude limits and the stellar-mass completeness limits described in Paper I. While Paper I focused on star-forming galaxies, in this paper we measure the distribution of AGN activity across the entire galaxy population as a function of stellar mass, redshift, and galaxy type. We classify galaxies as star-forming or quiescent based on their SFRs relative to the evolving “star-forming main sequence”. Figure 1 shows the distribution of SFR and $\mathcal{M}_{*}$ for our galaxy sample at different redshifts (pink crosses correspond to directly detected X-ray sources). A star-forming main sequence with a roughly constant slope and a normalization that increases toward higher redshifts is apparent in our samples out to at least $z\sim 3$ and is reasonably well described by the relation given by equation 6 in Paper I. We define a cut (shown by the green dotted lines in Figure 1) at 1.3 dex below the evolving main sequence from Paper I, which roughly corresponds to the minimum in the distribution of SFRs relative to the main sequence and is given by

[TABLE]

We classify galaxies that fall below this cut as quiescent; those above the green dotted line are classified as star-forming galaxies. The vertical grey dashed lines in Figure 1 indicate the stellar mass completeness limits we apply for the CANDELS and UltraVISTA samples. We note that the majority of quiescent galaxies in our sample are identified in the larger-area (but shallower) UltraVISTA imaging; the paucity of lower-mass quiescent galaxies at high-redshifts in our sample is due to a combination of the intrinsic rarity of such galaxies and the lack of data to the required (i.e. CANDELS) depth over a sufficiently large area.

Sources identified as “AGN-dominated” at optical to IR wavelengths based on our two-component SED fitting (i.e. where more than 50 per cent of the light at rest-frame 5000Å is from the AGN) are included in our sample of all galaxies. For the majority of these sources ($\sim 77$ per cent) the galaxy template dominates the emission at $\sim 1$µm (in our best fit), indicating that we can still estimate the stellar mass of the host galaxy to a reasonable degree of accuracy. However, we exclude all AGN-dominated sources (which constitute $\lesssim 20$ per cent of the X-ray detected samples and $\lesssim 0.3$ per cent of the overall galaxy sample) when considering the star-forming or quiescent samples separately as it is not possible to reliably measure the SFR and thus classify these sources into either group. In Appendix C we show that excluding the AGN-dominated sources from this analysis has a negligible impact on our overall results and conclusions.

Our final sample consists of 126,971 galaxies (above our stellar-mass completeness limits), of which 106,636 are classified as star-forming, 19,971 as quiescent, and 364 as AGN-dominated (precluding a classification of the host galaxy type). A total of 1,797 sources in our sample are directly detected in the 2–7 keV X-ray band (including 318 of the sources classified as AGN-dominated).

3 The distribution of AGN accretion rates across the galaxy population

In this section, we present measurements of the distribution of AGN accretion rates as a function of galaxy type, stellar mass and redshift. In Section 3.1 we give a brief description of our methodology (further details are given in Appendix B). Section 3.2 presents measurements of the distributions of specific black hole accretion rates for all galaxies as a function of stellar mass and redshift, while Section 3.3 presents measurements for the star-forming and quiescent populations separately. In Section 3.4 we calculate the AGN duty cycle and the average accretion rate (above a fixed limit), summarizing our results on how the accretion rate distributions depend on host galaxy properties.

3.1 Methodology

Following A12, throughout this paper we measure the distribution of specific black hole accretion rates ($\lambda_{\mathrm{sBHAR}}$), the rate of accretion onto the central supermassive black hole scaled relative to the stellar mass of the host galaxy. Defining AGN activity in this manner allows us to trace the wide distribution in accretion rates within galaxies of fixed stellar mass. We also account for the bias whereby a higher mass galaxy with a given $\lambda_{\mathrm{sBHAR}}$ appears more luminous than a lower mass galaxy that is growing its central black hole at a similar rate. Based on the same approach, A12 found that at $z\lesssim 1$ the incidence of AGN activity was defined by a single, power-law distribution of $\lambda_{\mathrm{sBHAR}}$ across a wide range of galaxy stellar masses ($\mathcal{M}_{*}\sim 10^{9.5-12}\mathcal{M}_{\odot}$), indicating that $\lambda_{\mathrm{sBHAR}}$ may be a more fundamental property of AGNs than the directly observed luminosity.

We define $\lambda_{\mathrm{sBHAR}}$ in dimensionless units, such that

[TABLE]

where $L_{\mathrm{X}}$ is the rest-frame 2–10 keV X-ray luminosity, $k_{\mathrm{bol}}$ is a bolometric correction factor (we adopt a constant $k_{\mathrm{bol}}=25$ in this paper) and $\mathcal{M}_{*}$ is the total stellar mass of the AGN host galaxy (estimated from our SED fitting). The additional scale factors are chosen such that $\lambda_{\mathrm{sBHAR}}\approx\lambda_{\mathrm{Edd}}$, the Eddington ratio, assuming that the mass of the central black hole scales directly with the total host stellar mass. Thus, the $\lambda_{\mathrm{sBHAR}}$ distribution can be regarded as a tracer of the distribution of Eddington ratios. However, these scaling relations are expected to have significant intrinsic scatter and may break down in certain parameter regimes (e.g. at lower stellar masses). Nonetheless, $\lambda_{\mathrm{sBHAR}}$ provides a meaningful measurement of the rate at which a given galaxy is growing its black hole, relative to the mass of the galaxy. With a fixed $k_{\mathrm{bol}}$, with our definition $\lambda_{\mathrm{sBHAR}}\propto L_{\mathrm{X}}/\mathcal{M}_{*}$ (see Bongiorno et al., 2012, 2016).

To trace the distribution of $\lambda_{\mathrm{sBHAR}}$  we measure the probability density function of $\log\lambda_{\mathrm{sBHAR}}$ at a given stellar mass and redshift, $p(\log\lambda_{\mathrm{sBHAR}}\;|\;\mathcal{M}_{*},z)$.

Our measurement of $p(\log\lambda_{\mathrm{sBHAR}}\;|\;\mathcal{M}_{*},z)$ describes the probability that a galaxy with a given stellar mass and redshift hosts an AGN with a given specific accretion rate. Equivalently, this quantity reflects the distribution of specific accretion rates within a sample of galaxies in a fixed range of stellar mass and redshift.

We adopt the Bayesian mixture modelling approach described in appendix B of Paper I, adapted to consider $\lambda_{\mathrm{sBHAR}}$ rather than $L_{\mathrm{X}}$ (see Appendix B below for more details). Our method provides a flexible method of recovering the underlying distribution that does not assume a particular functional form and has relatively few constraints, other than a prior that prefers a smoothly varying distribution and a requirement that the total probability distribution function must integrate to less than 1.

Our method accounts for uncertainties in individual X-ray flux measurements and the varying sensitivity of the Chandra X-ray imaging over the CANDELS and UltraVISTA fields.

The variation in the X-ray sensitivity of the deep Chandra surveys—varying both between and within our fields—is extremely well characterized and is driven by the vignetting of the telescope and degradation of the point spread function with off-axis angle, in addition to the overall exposure time. Our statistical approach uses the available X-ray information for every galaxy in our sample, including the observed counts and the sensitivity information at the position of galaxies that fall below the nominal X-ray detection threshold. By using the available X-ray information from every galaxy in our sample we are able to place much tighter constraints on $p(\log\lambda_{\mathrm{sBHAR}}\;|\;\mathcal{M}_{*},z)$. In appendix C of Paper I we performed a number of simulations to verify that our method is able to combine the variable-depth Chandra data and accurately recover the shape of different underlying distributions.

At the lowest X-ray luminosities probed by our data, galactic non-AGN processes (high- and low-mass X-ray binaries, hot gas, supernovae) can make a significant contribution to the observed luminosity, as explored in Paper I. We have thus updated our statistical method to allow for the expected contribution from galactic processes to the observed X-ray counts, based on our multiwavelength measurements of SFR and $\mathcal{M}_{*}$ for each galaxy in our sample. We correct for this potential contamination when recovering the distribution of AGN accretion rates, $p(\log\lambda_{\mathrm{sBHAR}}\;|\;\mathcal{M}_{*},z)$. The updates to our statistical method are described in Appendix B below. However, our estimates of $p(\log\lambda_{\mathrm{sBHAR}}\;|\;\mathcal{M}_{*},z)$ remain highly uncertain in regimes where the average galactic emission dominates over the X-ray AGN signal. In Figure 2 (described fully in Section 3.2 below), we use dashed and dotted lines to indicate regimes where $\lambda_{\mathrm{sBHAR}}$ corresponds to a luminosity within 0.5 dex of the “X-ray main sequence of star formation”, identified in Paper I. Our measurements of $p(\log\lambda_{\mathrm{sBHAR}}\;|\;\mathcal{M}_{*},z)$ should be treated with caution in these regimes. In subsequent plots we only consider results above these limits, where we can reliably trace the distribution of AGN emission.

3.2 The distribution of accretion rates within all galaxies as a function of stellar mass and redshift

Figure 2 presents measurements of $p(\log\lambda_{\mathrm{sBHAR}}\;|\;\mathcal{M}_{*},z)$ for our full galaxy sample as a function of stellar mass and redshift. The thick coloured lines and hatched regions show our best estimates of $p(\log\lambda_{\mathrm{sBHAR}}\;|\;\mathcal{M}_{*},z)$ and the 90 per cent confidence intervals, respectively, as recovered by our Bayesian modelling. The grey dashed histograms indicate the observed distribution of $\lambda_{\mathrm{sBHAR}}$ based on the X-ray detected population only without any corrections for the X-ray incompleteness, which differ significantly from our estimates of the intrinsic distributions at lower values of $\lambda_{\mathrm{sBHAR}}$ that use the available X-ray data and sensitivity information from all galaxies in a given sample (see Section 3.1 above).

The results shown in Figure 2 reveal the wide distribution of $\lambda_{\mathrm{sBHAR}}$, rising towards lower $\lambda_{\mathrm{sBHAR}}$ and spanning at least $2-3$ orders of magnitude, within galaxies of all stellar masses and redshift ranges probed here. Generally, the distributions have a relatively steep cutoff at the highest accretion rates ($\log\lambda_{\mathrm{sBHAR}}\gtrsim-1$, but dependent on stellar mass and redshift), a flatter distribution at moderate accretion rates ($-3\lesssim\log\lambda_{\mathrm{sBHAR}}\lesssim-1$), and in some cases a turnover towards lower $\lambda_{\mathrm{sBHAR}}$ resulting in a very broad peak. The orange lines in the low-redshift ($z<1.5$) panels of Figure 2 indicate the power-law relation measured by A12 (dashed lines indicate an extrapolation of the measured relation), which is independent of stellar mass but has a normalization that evolves with redshift. While this relation is broadly consistent with our measurements, our updated analysis reveals that the accretion rate distribution has a more complicated shape. There is also clear evidence of a stellar-mass dependence, at least at $z>0.5$, whereby the probability of hosting an AGN of $-3\lesssim\log\lambda_{\mathrm{sBHAR}}\lesssim-1$ is higher than the A12 relation at high stellar masses ($\log\mathcal{M}_{*}/\mathcal{M}_{\odot}\gtrsim 10.5$) and is below the extrapolation at low stellar masses ($\log\mathcal{M}_{*}/\mathcal{M}_{\odot}\lesssim 9.5$).

3.3 The distribution of accretion rates within the star-forming and quiescent galaxy populations

The overall galaxy population is built up by a mixture of star-forming and quiescent galaxies that varies with stellar mass and redshift. The incidence and distribution of AGN activity could be very different within the different galaxy populations given their different gas content, morphologies and star formation histories that potentially complicate a physical interpretation of the overall distributions traced in Figure 2. Thus, in Figure 3 we divide our galaxy sample into star-forming and quiescent populations based on their SFRs (see Section 2 for details) and present measurements of $p(\log\lambda_{\mathrm{sBHAR}}\;|\;\mathcal{M}_{*},z)$ within each population, as a function of stellar mass and redshift.

We find AGNs with a wide range of $\lambda_{\mathrm{sBHAR}}$ in both star-forming and quiescent galaxies at all redshifts and stellar masses (see also Azadi et al., 2015). In the majority of the stellar mass and redshift panels of Figure 3, we are unable to identify any significant differences in $p(\log\lambda_{\mathrm{sBHAR}}\;|\;\mathcal{M}_{*},z)$ between the two galaxy populations, given the uncertainties. However, at higher stellar masses ($\log\mathcal{M}_{*}/\mathcal{M}_{\odot}\gtrsim 10.0$) and lower redshifts ($z\lesssim 2.0$), where we have the best constraints, there are significant differences whereby $p(\log\lambda_{\mathrm{sBHAR}}\;|\;\mathcal{M}_{*},z)$ is higher at $\log\lambda_{\mathrm{sBHAR}}\gtrsim-2$ in star-forming galaxies. Thus, the probability of finding a rapidly accreting AGN is higher by up to a factor $\sim$5 in a star-forming galaxy than in a quiescent galaxiy of equivalent stellar mass and redshift.111 We note that quiescent galaxies dominate the overall number density of galaxies at high stellar masses ($\log\mathcal{M}_{*}/\mathcal{M}_{\odot}\gtrsim 10.5$) and low redshifts ($z\lesssim 1$) and thus there are more quiescent galaxies in our samples, as indicated in Figure 3, whereas star-forming galaxies dominate at low masses and higher redshifts. Our measurements, however, track the fraction of AGNs within a given galaxy sample.

The turnover in $p(\log\lambda_{\mathrm{sBHAR}}\;|\;\mathcal{M}_{*},z)$ at low $\lambda_{\mathrm{sBHAR}}$ is also seen more clearly for the star-forming galaxy samples. The entire distribution may be shifted to higher $\lambda_{\mathrm{sBHAR}}$ for star-forming galaxies compared to the quiescent galaxies and $p(\log\lambda_{\mathrm{sBHAR}}\;|\;\mathcal{M}_{*},z)$ may rise more steeply to lower $\lambda_{\mathrm{sBHAR}}$ in quiescent galaxies, although the uncertainties are large.

In Figure 4 we examine the redshift evolution of $p(\log\lambda_{\mathrm{sBHAR}}\;|\;\mathcal{M}_{*},z)$ at a given stellar mass within the entire galaxy population (top row), star-forming galaxies (middle row) and quiescent galaxies (bottom row), reproducing the results from Figures 2 and 3 to reveal the dependence on redshift. There is no evidence for evolution in $p(\log\lambda_{\mathrm{sBHAR}}\;|\;\mathcal{M}_{*},z)$ at low stellar masses ($\log\mathcal{M}_{*}/\mathcal{M}_{\odot}<10.0$) for any of the galaxy populations, although the uncertainties are large and we are only sensitive to a limited redshift range at these masses. However, we observe strong evolution of $p(\log\lambda_{\mathrm{sBHAR}}\;|\;\mathcal{M}_{*},z)$ in both star-forming and quiescent galaxies in all of the higher ($\log\mathcal{M}_{*}/\mathcal{M}_{\odot}>10.0$) mass bins, which combine to produce the observed evolution for all galaxies. The evolution appears to be characterised by an overall shift of $p(\log\lambda_{\mathrm{sBHAR}}\;|\;\mathcal{M}_{*},z)$ to higher $\lambda_{\mathrm{sBHAR}}$ at higher redshifts for both star-forming and quiescent galaxies. Interestingly, the evolution appears to be somewhat stronger in massive star-forming galaxies from $z\sim 0.1-1$, with a rapid shift in the distribution toward higher $\lambda_{\mathrm{sBHAR}}$, but with little or only mild evolution to higher redshifts. In massive quiescent galaxies, evolution is seen over the entire redshift range, although the distributions are generally at lower $\lambda_{\mathrm{sBHAR}}$, on average, compared to star-forming galaxies of equivalent mass and redshift. We attempt to quantify these trends using estimates of the AGN duty cycle in Section 3.4 below.

Figure 5 reproduces the results from Figures 2 and 3 to reveal the stellar mass dependence in fixed redshift bins. For the quiescent galaxy population (bottom row of Figure 5) our constraints on $p(\log\lambda_{\mathrm{sBHAR}}\;|\;\mathcal{M}_{*},z)$ are generally consistent with a single distribution with no stellar-mass dependence, although at lower redshifts ($z\lesssim 1.5$) we find that $p(\log\lambda_{\mathrm{sBHAR}}\;|\;\mathcal{M}_{*},z)$ may be shifted to lower $\lambda_{\mathrm{sBHAR}}$ at the highest stellar masses ($\mathcal{M}_{*}\gtrsim 10^{11}\mathcal{M}_{\odot}$). Furthermore, our constraints on $p(\log\lambda_{\mathrm{sBHAR}}\;|\;\mathcal{M}_{*},z)$ are poor at $\log\mathcal{M}_{*}/\mathcal{M}_{\odot}\lesssim 10.0$ as our samples of quiescent galaxies at these masses are relatively small and thus we cannot rule out a stellar-mass dependence. For star-forming galaxies (middle row of Figure 5), we find clear evidence of a stellar-mass dependence in $p(\log\lambda_{\mathrm{sBHAR}}\;|\;\mathcal{M}_{*},z)$ at a fixed redshift, at least at $0.5\lesssim z\lesssim 2$ where we have good constraints. At a given redshift, the maximum accretion rate appears to be limited by the same steep power law, while at lower accretion rates the distribution has a mass-dependent normalisation such that the probability of a star-forming galaxy hosting an AGN with $-3\lesssim\log\lambda_{\mathrm{sBHAR}}\lesssim-1$ is suppressed at lower stellar masses. The total galaxy population (top row of Figure 5) also shows evidence for a stellar-mass dependence at a fixed redshift (to at least $z\sim 2$). The mixing of the star-forming and quiescent populations leads to the complex structure and mass-dependence seen in $p(\log\lambda_{\mathrm{sBHAR}}\;|\;\mathcal{M}_{*},z)$ for the overall galaxy sample that is difficult to interpret without considering the constituent populations separately.

3.4 The AGN duty cycle and the average AGN accretion rate

In this section we attempt to summarize and quantify the information contained in Figures 2–5 by calculating the duty cycle of AGN activity as a function of stellar mass as well as the average accretion rate associated with these AGNs. We are thus able to condense the information contained in the full accretion rate distributions down to just two numbers, making it easier to identify and interpret the overall trends with redshift, stellar mass or galaxy type, albeit at the loss of some information on the overall shapes of the distributions as a function of $\lambda_{\mathrm{sBHAR}}$.

We define the AGN duty cycle as the fraction of galaxies (at a given range of stellar mass and redshift) with supermassive black holes growing at $\lambda_{\mathrm{sBHAR}}>0.01$. We can calculate the AGN duty cycle from our measurements of $p(\log\lambda_{\mathrm{sBHAR}}\;|\;\mathcal{M}_{*},z)$ described above,

[TABLE]

where $f(\lambda_{\mathrm{sBHAR}}>0.01)$ corresponds to our definition of the duty cycle. Uncertainties in the duty cycle are propagated from the posterior distributions of $p(\log\lambda_{\mathrm{sBHAR}}\;|\;\mathcal{M}_{*},z)$ and correspond to the 1$\sigma$-equivalent (i.e. 68 per cent) central confidence intervals.

Figure 6 presents measurements of the duty cycle, $f(\lambda_{\mathrm{sBHAR}}>0.01)$, as a function of redshift for different stellar masses. A number of trends, previously identified in the full accretion distributions in Sections 3.2 and 3.3 above, are apparent in this figure. Within star-forming galaxies, we observe a clear stellar mass dependence between $\log\mathcal{M}_{*}/\mathcal{M}_{\odot}\approx 10.0$ and $\log\mathcal{M}_{*}/\mathcal{M}_{\odot}\approx 11.5$, whereby the duty cycle is higher in higher stellar mass galaxies at any given redshift. At high stellar masses ($\log\mathcal{M}_{*}/\mathcal{M}_{\odot}\gtrsim 10.0$), we find that the duty cycle evolves with redshift, increasing up to $z\sim 2$ and flattening off or potentially declining to higher redshifts. The redshift evolution appears to get weaker at lower stellar masses and in our lowest stellar mass bins ($\log\mathcal{M}_{*}/\mathcal{M}_{\odot}<9.5$) our results are consistent with no evolution (although we only probe a limited range of redshifts and the uncertainties are large). For quiescent galaxies, conversely, the duty cycle is generally consistent with no stellar mass dependence between $\log\mathcal{M}_{*}/\mathcal{M}_{\odot}\approx 10$ and $\log\mathcal{M}_{*}/\mathcal{M}_{\odot}\approx 11$ and a drop at the highest stellar mass ($11<\log\mathcal{M}_{*}/\mathcal{M}_{\odot}<11.5$), at least out to $z\sim 2$. The data for quiescent galaxies are consistent with a strong evolution with redshift for all stellar masses between $\log\mathcal{M}_{*}/\mathcal{M}_{\odot}=10.0$ and $\log\mathcal{M}_{*}/\mathcal{M}_{\odot}=11.5$; we do not show data for the lower stellar mass bins where there are very few quiescent galaxies and the duty cycle is poorly constrained.

In Figure 7 we directly compare the redshift evolution of the duty cycle in star-forming and quiescent galaxies at a range of stellar masses. For the lowest mass bin we include all galaxies (of a given type) with $8.5<\log\mathcal{M}_{*}/\mathcal{M}_{\odot}<10.0$ and above our stellar mass completeness limits. At these low masses there is no evidence of a difference in the duty cycle between star-forming or quiescent galaxies or any evolution with redshift (although we cannot rule out an increase of the duty cycle with redshift for quiescent galaxies, given the large uncertainties). At $10.0<\log\mathcal{M}_{*}/\mathcal{M}_{\odot}<10.5$ we find that the duty cycle increases between $z\sim 0$ and $z\sim 1-2$ for both galaxy types and is roughly constant to higher redshifts. The duty cycle is consistent between the the two galaxy types at $z>1.0$ (within the uncertainties) but at lower redshifts we measure a significantly lower duty cycle (by a factor $\sim 2-6$) in quiescent galaxies compared to star-forming galaxies. In the two highest mass bins in Figure 7 we find that the duty cycle is substantially lower in quiescent galaxies at $z<2$ (by a factor $\sim 10$) compared to star-forming galaxies but evolves more rapidly with redshift. The full distributions shown in Figure 3 confirm this trend. At $z>2$ the duty cycles in star-forming and quiescent galaxies are fairly similar in the two high mass bins, although we note that at these redshifts there are substantially fewer quiescent galaxies than star-forming galaxies (even at these high stellar masses). Thus, the bulk of AGN at high redshifts are hosted by star-forming galaxies, despite the similar duty cycles in both galaxy types.

To supplement the AGN duty cycle and gain further insights into the accretion rates of the AGNs contained within our galaxy samples, we calculate the average accretion rate, $\langle\lambda_{\mathrm{sBHAR}}\rangle$, for galaxies that contain AGNs accreting above a fixed limit of $\lambda_{\mathrm{sBHAR}}>0.01$, given by

[TABLE]

While the duty cycle represents the fraction of galaxies with AGNs, the average accretion rate defined in this way provides information on the typical accretion rate of those AGNs, where in both cases we are effectively defining an “AGN” as a supermassive black hole that is accreting above a rate of $\lambda_{\mathrm{sBHAR}}=0.01$. This quantity thus allows us to distinguish between galaxies with similar overall duty cycles but where the shapes of $p(\log\lambda_{\mathrm{sBHAR}}\;|\;\mathcal{M}_{*},z)$(above the $\lambda_{\mathrm{sBHAR}}=0.01$ limit) are very different, e.g. where the bulk of those AGNs are close to the limit versus at much higher $\lambda_{\mathrm{sBHAR}}$.

Our estimates of $\langle\lambda_{\mathrm{sBHAR}}\rangle$ are presented in Figure 8. Our measurements reveal a clear pattern within both star-forming galaxy and all galaxy samples (dominated by star-forming galaxies at all stellar masses for $z\gtrsim 1$), whereby the average accretion rate is increasing with redshift at all stellar masses (increasing by $\sim 0.5$ dex between $z\sim 0.5$ and $z\sim 3$) and becomes progressively higher with decreasing stellar mass at a fixed redshift.

This behaviour contrasts with the evolution of the duty cycle (see Figure 6), where the duty cycle becomes progressively lower with decreasing stellar mass. Thus, while the fraction of galaxies hosting an AGN is lower in lower massive galaxies, the typical accretion rates of that small fraction of AGNs are generally higher. The underlying cause of this pattern is more apparent in the full accretion rate distributions (e.g. central row of Figure 5), where we can see at a fixed redshift (e.g. $0.5<z<1.0$) that $p(\log\lambda_{\mathrm{sBHAR}}\;|\;\mathcal{M}_{*},z)$ is suppressed at lower accretion rates in lower mass galaxies (reducing the overall duty cycle) but that the distribution extends to higher values (the break in the distribution is at higher $\lambda_{\mathrm{sBHAR}}$), resulting in the increase in $\langle\lambda_{\mathrm{sBHAR}}\rangle$. In the discussion below (Section 4) we present a simplified sketch that attempts to summarize these overall trends (see Figure 9).

Our measurements of $\langle\lambda_{\mathrm{sBHAR}}\rangle$ within quiescent galaxies (right panel of Figure 8) are fairly noisy, although at $10<\log\mathcal{M}_{*}/\mathcal{M}_{\odot}<10.5$ and $10.5<\log\mathcal{M}_{*}/\mathcal{M}_{\odot}<11$ (yellow and brown points) we do see a rise of the average accretion rate between $z\sim 0.5$ and $=z\sim 3$, similar to that seen in the star-forming galaxy sample. We note that the average accretion rates are generally lower for the quiescent galaxies than for star-forming galaxies at the same redshift and stellar mass, consistent with the differences seen in Figure 3. We also caution that our measurements are most uncertain in our lowest redshift bin ($0.1<z<0.5$) for all galaxy types, where the duty cycle is generally lowest and it is hard to identify clear trends.

In summary, our measurements show clear evidence that the duty cycle of AGN, defined above a fixed $\lambda_{\mathrm{sBHAR}}$ threshold, increases with stellar mass in star-forming galaxies, while the typical accretion rates of these AGNs are lower at higher stellar masses. The duty cycle in quiescent galaxies is generally lower than in star-forming galaxies and decreases at the highest stellar masses ($11<\log\mathcal{M}_{*}/\mathcal{M}_{\odot}<11.5$ at $z\lesssim 2$). Both the duty cycle and the average accretion rates increase with redshift in star-forming and quiescent galaxies. The strongest evolution of the duty cycle is observed in high stellar mass quiescent galaxies ($11<\log\mathcal{M}_{*}/\mathcal{M}_{\odot}<11.5$), where the duty cycle increases from $f(\lambda_{\mathrm{sBHAR}}>0.01)\approx 0.1$ per cent at $z\approx 0.3$ to $f(\lambda_{\mathrm{sBHAR}}>0.01)\approx 10$ per cent at $z\approx 2$. In Figure 9 below we attempt to encapsulate all these trends into a simplified sketch summarizing our measurements of the full accretion rate probability distributions, $p(\log\lambda_{\mathrm{sBHAR}}\;|\;\mathcal{M}_{*},z)$, as a function of stellar mass and redshift in both the star-forming and quiescent galaxy populations (see Section 4 for a full discussion).

4 Discussion

In this paper we have presented new measurements of the probability distribution function of AGN accretion rates, $p(\log\lambda_{\mathrm{sBHAR}}\;|\;\mathcal{M}_{*},z)$222For simplicity, we denote $p(\log\lambda_{\mathrm{sBHAR}}\;|\;\mathcal{M}_{*},z)$ as $p(\lambda_{\mathrm{sBHAR}})$ throughout the remainder of this discussion section, implicitly assuming logarithmic units and that the probability distribution is defined within a sample of galaxies of given $\mathcal{M}_{*}$ and $z$.,

as a function of stellar mass and redshift within both the star-forming and quiescent galaxy populations. Our method (developed in the companion Paper I) starts with near-infrared selected samples of galaxies and extracts X-ray information for every galaxy from deep Chandra imaging. We are thus able to measure the probability of hosting an AGN directly (as a function of $\lambda_{\mathrm{sBHAR}}$) for different galaxy types and place robust constraints on the shape of the accretion rate distribution (using a non-parametric model with minimal assumptions) over a wide range in both stellar mass and redshift. Our measurements reveal a somewhat complex picture where $p(\lambda_{\mathrm{sBHAR}})$ differs between galaxy types and depends on both stellar mass and redshift. Here, we present a simplified summary of our findings, discuss the possible physical underpinnings of the observed behaviour, and compare our measurements with previous studies.

4.1 The distribution and evolution of AGN activity in star-forming galaxies

In the top row of Figure 9 we attempt to summarize our main findings for the shape, redshift evolution and stellar mass dependence of $p(\lambda_{\mathrm{sBHAR}})$ in star-forming galaxies. We note that this figure is intended as an illustration to guide the discussion and interpretation of our results, as recovered by our non-parametric method. A full parametric modelling of our data, informed by these findings, is deferred to future work.

In low-redshift ($z\sim 0.5$) star-forming galaxies (illustrated in the top-left panel of Figure 9) our results are consistent with a broad distribution of accretion rates with a relatively shallow slope at lower accretion rates (with exponent $\sim 0.2-0.8$, consistent with prior findings e.g. A12; Bongiorno et al. 2012; Azadi et al. 2015), a break at $\lambda_{\mathrm{sBHAR}}\approx 0.1$, and a steeper power law extending to higher redshifts (as inferred by Aird et al. 2013; Caplar et al. 2015; Bongiorno et al. 2016; see also Georgakakis et al. 2017). This broad distribution likely reflects long-term variability in the levels of AGN accretion within a given star-forming galaxy over timescales $\lesssim 100$ Myr (Hickox et al., 2014; Schawinski et al., 2015) due to the stochastic nature of the supply of gas to the nuclear regions of galaxies, in addition to localised heating and feedback from either the AGN itself or stellar processes (see e.g. Hopkins & Hernquist, 2006; Novak et al., 2011; Gaspari et al., 2017). At these redshifts ($z\sim 0.5$), the break at $\lambda_{\mathrm{sBHAR}}\approx 0.1$ indicates a lack of high accretion rate events, possibly related to the distribution of cold molecular gas within low-redshift star-forming galaxies (and a lack of massive gas clouds) or the rarity of violent triggering events (e.g. mergers) that can drive sufficiently high gas densities into the central regions to fuel the highest accretion rate AGN (e.g. Kauffmann & Haehnelt, 2000; Hopkins et al., 2008)

We also identify a possible turnover or flattening of $p(\lambda_{\mathrm{sBHAR}})$ at the lowest accretion rates ($\lambda_{\mathrm{sBHAR}}\lesssim 10^{-2}-10^{-3}$), at least at moderate-to-high stellar masses ($\log\mathcal{M}_{*}/\mathcal{M}_{\odot}\gtrsim 10$) where we are able to probe sufficiently low accretion rates out to $z\sim 2$.

We note that the position of the low-$\lambda_{\mathrm{sBHAR}}$ turnover is close to the regime where star formation processes dominate (see e.g. Figure 2), making it difficult to robustly determine the full extent of any turnover. An eventual flattening or turnover at low accretion rates is not surprising. Indeed, a flattening is required to ensure that the probability distribution function integrates to less than unity – a reasonable prior that is implicit in our Bayesian modelling. However, the position of the observed turnover is not driven by this prior as the integral of $p(\lambda_{\mathrm{sBHAR}})$ is generally much less than 1 over the range of $\lambda_{\mathrm{sBHAR}}$ that we probe (i.e. we do not reach such low accretion rates that every galaxy effectively hosts an AGN). The low-$\lambda_{\mathrm{sBHAR}}$ turnover may therefore represent a natural lower limit to the range of accretion rates in galaxies, indicating that a minimal fuelling level produces a short-lived radiatively efficient AGN (observable in the X-ray) that consumes the gas in the central region and rapidly fades, with the black hole then remaining in a state of relative quiescence for a period of time333Very low accretion rate, radiatively inefficient AGN activity, seen in the radio, may be present during this time but is not revealed by our X-ray–based measurements. before a subsequent fuelling event triggers the next AGN episode.

Our measurements also reveal a stellar mass dependence whereby the probability of a lower mass galaxy hosting a low-accretion-rate AGN ($\lambda_{\mathrm{sBHAR}}\lesssim 0.1$) is suppressed compared to the probability within higher mass galaxies (as indicated in Figure 9), although a broad range of accretion rates—spanning the range observed in higher mass galaxies—are still observed. This suppression of low $\lambda_{\mathrm{sBHAR}}$ fuelling events could indicate that it is more difficult to drive smaller pockets of gas into the central regions of lower mass galaxies, possibly due to the difficulty of releasing angular momentum from such gas clouds in lower mass galaxies (e.g. Rosas-Guevara et al., 2015) or an increase in the effects of stellar feedback at low galaxy masses (e.g. Dekel & Silk, 1986; White & Frenk, 1991; Cole et al., 2000; Hopkins et al., 2014). The nature of this low-$\lambda_{\mathrm{sBHAR}}$ suppression reduces the overall probability of hosting an AGN at lower $\lambda_{\mathrm{sBHAR}}$ (resulting in a lower duty cycle in low-mass star-forming galaxies, see Figure 6), yet the observed range of accretion rates remains the same and there is thus an increase in the average accretion rate (see Figure 8).

As redshift increases (from $z\sim 0.5$ to $z\sim 1$), the accretion rate distributions in star-forming galaxies appear to shift systematically toward higher $\lambda_{\mathrm{sBHAR}}$ (as illustrated in Figure 9; see also the measurements in Figure 4). We propose that this shift is related to the substantial increase in the cold gas fraction toward higher redshifts within star-forming galaxies at all stellar masses. Such an increase is also thought to drive the evolution of the main sequence of star formation and the associated increase in the typical SFRs of star-forming galaxies to higher redshifts (e.g. Lagos et al., 2011; Combes et al., 2013; Saintonge et al., 2013; Santini et al., 2014). This evolution of $p(\lambda_{\mathrm{sBHAR}})$ toward higher $\lambda_{\mathrm{sBHAR}}$ drives the rapid evolution of the duty cycle (above our $\lambda_{\mathrm{sBHAR}}>0.01$ threshold) in higher mass galaxies due to the transition from the steeper high-$\lambda_{\mathrm{sBHAR}}$ slope to the flatter low-$\lambda_{\mathrm{sBHAR}}$ slope in this regime. The average accretion rate, $\langle\lambda_{\mathrm{sBHAR}}\rangle$, on the other hand, increases in a consistent manner across all stellar masses.

Our results indicate that the accretion rate distributions are ultimately limited by the self-regulation of black hole growth, at a point that roughly corresponds to the Eddington limit ($\lambda_{\mathrm{sBHAR}}\approx 1$ in our adopted units), indicated by the grey shaded region in Figure 9. Rather than a hard limit, however, the data suggest a steep power-law relation defines this limit as rare and short-lived periods of super-Eddington accretion may still take place and, more importantly, scatter in the black hole mass-to-stellar mass scaling relation will soften a hard cut-off in our observations of $\lambda_{\mathrm{sBHAR}}$ (we do not directly measure the Eddington ratio, see Aird et al. 2013). Effectively, galaxies with over-massive black holes (compared to the baseline relation) can grow at specific accretion rates (relative to the host mass) that exceed $\lambda_{\mathrm{sBHAR}}\approx 1$ without violating the true Eddington limit of the black hole. Such specific accretion rates ($\lambda_{\mathrm{sBHAR}}>1$) would also be observed if black hole growth precedes the bulk of galaxy growth (e.g. Bennert et al., 2011; Volonteri et al., 2015; Trakhtenbrot et al., 2015).

We suggest that the self-regulation of black hole growth prevents further positive evolution of $p(\lambda_{\mathrm{sBHAR}})$ at high $\lambda_{\mathrm{sBHAR}}$ above $z\sim 1$ in star-forming galaxies, despite further increases in the gas fraction. While more gas is supplied to the central regions of the galaxy, self-regulation of the accretion process prevents $p(\lambda_{\mathrm{sBHAR}})$ from exceeding this limit. Nevertheless, the copious supply of gas may continue to push up the average accretion rates (see Figure 8), increase the minumum accretion rates (i.e. the low-$\lambda_{\mathrm{sBHAR}}$ turnover) and lead to a narrowing of the distribution close to the self-regulation limit at the highest redshifts (as illustrated in the top-right panel of Figure 9).

4.2 The distribution and evolution of AGN activity in quiescent galaxies

We now turn our attention to the AGNs within quiescent galaxies. The bottom row of Figure 9 attempts to summarize our main results and the overall trends that we identify with redshift and stellar mass, focusing on moderate-to-high stellar masses ($\log\mathcal{M}_{*}/\mathcal{M}_{\odot}\sim 10-11.5$), where our measurements of $p(\lambda_{\mathrm{sBHAR}})$ are most reliable.

As in star-forming galaxies, we measure a broad distribution of $\lambda_{\mathrm{sBHAR}}$ within quiescent galaxies across all stellar masses and redshifts that we probe, likely indicating that stochastic fuelling mechanisms are at work. At lower redshifts (e.g. $z\sim 0.5-1$), the distribution of $\lambda_{\mathrm{sBHAR}}$ in quiescent galaxies is shifted to lower values than in star-forming galaxies, possibly reflecting the much lower gas fraction in the quiescent galaxy population. There are also indications of a mild stellar mass dependence at $z<2$, whereby at the highest stellar masses ($11<\log\mathcal{M}_{*}/\mathcal{M}_{\odot}<11.5$) the duty cycle decreases (see right panel of Figure 6) and the accretion distribution may be shifted to lower $\lambda_{\mathrm{sBHAR}}$ on average (see bottom row of Figure 5 and the illustration in Figure 9). Such a shift to lower $\lambda_{\mathrm{sBHAR}}$ could be due to a decrease in the gas fraction in the highest mass quiescent galaxies (e.g. Sargent et al., 2015), potentially due to the denser environments of higher mass quiescent galaxies and an associated decrease in the external gas supply. Alternatively, the primary fuel supply for AGN activity in quiescent galaxies may be mass-loss from the stars in the galaxy itself (e.g. Norman & Scoville, 1988; Ciotti & Ostriker, 2007; Kauffmann & Heckman, 2009), which can provide a sustained supply of low angular momentum gas that could drive low-level AGN activity. The shift to lower $\lambda_{\mathrm{sBHAR}}$ at the highest stellar masses would then reflect the older stellar ages of the most massive galaxies and the decreased mass loss rate from a more evolved stellar population.

Our measurements reveal a strong evolution of $p(\lambda_{\mathrm{sBHAR}})$ within quiescent galaxies for all masses $\log\mathcal{M}_{*}/\mathcal{M}_{\odot}\gtrsim 10$ (see bottom row of Figure 4), which is reflected by the rapid increase in the AGN duty cycle (see right panel of Figure 6). As in the star-forming galaxy population, this evolution may simply be driven by the overall increase in the gas fraction; at $z\gtrsim 1$, the “quiescent” galaxy population, defined relative to the evolving star-forming main sequence, can have relatively high SFRs (compared to quiescent galaxies at lower $z$) and may have correspondingly higher gas fractions (e.g. Gobat et al., 2017). Alternatively, the evolution may reflect the younger stellar ages of some higher redshift quiescent galaxies and the associated increase in stellar mass loss. At high redshifts ($z\sim 2$), a significant fraction of massive quiescent galaxies have relatively young stellar populations, indicating they have undergone recent and rapid quenching of their star formation activity (e.g. Whitaker et al., 2012), which could explain the rapid evolution of $p(\lambda_{\mathrm{sBHAR}})$ from $z\sim 0.5$ to 2 and the high duty cycle in quiescent galaxies at $z\gtrsim 2$ (see also Wang et al., 2017).

We also suggest that AGN activity in quiescent galaxies is subject to the same self-regulation—associated with Eddington-limited growth—that restricts the evolution of $p(\lambda_{\mathrm{sBHAR}})$ in star-forming galaxies. However, quiescent galaxies appear to reach this limit at higher redshifts ($z\gtrsim 2$, compared to $z\sim 1$ for star-forming galaxies), with high $\mathcal{M}_{*}$ quiescent galaxies reaching this limit at higher $z$ than in more moderate $\mathcal{M}_{*}$ galaxies (as illustrated in Figure 9).

4.3 Comparison with prior studies of AGN accretion rate distributions

Here we compare our findings with prior studies of the distribution of AGN accretion rates (or Eddington ratios) and their relation to the galaxy population.

Early work by Kauffmann & Heckman (2009) used optical spectroscopic diagnostics to determine the distribution of Eddington ratios within nearby ($z\lesssim 0.3$) galaxies from the Sloan Digital Sky Survey (SDSS: York et al., 2000) and identified two, distinct Eddington ratio distributions: a relatively narrow, lognormal distribution in young (i.e. star-forming) galaxies; and a broader, power-law distribution in older (i.e. quiescent) galaxies. However, more recent work by Jones et al. (2016) showed that the intrinsic distribution of Eddington ratios in star-forming galaxies—once contamination of the optical emission lines due to star formation processes is accounted for—is much broader and can be described by a Schechter function (i.e. a power-law like dependence at low Eddington ratios with an exponential cutoff at the highest accretion rates). Our measurements of $p(\lambda_{\mathrm{sBHAR}})$ based on the X-ray emission thus appear consistent with lower redshift optical studies: we also find a broad distribution of $\lambda_{\mathrm{sBHAR}}$ in both star-forming and quiescent galaxies with a relatively flat, power-law like slope at lower $\lambda_{\mathrm{sBHAR}}$ and a steeper tail to the highest $\lambda_{\mathrm{sBHAR}}$. In addition, we identify a turnover or flattening of $p(\lambda_{\mathrm{sBHAR}})$ at the lowest accretion rates ($\lambda_{\mathrm{sBHAR}}\lesssim 10^{-2}-10^{-3}$) that may indicate a rapid consumption of gas in the nuclear regions by the central black hole, somewhat akin to the “feast” mode introduced by Kauffmann & Heckman (2009), although our distributions are broader than the lognormal originally identified by those authors and are in better agreement (at higher $\lambda_{\mathrm{sBHAR}}$) with the Schechter function proposed by Jones et al. (2016). Nevertheless, the differences between $p(\lambda_{\mathrm{sBHAR}})$ for star-forming and quiescent galaxies that we identify are indicative of different modes of black hole growth in the different galaxy types and are thus consistent with the general picture of a “feast” mode (associated with the rapid consumption of cold gas) dominating in star-forming galaxies and a “famine” mode (fuelled by stellar mass loss) in quiescent galaxies, as proposed by Kauffmann & Heckman (2009).

Prior work at higher redshifts have primarily used X-ray emission as the tracer of the AGN accretion rate and are thus more directly comparable with our study. A12 measured $p(\lambda_{\mathrm{sBHAR}})$ out to $z\sim 1$ based on X-ray selection of AGNs within a sample of galaxies from the PRIMUS survey, and concluded that the distribution of accretion rates could be described by a mass-independent, power-law distribution (see also Bongiorno et al., 2012). In contrast to A12, our new results presented here find evidence for a stellar-mass-dependent $p(\lambda_{\mathrm{sBHAR}})$, both within the entire galaxy population and the star-forming and quiescent galaxy samples (see Figures 2 and 3). However, it is important to note our findings are consistent (within the uncertainties) with these prior studies. In fact, the strong, observational bias identified by A12—whereby lower specific accretion rate AGNs are harder to identify in lower mass galaxies—remains a significant effect in our study. By measuring the distribution of specific accretion rates, normalizing by the stellar mass, we account for the observational bias identified by A12. The stellar mass dependencies revealed by our measurements corresponds to a second-order effect that is in addition to the primary observational effect identified by A12. Our order of magnitude increase in galaxy sample size, the increase in the range of stellar masses and redshifts probed by our study, and our advanced techniques to exploit the deep Chandra imaging have allowed us to identify this additional stellar mass dependence that was not apparent in the earlier work of A12.

Our flexible, non-parametric approach also removes the need to assume an underlying functional form (e.g. the power-law function assumed by A12; Azadi et al., 2015; Wang et al., 2017) and thus enables us to recover the shape of the underlying distribution, informed directly by the data, as a function of stellar mass and redshift. We find that $p(\lambda_{\mathrm{sBHAR}})$ generally has a double power-law shape at $\lambda_{\mathrm{sBHAR}}\gtrsim 10^{-3}$, with a sharp turnover at the highest accretion rates and evolving toward higher $\lambda_{\mathrm{sBHAR}}$ at higher $z$. The earlier observational study of Bongiorno et al. (2012) found indications of a high $\lambda_{\mathrm{sBHAR}}$ turnover in $p(\lambda_{\mathrm{sBHAR}})$, while Aird et al. (2013) suggested a double power-law shape (with a steep high-$\lambda_{\mathrm{sBHAR}}$ slope) was required to reconcile the measurements of $p(\lambda_{\mathrm{sBHAR}})$ and the galaxy stellar mass function with the observed shape of the AGN luminosity function (see also Caplar et al., 2015; Bongiorno et al., 2016). Our study confirms these findings and provides new observational constraints on the overall shape and evolution of $p(\lambda_{\mathrm{sBHAR}})$, independent of any assumed functional form. Our measurements also indicate a flattening or turnover of $p(\lambda_{\mathrm{sBHAR}})$ at the lowest accretion rates ($\lambda_{\mathrm{sBHAR}}\lesssim 10^{-2}-10^{-3}$), as discussed in Section 4.1 above.

By dividing our galaxy sample into star-forming and quiescent galaxies (based on their position in the SFR-$\mathcal{M}_{*}$ plane) and measuring $p(\lambda_{\mathrm{sBHAR}})$ within each population as a function of stellar mass and redshift, we are able to identify potential differences in the triggering and fuelling mechanisms of AGNs across the galaxy population. Earlier studies also found evidence that the normalization of $p(\lambda_{\mathrm{sBHAR}})$ is higher in galaxies with blue colours (e.g. A12) or with higher SFRs (e.g. Azadi et al., 2015). Georgakakis et al. (2014) also showed that the absolute space density as a function of $\lambda_{\mathrm{sBHAR}}$ is higher for AGNs with star-forming host galaxies than for those with quiescent host galaxies (separated based on the rest-frame UVJ colour criteria, e.g. Williams et al. 2009), and found initial indications of differences in the shapes of the $\lambda_{\mathrm{sBHAR}}$ distributions. More recently, Wang et al. (2017) measured $p(\lambda_{\mathrm{sBHAR}})$ directly in samples of red, green and blue galaxies (separated using dust-corrected rest-frame colours) in the GOODS-North and -South fields. Consistent with our study, they found that the normalization of $p(\lambda_{\mathrm{sBHAR}})$ is higher in blue (i.e. star-forming) galaxies than in red (i.e. quiescent) galaxies, has a different slope (fitting with a simple power-law function), and may evolve differently with redshift. In particular, Wang et al. (2017) identified a strong evolution of $p(\lambda_{\mathrm{sBHAR}})$ with redshift in red galaxies and a correspondingly high duty cycle in high-redshift ($z\sim 2$) quiescent galaxies that is confirmed by our measurements. However, by adopting a larger sample of galaxies—covering a broader range in redshift and stellar mass—and exploiting new techniques to push the limits of the Chandra survey data, we are able to identify a stellar mass dependence (that was not found by Wang et al., 2017) and place tighter constraints on the shape and evolution of $p(\lambda_{\mathrm{sBHAR}})$ within both star-forming and quiescent galaxies.

Bongiorno et al. (2016) and Georgakakis et al. (2017) both use alternative approaches to probe the distribution of $\lambda_{\mathrm{sBHAR}}$ as a function of stellar mass and redshift. We start with a sample of galaxies and directly measure the AGN content of subsets of these galaxies using deep X-ray imaging (see also A12; Azadi et al., 2015; Wang et al., 2017), whereas Bongiorno et al. (2016) and Georgakakis et al. (2017) start from a sample of X-ray selected AGNs and measure the stellar masses of the host galaxies of this sample. Bongiorno et al. (2016) take the X-ray AGN sample from the XMM-COSMOS survey (Hasinger et al., 2007; Brusa et al., 2010) and attempt to find parameterizations for the distribution of $\lambda_{\mathrm{sBHAR}}$ and the host galaxy stellar mass function that, after accounting for the X-ray selection effects, are consistent with the observed X-ray AGN sample. Direct comparisons with our measurements are difficult as this method does not constrain the absolute normalization of $p(\lambda_{\mathrm{sBHAR}})$. Nonetheless, the basic patterns identified by Bongiorno et al. (2016)—a double power-law shaped distribution of $\lambda_{\mathrm{sBHAR}}$that shifts to higher $\lambda_{\mathrm{sBHAR}}$ at both higher redshifts and lower stellar masses—are consistent with our measurements (see e.g. Figure 8), although the parametrization adopted by Bongiorno et al. (2016) lacks the low-$\lambda_{\mathrm{sBHAR}}$ turnover identified using our data-driven method.

Georgakakis et al. (2017) adopt a larger sample of X-ray AGNs (combining data from a wide range of Chandra and XMM-Newton surveys) and use a non-parametric approach to find a $p(\lambda_{\mathrm{sBHAR}})$ that, when combined with independent measurements of the overall galaxy stellar mass function (taken from Ilbert et al., 2013), are consistent with the observed redshifts, luminosities, and host stellar masses of the X-ray selected sample. Despite the substantial differences between these studies,444As well as a distinct, independently developed methodology, Georgakakis et al. 2017 use data from a number of wide XMM-Newton surveys (in addition to the Chandra surveys used here), different multiwavelength imaging, and independent measurements of host galaxy stellar masses. the measurements of $p(\lambda_{\mathrm{sBHAR}})$ are in remarkably good agreement (see appendix B of Georgakakis et al. 2017 for a direct comparison). Both of our studies find a broad distribution of $\lambda_{\mathrm{sBHAR}}$ with a sharp cutoff at the highest accretion rates ($\lambda_{\mathrm{sBHAR}}\gtrsim 0.1-1$) and indications of a flattening or turnover at very low accretion rates ($\lambda_{\mathrm{sBHAR}}\lesssim 10^{-3}-10^{-2}$). As in Georgakakis et al. (2017), we identify a significant increase in the AGN duty with both increasing redshift and stellar mass. While Georgakakis et al. (2017) attribute the redshift evolution to an increase in the normalization of $p(\lambda_{\mathrm{sBHAR}})$, our measurements reveal a shift of the distributions toward higher $\lambda_{\mathrm{sBHAR}}$ at higher $z$ (although some level of evolution in normalization cannot be ruled out). Our approach, extracting X-ray data from all galaxies in a sample to probe below the nominal sensitivity limits of Chandra, allows us to place improved constraints at lower $\lambda_{\mathrm{sBHAR}}$ and across a wider range of galaxy stellar mass. In addition, we measure $p(\lambda_{\mathrm{sBHAR}})$ independently as a function of stellar mass and redshift in star-forming and quiescent galaxies. Separately probing $p(\lambda_{\mathrm{sBHAR}})$ in star-forming and quiescent galaxies, as in our study, provides crucial insights into the different AGN fuelling mechanisms that appear to drive black hole growth in these different galaxy types.

4.4 Average AGN accretion rates and the AGN main sequence

A number of recent studies have used X-ray stacking techniques to measure average AGN accretion rates for samples of galaxies as a function of redshift, stellar mass or SFR (e.g. Mullaney et al., 2012; Chen et al., 2013; Delvecchio et al., 2015; Yang et al., 2017). As discussed by Hickox et al. (2014), stacking can account for the long-term variability in AGN accretion rates and reveal relationships between the average black hole growth rates and the relatively slowly varying properties of galaxies. In contrast, our measurements of $p(\lambda_{\mathrm{sBHAR}})$ track the distribution of AGN accretion rates, directly tracing the variability of AGN activity and revealing how these overall distributions depend on galaxy properties. X-ray stacking is equivalent to averaging over these distributions, losing much of the detailed information (such as the fraction of galaxies with black holes growing above a given accretion rate i.e. the duty cycle). Here, we briefly discuss how our findings relate to these prior studies.

Mullaney et al. (2012, hereafter M12) stacked the X-ray emission from samples of star-forming galaxies at $z\sim 1$ and $z\sim 2$ and identified a relation between the average X-ray luminosity and stellar mass with slope of $\sim$unity, designating this relation as a hidden “AGN main sequence”.

As noted by M12, this relation is primarily due to higher stellar mass galaxies that are accreting at the same specific accretion rates producing higher observed X-ray luminosities – effectively the stellar mass selection bias identified by A12 and accounted for in this paper by measuring the distribution of $\lambda_{\mathrm{sBHAR}}$ directly.

In Figure 10 we estimate the mean X-ray luminosity of star-forming galaxy samples, $\langle L_{\mathrm{X}}\rangle$, as a function of stellar mass and redshift using our measurements of $p(\lambda_{\mathrm{sBHAR}})$, integrating over the full range of $\lambda_{\mathrm{sBHAR}}$ and converting back to $L_{\mathrm{X}}$ using Equation 2. Our results (shaded regions) are consistent with a rising $\langle L_{\mathrm{X}}\rangle$ with $\mathcal{M}_{*}$ i.e. an “AGN main sequence”. Our results agree with the M12 measurements at $z\sim 1$, although our estimates lie below the extrapolation of a linear fit to the M12 measurements (dashed blue line) at lower stellar masses (see also Yang et al., 2017). This suppression of AGN activity at lower stellar masses is related to the reduction in the AGN duty cycle that we identify above (see Figure 6), although the strength of the suppression of $\langle L_{\mathrm{X}}\rangle$ is somewhat counteracted by the increase in the average specific accretion rate, $\langle\lambda_{\mathrm{sBHAR}}\ \rangle$, of those AGNs (see Figure 8).

At $z\sim 2$ our estimates of $\langle L_{\mathrm{X}}\rangle$ in Figure 10 lie systematically below the M12 measurements of $\langle L_{\mathrm{X}}\rangle$, most likely due to differences in the corrections for the contribution from star formation to the X-ray luminosity (primarily from high-mass X-ray binaries). M12 apply corrections based on locally-calibrated scaling relations (Ranalli et al., 2003), whereas we apply a larger correction based on our work in Paper I that indicates a larger X-ray luminosity per unit SFR at these redshifts (see also Lehmer et al., 2016). Our measurements still reveal a clear “AGN main sequence” at $z\sim 2$, although there are indications that the slope may be steeper than unity due to a suppression of AGN activity at lower stellar masses.

Georgakakis et al. (2017) also present estimates of $\langle L_{\mathrm{X}}\rangle$ based on their independent measurements of $p(\lambda_{\mathrm{sBHAR}})$. In contrast to our work, Georgakakis et al. (2017) find a significantly flatter slope in the relation between $\log\langle L_{\mathrm{X}}\rangle$ and $\log\mathcal{M}_{*}$, associated with a suppression in AGN activity at high stellar masses (cf. our suppression at low stellar mass seen in Figure 10 and steeper slope). However, Georgakakis et al. (2017) do not separate the star-forming and quiescent galaxy populations and thus this flatter relation is likely associated with the lower accretion rates of AGN in quiescent galaxies (that start to dominate galaxy samples at $\log\mathcal{M}_{*}/\mathcal{M}_{\odot}\gtrsim 10.5$ and $z\lesssim 1$), as revealed in our direct measurements of $p(\lambda_{\mathrm{sBHAR}})$ in the separated star-forming and quiescent galaxy populations (e.g. Figure 3).

Other studies have measured the stacked X-ray emission from samples of galaxies as a function of SFR (rather than stellar mass), finding a positive correlation that could indicate a global connection between the build of galaxies and the growth of their central black holes (e.g. Chen et al. 2013; see also Hickox et al. 2014). However, whether this correlation is due to the underlying relationship between stellar mass and SFR in galaxies (i.e. the galaxy main sequence) and thus whether stellar mass or SFR are the primary driver of AGN activity remains an area of active study (e.g. Vito et al., 2014; Delvecchio et al., 2015; Azadi et al., 2015; Yang et al., 2017). Such work is further complicated by the presence of AGNs in quiescent galaxies (with very low SFRs) that may be fuelled by different physical mechanisms, as suggested in this paper (see also Kauffmann & Heckman, 2009). Detailed measurements of $p(\lambda_{\mathrm{sBHAR}})$ as a function of both stellar mass and SFR are deferred to Paper III.

4.5 AGN activity in dwarf galaxies

At the lower redshifts probed by our study ($z<1.5$), our galaxy sample covers very low stellar masses corresponding to the dwarf galaxy regime ($\mathcal{M}_{*}/\mathcal{M}_{\odot}\lesssim 3\times 10^{9}$). The existence, masses and growth rates of massive black holes at the centres of such low-mass galaxies have been the focus of extensive study and interest (e.g. Greene & Ho, 2007; Reines et al., 2013; Mezcua et al., 2016; Pardo et al., 2016; Baldassare et al., 2016, 2017; Hood et al., 2017; Chen et al., 2017). Such galaxies may provide the conditions under which the seeds of supermassive black holes can form and thus may provide local analogues that can help shed light on how supermassive black holes formed and subsequently grew in the very early ($z\gtrsim 6$) universe (e.g. Volonteri, 2010; Reines & Comastri, 2016).

Our measurements place important constraints on both the black hole occupation fraction (the fraction of galaxies with a central, massive black hole) and the extent of AGN activity within dwarf galaxies. Figure 6 indicates an AGN duty cycle of $f(\lambda_{\mathrm{sBHAR}}>0.01)\approx 0.1-0.9$ per cent in galaxies with $8.5<\log\mathcal{M}_{*}/\mathcal{M}_{\odot}<9.5$, which is in agreement with prior studies of X-ray AGNs in dwarf galaxy samples (e.g. Pardo et al. 2016 find an AGN fraction of $\sim$0.6 – 3 per cent in galaxies with $\mathcal{M}_{*}\leq 3\times 10^{9}\mathcal{M}_{\odot}$ to a luminosity threshold of $L_{\mathrm{X}}>10^{41}$ $\mathrm{erg\;s^{-1}}$, comparable to our $\lambda_{\mathrm{sBHAR}}$ limit). Our measurements are also consistent with no evolution with redshift in the AGN duty cycle for dwarf galaxies out to at least $z\sim 1$, in contrast to higher mass galaxies where there is significant evolution. Despite this low overall AGN duty cycle, we do identify a population of high accretion rate ($\lambda_{\mathrm{sBHAR}}\gtrsim 0.1$) AGNs within our dwarf galaxy samples, as indicated by our measurements of $p(\lambda_{\mathrm{sBHAR}})$ (see e.g. Figure 2). Thus, a small fraction of dwarf galaxies ($\sim 0.1$ per cent) must contain massive black holes that are growing at high accretion rates relative to their stellar mass. Indeed, our measurements of $p(\lambda_{\mathrm{sBHAR}})$ at high $\lambda_{\mathrm{sBHAR}}$ are consistent across the full range of galaxy masses probed in this study, with significant differences only becoming apparent at lower $\lambda_{\mathrm{sBHAR}}$ (see e.g. the two lowest redshift panels of Figure 5 for all galaxies or star-forming galaxies).

It is also important to note, especially for this low-mass regime, that our measurements place constraints on the distribution of specific accretion rates: the levels of black hole growth relative to the total galaxy stellar mass. The scaling between the black hole mass and the host galaxy mass may have significantly larger scatter or break down completely in these very low-mass galaxies (e.g. Greene et al., 2008, 2010; Kormendy & Ho, 2013; Graham & Scott, 2015; Reines & Volonteri, 2015). Thus, our measurements of $p(\lambda_{\mathrm{sBHAR}})$ may not reflect the underlying distribution of Eddington ratios; in fact, our measurements track the convolution of the probability distribution of Eddington ratios and the distribution of black hole masses at a given galaxy stellar mass. Our measurements are also degenerate with the overall black hole occupation fraction: the low AGN duty cycle that we find could indicate that a large fraction of galaxies in this mass regime simply lack central massive black holes and thus would never be observed as an AGN (see also Trump et al., 2015). It is not possible with the current data to break such degeneracies, although our measurements do provide a lower limit on the occupation fraction and show that actively accreting massive black holes can exist in such low-mass galaxies. Thus, despite these limitations, our measurements of $p(\lambda_{\mathrm{sBHAR}})$ provide important observational constraints on the incidence of AGNs and the extent of black hole growth within the dwarf galaxy population.

4.6 AGN accretion rate distributions in cosmological simulations

Recent large-volume cosmological hydrodynamic simulations (e.g. Dubois et al., 2014; Vogelsberger et al., 2014; Schaye et al., 2015) are able to model the flow of gas into the very central regions of galaxies and—via sub-grid modelling—predict the accretion rates onto the central black holes. Here we briefly discuss how the results of recent cosmological simulations compare with our measurements of $p(\lambda_{\mathrm{sBHAR}})$ and AGN duty cycles.

Volonteri et al. (2016) presented predictions for both Eddington ratio distribution functions and distributions of $\lambda_{\mathrm{sBHAR}}$ (accretion rate scaled relative to the total predicted galaxy stellar mass) based on the Horizon-AGN simulation. The predictions are in good agreement with our observational results for galaxies with moderate stellar masses, $\mathcal{M}_{*}\gtrsim 10^{10}\mathcal{M}_{\odot}$, revealing a broad, roughly power-law distribution at moderate accretion rates ($\lambda_{\mathrm{sBHAR}}\sim 10^{-3}-0.1$), a break at higher accretion rates ($\lambda_{\mathrm{sBHAR}}\gtrsim 0.1-1$), and a shift toward higher $\lambda_{\mathrm{sBHAR}}$ at higher redshifts. The simulation also predicts the low $\lambda_{\mathrm{sBHAR}}$ flattening or turnover that is seen in our measurements. At lower stellar masses, the simulation tends to over-predict $p(\lambda_{\mathrm{sBHAR}})$, possibly indicating the need for stronger supernova-driven feedback in the simulations that would limit the fuelling AGN activity in lower mass galaxies (see discussion in Volonteri et al. 2016 and Dubois et al. 2015).

Rosas-Guevara et al. (2016) presented predictions for the observed properties of AGNs in the EAGLE simulation (see also McAlpine et al., 2017). They find that the average Eddington ratio for the simulated galaxy sample increases with redshift, consistent with our observed trends for $\langle\lambda_{\mathrm{sBHAR}}\rangle$ (see Figure 8). The average Eddington ratio is found to decrease at higher black hole masses at a fixed redshift, which appears broadly consistent with our measurements of $\langle\lambda_{\mathrm{sBHAR}}\rangle$ seen in Figure 8 (left panel) but may be in conflict with the increase in the overall duty cycle at higher $\mathcal{M}_{*}$ that we observe (see Figure 6). Further, detailed comparisons of EAGLE predictions of $p(\lambda_{\mathrm{sBHAR}})$ that can be directly compared to our measurements are required to reconcile these issues and assess the modelling of AGN accretion in EAGLE.

Sijacki et al. (2015) presented predicted Eddington ratio distributions as a function of redshift based on the Illustris simulation. Consistent with our observations of $p(\lambda_{\mathrm{sBHAR}})$, the predicted distributions exhibit a broad range of accretion rates at any epoch. The distributions shift toward higher accretion rates at higher redshifts and narrow as they reach the Eddington limit, as seen in our measurements. Sijacki et al. (2015) also predict that average Eddington ratios are roughly constant over a broad range of black hole mass (at a given redshift) but drop very rapidly for the very highest mass black holes ($\mathcal{M}_{\mathrm{BH}}\gtrsim 10^{9}-10^{10}\mathcal{M}_{\odot}$, depending on redshift). This mass dependence differs from the much more gradual decrease in average Eddington ratio with increasing black hole mass found by Rosas-Guevara et al. (2016) in the EAGLE simulation and the changes in the duty cycle and $\langle\lambda_{\mathrm{sBHAR}}\rangle$ with stellar mass that we measure. Direct comparisons of the Illustis predictions of $p(\lambda_{\mathrm{sBHAR}})$ as a function of $\mathcal{M}_{*}$ are required to explore this issue.

In conclusion, our measurements of $p(\lambda_{\mathrm{sBHAR}})$ are generally consistent with the broad range of accretion rates predicted by recent cosmological simulations and the overall evolution of these distributions, shifting to higher accretion rates at higher redshifts. However, none of the simulations appear to produce the exact stellar mass dependence of $p(\lambda_{\mathrm{sBHAR}})$ that we measure. Directly comparing our measurements with predications of $\lambda_{\mathrm{sBHAR}}$ (i.e. relative to host stellar mass) as a function of galaxy type, stellar mass and redshift will enable crucial tests of new and ongoing cosmological simulations, providing a more refined observational comparison than generic characterizations of the AGN population (e.g. luminosity functions).

5 Conclusions

In this paper, we use deep Chandra X-ray data to probe the distribution of AGN accretion rates within large samples of stellar-mass-selected star-forming and quiescent galaxies out to $z\sim 4$. Our main conclusions are as follows:

Our measurements show that $p(\log\lambda_{\mathrm{sBHAR}}\;|\;\mathcal{M}_{*},z)$ has a broad distribution that is roughly consistent with a power-law with a steep cutoff at the highest accretion rates ($\lambda_{\mathrm{sBHAR}}\gtrsim 0.1-1$) and a flattening or possibly a turnover at low accretion rates ($\lambda_{\mathrm{sBHAR}}\lesssim 10^{-3}-10^{-2}$). These broad distributions are consistent with stochastic processes being primarily responsible for the fuelling of AGNs and substantial variability in the level of black hole growth over short timescales relative to the host galaxy. 2. 2.

This broad distribution of accretion rates is found for both the star-forming and quiescent galaxy populations, across a wide range in stellar mass and out to $z\sim 4$. At higher stellar masses ($\mathcal{M}_{*}\gtrsim 10^{10}\mathcal{M}_{\odot}$) and lower redshifts ($z\lesssim 2$) we are able to identify significant differences in $p(\log\lambda_{\mathrm{sBHAR}}\;|\;\mathcal{M}_{*},z)$ in star-forming versus quiescent galaxies (at fixed $\mathcal{M}_{*}$ and $z$), whereby the distribution in quiescent galaxies appears shifted to lower $\lambda_{\mathrm{sBHAR}}$. Thus, quiescent galaxies tend to host weaker (lower $\lambda_{\mathrm{sBHAR}}$) AGNs and the AGN duty cycle (the fraction of galaxies with AGN above a fixed limit of $\lambda_{\mathrm{sBHAR}}>0.01$) is generally lower. 3. 3.

We find that the $p(\log\lambda_{\mathrm{sBHAR}}\;|\;\mathcal{M}_{*},z)$ evolves strongly with redshift for star-forming galaxies with $\mathcal{M}_{*}\gtrsim 10^{10}\mathcal{M}_{\odot}$, shifting toward higher $\lambda_{\mathrm{sBHAR}}$ at higher redshifts and leading to a significant increase in the AGN duty cycle. This evolution may be due to the increased availability of cold gas at higher redshifts that increases both the frequency and luminosities of accretion events at earlier cosmic times. However, the distributions appear to be truncated at a point roughly corresponding to the Eddington limit, indicating that black holes may ultimately self-regulate their own growth at high redshifts when copious gas is available. 4. 4.

Our results indicate that $p(\log\lambda_{\mathrm{sBHAR}}\;|\;\mathcal{M}_{*},z)$ in star-forming galaxies depends on stellar mass. The probability of a lower mass star-forming galaxy hosting a low-accretion-rate AGN ($\lambda_{\mathrm{sBHAR}}\lesssim 0.1$) is suppressed compared to within higher mass galaxies. The evolution of the AGN duty cycle with redshift is thus weaker at lower stellar masses. We measure a low duty cycle ($\lesssim 0.9$ per cent) in the lowest mass, dwarf galaxy regime ($\mathcal{M}_{*}\lesssim 3\times 10^{9}\mathcal{M}_{\odot}$) and find no evidence for evolution out to $z\sim 1$. Nonetheless, high accretion rate AGNs ($\lambda_{\mathrm{sBHAR}}\gtrsim 0.1$) are identified in a small fraction ($\sim 0.1$ per cent) of dwarf galaxies. 5. 5.

In quiescent galaxies, we also find strong evolution whereby $p(\log\lambda_{\mathrm{sBHAR}}\;|\;\mathcal{M}_{*},z)$ shifts to higher $\lambda_{\mathrm{sBHAR}}$ at higher $z$ for $\mathcal{M}_{*}\gtrsim 10^{10}\mathcal{M}_{\odot}$. The duty cycle appears to decrease at the highest stellar masses, $\mathcal{M}_{*}\gtrsim 10^{11}\mathcal{M}_{\odot}$ (in contrast to the increase with stellar mass found for star-forming galaxies). We also find that the duty cycle for quiescent galaxies evolves strongly with redshift and is comparable to star-forming galaxies at $z\sim 2$. The rapid evolution of the duty cycle may also be related to the increased availability of cold gas at higher redshifts or could be due to younger (recently quenched) stellar populations with increased stellar mass loss.

Our new measurements have enabled a detailed mapping of the distribution of black hole growth over the galaxy population. While our results support a picture of stochastically driven AGN activity, our measurements show that such activity is enhanced at higher redshifts, in high-mass star-forming galaxies, and in high-redshift quiescent galaxies.

These findings indicate that black hole growth (via AGN activity) does not simply track galaxy growth (via star formation) and that the triggering and fuelling of AGNs may be driven by different physical mechanisms across the galaxy population (e.g. in star-forming versus quiescent galaxies). Our detailed measurements place strong constraints on the extent of these different AGN fuelling mechanisms over a wide range of stellar mass and redshift.

The measurements of accretion rate probability distributions, $p(\log\lambda_{\mathrm{sBHAR}}\;|\;\mathcal{M}_{*},z)$, and the AGN duty cycle, $f(\lambda_{\mathrm{sBHAR}}>0.01)$, presented in this paper are made available in machine-readable format at http://doi.org/10.5281/zenodo.1009605.

acknowledgements

We thank the anonymous referee for helpful comments that contributed to this paper. We acknowledge helpful conversations with Dave Alexander, Francesca Civano, Andy Fabian, Ryan Hickox, James Mullaney, Kirpal Nandra and Marta Volonteri. JA acknowledges support from ERC Advanced Grant FEEDBACK 340442. AG acknowledges the thales project 383549 that is jointly funded by the European Union and the Greek Government in the framework of the programme “Education and lifelong learning”. This work is based in part on observations taken by the CANDELS Multi-Cycle Treasury Program and the 3D-HST Treasury Program (GO 12177 and 12328) with the NASA/ESA HST, which is operated by the Association of Universities for Research in Astronomy, Inc., under NASA contract NAS5-26555. The scientific results reported in this article are based to a significant degree on observations made by the Chandra X-ray Observatory.

Appendix A Two-component (galaxy+AGN) fitting of spectral energy distributions and AGN-corrected estimates of star formation rates

In appendix A of Paper I we describe our approach to estimate the SFRs of galaxies, based on either summing the UV+IR emission or from fitting the UV-to-MIR spectral energy distributions (SEDs) with stellar population synthesis (SPS) models using a modified version of the FAST code (Kriek et al., 2009). The SED fitting is also used to estimate the total galaxy stellar mass, $\mathcal{M}_{*}$. However, for sources with an AGN, the host galaxy light can be severely contaminated by the AGN emission, especially at UV and MIR wavelengths. In this appendix, we describe how we correct for potential AGN contamination by including AGN templates in our SED fitting and SFR estimates, implementing a two-component (galaxy+AGN) fitting approach (see also e.g. Bongiorno et al., 2007, 2012; Lusso et al., 2012; Rovilos et al., 2014; Ciesla et al., 2015).

We allow for a potential AGN contribution for any galaxy with a significant X-ray detection (i.e. detected with a false probability $<4\times 10^{-6}$ in the full, soft, hard or ultra-hard Chandra energy bands, see Section 2 above and Laird et al., 2009). Allowing for a potential AGN contribution for every galaxy in our overall sample, regardless of X-ray properties, is computationally prohibitive and risks systematically altering our SFR and $\mathcal{M}_{*}$ estimates when we know that the majority of galaxies are not expected to have a significant AGN contribution. There will be AGN that are missed by our X-ray selection (in particular, luminous but obscured sources that may be revealed by their MIR colours but are not seen at X-ray wavelengths e.g. Donley et al., 2012; Mendez et al., 2013). However, such sources are expected to constitute a small fraction of our overall galaxy sample and thus, while the SFR and $\mathcal{M}_{*}$ esitimates may be inaccurate for individual sources, they will not have a significant impact on our overall results.

We adopt a library of eight, empirically determined templates for the SEDs of AGNs (that exclude any host galaxy contribution) by combining the five AGN-dominated templates from the Polletta et al. (2007) SWIRE template library (namely, the Torus, TQSO1, BQSO1, QSO1 and QSO2 templates) with the three composite SEDs of X-ray selected AGNs with absorption column densities $\mathrm{N_{H}}=10^{22-23},10^{23-24}$ and $10^{24-25}$ cm*-2* from Silva et al. (2004). These templates are chosen to provide a reasonable sampling of the range of possible AGN SEDs. The templates from Silva et al. (2004) are adopted as they correspond to optically obscured AGN with very red, power-law emission that can dominate at IR wavelengths (with negligible UV-to-optical emission). We do not attempt to estimate $\mathrm{N_{H}}$ for our X-ray selected AGN sample and thus do not apply any prior preference for a particular template based on the X-ray absorption or any other data when fitting the SED of a given source. Our AGN template library is shown in Figure 11.

We allow for linear combinations of each of our eight AGN templates with each of our grid of galaxy templates (described in Paper I), such that

[TABLE]

where $f_{\mathrm{gal}}(\lambda)$ and $f_{\mathrm{AGN}}(\lambda)$ represent a given galaxy and AGN SED template, respectively, and $A_{1}$ and $A_{2}$ are unknown, free parameters corresponding to the relative scalings of the two templates. For each possible combination (i.e. for every galaxy template combined with each of the eight AGN templates), we calculate $A_{1}$ and $A_{2}$ using $\chi^{2}$ minimization to fit to the observed photometry. We also determine the $\chi^{2}$ for every galaxy template without any AGN component (i.e. fixing $A_{2}=0$). We retain the standard FAST template error function when calculating $\chi^{2}$, although we truncate the uncertainty to 50 per cent (in flux) at long wavelengths (where the standard template error function becomes large, mostly to allow for uncertainties related to an AGN contribution that we now model directly). We idenitfy the minimum $\chi^{2}$ over the entire grid of galaxy and AGN combinations, which we retain as our best-fitting estimate of the overall SED. Following our approach in Paper I, we also calculate a posteriori estimates of galaxy properties (i.e. $\mathrm{SFR_{SED}}$, $\mathcal{M}_{*}$) by marginalizing over the full grid of possible galaxy and AGN combinations (as well as the “no-AGN” possibility), although we only take the single best $\chi^{2}$ (and corresponding $A_{1}$ and $A_{2}$) for a given combination rather than performing a fully Bayesian marginalization. We use the same approach to provide a best a posteriori estimate of the fractional AGN contribution, $F_{AGN,5000}$, defined as the fraction of the light at 5000Å that is attributed to the AGN template component. We flag sources where $F_{AGN,5000}$ is greater then 50 per cent as AGN-dominated in the UV-optical.

The next step of our analysis is to determine our best estimate of the SFR for each galaxy in the X-ray detected sample. In Paper I, we describe our “SFR ladder” for galaxies: when a galaxy is detected at 24µm we adopt $\mathrm{SFR_{UV+IR}}$; for galaxies without 24µm detections we adopt $\mathrm{SFR_{SED}}$. Here, we adopt a similar approach but apply a correction for any AGN contribution. First, we extrapolate the best-fitting (minimum $\chi^{2}$) AGN template to the infrared and estimate the expected flux in the MIPS 24µm band due to the AGN component. We subtract this expected AGN contribution from any observed flux at 24µm. If the source is still detected at $>3\sigma$ after subtracting the AGN contribution, we calculate $\mathrm{SFR_{UV+IR}}$ using equation A1 of Paper I but subtracting the expected AGN contribution to both the 24µm flux and $L_{\nu 2800}$ based on our two-component SED fitting. For sources where all of the 24µm flux is associated with the AGN component, or sources without 24µm detections, we adopt $\mathrm{SFR_{SED}}$ as the best estimate of the SFR. Figure 12 shows several examples of our two-component SED fits for a range of sources with differing levels of AGN contribution. In panels i and iii all of the observed 24µm flux can be accounted for with the AGN template and thus we estimate the SFR based on the best-fitting SPS model for the galaxy component (the SFR estimate will be highly uncertain in the case of panel i, where the observed SED is dominated by the AGN component). In panels ii and iv the extrapolation of the AGN component is not suffiicient to account for the 24µm flux (the galaxy models do not include the far-IR emission from star-formation-heated dust that is needed to account for this emission). In such cases we estimate the SFR from the UV+IR emission (via equation A1 of Paper I), subtracting any AGN contribution.

In Figure 13 we compare our estimates of $\mathcal{M}_{*}$ and SFR based on our two-component fitting to the previous (galaxy only, single-component) approach. The stellar masses (top panel) generally change by $\lesssim 0.3$ dex for galaxy-dominated sources ($F_{\mathrm{AGN,5000}}<50$ per cent, shown in orange) and there is no overall systematic offset. For AGN-dominated sources (blue points, $\sim$ 20 per cent of the X-ray selected sample) there is a larger scatter, although our two-component fitting should provide an improved estimate of the stellar mass for such sources. For the SFRs (bottom panel) there is a clearer systematic shift for both AGN- and galaxy-dominated sources that reduces our estimates by $\sim 0.1-0.3$ dex on average and by up to an order of magnitude in extreme cases, demonstrating the importance of applying corrections for the AGN contribution to the UV, optical and IR light.

An updated version of the FAST code, implementing the changes described here and in appendix A of Paper I, is made available at https://github.com/jamesaird/FAST.

Appendix B Bayesian mixture modelling of the distribution of AGN accretion rates and corrections for the contribution of the host galaxy

In this paper we adopt the flexible Bayesian mixture modelling approach—described in detail in appendix B of Paper I—to estimate $p(\log\lambda_{\mathrm{sBHAR}}\;|\;\mathcal{M}_{*},z)$, the intrinsic distribution of $\lambda_{\mathrm{sBHAR}}$ for galaxies of a given stellar mass, redshift, and galaxy classification (e.g. star-forming or quiescent). Here, we describe two further refinements of the methodology for the present study:

1. adapting the method to recover the distribution of specific black hole accretion rates, $\lambda_{\mathrm{sBHAR}}$, rather than X-ray luminosities; and
2. accounting for (and thus removing) the expected contribution to the observed X-ray flux from galactic (non-AGN) processes.

The first amendment—converting from X-ray luminosity to $\lambda_{\mathrm{sBHAR}}$—simply requires us to replace all $L_{\mathrm{X}}$ terms in appendix B of Paper I with $\lambda_{\mathrm{sBHAR}}$ and adjust the scaling factor for each galaxy, $k_{i}$, that converts this quantity to a predicted number of X-ray counts. Thus, following Equation 2 above that defines $\lambda_{\mathrm{sBHAR}}$, equation B3 of Paper I is adapted to

[TABLE]

where $\kappa=1.04\times 10^{34}$ is a constant factor defined by our assumed (constant) bolometric correction and our scaling of $\lambda_{\mathrm{sBHAR}}$ to “Eddington-ratio equivalent” units. The scaling for an individual galaxy, $i$, also depends on the conversion factor, $\eta(z_{i})$,555The $\eta(z_{i})$ factor converts between rest-frame 2–10 keV luminosity and the observed count rate in the 2–7 keV Chandra band. We assume an unabsorbed spectrum with photon index $\Gamma=1.9$ for every source. Thus, the conversion factor depends only on the redshift of the source, $z_{i}$. the effective X-ray exposure at the galaxy position, $t_{i}$, and the total stellar mass of the galaxy, $\mathcal{M}_{i}$. We model $p(\log\lambda_{\mathrm{sBHAR}}\;|\;\mathcal{M}_{*},z)$ as a sum of Gamma distributions, as described by equations B7-B10 of Paper I (substituting $L_{\mathrm{X}}$ with $\lambda_{\mathrm{sBHAR}}$). Similar to Paper I, we fix the positions of the Gamma distributions in $\lambda_{\mathrm{sBHAR}}$, adopting a logarithmically spaced grid over the range $-5\leq\log\lambda_{j}\leq 2$ in steps of 0.2 dex, where $\lambda_{j}$ is the scale parameter of each mixture component. We also fix the shape parameter to $\alpha_{j}=3.0$ for every component and apply a prior that prefers a smooth variation between adjacent components. To ensure sufficient flexibility, we allow for an additional component at $\log\lambda_{j}=-7$ that accounts for galaxies that do not contribute to the observed distribution of $\lambda_{\mathrm{sBHAR}}$ and is not subjected to our smoothness prior. The integral over all mixture components is required to sum to 1. Our overall model for $p(\log\lambda_{\mathrm{sBHAR}}\;|\;\mathcal{M}_{*},z)$ for a sample of galaxies in a given range of $\mathcal{M}_{*}$ and $z$ thus has 36 free parameters (linked via the smoothness prior), corresponding to the normalizations, $A_{j}$, of each mixture component.

Our second refinement is to apply a correction for the contribution from galactic (non-AGN) X-ray processes. In Paper I, we found that the observed distributions of $L_{\mathrm{X}}$ for galaxies exhibit a narrow peak at $L_{\mathrm{X}}\lesssim 10^{42}$ $\mathrm{erg\;s^{-1}}$, associated with the emission from (predominantly) high- and low-mass X-ray binaries. We used this peak to trace the star-forming “X-ray main sequence” and re-calibrate the scaling between X-ray luminosity and galaxy properties (SFR, $\mathcal{M}_{*}$). In the present paper, we aim to measure the distribution of AGN activity, as traced by the X-ray emission, and thus need to remove this non-AGN contribution. We adjust the likelihood function from Paper I for the data from a single galaxy to allow for three possible origins of the observed X-ray counts, each described by an independent Poisson process: 1) emission from the AGN, 2) emission from other processes in the galaxy, and 3) a background component. Thus, equation B4 of Paper I can be re-written as

[TABLE]

where $C_{\mathrm{A}}$, $C_{\mathrm{G}}$ and $C_{\mathrm{B}}$ are unknown, nuisance parameters that represent the integer X-ray counts due to the AGN, galactic processes, and the background, respectively. The summations in Equation 7 are equivalent to marginalizing over the possible values of the nuisance parameters: $C_{\mathrm{A}}$, $C_{\mathrm{G}}$ and $C_{\mathrm{B}}$. We require

[TABLE]

so that these three unknown parameters sum to the total observed counts, $N_{i}$, from galaxy $i$. The $k_{i}\lambda_{\mathrm{sBHAR}}$, $l_{i}L_{\mathrm{G}}$, and $b_{i}$ terms correspond to the (non-integer) underlying count rates related to the AGN, galactic, and background components, respectively, which produce the integer counts from each component ($C_{\mathrm{A}}$, $C_{\mathrm{G}}$ and $C_{\mathrm{B}}$) via a Poisson process. The background count rate, $b_{i}$, is well-determined and is taken from our background maps (see Paper I; Georgakakis et al., 2008). The likelihood function of Equation 7 describes our knowlege of the AGN count rate ($k_{i}\lambda_{\mathrm{sBHAR}}$) for galaxy $i$, which is ultimately combined with the data from other galaxies in a sample to place posterior constraints on $p(\log\lambda_{\mathrm{sBHAR}}\;|\;\mathcal{M}_{*},z)$.

The expected galaxy count rate, $l_{i}L_{\mathrm{G}}$, depends on a conversion factor,666The $l_{i}$ factor converts from a rest-frame 2–10 keV luminosity to observed 2–7 keV counts, assuming an unabsorbed, $\Gamma=1.9$ X-ray spectrum appropriate for the galactic emission. Thus $l_{i}=\eta(z_{i})\times t_{i}$. $l_{i}$, and the X-ray luminosity due to galactic processes, $L_{\mathrm{G}}$. We can apply a prior for $L_{\mathrm{G}}$ for a given galaxy based on information at other wavelengths: the SFR and $\mathcal{M}_{*}$ measured from the UV-to-IR SED (see Appendix A above). To scale from SFR and $\mathcal{M}_{*}$ to an expected X-ray luminosity, we adopt Model 5 from Paper I, which accounts for contributions from both high- and low-mass X-ray binaries that should roughly scale with SFR and $\mathcal{M}_{*}$, respectively. To allow for uncertainty in the measured SFR and $\mathcal{M}_{*}$, as well as the large uncertainty in the conversion to X-ray luminosity, we describe our prior using a single Gamma distribution for each galaxy,

[TABLE]

with a shape parameter ${\alpha=2}$ and a scale parameter ${\theta=\tilde{L_{\mathrm{G}}}/(\alpha-1)}$. A Gamma distribution with $\alpha=2$ is roughly equivalent to a logarthmic “1$\sigma$” uncertainty of $\sim 0.8$ dex, thus allowing for a large uncertainty in the galactic X-ray luminosity, $L_{\mathrm{G}}$, for an individual galaxy. We choose to model this uncertaintiy with a Gamma distribution (rather than a normal or lognormal distribution) as the Gamma distribution is conjugate to our Poisson likelihood function and thus simplifies our calculations. The scale parameter, $\theta$, is chosen such that the mode of the prior distribution corresponds to $\tilde{L_{\mathrm{G}}}$, our best estimate of the X-ray luminosity from galactic processes based on our measurement of the stellar mass ($\tilde{\mathcal{M}_{i}}$) and SFR ($\tilde{\mathrm{SFR}_{i}}$) for galaxy $i$.

Our overall likelihood function for the data from all galaxies in a given stellar mass–redshift bin ($\mathbf{D}_{\mathrm{bin}}$) is then given by

[TABLE]

where $\lambda\equiv\lambda_{\mathrm{sBHAR}}$ and

[TABLE]

which is the distribution of sBHAR that we model as a mixture of Gamma distributions, as discussed above. As in Paper I, the overall likelihood function described by Equation 10 can be reduced to

[TABLE]

where $A_{j}$ are the normalizations of each mixture component in our model of $p(\log\lambda_{\mathrm{sBHAR}}\;|\;\mathcal{M}_{*},z)$ and $\mathcal{W}_{ij}$ are pre-computed weights for each component $j$ and source $i$. The sum is taken over all $K$ mixture components and the product is taken over all $n_{\mathrm{source}}$ galaxies with a given $\mathcal{M}_{*}$, $z$, and classification.

This revised scheme effectively allows for an additional “background” X-ray emission for every galaxy in our sample based on the expected galactic emission (allowing for an uncertainty in the luminosity of such emission). By modelling the galactic emission separately from the AGN X-ray luminosity we are able to remove this component and recover the underlying distribution of $\lambda_{\mathrm{sBHAR}}$ related to the AGN population within our galaxy samples (removing the low-luminosity peaks found in Paper I). In practice, for galaxies where we directly detect a luminous X-ray AGN ($L_{\mathrm{X}}\gtrsim 10^{42}$ $\mathrm{erg\;s^{-1}}$) the expected galactic contribution is usually negligible. For the large number of galaxies in our samples without direct X-ray detections, but where we still use the available X-ray information, this correction is more important. Furthermore, for high-mass, high-redshift, or high-SFR galaxies the galactic X-ray emission can be comparable to or exceed a weak AGN. Nevertheless, the results of this paper are mostly based on the broad, high-luminosity tail of the luminosity distribution (related to AGN). We generally do not consider ranges of $\lambda_{\mathrm{sBHAR}}$ where the galactic X-ray emission is expected to be the dominant component and our measurements of $p(\log\lambda_{\mathrm{sBHAR}}\;|\;\mathcal{M}_{*},z)$, even after applying the correction described in this appendix, are uncertain.

Appendix C The impact of AGN-dominated sources on measurements of the accretion rate probability distribution

Via our SED-fitting (see Appendix A), we are able to identify roughly 20 per cent of our X-ray detected sample ($\lesssim 0.3$ per cent of the overall galaxy sample) as “AGN-dominated”, where the light associated with the AGN component dominates over any galaxy component at optical wavelengths (specifically, $>$50 per cent of the light at rest-frame 5000 Å is associated with the AGN component in the SED fit). For such sources, it is extremely difficult to accurately measure the SFR of the galaxy (the UV-to-MIR SED template is generally poorly constrained and in many cases the 24µm light is also dominated by the AGN component), precluding a classification of the host galaxy as star-forming or quiescent (see Figure 1). However, as the galaxy component tends to constitute a significant fraction of the light at rest-frame $\sim$1µm ($>$ 50 per cent of the light at 1µm is associated with the galaxy component in 77 per cent of the AGN-dominated sources), we can still obtain reasonably accurate estimates of the host stellar mass (see also Georgakakis et al., 2017). We thus choose to include AGN-dominated sources in our measurements of $p(\log\lambda_{\mathrm{sBHAR}}\;|\;\mathcal{M}_{*},z)$ for all galaxies, but exclude these sources when dividing our galaxy into star-forming and quiescent (see Sections 3.2–3.3). In this appendix, we examine the impact of excluding or including the AGN-dominated sources on our measurements.

In Figure 14 we compare estimates of $p(\log\lambda_{\mathrm{sBHAR}}\;|\;\mathcal{M}_{*},z)$ for star-forming galaxies as a function of stellar mass and redshift when excluding AGN-dominated sources (black hatched regions, corresponding to our main results presented in Figure 3) and when all sources identified as AGN-dominated are included in the star-forming galaxy sample (orange hatched regions). The grey histograms correspond to the observed distribution of $\lambda_{\mathrm{sBHAR}}$ for X-ray detected sources in star-forming galaxies (neglecting any correction for X-ray completeness and excluding AGN-dominated sources), whereas the green histograms indicate the observed distributions for the AGN-dominated sample. We note that including all of the AGN-dominated sources in the star-forming galaxy sample is conservative as some of these sources may be hosted by quiescent galaxies. However, for the majority of the stellar masses and redshifts probed in this study the star-forming galaxies dominate the total number density. Our results also indicate an increase in the probability of hosting an AGN for star-forming galaxies (compared to quiescent galaxies) and independent studies have also found that AGN (both Type-1 and Type-2 sources) are primarily hosted by star-forming galaxies (e.g. Rosario et al., 2012, 2013; Stanley et al., 2015). Thus, it is reasonable to assume that the majority of the AGN-dominated sources in our sample are in fact hosted by star-forming galaxies.

It is immediately apparent in Figure 14 that AGN-dominated sources are preferentially associated with the highest $\lambda_{\mathrm{sBHAR}}$ sources at all stellar masses and redshifts. Assuming that AGN-dominated sources correspond to Type-1, unobscured AGNs, this trend could indicate an association between $\lambda_{\mathrm{sBHAR}}$ (and thus, Eddington ratio) and the levels of obscuration (see also e.g. Raimundo et al., 2010; Merloni et al., 2014). However, to be identified as AGN-dominated with our method requires a high AGN luminosity relative to the host and thus a high $\lambda_{\mathrm{sBHAR}}$. An unobscured, “Type-1” source with a low $\lambda_{\mathrm{sBHAR}}$ would not be classified as AGN-dominated. The association between higher $\lambda_{\mathrm{sBHAR}}$ and identification as an AGN-dominated source is a natural consequence of our method. A full investigation of the relationship between obscuration (at optical or X-ray wavelengths) and $\lambda_{\mathrm{sBHAR}}$ is beyond the scope of this work.

The prevalence of AGN-dominated sources at higher $\lambda_{\mathrm{sBHAR}}$ introduces a small, systematic shift in our estimates of $p(\log\lambda_{\mathrm{sBHAR}}\;|\;\mathcal{M}_{*},z)$, with the orange hatched regions in Figure 14 being systematically higher at $\lambda_{\mathrm{sBHAR}}\gtrsim 0.1$ than when the AGN-dominated sources are excluded (black hatched regions). However, in general the 90 per cent confidence intervals overlap, indicating that the difference is not significant. Even at high $\lambda_{\mathrm{sBHAR}}$, the AGN-dominated sources tend to constitute at most $\sim$50 per cent of our X-ray AGN sample. Our measurements track the broad shape of the probability distribution over logarithmic scales; thus, excluding the AGN-dominated population does not have a significant effect on these measurements. Significant discrepancies (i.e. no overlap between the 90 per cent confidence regions) are seen at the highest accretion rates for the $10<\log\mathcal{M}_{*}/\mathcal{M}_{\odot}<10.5$ stellar mass bin and both the $1.0<z<1.5$ and $1.5<z<2.0$ redshift bins, although the shape of the distribution remains the same and these differences do not alter any of the conclusions we draw in this paper. We thus conclude that our measurements of $p(\log\lambda_{\mathrm{sBHAR}}\;|\;\mathcal{M}_{*},z)$ for the star-forming galaxy sample (as a function of $\mathcal{M}_{*}$ and $z$) are not significantly affected by the exclusion of the AGN-dominated sources and retain the measurements that exclude these sources as our best estimates presented in Figure 3 and used for subsequent analysis.

We note that a similar result (i.e. no significant impact from excluding AGN-dominated sources) is found when considering all galaxies. However, as we are still able to constrain the stellar mass of these sources we choose to adopt the measurements that include the AGN-dominated sources in Figure 2 and subsequent analysis of the “All galaxies” sample.

In Figure 15 we perform a similar analysis but considering the quiescent galaxy samples. We assign all of the AGN-dominated sources (that cannot be classified individually) to the quiescent galaxy sample for the “Including AGN-dominated” estimates (orange hatched regions). As discussed above, we would expect the vast majority of AGN-dominated sources to actually lie in star-forming host galaxies, thus this approach is extremely conservative and provides strict upper limits on $p(\log\lambda_{\mathrm{sBHAR}}\;|\;\mathcal{M}_{*},z)$ for the quiescent galaxy population. Figure 15 shows that including all AGN-dominated sources in the quiescent galaxy sample can have a significant impact on our measurements, mostly due to the relatively small size of our quiescent galaxy samples and the lack of the high $\lambda_{\mathrm{sBHAR}}$ sources within the original sample. Nonetheless, in some bins—given the large uncertainties on the original “Excluding AGN-dominated” measurements—the differences are not significant. In a number of bins (e.g. $1.0<z<1.5$ and $10.0<\log\mathcal{M}_{*}/\mathcal{M}_{\odot}<10.5$) including AGN-dominated sources substantially warps the recovered shape of the probability distribution, possibly indicating that the AGN-dominated sources are drawn from a distinct population. The overall impact of including the AGN-dominated sources is to bring the probability distributions into closer agreement with the measurements for star-forming galaxies.

We thus conclude that excluding AGN-dominated sources from our measurements of $p(\log\lambda_{\mathrm{sBHAR}}\;|\;\mathcal{M}_{*},z)$ for star-forming galaxies does not have a significant impact on our results or conclusions. The impact of AGN-dominated sources (primarily, but not exclusively, expected to be Type-1 AGN) on our measurements of $p(\log\lambda_{\mathrm{sBHAR}}\;|\;\mathcal{M}_{*},z)$ for quiescent galaxies is less clear, but we believe our measurements are robust given the reasonable assumption that the majority of the AGN-dominated sources are in fact hosted by star-forming galaxies.