The Evolving AGN Duty Cycle in Galaxies Since z ∼ 3 as Encoded in the X-Ray Luminosity Function

Our analysis relies on three prior assumptions, which are listed below.

(1) We assume that the X-ray AGN LF is predominantly made by MS and SB galaxies. Passive systems, meant to be galaxies well below the MS relation, are assumed to have a negligible contribution (<10%) at all redshifts and at all L_X. Thus, hereafter we refer to the combined (MS+SB) XLF as the "total XLF." Below we report a number of evidence and quantitative arguments supporting our hypothesis.

The lesser role of quiescent galaxies in the XLF is suggested by studies of the nuclear properties of early-type galaxies, both at z ∼ 0 (e.g., Pellegrini 2010) and at z ∼ 2 (Olsen et al. 2013; Civano et al. 2014). These works generally found low-level X-ray AGN activity, with L_X < 10⁴³ erg s⁻¹, predominantly attributed to free–free emission from hot (T ∼ 10^6–7 K) virialized gas in the galaxy halo (Kim & Fabbiano 2013). Our assumption is also supported by the significantly smaller reservoirs of cold gas measured in passive galaxies compared to those observed in typical galaxies on the MS relation, despite an important redshift increase at least up to z ∼ 1.5 (see Gobat et al. 2018). Another justification comes from the prevalence of radio AGNs within massive and passive galaxies at z < 1.4 (Hickox et al. 2009; Goulding et al. 2014), which display systematically lower λ_EDD (<10⁻³) than X-ray and MIR-selected AGNs (>10⁻²). This is also supported by studies on the intrinsic λ_EDD distribution in quiescent versus star-forming galaxies, reporting systematically lower mean λ_EDD values in quiescent systems (e.g., Wang et al. 2017; Aird et al. 2019). Finally, it is worth noticing that the number density of passive galaxies notably drops at z > 1 (Davidzon et al. 2017), therefore mitigating the incidence of this population at high redshift. Though we acknowledge that passive galaxies might display substantial X-ray emission from hot ionized gas, in this paper we focus on the X-ray emission directly attributed to SMBH accretion.

A quantitative estimate of the subdominant role of quiescent galaxies in the global XLF is presented in Section 3.1.1. Briefly, we conservatively assumed that quiescent galaxies follow the same intrinsic λ_EDD distribution of MS galaxies at each redshift. Then we rescaled the λ_EDD distribution to match empirical mean L_X/M_⋆ measurements for the quiescent population (Carraro et al. 2020). This enabled us to quantify upper limits on the space density and luminosity density of quiescent galaxies, confirming their negligible contribution across the L_X and redshift range explored in this study.

(2) We assume that the intrinsic λ_EDD distribution of AGNs follows a broken-power-law profile, as parameterized in a number of recent studies (e.g., Caplar et al. 2015, 2018; Weigel et al. 2017; Bernhard et al. 2018). This will be further motivated in Section 2.6.

(3) Lastly, we assume that the faint-end (α) and bright-end (β) slopes of the λ_EDD distributions do not differ between MS and SB galaxies (Section 2.6), with only the corresponding break values (${\lambda }_{\mathrm{EDD}}^{* }$) and normalizations being allowed to vary. As a consequence of this assumption, the only free parameter allowed to vary independently among the two populations is ${\lambda }_{\mathrm{EDD}}^{* }$. The main reason is that a simple shift in ${\lambda }_{\mathrm{EDD}}^{* }$ between MS and SB galaxies resembles the well-known double-Gaussian sSFR profile seen in the two populations (e.g., Rodighiero et al. 2011; Sargent et al. 2012). More specifically, since SB-driven star formation is vigorous but short-lived (relative to the galaxy life cycle), it does not primarily drive the growth of galaxy M_⋆. Similarly, we can easily parameterize BH accretion in SBs as having much larger BHAR fluctuations compared to the variation of the cumulative BH mass. In addition, our simplistic treatment of SBs is motivated by the unnecessarily high number of free parameters that would otherwise be allowed to vary simultaneously, leading to large degeneracies and poor constraints on the overall behavior of the λ_EDD function. We briefly discuss the effect of relaxing this condition in Section 2.7. More details on the λ_EDD profiles for the MS and SB populations are given in Section 2.6.

Our approach is schematically summarized in Figure 1 and consists of five steps. We briefly summarize each of them, while a detailed description is presented in the corresponding sections.

• 1.  
  We parameterize the galaxy M_⋆ function of SF galaxies at each redshift (0.5 < z < 3).
• 2.  
  At each M_⋆ and redshift, we assign the corresponding SFR by randomly extracting each value from a lognormal SFR kernel, for both MS and SB galaxies, in order to account for the dispersion of the corresponding locus.
• 3.  
  We derive the expected mean $\langle \mathrm{BHAR}\rangle$ (or mean $\langle {L}_{{\rm{X}}}\rangle$) by multiplying the SFR by various M_⋆-dependent BHAR/SFR relations from the literature.
• 4.  
  We assume a large set of ${\lambda }_{\mathrm{EDD}}$ distributions, each normalized to have a mean value equal to the corresponding $\langle {L}_{{\rm{X}}}\rangle$ expected from (3), at a given M_⋆ and redshift.
• 5.  
  Each simulated λ_EDD distribution is convolved with the M_⋆ function, yielding the predicted XLF. These five steps are run separately for MS and SB galaxies. Each predicted XLF (after combining both MS and SB galaxies) is then compared with the observed XLF of Aird et al. (2015), as detailed in Section 3.

Zoom In Zoom Out Reset image size

In the following sections, we expand each of the above steps in more detail. A comprehensive list of the free parameters adopted in this work is given in Section 2.7 and Table 1, in which we also discuss the effect induced by each assumption.

Table 1.  List of Free Parameters, Ranges, and Relative Assumptions Made in This Work (See Also Section 2.7)

Parameter	Range	Assumptions	Effects
(1)	(2)	(3)	(4)
α	[0.01; 0.55]	α[z = 0.5] = 0.55	Reduce the parameter space
		and evolves as (1 + z)γ	in line with empirical studies
		with γ = [−4.22; 0]	(Section 2.7)
		(Section 2.6.1 and Figure 5)
		α(MS) = α(SB)	Reduce the number of free parameters
		(Section 2.1 and Figure 6)	that could not be constrained
			Link sSFR and sBHAR variations (Section 2.6.1)
		independent of M⋆	Simplify the shape of p(log λEDD)
		(Section 2.6.1)
β	[2, 3, 4, 5]	β(MS) = β(SB)	Reduce the number of free parameters
		(Section 2.1 and Figure 6)	that could not be constrained
			Link sSFR and sBHAR variations (Section 2.6.1)
		independent of M⋆	Simplify the shape of p(log λEDD)
		(Section 2.6.1)
$\mathrm{log}{\lambda }_{\mathrm{MS}}^{* }$	[−1; +0.5]	full range explored at each z	The positive shift with redshift
		(Section 2.6.1)	is genuine (Section 3.4)
		independent of M⋆	The z-evolution of λ* is mirrored in ${L}_{{\rm{X}}}^{* }$
		(Section 2.6.1)	at each ${M}_{\star }$ (Section 3.2)
$\mathrm{log}{\lambda }_{\mathrm{SB}}^{* }$	[−1; +0.5]	${\lambda }_{\mathrm{SB}}^{* }\gt {\lambda }_{\mathrm{MS}}^{* }$	The positive shift with redshift
		(Section 2.6.1)	is induced by ${\lambda }_{\mathrm{MS}}^{* }$ (Section 3.4)
		independent of M⋆	The SB–MS offset in λ* is mirrored in ${L}_{{\rm{X}}}^{* }$
		(Section 2.6.1)	at each M⋆ (Section 3.2)
BHAR/SFR	[0; 1.05]	18 values: 15 around A19	M⋆-dependent BHAR/SFR ratios
slope(**) with log M⋆		+3 to match M12, R15	are favored, but a flat trend is rejected
		and a flat trend at 10−3	at ∼3σ (Section 4.2)
		(Section 2.5 and Figure 4)
		The mean $\langle \mathrm{BHAR}\rangle$ anchors	The minimum ${\lambda }_{\mathrm{MIN}}$ changes with M⋆
		the mean p(log λEDD) value at each M⋆	to accommodate a M⋆-independent
		(Section 2.6.2)	p(log λEDD) shape (Section 2.6.2)
		same relation for MS and SB	The constant $\tfrac{{\mathrm{SFR}}_{\mathrm{SB}}}{{\mathrm{SFR}}_{\mathrm{MS}}}$ induces a constant
		(Section 2.5)	mean $\tfrac{\langle {\mathrm{BHAR}}_{\mathrm{SB}}\rangle }{\langle {\mathrm{BHAR}}_{\mathrm{MS}}\rangle }\,\approx$ 0.8 dex, at each M⋆ and z
			(Section 3.3 and Figure 9)
		full range explored at each z	The nonevolution of this relation
		(Section 2.5)	with redshift is genuine (Section 3.4)

Note. The motivation behind each assumption is described in the corresponding sections. We briefly summarize (column 4) the effect produced by each assumption to help the reader distinguish the genuine trends from those induced by our prior hypotheses. The reference BHAR/SFR trends are taken from Mullaney et al. (2012, M12), Rodighiero et al. (2015, R15), and Aird et al. (2019, A19). (∗∗) The relative normalization is chosen to fit the corresponding data points of Figure 4.

Download table as:  ASCIITypeset image

The first step displayed in Figure 1 consists in setting the input galaxy M_⋆ function at different redshifts. The prescription is taken from Davidzon et al. (2017), who exploited the latest optical to infrared photometry collected in the COSMOS field over the UltraVISTA area (1.5 deg², see Laigle et al. 2016). They provide the M_⋆ function separately between SF and passive galaxies, based on the [NUV − r]/[r − J] colors (Ilbert et al. 2013). Throughout this work, we consider only the M_⋆ function relative to the SF galaxy population (i.e., MS and SB galaxies).

Next, we split the M_⋆ function of SF galaxies among the MS and SB populations. Given that the relative fraction of the two populations has been shown not to vary with M_⋆ (e.g., Rodighiero et al. 2011; Sargent et al. 2012), we only consider how their relative fraction evolves with redshift by following the prescription of Béthermin et al. (2012). The fraction of SB galaxies (f_SB) appears to evolve linearly, from f_SB = 0.015 (z = 0) to f_SB = 0.03 (z = 1), while it stays flat at higher redshifts. By simply scaling the SF galaxy M_⋆ function down by f_SB (or 1 − f_SB), we end up with the M_⋆ function of MS and SB galaxies at each redshift (Figure 2). We note that the input M_⋆ function of Davidzon et al. (2017) is already corrected for the Eddington bias, which might have some impact at the high-M_⋆ end. After dissecting among MS and SB galaxies, we interpolate the corresponding M_⋆ function at redshift z = 0.5, z = 1, z = 2, and z = 3 across a M_⋆ range of 10⁸ < M_⋆ < 10¹² M_⊙.

Zoom In Zoom Out Reset image size

We use the MS prescription presented in Schreiber et al. (2015), which incorporates a redshift evolution and a bending toward higher M_⋆.

Schreiber et al. (2015) studied a sample of Herschel-selected galaxies out to z ∼ 4, dissecting the observed distribution of MS offsets (=SFR/SFR_MS, see their Equation (10) and Figure 19) among the MS and SB populations. By fitting that distribution via a double lognormal function, MS galaxies are centered at 0.87 × SFR_MS, while SB galaxies are centered at 5.3 × SFR_MS (Schreiber et al. 2015). Both relations were rescaled to a Chabrier (2003) IMF. The 1σ dispersion for both the MS and SB loci is assumed to be 0.3 dex (e.g., Speagle et al. 2014). We note that such an MS relation displays a bending toward the highest M_⋆, which makes the transition from MS to SB galaxies not a linear function of sSFR. The adopted MS is qualitatively similar to other recent M_⋆-dependent prescriptions (e.g., Lee et al. 2015; Scoville et al. 2017), while a single-power-law MS would deliver slightly higher SFR estimates (e.g., Rodighiero et al. 2015), yet consistent results within the uncertainties (e.g., Yang et al. 2018).

At fixed M_⋆ and redshift, we account for the MS dispersion by randomly extracting each SFR from a lognormal SFR kernel centered as described above. This is shown in Figure 3 at various redshifts, and separately for MS (solid lines) and SB (dashed lines) galaxies.

Zoom In Zoom Out Reset image size

As shown in Figure 1, the third step consists of converting the derived SFR into BHAR. A number of BHAR/SFR relations have been proposed in the literature, mostly relying on X-ray and IR observations of AGN samples (e.g., Shao et al. 2010; Mullaney et al. 2012; Page et al. 2012; Rosario et al. 2012; Chen et al. 2013; Delvecchio et al. 2015; Stanley et al. 2015).

In order to mitigate possible selection biases induced by short-term (<1 Myr, Schawinski et al. 2015) AGN variability, an effective approach is starting from large M_⋆-selected samples and averaging AGN activity over galaxy timescales (>100 Myr) to unveil the "typical" SMBH accretion rate across the full galaxy life cycle (Hickox et al. 2014).

This approach was first pioneered in Mullaney et al. (2012), who used a Ks-selected galaxy sample to investigate the BHAR/SFR relationship in the GOODS-South field. Interestingly, they found a roughly constant BHAR/SFR ∼ 10⁻³ with redshift (at 0.5 < z < 2.5), which nicely reproduced the local M_BH–M_bulge correlation as the consequence of steady SMBH accretion and SF activity over cosmic time (Kormendy & Ho 2013).

Moreover, Rodighiero et al. (2015) analyzed BzK-selected galaxies at z ∼ 2, split between MS, SB, and passive systems (Daddi et al. 2004), in the COSMOS field (Scoville et al. 2007). The authors found that MS galaxies display a M_⋆-dependent BHAR/SFR relation, as $\mathrm{BHAR}/\mathrm{SFR}\propto {M}_{\star }^{0.44}$. In addition, they argued that SB galaxies at z ∼ 2 show 2× lower BHAR/SFR ratios relative to MS analogs at the same M_⋆.

More recently, Aird et al. (2019) adopted a Bayesian approach to reconstruct the intrinsic λ_EDD distribution across the full galaxy population: they corroborated the need for a linearly M_⋆-dependent BHAR/SFR at 0.5 < z < 2.5, roughly independent of redshift and of galaxy sSFR. Delvecchio et al. (2015) explored the average BHAR/SFR in a sample of Herschel-selected galaxies at z < 0.5, finding no obvious difference when moving above the MS relation. From those studies, it is therefore still unclear whether such a BHAR/SFR relation evolves with M_⋆ and redshift, and whether SB galaxies are truly offset from this trend.

Given these open questions, we prefer to explore a wide set of BHAR/SFR correlations at each redshift, spanning the full parameter space of slopes and normalizations probed by previous studies. Specifically, 15 different slopes have been explored around the most recent derivation of Aird et al. (2019), plus three additional ones to account for different results (Figure 4, green dotted–dashed lines). These 18 slopes are chosen as follows: 15 are uniformly extracted within 2σ from the most recent derivation by Aird et al. (2019); the remaining 3 are taken to match the relationships found by Mullaney et al. (2012) and Rodighiero et al. (2015) and a flat BHAR/SFR = 10⁻³ as the most extreme case. For each slope around the best fit of Aird et al. (2019), the relative normalization is set accordingly to fit the corresponding data points of Figure 4. Therefore, the slope and normalization of each relation are covariant and count as a single free parameter.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

We explore the full set of BHAR/SFR relations at each redshift, assuming that MS and SB galaxies share the same trend, since no stringent constraints on a potential deviation are clearly found in the literature (e.g., Delvecchio et al. 2015; Yang et al. 2018; Aird et al. 2019). Although we acknowledge that some studies argued in favor of a 2× lower BHAR/SFR in SB galaxies at z = 2 relative to MS galaxies, we caution that a substantial BHAR contribution might be highly obscured, especially in compact SB galaxies at high redshift, and unaccounted for via a simple hardness ratio technique (Aird et al. 2015; Bongiorno et al. 2016). We note that a positive redshift dependence was claimed in Yang et al. (2018), who assumed a single-power-law MS relation (from Behroozi et al. 2013) at each redshift. However, if taking a bending MS toward high M_⋆, especially at lower redshifts (e.g., Schreiber et al. 2015; Scoville et al. 2017), we remark that all previous studies are consistent with a redshift-invariant BHAR/SFR ratio.

For each (M_⋆, z, SFR), the resulting BHAR is simply calculated by multiplying the BHAR/SFR ratio by the corresponding SFR obtained from Section 2.4. We stress that such a BHAR is meant to be the "mean" linear BHAR ($\langle \mathrm{BHAR}\rangle$ hereafter). This is connected to the mean X-ray luminosity $\langle {L}_{{\rm{X}}}\rangle$ as follows:

$$\begin{eqnarray}&&\langle {L}_{{\rm{X}}}\rangle =\displaystyle \frac{\epsilon \,{c}^{2}}{1-\epsilon }\cdot \displaystyle \frac{\langle \mathrm{BHAR}\rangle }{{k}_{\mathrm{BOL}}}\end{eqnarray} \tag{ 1 }$$

where c is the speed of light in the vacuum, is the matter-to-radiation conversion efficiency, and k_BOL is the [2–10] keV bolometric correction. If we assume  = 0.1 (e.g., Marconi et al. 2004; Hopkins et al. 2007) and a single k_BOL = 22.4 (median value found by Vasudevan & Fabian 2007 in local AGN samples), Equation (1) reduces to $\langle {L}_{{\rm{X}}}/(\mathrm{erg}\,{{\rm{s}}}^{-1})\rangle =2.8\,\times {10}^{44}\langle \mathrm{BHAR}/[{M}_{\odot }\,{{yr}}^{-1}]\rangle$.

We acknowledge that the k_BOL is known to exhibit a positive L_X-dependence (e.g., Marconi et al. 2004; Lusso et al. 2012). However, in this study we are not assuming a k_BOL. More simply, in order to scale the average BHAR back to L_X, we need to adopt the same k_BOL value (e.g., 22.4) used in previous works (Mullaney et al. 2012; Rodighiero et al. 2015; Aird et al. 2019), otherwise we would obtain inconsistent L_X from what they started. In other words, we used the same k_BOL as previous studies to get rid of the k_BOL dependence when computing the XLF.

In this work, we express the Eddington ratio λ_EDD as a proxy for L_X/M_⋆ (or BHAR/M_⋆), traditionally named "specific L_X" (or "specific BHAR," sBHAR). This formalism has been used by many authors to quantify how fast the SMBH is accreting relative to the M_⋆ of the host galaxy (e.g., Aird et al. 2012, 2018). This quantity is likely more physically meaningful than the absolute L_X, since it accounts for the bias that a more massive galaxy with a given BH λ_EDD would appear more luminous than a less massive galaxy at the same λ_EDD. In particular, assuming a fixed k_BOL (=22.4, see Section 2.5) and a fixed M_BH/M_⋆ ratio of 1/500 (Häring & Rix 2004), the λ_EDD can be linked to L_X/M_⋆ via

$$\begin{eqnarray}&&{\lambda }_{\mathrm{EDD}}=\displaystyle \frac{{k}_{\mathrm{BOL}}\,{L}_{{\rm{X}}}}{1.3\times {10}^{38}\times 0.002\,{M}_{\star }/{M}_{\odot }}.\end{eqnarray} \tag{ 2 }$$

We will briefly discuss in Section 4.5 how a M_⋆-dependent M_BH/M_⋆ ratio (see e.g., Delvecchio et al. 2019) would affect the derived λ_EDD distributions and our global picture of cosmic BH growth.

A single, universal power-law shape was first proposed by Aird et al. (2012) by analyzing the incidence of X-ray AGN activity in galaxies at 0.2 < z < 1.0. The assumption of a broken power law, with a break close to the Eddington limit, has been found to better reproduce the observed shape of the XLF (Aird et al. 2013). Despite the M_⋆-invariant distribution, they observed a steep redshift increase of its normalization, as ∝(1 + z)^3.8 (see also Bongiorno et al. 2012).

In order to reproduce the XLF at z ≳ 1, a M_⋆-dependent λ_EDD distribution has been implemented in a number of studies (e.g., Bongiorno et al. 2016; Georgakakis et al. 2017; Jones et al. 2017; Aird et al. 2018; Bernhard et al. 2018). Building on those findings, we attempt to keep the λ_EDD distribution as simple as possible while being consistent with the BHAR/SFR trends reported in the recent literature.

2.6.1. Broken Power Law and Flattening with Redshift

The fourth step of Figure 1 shows the shape of the assumed λ_EDD distribution. As mentioned in Section 2, we assume a broken power law with faint- and bright-end slopes α and β, respectively, that meet at the break λ*. In order to mitigate the parameter degeneracy, we assume that MS and SB galaxies share the same slopes (α, β), while the corresponding breaks ${\lambda }_{\mathrm{MS}}^{* }$ and ${\lambda }_{\mathrm{SB}}^{* }$ are allowed to vary independently within the range [−1; +0.5] in log space (with steps of 0.1) at each redshift. This λ* range was chosen to be consistent with the typical position of the knee found in recent determinations of the λ_EDD distribution at z ≲ 2 (e.g., Bernhard et al. 2018; Caplar et al. 2018; Aird et al. 2019).

As mentioned in Section 2.1, while a double-Gaussian profile describes the sSFR variation between MS and SB galaxies (e.g., Rodighiero et al. 2011; Sargent et al. 2012), similarly a shift in ${\lambda }_{\mathrm{EDD}}^{* }$ enables us to describe the difference in sBHAR (or Eddington ratio) among those populations. In particular, SB galaxies show intense and short-lasting SFR variations relative to MS analogs, which thus are not important to explaining the growth of galaxy M_⋆ (e.g., Rodighiero et al. 2011; Sargent et al. 2012). With such a formalism, we can parameterize BH accretion in SBs as an intense and short-lasting phenomenon too, characterized by much larger BHAR fluctuations compared to the variation of the cumulative BH mass.

We further reduce the parameter space by imposing that ${\lambda }_{\mathrm{SB}}^{* }\gt {\lambda }_{\mathrm{MS}}^{* }$ in order to reproduce the systematically higher BHAR found in SB galaxies constrained by previous studies (Delvecchio et al. 2015; Rodighiero et al. 2015; Aird et al. 2019; Grimmett et al. 2019). Moreover, both slopes and λ* values are assumed to be M_⋆-invariant. Although we acknowledge that the intrinsic λ_EDD shape might be more complex (M_⋆-dependent, see, e.g., Bernhard et al. 2018; Aird et al. 2019; Grimmett et al. 2019), this simplistic prescription allows us to link the evolution of the mean expected λ_EDD to a rigid shift in λ*. The only foreseen M_⋆ dependence comes from the adopted BHAR/SFR relation (Section 2.5), as described in Section 2.6.2.

In addition, we implement a flattening of the faint-end slope α with redshift. This is supported by recent studies attempting to reproduce the AGN bolometric LF (Caplar et al. 2015; Weigel et al. 2017; Jones et al. 2019) by convolving the galaxy M_⋆ function with some M_⋆-dependent p(log λ_EDD) distribution of AGNs. Similarly for our XLF, if the faint end of the M_⋆ function steepens with redshift (Figure 2) while the faint-end XLF flattens with redshift, this latter feature can be reproduced if the p(log λ_EDD) at low λ_EDD intrinsically flattens with redshift (Bongiorno et al. 2016; Bernhard et al. 2018; Caplar et al. 2018).

We parameterize the redshift flattening of α as follows:

$$\begin{eqnarray}&&\alpha (z)\,=\,\alpha ({z}_{0})\cdot {\left(\displaystyle \frac{1+z}{1+{z}_{0}}\right)}^{\gamma }.\end{eqnarray} \tag{ 3 }$$

For a given λ_EDD distribution, the bright-end slope β is directly linked to the bright-end slope of the AGN bolometric LF (Caplar et al. 2015) because it is much flatter than the exponential decline of the galaxy M_⋆ function at the high-M_⋆ end. Previous studies found that β ranges between 1.8 and 2.5 (Hopkins et al. 2007; Caplar et al. 2015; Caplar et al. 2018). However, steeper slopes might be still accommodated in case multiple contributions are superimposed on one another (i.e., MS and SB). In order to account for this and for a possible redshift evolution of β, we assume β takes the values [2, 3, 4, 5]. As pointed out in Caplar et al. (2018), we stress that changing β within those values has a negligible impact on the integrated X-ray luminosity density. We therefore anticipate that our analysis is not able to tightly constrain this parameter (see Section 3.4). We refer the reader to Table 1 for an exhaustive list of the aforementioned assumptions.

2.6.2. Probability Density Function

To trace the distribution of λ_EDD, we measure the probability density function p($\mathrm{log}{\lambda }_{\mathrm{EDD}}|$M_⋆, z) as a function of (M_⋆, z), which is defined as follows (see Aird et al. 2019):

$$\begin{eqnarray}&&{\int }_{-\infty }^{\infty }p(\mathrm{log}{\lambda }_{\mathrm{EDD}}| {M}_{\star },z)\,d\mathrm{log}{\lambda }_{\mathrm{EDD}}\,\,=\,\,1.\end{eqnarray} \tag{ 4 }$$

This approach assumes that all SMBHs are accreting, however weak their accretion rate might be. Therefore, p(log λ_EDD) reflects the entire distribution of specific L_X/M_⋆ encompassed by SMBHs during their life cycle. According to this formalism, the mean λ_EDD of the model ($\langle {\lambda }_{\mathrm{mod}}\rangle$) defines the "typical" $\langle {L}_{{\rm{X}}}/{M}_{\star }\rangle$ averaged over the entire SMBH life cycle. This quantity can be written as

$$\begin{eqnarray}&&\langle {\lambda }_{\mathrm{mod}}\rangle ={\int }_{\mathrm{log}{\lambda }_{\mathrm{MIN}}}^{\mathrm{log}{\lambda }_{\mathrm{MAX}}}{\lambda }_{\mathrm{EDD}}\cdot p(\mathrm{log}{\lambda }_{\mathrm{EDD}})\,d\mathrm{log}{\lambda }_{\mathrm{EDD}}.\end{eqnarray} \tag{ 5 }$$

For simplicity, we do not assume a M_⋆-dependent shape of the λ_EDD distribution. Instead, at each (M_⋆, z) we tailor the minimum λ_EDD (λ_MIN) in order to normalize our p(log λ_EDD) to 1 (Equation (4)), while anchoring the mean $\langle {\lambda }_{\mathrm{mod}}\rangle$ (Equation (5)) to match empirical BHAR/SFR trends, as explained below.

First, at fixed (M_⋆, z), we can set the expected average SFR (Section 2.4) and the average BHAR (Section 2.5), which yield a mean expected Eddington ratio (or $\langle {L}_{{\rm{X}}}\rangle$), namely $\langle {\lambda }_{\exp }\rangle$.

Second, in order to match $\langle {\lambda }_{\mathrm{mod}}\rangle$ and $\langle {\lambda }_{\exp }\rangle$, we scan each simulated λ_EDD distribution backward from the maximum (log λ_MAX = 2), assuming a logarithmic step Δ(log λ_EDD) = 0.02. At each iteration, we calculate the corresponding $\langle {\lambda }_{\mathrm{mod}}\rangle$ (Equation (5)) and compare it with the expected $\langle {\lambda }_{\exp }\rangle$ taken from Section 2.5. We stop when $\langle {\lambda }_{\mathrm{mod}}\rangle$ equals $\langle {\lambda }_{\exp }\rangle$ within 0.02 dex, which sets λ_MIN. Below this value, we impose p(log λ_EDD) = 0. When log λ_MIN < −6 (i.e., L_X ≲ 10⁴⁰ erg s⁻¹), we truncate the λ_EDD distribution at that value, since current observational data do not probe down to the corresponding L_X (Figure 7). Our arbitrary choice of log λ_MAX = 2 does not impact our procedure because the distribution drops steeply above the Eddington limit (Section 2.7).

We iterate the procedure described above at each M_⋆, redshift, and BHAR/SFR trend, and for every combination of the p(log λ_EDD) parameters: α (or equivalently γ), β, ${\lambda }_{\mathrm{MS}}^{* }$, and ${\lambda }_{\mathrm{SB}}^{* }$.

Figure 6 shows the set of p($\mathrm{log}{\lambda }_{\mathrm{EDD}}$) (for galaxies at M_⋆ = 10^10.5 M_⊙) that best reproduce the observed XLF at different redshifts (see Section 3.1). Each function is defined down to its corresponding λ_MIN and normalized to unity.

Zoom In Zoom Out Reset image size

In this section, we summarize the five free parameters introduced in our analysis: (α, β, ${\lambda }_{\mathrm{MS}}^{* }$, ${\lambda }_{\mathrm{SB}}^{* }$, and the BHAR/SFR relation). A comprehensive list of all prior assumptions made for these parameters is detailed in Table 1. Next to each assumption, we report the effect produced in this work in order to help the reader distinguish between genuine trends and the behaviors obtained by construction.

The faint- and bright-end slopes (α, β) of the λ_EDD distribution are assumed for simplicity to be the same between MS and SB galaxies. Relaxing this condition would increase the parameter degeneracy without adding constraints on the intrinsic α(SB) and β(MS) at low and high L_X, respectively. Specifically, our prior assumption on α consists in a progressive flattening of p(λ_EDD) with redshift, in order to reproduce the faint-end flattening of the XLF toward higher redshifts (Section 2.6.1). We start from α = 0.55 at z = 0.5, which is consistent with the faint-end λ_EDD slope presented in previous studies at z < 1 (α = 0.65 ± 0.05, Aird et al. 2012; α = 0.45, Caplar et al. 2018).

The bright-end slope β is instead assumed to take the values [2, 3, 4, 5] in order to cover the typical range of bright-end slopes found in the AGN bolometric LF (Hopkins et al. 2007; Caplar et al. 2015; Caplar et al. 2018).

The break Eddington ratios of MS and SB galaxies (${\lambda }_{\mathrm{MS}}^{* }$, ${\lambda }_{\mathrm{SB}}^{* }$) are instead allowed to freely vary over the range [−1; +0.5] with a uniform logarithmic step of 0.1. In order to be consistent with recent papers finding systematically higher mean BHAR in SB relative to MS galaxies (Rodighiero et al. 2015; Yang et al. 2018; Aird et al. 2019; Grimmett et al. 2019; Carraro et al. 2020), we accordingly impose that ${\lambda }_{\mathrm{SB}}^{* }\gt {\lambda }_{\mathrm{MS}}^{* }$.

Finally, the BHAR/SFR slope ranges between 0 and 1.05, covering various empirical trends with M_⋆ reported in the recent literature (Section 2.5 and Figure 4). Each normalization is set to minimize the corresponding χ².

By combining the free parameters listed in Table 1, we generate 129,600 predicted XLFs in total. This comes by multiplying the following numbers: 15 (α values), 4 (β values), 18 (BHAR/SFR trends with M_⋆), and (16 × 15)/2 combinations of λ* (accounting for the condition ${\lambda }_{\mathrm{SB}}^{* }\gt {\lambda }_{\mathrm{MS}}^{* }$).

Following Equation (6), we calculate the total XLF and compare each derivation with the latest observed XLF presented in Aird et al. (2015). The authors separately calculated the XLF both in the soft (0.5–2 keV) and hard (2–10 keV) X-ray bands and combined them consistently in a single data set at 2–10 keV. They also subtracted the X-ray emission expected from star formation (Aird et al. 2017), and further corrected for incompleteness and AGN obscuration. Therefore, this compilation is the most complete XLF constrained by X-ray observations over such a luminosity and redshift range.

The observed data points of Aird et al. (2015) are given both in the soft (magenta) and hard (blue) X-ray bands. Some redshift bins in Figure 7 display two data sets from Aird et al. (e.g., at both 0.8 < z < 1.0 and 1.0 < z < 1.2 in our z = 1 bin), which were taken in order to match the mean redshift between the data and our model predictions. Among the Aird et al. (2015) published data points, we exploited only those lying at high enough L_X where the contamination from galaxy star formation is negligible (see Figure 8 in Aird et al. 2015), namely >10^41.3 erg s⁻¹ at z = 0.5 and >10⁴² erg s⁻¹ in the other bins.

Zoom In Zoom Out Reset image size

We select the best XLF model prediction via a simple χ² minimization. Starting from α(z = 0.5) = 0.55, we first identify the γ value that best describes the observed XLF across all redshift bins (i.e., minimizing the global χ² at 0.5 ≤ z ≤ 3). This led us to $\gamma =-{3.16}_{-0.00}^{+0.79}$ ¹⁴ (at 1σ level), which defines the flattening trend with redshift. Second, among the pool of models within 1σ from the best γ, we searched for the best XLF at each redshift, based on χ² minimization.

Although the comparison with Aird et al. (2015) does not allow us to test the separate contributions of MS and SB galaxies to the observed XLF, it is important to verify that our combined (MS+SB) XLF agrees with current data. This is not obvious, as we stress again that our best XLF is not actually a fit but the model prediction that best agrees with the observed XLF of Aird et al. (2015). Figure 7 shows the best XLF (red solid lines) at each redshift and split between the MS (green dashed lines) and SB (blue dashed lines) populations.

The range of XLFs corresponding to ±1σ confidence intervals is enclosed within the corresponding dotted–dashed lines. Such a range is delimited by all the predicted XLFs within a certain Δχ² threshold with N_dof degrees of freedom from the best XLF¹⁵ . The confidence interval around the best XLF also incorporates the propagation of the uncertainties on γ.

The agreement with the XLF of Aird et al. (2015) is striking in all redshift bins, suggesting that our simple statistical approach, constrained by empirical grounds, is able to successfully reproduce the XLF since z ∼ 3 without invoking complex λ_EDD shapes or large numbers of free parameters.

As reported in Aird et al. (2015), the observed XLF is best reproduced with a flexible double-power-law (FDPL) model, incorporating both an L_X-dependent flattening at the faint end and a positive L_X shift with redshift. In our modeling, we also assumed that MS and SB galaxies follow the same intrinsic shape in λ_EDD, independent of M_⋆ (Sections 2.1 and 2.6.1). The only difference in λ_EDD is driven by the corresponding break λ*. With this simple formalism, the flattening of the XLF with redshift is reproduced through a significant flattening of the α slope (Figure 5), whereas the L_X shift is obtained through a gradual predominance of SB galaxies toward higher L_X. This feature comes naturally from our initial assumptions that ${\lambda }_{\mathrm{EDD},\mathrm{SB}}\gt {\lambda }_{\mathrm{EDD},\mathrm{MS}}$, in accordance with the higher specific BHARs found in SB relative to MS galaxies (Delvecchio et al. 2015; Rodighiero et al. 2015; Bernhard et al. 2016; Yang et al. 2018; Aird et al. 2019). This trend also agrees with previous studies (e.g., Caplar et al. 2015, 2018) finding that λ* ∝ 0.048(1 + z)^2.5 at z < 2 and constant at z ≳ 2 (see Figure 10).

Figure 7 clearly shows that MS galaxies dominate Φ(L_X) at L_X ≲ 10^44.5 erg s⁻¹, while SB galaxies tend to take over at higher luminosities, with a crossover ${L}_{{\rm{X}}}^{\mathrm{cross}}$ that slightly evolves with redshift (see Section 4.2 and Figure 13).

3.1.1. Subdominant Contribution of Quiescent Galaxies to the Global XLF

We further explore the differential contribution of galaxies of different M_⋆ to the observed XLF. To do this, we dissect our best model prediction into three M_⋆ bins (M_⋆ < 10¹⁰, 10¹⁰ < M_⋆ < 10¹¹, and M_⋆ > 10¹¹ M_⊙) and separately between MS and SB galaxies, as shown in Figure 8. We remind the reader that the XLF comes from the convolution of the galaxy M_⋆ function and the λ_EDD distribution: because we assumed a M_⋆-independent shape of the λ_EDD and a M_⋆-dependent normalization, the differential contribution in M_⋆ is mostly driven by the M_⋆-dependent BHAR/SFR ratio. This dictates a shift of the mean λ_EDD with M_⋆, which translates to a different break L_X with M_⋆.

Zoom In Zoom Out Reset image size

Given these considerations, unsurprisingly we see that galaxies of different M_⋆ dominate at different L_X. Particularly, M_⋆ < 10¹⁰ galaxies make a negligible contribution to the XLF, in both MS and SB populations, out to z ∼ 2. Instead, at higher redshifts the galaxy M_⋆ function steepens, meaning the number density of less massive to more massive systems increases, thus strengthening the contribution of M_⋆ < 10¹⁰ galaxies, especially at faint L_X. As for galaxies at 10¹⁰ < M_⋆ < 10¹¹ M_⊙, they tend to dominate the XLF up to the knee L_X of the corresponding population, thus contributing to the bulk X-ray emission at all redshifts. Lastly, galaxies with M_⋆ > 10¹¹ M_⊙ contribute to the bright-end XLF at all cosmic epochs, both for MS and SB systems. Starbursts populate higher L_X than MS galaxies matched in M_⋆ and redshift, given their λ_EDD distribution is shifted toward higher values.

Given the derived XLF, we are able to derive the black hole accretion rate density (BHARD or Ψ_BHAR) since z ∼ 3, separately between MS and SB galaxies. This quantity is fundamental for characterizing the overall luminosity-weighted SMBH growth and is defined by the following expression:

$$\begin{eqnarray}&&{{\rm{\Psi }}}_{\mathrm{BHAR}}(z)={\int }_{0}^{\infty }\displaystyle \frac{1-{\epsilon }_{\mathrm{rad}}}{{\epsilon }_{\mathrm{rad}}\,{c}^{2}}\,{L}_{\mathrm{AGN}}\,\phi ({L}_{\mathrm{AGN}})\,d\mathrm{log}{L}_{\mathrm{AGN}}\end{eqnarray} \tag{ 7 }$$

where _rad is the matter-to-radiation conversion efficiency, which is assumed for simplicity to take the constant value _rad = 0.1, in line with previous studies (Marconi et al. 2004; Hopkins et al. 2007; Merloni & Heinz 2008; Delvecchio et al. 2014; Ueda et al. 2014; Aird et al. 2015).

The AGN bolometric luminosity L_AGN is simply scaled from L_X via a set of k_BOL. Once we remove the dependence on the single k_BOL = 22.4 used in previous studies (Section 2.5), we choose to adopt a set of L_X-dependent k_BOL from Lusso et al. (2012) when calculating our Ψ_BHAR(z) estimates. Nevertheless, we will discuss below the effect that a constant k_BOL = 22.4 would have on the estimated Ψ_BHAR(z).

Figure 9 shows the Ψ_BHAR(z) derived from Equation (7) in our four redshift bins (red points). The relative contributions from MS and SB galaxies are the green and blue points, respectively. The hatched areas enclose the ±1σ uncertainties by simply propagating the ±1σ confidence interval of the XLF (see Figure 7). We integrate the upper limit XLF of quiescent galaxies at each redshift (Section 3.1.1) and report the corresponding upper limits on Ψ_BHAR(z) (orange downward arrows). As displayed in Figure 9, the MS population makes the bulk Ψ_BHAR at all cosmic epochs, while SBs are subdominant and display a similar redshift evolution. The integrated BHARD shown in Figure 9 agrees by design with the derivation by Aird et al. (2015, purple stars) and displays a broadly similar shape to that of the star formation rate density (SFRD; Madau & Dickinson 2014), here scaled down by 3300 for illustrative purposes (gray dashed line). This similarity is a natural consequence of the redshift-invariant BHAR/SFR ratio constrained from our analysis (Section 3.4), which is in turn a genuine result of our modeling. The bottom panel of Figure 9 displays the fractional contribution of the SB population (f_SB) relative to the total Ψ_BHAR(z) at each redshift: f_SB ranges between 20% and 30% and stays roughly constant with redshift. The redshift-invariant f_SB is, instead, artificially induced by our assumption that MS and SB galaxies share the same BHAR/SFR ratio at each M_⋆ and redshift. Indeed, the fraction of SFRD made by SB galaxies (Sargent et al. 2012, magenta squares) ranges between 10% and 15%, consistent with the two galaxy populations contributing in similar proportions to both SMBH accretion and star formation at all cosmic epochs (see also Gruppioni et al. 2013; Magnelli et al. 2013).

Zoom In Zoom Out Reset image size

We note that adopting an L_X-dependent k_BOL does change slightly the resulting BHARD relative to the case of a single k_BOL. Indeed, at z = 0.5, the assumption of k_BOL = 22.4 is consistent with the mean L_X-evolving k_BOL at the break L_X. Nevertheless, at z > 0.5, the break L_X shifts to higher values, corresponding to about 2× higher k_BOL, implying a differential boost of the integrated BHARD. Assuming a single k_BOL = 22.4 would instead lower f_SB down to ≈15%–20% at all redshifts.

Our choice of neglecting the quiescent galaxy population might lead us to slightly overestimate the relative contribution of MS and SB galaxies to the full BHARD. We quantify these fractions based on the derived upper limits on the quiescent Ψ_BHAR, being <4.3% (z = 0.5), <6.5% (z = 1), <1.8% (z = 2), and <1.6% (z = 3). If the total Ψ_BHAR was rescaled to accommodate the quiescent population, the relative f_SB would change by the following amount: from 20% to 19% at z = 0.5, from 23% to 21% at z = 1, and unchanged at higher redshifts. We note that these small differences represent upper limits (see Section 3.1). To this end, disentangling the small contribution of the quiescent population is beyond the scope of this paper (but see Bernhard et al. 2018).

In this section, we discuss the uncertainties on each free parameter assumed in this work, as well as their redshift evolution. As listed in Table 1, five free parameters are adopted in this work: the faint- and bright-end slopes of the ${\lambda }_{\mathrm{EDD}}$ distribution (α, β), the break Eddington ratios of MS and SB galaxies (${\lambda }_{\mathrm{MS}}^{* }$, ${\lambda }_{\mathrm{SB}}^{* }$), and the BHAR/SFR relation with ${M}_{\star }$. The combination of these parameters identifies the ±1σ confidence range around the best XLF shown in Figure 7 (crossed lines). We checked that adding the SB component to the XLF improves the reduced ${\chi }^{2}$ at >90% significance level in all redshift bins. Furthermore, from the inspection of the covariance matrices between all free parameters, we verified that their uncertainties seem to be unrelated to each other. Below we discuss the confidence range of each free parameter and refer the reader to Table 2 for a comprehensive list of uncertainties.

The faint-end slope α flattens with redshift in the form $\alpha \propto {(1+z)}^{\gamma }$ (Equation (3)), which yields a best value of $\gamma =-{3.16}_{-0.00}^{+0.79}$. This implies a steep flattening (though not the steepest trend assumed in our input grid), corresponding to a nearly flat ${\lambda }_{\mathrm{EDD}}$ distribution at z = 3 (top panel of Figure 10).

Zoom In Zoom Out Reset image size

The bright-end slope β is unconstrained from our analysis. Among the β values initially allowed by our grid (2, 3, 4, and 5), the resulting ±1σ confidence intervals range between 3 and 5 (see Table 2). Conversely, the flatter value of β = 2 is marginally (at >1σ) disfavored, despite being the closest to what was found in previous studies ($\beta \sim$ 1.8–2.5; e.g., Hopkins et al. 2007; Caplar et al. 2015, Caplar et al. 2018). This apparent discrepancy is likely driven by the fact that two ${\lambda }_{\mathrm{EDD}}$ distributions are here assumed to contribute to the bright-end XLF, with their relative breaks being placed at different ${L}_{{\rm{X}}}$. This overlap generates an overall flatter slope (e.g., $\beta \sim 2$), which in our study is parameterized with two steeper slopes offset from each other.

The break ${\lambda }_{\mathrm{EDD}}$ of MS galaxies, ${\lambda }_{\mathrm{MS}}^{* }$, is allowed to freely vary across the logarithmic range [−1; +0.5] in all redshift bins. Our best solution prefers a progressive shift of ${\lambda }_{\mathrm{MS}}^{* }$ with redshift, from ${\lambda }_{\mathrm{MS}}^{* }=-0.7$ at z = 0.5 to ${\lambda }_{\mathrm{MS}}^{* }=-0.2$ at z = 3 (see Table 2). We parameterize this redshift trend as follows (see green dashed line in Figure 10, bottom panel):

$$\begin{eqnarray}&&{\lambda }_{\mathrm{MS}}^{* }={10}^{-0.9\pm 0.4}\cdot {\left(1+z\right)}^{1.12\pm 0.19}.\end{eqnarray} \tag{ 8 }$$

Not only is ${\lambda }_{\mathrm{MS}}^{* }$ constrained fairly well by our analysis, but its evolution with redshift closely resembles the trend presented in Caplar et al. (2018, hatched area). The authors studied the characteristic ${\lambda }^{* }$ of AGN host galaxies based on the ratio between BHAR and ${M}_{\star }$ density (from Aird et al. 2015 and Ilbert et al. 2013, respectively), from which they identified the corresponding knee ${L}_{\mathrm{AGN}}$ and ${M}_{\star }$. For the whole star-forming galaxy population (i.e., MS and SB), they obtained ${\lambda }^{* }\propto {(1+z)}^{2.5}$ out to z ∼ 2, and constant at z > 2.

In addition, Figure 10 (bottom panel) shows the evolution of the break ${\lambda }_{\mathrm{EDD}}$ of SB galaxies (${\lambda }_{\mathrm{SB}}^{* }$, blue dashed line):

$$\begin{eqnarray}&&{\lambda }_{\mathrm{SB}}^{* }={10}^{0.0\pm 1.0}\cdot {\left(1+z\right)}^{0.80\pm 0.57}.\end{eqnarray} \tag{ 9 }$$

As previously mentioned, our prior assumption that MS and SB galaxies share the same BHAR/SFR ratio, at each ${M}_{\star }$ and redshift, produces a constant gap between ${\lambda }_{\mathrm{SB}}^{* }$ and ${\lambda }_{\mathrm{MS}}^{* }$ of the order of 0.8 dex (i.e., a factor of 6, Section 3.3). This gap mimics the ×6 higher SFR between SB and MS galaxies (Schreiber et al. 2015; Béthermin et al. 2017). Notwithstanding our prior assumptions to induce the redshift trend of ${\lambda }_{\mathrm{SB}}^{* }$, empirical studies do support the condition ${\lambda }_{\mathrm{SB}}^{* }\gt {\lambda }_{\mathrm{MS}}^{* }$ (e.g., Delvecchio et al. 2015; Rodighiero et al. 2015; Aird et al. 2019; Bernhard et al. 2019; Grimmett et al. 2019). An interesting implication is that SB galaxies are expected to host AGNs accreting close to or slightly above the Eddington limit, especially toward higher redshifts. This will be further discussed in Section 4.2.

The BHAR/SFR versus M_⋆ relations assumed in this work include 18 trends taken from the recent literature, as displayed in Figure 4. On the one hand, the best slope and normalization listed in Table 2 at each redshift are consistent with a redshift-invariant relation. On the other hand, our analysis strongly supports a M_⋆-dependent BHAR/SFR relation, with an average slope of ${0.73}_{-0.29}^{+0.22}$, consistent with Aird et al. (2019) within the uncertainties. We verified that flatter relations, like that proposed by Mullaney et al. (2012) or the flat trend BHAR/SFR = 10⁻³, would lead to a higher than observed faint-end XLF, due to excessive mean BHAR values arising from low-mass galaxies.

A number of studies have recently reported an apparently flat, or slightly positive, relationship between average SFR and L_X for X-ray-selected AGNs at various redshifts. The origin of this trend is still debated. On the one hand, at moderate (<10⁴⁴ erg s⁻¹) X-ray luminosities, a flat trend is commonly observed, which argues in favor of a weak dependence between SMBH accretion and star formation in galaxies (Hickox et al. 2014), possibly driven by stochastic fueling mechanisms that wash out a potential correlation with the instantaneous AGN accretion rate traced by X-ray emission (e.g., Rosario et al. 2012; Stanley et al. 2015).

On the other hand, at the highest (quasar-like, >10⁴⁴ erg s⁻¹) X-ray luminosities, other studies argue in favor of a slightly positive trend (e.g., Netzer et al. 2016; Duras et al. 2017; Schulze et al. 2019), suggesting concomitant star formation and AGN activity, possibly driven by major mergers. The transition between these two modes is still poorly understood, as it depends not only on L_X but also on sample selection and redshift. Testing our modeling with a number of observed L_X–SFR trends is therefore a useful test case for double-checking the validity of our approach based on solid empirical grounds.

We collect a compilation of L_X–SFR trends across our full redshift range 0.5 < z < 3, as shown in Figure 11. Samples were taken from Stanley et al. (2015, filled stars), Rosario et al. (2012, filled triangles), Rosario et al. (2013, empty triangles), and Bernhard et al. (2016, filled circles). At z ∼ 2 we collect data from Stanley et al. (2017, empty stars), Scholtz et al. (2018, downward triangles), Duras et al. (2017, empty squares), Schulze et al. (2019, filled squares), and Netzer et al. (2016, empty circles), this latter extending out to z ∼ 3.5. The color bar indicates the redshift of each data set.

Zoom In Zoom Out Reset image size

Solid lines highlight our linear mean SFR estimates, derived in bins of L_X and redshift, by weighting the total (MS and SB combined) XLF by the SFR of each galaxy population contributing at each L_X. We note that each data set was properly rescaled in SFR to match the closest redshift assumed in this work, by a factor corresponding to the evolution of the MS normalization between the two redshifts (at M_⋆ = 10^10.8 M_⊙). Whenever necessary, we converted the SFRs taken from the literature to a Chabrier (2003) IMF. Our modeling displays a good agreement with the observed L_X–SFR trends at all redshifts. The high-L_X "bump" predicted by our curves is likely driven by the gradual predominance of SB galaxies, with the position of the bump shifting with redshift (see top panel of Figure 13).

As discussed in Stanley et al. (2015), the predicted mean SFRs are strongly dependent on the assumed λ_EDD distribution, which incorporates the stochasticity of SMBH accretion. However, we are able to circumvent this issue by constraining the main λ_EDD distribution slope and break via empirical BHAR/SFR relations and comparison with the observed XLF. This test is encouraging, as it demonstrates that our simple, empirically motivated modeling is successful in predicting the average SFR observed across a wide range of L_X and redshift.

Our study supports a nonnegligible contribution of SB galaxies to the integrated BHARD (20%–30%, see Figure 9). While the bulk SMBH accretion history is made by massive (10¹⁰ < M_⋆ < 10¹¹ M_⊙) MS galaxies, the SB population appears to take over at relatively high L_X, enabling us to reproduce the bright-end XLF since z ∼ 3. Specifically, Figure 12 displays the fractional contribution of SB galaxies to the XLF (f_SB) as a function of L_X and redshift. Solid lines mark the best f_SB(L_X), with colors used to indicate our four redshift values. Dashed lines indicate the corresponding ±1σ interval propagated from the XLF of SBs. We observe that f_SB increases with L_X from a few percent (L_X < 10^43.5 erg s⁻¹) to nearly 100% (L_X > 10⁴⁵ erg s⁻¹). As expected, the scatter becomes narrower toward higher L_X, where SBs increasingly dominate. The black horizontal line marks f_SB = 0.5, which identifies the crossover L_X^cross (see top panel of Figure 13). This value scales as

$$\begin{eqnarray}&&{L}_{{\rm{X}}}^{\mathrm{cross}}(\mathrm{erg}\,{{\rm{s}}}^{-1})={10}^{44.36\pm 0.20}\cdot {\left(1+z\right)}^{1.28\pm 0.33}.\end{eqnarray} \tag{ 10 }$$

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Interestingly, the total IR (rest-frame 8–1000 μm) galaxy LF derived with deep Herschel data also displays a strong luminosity (and density) evolution with redshift, as well as a constant ${{\rm{f}}}_{\mathrm{SB}}\sim 20 \%$ (Gruppioni et al. 2013; Magnelli et al. 2013). This similarity is consistent with star formation and AGN accretion in galaxies evolving in a similar fashion over cosmic time, and similarly between MS and SB galaxies.

From our analysis, the best BHAR/SFR slope with M_⋆ is found to range between 0.73 and 0.95, without a significant redshift dependence (see Table 2). These values agree with the linear fit obtained by Aird et al. (2019), while they seem to reject at ∼3σ significance a flat BHAR/SFR trend with M_⋆. This positive M_⋆ dependence introduces a nonlinearity in the cosmic buildup of galaxies and their central SMBHs: while at low M_⋆ the galaxy grows in mass faster than the SMBH (i.e., low BHAR/SFR ratio), when the galaxy reaches high enough M_⋆ the SMBH gradually catches up (high BHAR/SFR ratio). This twofold behavior might be primarily driven by the ability of the dark matter halo mass to set the usable amount of cold gas for the host (Delvecchio et al. 2019).

As mentioned in Section 3.4, another genuine trend is the observed shift of the λ_EDD distribution with redshift, which comes directly from the comparison with the observed XLF of Aird et al. (2015). As a consequence, our results suggest that SB galaxies have a characteristic λ* close to or slightly above the Eddington limit (Section 3.4). Although this might sound unlikely, we note that (1) the uncertainties are broadly consistent with the Eddington limit; (2) the Eddington limit is not a physical boundary, but it can be exceeded in the case of nonspherical gas accretion; (3) theoretical predictions of the λ_EDD distribution support a progressive flattening at low λ_EDD, as well as an increasing fraction of super-Eddington accretion with redshift (e.g., Kawaguchi et al. 2004; Shirakata et al. 2019); (4) we verified that imposing a maximum λ_EDD equal to the Eddington limit fails to reproduce the bright-end XLF at z > 1, thus supporting the need for a minor fraction (1%–5%) of super-Eddington accretion in SB (while only <0.1% for MS) galaxies at >3σ significance.

Our results reveal a significant evolution of λ*_MS with redshift, which induces the offset trend for SB galaxies (bottom panel of Figure 10). This close-to-linear redshift dependence suggests that the probability of finding SMBHs accreting above a certain λ_EDD is higher toward earlier times, for both populations. A qualitatively similar cosmic evolution of the active AGN fraction and λ_EDD distribution is also independently seen in optically selected quasars at 1 < z < 2 (Schulze et al. 2015) and ascribed to an increasing intensity of SMBH growth toward earlier times (see their Figures 18 and 23).

It is well established that galaxies were more gas rich at earlier epochs, with the cold gas fraction f_gas increasing out to z ∼ 2–3 (as f_gas ∝ (1 + z)², e.g., Saintonge et al. 2013). Larger (molecular) gas reservoirs coincide with higher SFR densities via the SK relation, but also with more probable triggering of the central SMBH and the onset of radiative AGN activity (e.g., Alexander & Hickox 2012; Vito et al. 2014). In addition, higher redshift galaxies tend to be more compact and to show more disturbed and irregular morphologies (Förster Schreiber et al. 2009; Kocevski et al. 2012). This profound structural transformation of galaxies over cosmic time might explain the concomitant evolution of the typical λ_EDD of their central SMBHs. A good place to witness the enhancement of both phenomena is given by high-redshift (z > 1) SB galaxies, which are characterized by higher SFE and denser gas reservoirs (Daddi et al. 2010; Genzel et al. 2010) relative to MS analogs at the same redshift. In these systems, several mechanisms, such as major mergers, violent disk instabilities (Bournaud et al. 2011), and cold gas inflows (Di Matteo et al. 2012), might be at play, triggering starbursting star formation, mostly enshrouded in compact and highly obscured molecular clouds. This is observed in numerical simulations to be the case at high redshift, where the typical M_gas exceeds M_⋆ (Dubois et al. 2014, 2016), which makes them ideal environments for triggering highly accreting SMBHs (λ_EDD > 10%).

In the light of the above considerations, our results support a picture in which the cosmic evolution of galaxies' cold gas content might be the main driver of the redshift-invariant BHAR/SFR relation and the positive shift of the λ_EDD distribution.

The evolution of the characteristic λ* with redshift in both MS and SB galaxies links to the question of whether the AGN duty cycle changes over cosmic time. Putting constraints on the relative time spent by AGNs above a certain λ_EDD is crucial for understanding the global incidence of AGN activity across the galaxy population.

Following previous studies on this topic (Aird et al. 2019; Grimmett et al. 2019), we explore the fraction of AGNs accreting above 10% Eddington (f[λ_EDD > 0.1]) and how this evolves with redshift across MS and SB galaxies. We consider the λ_EDD distribution corresponding to the best XLF at each redshift and M_⋆, separately for MS and SB, and integrate each function from λ_MAX = 100 (our maximum value) down to λ_EDD = 0.1. Because the λ_EDD is normalized to unity (Equation (4)), it describes the stochasticity of SMBHs across the entire galaxy life cycle. Therefore, we interpret the fraction above a certain λ_EDD as being proportional to the time spent above that λ_EDD. At each redshift, we weight all the f[λ_EDD > 0.1] estimates obtained at various M_⋆ by the contribution of each M_⋆ bin to the total BHARD at that redshift (Section 3.2). In this way, we infer the luminosity-weighted f[λ_EDD > 0.1] at each redshift, separately for MS and SB galaxies. The bottom panel of Figure 13 shows f[λ_EDD > 0.1] for both MS (green points) and SB (blue points) galaxies. Error bars are given at the 1σ level and incorporate the propagation of the XLF uncertainties. We fit the trend for MS and SB galaxies with a power-law function in (1 + z), obtaining the following expressions:

$$\begin{eqnarray}&&f[{\lambda }_{\mathrm{EDD},\mathrm{MS}}\gt 0.1]={10}^{-2.54\pm 0.19}\cdot {\left(1+z\right)}^{2.26\pm 0.40}\end{eqnarray} \tag{ 11 }$$

$$\begin{eqnarray}&&f[{\lambda }_{\mathrm{EDD},\mathrm{SB}}\gt 0.1]={10}^{-1.82\pm 0.57}\cdot {\left(1+z\right)}^{1.68\pm 1.02},\end{eqnarray} \tag{ 12 }$$

which imply a steep rising toward earlier cosmic epochs, from 0.4% (3.0%) at z = 0.5%–6.5% (15.3%) at z = 3 in MS (SB) galaxies. Data taken from Grimmett et al. (2019, stars) and Aird et al. (2019, triangles) are reported for comparison, based on more complex λ_EDD profiles. Though our f[λ_EDD > 0.1] estimates stand slightly higher than their data, we consistently observe a monotonic increase with redshift, suggesting that SMBHs spend longer at high accretion rates at earlier epochs, with SMBHs in SB galaxies being the most active.

We stress that a longer AGN duty cycle does not necessarily only imply a longer duration of single episodes of AGN activity, but also similar timescales of AGN activity repeated more frequently across the galaxy life cycle. Therefore, the AGN duty cycle we refer to coincides with the average fraction of the SMBH lifetime spent above a given ${\lambda }_{\mathrm{EDD}}$.

Interestingly, our trend broadly follows the evolution of the molecular gas fraction f_gas with redshift, namely f_gas ∝ (1 + z)² (e.g., Saintonge et al. 2013; Tacconi et al. 2013). We speculate that a link between f_gas and typical AGN Eddington ratio (∝L_X/M_BH ∼ L_X/M_⋆ if assuming a fixed M_BH/M_⋆ ratio) could be foreseen if the available cold gas supply regulates the triggering and duration of both stellar and BH growth. At higher redshift, larger gas reservoirs are condensed in more compact regions than in local galaxies, so the gas depletion timescale (t_dep = M_gas/SFR) is shorter and star formation is more likely and takes place more efficiently (Daddi et al. 2010; Genzel et al. 2010). Moreover, in the early universe, the merger rate was significantly higher than today (Le Fèvre et al. 2013), providing a further mechanism for fueling and sustaining SMBH–galaxy growth. In this scenario, if SB-driven star formation and BH accretion are mostly triggered by major mergers (e.g., Calabrò et al. 2019; Cibinel et al. 2019), it might be plausible to expect higher average AGN accretion rates and longer AGN duty cycles (e.g., Glikman et al. 2015; Ricci et al. 2017), as we observe in this work.

Recent studies (e.g., Vito et al. 2018; Cowie et al. 2020) tried to constrain the XLF at 3 < z < 6 based on the currently deepest Chandra data in the 7 Ms Chandra Deep Field South (Luo et al. 2017). An interesting outcome of those observational studies is that the global BHARD seems to display a steep decline at z > 3, much stronger than that observed at z < 1. Since these findings link directly to the integral of the XLF, here we compare our modeling against the most recent data at z > 3, to further test whether our extrapolated XLF at z > 3 yields a similar behavior (Figure 14). To this end, we collect the latest observed XLF at 3 < z < 6 from Vito et al. (2018), centered at z ∼ 3.3 (3 < z < 3.6, yellow circles) and z ∼ 4.1 (3.6 < z < 6, red circles). Then we extrapolate our best XLF (obtained at 0.5 < z < 3) as follows. We convolve the galaxy M_⋆ function (Davidzon et al. 2017) of MS and SB galaxies, at z ∼ 3.3 and z ∼ 4.1, with the corresponding λ_EDD distributions. These latter are derived by extrapolating the redshift evolution of the slopes (α, β) and break λ* of each corresponding population (see Table 2 and Equations (8)–(9)). As shown in Figure 14, our predicted XLF shows a very good agreement with current observational data (Vito et al. 2018) at z ∼ 3.3 and z ∼ 4.1, which we think further proves our modeling solid. We note that the CDF-S is a pencil-beam field, so the brightest data points of Vito et al. (2018) are more uncertain.

Zoom In Zoom Out Reset image size

As discussed in Section 2.6, the behavior of the XLF out to z ∼ 5 is driven by a progressive flattening of the faint-end and a positive shift of the bright-end λ_EDD distribution. As confirmed by previous studies (e.g., Aird et al. 2012; Shankar et al. 2013; Aversa et al. 2015; Weigel et al. 2017; Caplar et al. 2018), such a trend is consistent with an antihierarchical growth of BHs, possibly linked to an intrinsically longer AGN duty cycle (Section 4.3).

We publicly release our best XLF (and its ±1σ confidence intervals) out to z ∼ 5 in the online supplementary material, for both MS and SB galaxies.

We briefly explore the implications of relaxing the assumption of a constant M_BH/M_⋆ = 1/500 (e.g., Häring & Rix 2004, see Figure 6). We stress again that such an assumption does not enter our modeling, which is fully based on the observed L_X/M_⋆ (i.e., sBHAR), while it comes into play when conceptually linking sBHAR to Eddington ratio (e.g., Aird et al. 2012). In this respect, recent studies support a M_⋆-increasing M_BH/M_⋆ ratio in star-forming galaxies, from both observational (Reines & Volonteri 2015; Shankar et al. 2016; Shankar et al. 2019) and theoretical (e.g., Bower et al. 2017; Habouzit et al. 2017; Lupi et al. 2019) arguments. An intriguing implication of this behavior is presented in Delvecchio et al. (2019), where we integrate over time the BHAR/SFR trend obtained in this work (Section 3.4) to track the cosmic assembly of the M_BH–M_⋆ scaling relation. In agreement with the above-mentioned literature, we found that in galaxies with M_⋆ ≳ 10¹⁰ M_⊙ the M_BH/M_⋆ ratio increases with M_⋆ as log(M_BH/M_⋆) = −11.14 + 0.70 × logM_⋆. In less massive galaxies, instead, the relation may flatten out depending on the assumed BH seed mass, which typically ranges between 10² and 10⁶ M_⊙ (e.g., Begelman & Rees 1978).

In this work, we have shown that the bulk of SMBH growth occurs in MS galaxies with 10¹⁰ < M_⋆ < 10¹¹ M_⊙ (Section 3.2 and Figure 8). In that range, the relation obtained in Delvecchio et al. (2019) yields M_BH/M_⋆ ≈ 1/5000. From Equation (2), this implies that the intrinsic Eddington ratio (${\lambda }_{\mathrm{EDD},\mathrm{INTR}}$) would shift up by a factor of 10. Figure 15 shows the comparison between the standard λ_EDD distributions (M_BH/M_⋆ = 1/500, solid lines) obtained in this work against the λ_EDD,INTR distributions (dashed lines) for MS galaxies at different redshifts. The effect of adopting empirical M_BH/M_⋆ ratios is that the break ${\lambda }_{\mathrm{EDD},\mathrm{INTR}}$ shifts above the Eddington limit. This trend is significant, given the relatively small uncertainties on our break λ_EDD (see Table 2). Although we remind the reader that the Eddington limit is physically binding only for idealized conditions of BH accretion (see Section 4.2), our framework predicts between 0.5% (z = 0.5) and 5% (z = 3) super-Eddington BH growth in massive MS galaxies. The mean ${\lambda }_{\mathrm{EDD},\mathrm{INTR}}$ rises from 0.03 (z = 0.5) to 0.19 (z = 3), thus well below the Eddington limit.

Zoom In Zoom Out Reset image size

This effect would be further amplified in M_⋆ < 10¹⁰ M_⊙ galaxies, albeit it could be partly mitigated by assuming quite massive BH seeds (M_BH,SEED ≳ 10⁵ M_⊙) and accounting for the minor contribution of low-M_⋆ galaxies to the integrated SMBH accretion density.

In order to alleviate the deviation from the canonical M_BH/M_⋆ = 1/500, a lower radiative efficiency  < 0.1 could be postulated, though it is insufficient to reconcile the two trends (see Delvecchio et al. 2019 for a detailed discussion). Alternatively, avoiding super-Eddington BH accretion would require the sBHAR distributions to peak at much lower L_X/M_⋆ with decreasing M_⋆, inconsistent with current observables (e.g., Aird et al. 2019). Shedding light on these issues would be relevant for a number of studies focusing on the evolution of λ_EDD distributions (e.g., Aird et al. 2012, 2019; Bongiorno et al. 2012; Weigel et al. 2017; Bernhard et al. 2018, 2019; Caplar et al. 2018; Grimmett et al. 2019).

We thus argue that testing this empirical prediction in the intrinsic Eddington-ratio space is essential to constraining the buildup of the M_BH/M_⋆ relation and standard prescriptions for BH accretion distributions. We propose that a highly complete sample of AGNs above a certain Eddington ratio (obtained via reliable M_BH measurements) would be useful for observationally tying down the average break λ_EDD,INTR in a statistical manner.