# The Evolving AGN Duty Cycle in Galaxies Since $z \sim 3$ as Encoded in the X-Ray Luminosity Function

I. Delvecchio, E. Daddi, J. Aird, J. R. Mullaney, E. Bernhard, L. P. Grimmett, R. Carraro, A. Cimatti, G. Zamorani, N. Caplar, et al.

► To cite this version:

I. Delvecchio, E. Daddi, J. Aird, J. R. Mullaney, E. Bernhard, et al., The Evolving AGN Duty Cycle in Galaxies Since $z \sim 3$ as Encoded in the X-Ray Luminosity Function. The Astrophysical Journal, 2020, 892, 10.3847/1538-4357/ab789c. insu-03746185

HAL Id: insu-03746185

https://insu.hal.science/insu-03746185v1

Submitted on 6 Aug 2025

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
<figure><img src="image_2.png" /></figure>

# The Evolving AGN Duty Cycle in Galaxies Since $z \sim 3$ as Encoded in the X-Ray Luminosity Function

I. Delvecchio<sup>1,2,13</sup>, E. Daddi<sup>1</sup>, J. Aird<sup>3</sup>, J. R. Mullaney<sup>4</sup>, E. Bernhard<sup>4</sup>, L. P. Grimmett<sup>4</sup>, R. Carraro<sup>5</sup>, A. Cimatti<sup>6,7</sup>,
G. Zamorani<sup>8</sup>, N. Caplar<sup>9</sup>, F. Vito<sup>10,11</sup>, D. Elbaz<sup>1</sup>, and G. Rodighiero<sup>12</sup>

<sup>1</sup>CEA, IRFU, DAp, AIM, Université Paris-Saclay, Université Paris Diderot, Sorbonne Paris Cité, CNRS, F-91191 Gif-sur-Yvette, France; ivan.delvecchio@cea.fr

<sup>2</sup>INAF—Osservatorio Astronomico di Brera, via Brera 28, I-20121, Milano, Italy

<sup>3</sup>Department of Physics & Astronomy, University of Leicester, University Road, Leicester LE1 7RJ, UK

<sup>4</sup>Department of Physics and Astronomy, The University of Sheffield, Hounsfield Road, Sheffield S3 7RH, UK

<sup>5</sup>Instituto de Física y Astronomía, Universidad de Valparaíso, Gran Bretaña 1111, Playa Ancha, Valparaíso, Chile

<sup>6</sup>University of Bologna, Department of Physics and Astronomy (DIFA), Via Gobetti 93/2, I-40129, Bologna, Italy

<sup>7</sup>INAF—Osservatorio Astrofisico di Arcetri, Largo E. Fermi 5, I-50125, Firenze, Italy

<sup>8</sup>INAF—Osservatorio di Astrofisica e Scienza dello Spazio, via Gobetti 93/3, I-40129, Bologna, Italy

<sup>9</sup>Department of Astrophysical Sciences, Princeton University, 4 Ivy Ln., Princeton, NJ 08544, USA

<sup>10</sup>Instituto de Astrofísica and Centro de Astroingenieria, Facultad de Física, Pontificia Universidad Católica de Chile, Casilla 306, Santiago 22, Chile

<sup>11</sup>Chinese Academy of Sciences South America Center for Astronomy, National Astronomical Observatories, CAS, Beijing 100012, People's Republic of China

<sup>12</sup>University of Padova, Department of Physics and Astronomy, Vicolo Osservatorio 3, I-35122, Padova, Italy

Received 2019 December 1; revised 2020 February 12; accepted 2020 February 19; published 2020 March 23

## Abstract

We present a new modeling of the X-ray luminosity function (XLF) of active galactic nuclei (AGNs) out to $z \sim 3$, dissecting the contributions of main-sequence (MS) and starburst (SB) galaxies. For each galaxy population, we convolved the observed galaxy stellar mass ($M_*$) function with a grid of $M_*$-independent Eddington ratio ($\lambda_{EDD}$) distributions, normalized via empirical black hole accretion rate (BHAR) to star formation rate (SFR) relations. Our simple approach yields an excellent agreement with the observed XLF since $z \sim 3$. We find that the redshift evolution of the observed XLF can only be reproduced through an intrinsic flattening of the $\lambda_{EDD}$ distribution and with a positive shift of the break $\lambda^*$, consistent with an antihierarchical behavior. The AGN accretion history is predominantly made by massive ($10^{10} < M_* < 10^{11} M_\odot$) MS galaxies, while SB-driven BH accretion, possibly associated with galaxy mergers, becomes dominant only in bright quasars, at $\log(L_X/\text{erg s}^{-1}) > 44.36 + 1.28 \times (1+z)$. We infer that the probability of finding highly accreting ($\lambda_{EDD} > 10\%$) AGNs significantly increases with redshift, from 0.4% (3.0%) at $z = 0.5\%–6.5\%$ (15.3%) at $z = 3$ for MS (SB) galaxies, implying a longer AGN duty cycle in the early universe. Our results strongly favor a $M_*$-dependent ratio between BHAR and SFR, as BHAR/SFR $\propto M_*^{0.73[+0.22, -0.29]}$, supporting a nonlinear BH buildup relative to the host. Finally, this framework opens potential questions on super-Eddington BH accretion and different $\lambda_{EDD}$ prescriptions for understanding the cosmic BH mass assembly.

*Unified Astronomy Thesaurus concepts: AGN host galaxies (2017); Starburst galaxies (1570); Active galactic nuclei (16); Active galaxies (17); Galaxy evolution (594)*

## 1. Introduction

One of the most outstanding achievements of modern astrophysics is the discovery that nearly every galaxy hosts a central supermassive black hole (SMBH) with mass $M_{BH} \sim 10^{6-10} M_\odot$ (e.g., Schmidt 1963; Lynden-Bell 1969). SMBHs are believed to grow in mass via accretion of cold gas within the galaxy, occasionally shining as active galactic nuclei (AGNs; Soltan 1982). Although almost all of today's SMBHs are quiescent, several empirical correlations have been found between $M_{BHs}$ and the properties of local galaxy bulges (e.g., Kormendy & Ho 2013), interpreted as the outcome of a long-lasting interplay between SMBH and galaxy growth (e.g., Magorrian et al. 1998; Ferrarese & Merritt 2000; Gültekin et al. 2009).

To explain this, state-of-the-art numerical simulations advocate for two phases of AGN feedback characterized by high radiative (“quasar mode”) and high kinetic (“jet mode”) luminosities, that combined are able to remove or heat up the gas within the galaxy, via outflows and relativistic jets (Sanders et al. 1988; Fabian 2012). Both types of AGN feedback are invoked for gradually hampering the star-forming (SF) content of massive (stellar mass $M_* > 10^{10} M_\odot$) galaxies,

thus preventing their runaway mass growth (e.g., Hopkins et al. 2008). While observations and models support this AGN-driven “quenching” paradigm to explain the color bimodality and $M_*$ function of local massive systems (e.g., Benson et al. 2003; Morganti et al. 2003, 2005; Croton et al. 2006; Fabian 2012; Heckman & Best 2014), other studies argue in favor of an AGN-driven enhancement of galaxy star formation rate (SFR; Santini et al. 2012; Rosario et al. 2013; Cresci et al. 2015).

Though the origin of the SMBH–galaxy coevolution is not yet fully understood, it is widely accepted that the gas content plays a crucial role in triggering both AGN and star formation activity. Indeed, the SFR is tightly linked to the (molecular) gas content through the Schmidt–Kennicutt relation (hereafter SK relation; Schmidt 1959; Kennicutt 1998). In parallel, radiative AGN activity (i.e., in the X-rays) is observed to be more prevalent in gas-rich, SF galaxies (e.g., Vito et al. 2014), which might explain the observed positive correlations between SFR and average black hole accretion rate (BHAR; e.g., Mullaney et al. 2012). However, still unclear is whether major mergers or secular processes (e.g., violent disk instabilities, minor mergers) are the leading actors in regulating the growth of SMBHs at different luminosities.

<footnote n="13"><sup>13</sup>Marie Curie Fellow.</footnote>
Two main modes of star formation are known to control the growth of galaxies: a relatively steady, secular mode in disk-like galaxies, defining a tight star-forming “main sequence” (MS; Noeske et al. 2007; Elbaz et al. 2011; Speagle et al. 2014; Schreiber et al. 2015) between SFR and $M_*$ (1$\sigma$ dispersion of 0.3 dex); and a “starbursting” mode above the MS, which is interpreted as being driven by mergers (Cibinel et al. 2019). This latter class of starburst (SB) galaxies is usually defined as showing SFR at least 4× above the MS, at fixed $M_*$ (e.g., Rodighiero et al. 2011).

Furthermore, multiple studies corroborated the idea that the cold gas fraction $f_{gas}$ (i.e., the ratio between cold gas mass and total baryonic mass, $M_{gas}/[M_{gas} + M_{*}]$) undergoes a strong redshift evolution ($f_{gas} \propto (1+z)^2$) in MS galaxies from the local universe to $z \sim 2$ (Leroy et al. 2008; Daddi et al. 2010; Tacconi et al. 2010; Geach et al. 2011; Saintonge et al. 2013), with a plateau at higher redshift ($z \sim 3$, Magdis et al. 2013). At fixed $M_{gas}$, SB galaxies are characterized by higher SFRs compared to MS galaxies, implying higher star formation efficiencies (SFE = SFR/$M_{gas}$, Daddi et al. 2010; Genzel et al. 2010).

In this context, Sargent et al. (2012) found that MS and SB galaxies display a bimodal distribution in their specific-SFR (sSFR = SFR/$M_*$), with SB systems contributing to 8%–14% of the total SFR density, up to $z \sim 2$. The luminosity threshold above which SB activity dominates the infrared (IR) luminosity function (LF) evolves with redshift in a fashion similar to the sSFR of MS galaxies (as $\propto (1+z)^{2.9-3.8}$ with the slope depending on $M_*$), suggesting a roughly constant bimodality at least up to $z \sim 2$.

While both galaxy populations are required to reproduce the total IR (8–1000 µm) luminosity function, several studies pointed out intrinsic differences between MS and SB galaxies in terms of structural and physical properties. At $z \sim 0$, MS galaxies are preferentially regular disks and less disturbed compared to SB galaxies, which are instead more compact and mostly identified as merging systems, particularly ultraluminous IR galaxies (ULIRGs; i.e., having IR luminosity $L_{IR} > 10^{12} L_{\odot}$, e.g., Veilleux et al. 2002). At intermediate redshifts ($z \sim 0.7$), Calabrò et al. (2019) observed an increasing incidence of SF clumps when moving above the MS relation, which might indicate a prevalence of merger-induced clumpy star formation toward higher sSFRs. At $z \sim 2$, the morphological dichotomy seen in the local universe becomes much less pronounced, since the fraction of irregular and disturbed morphologies is generally high and spread out quite uniformly across the SFR–$M_*$ plane (e.g., Elmegreen et al. 2007; Förster Schreiber et al. 2009; Kocevski et al. 2012).

Despite much progress having been made in characterizing the star formation, gas content, size, and morphology between MS and SB galaxies, still unclear are their separate contributions to the global SMBH accretion history.

Whether AGN activity and star formation evolve in a similar fashion between MS and SB galaxies is still a matter of debate (see Rodighiero et al. 2019). A seminal study of Mullaney et al. (2012) put forward the idea that the BHAR/SFR ratio is both redshift and $M_*$ invariant at $M > 10^{10} M_{\odot}$ and $0.5 < z < 2.5$. This “hidden AGN main sequence” lies at BHAR/SFR $\sim 10^{-3}$, thus calling for a constant $M_{BH}/M_*$ ratio over cosmic time, which would naturally explain the observed $M_{BH}-M_{bulge}$ relation at $z \sim 0$ (Kormendy & Ho 2013). Lately, other studies have argued in favor of a $M_*$-dependent BHAR/SFR ratio (Rodighiero et al. 2015; Yang et al. 2018; Aird et al. 2019; Bernhard et al. 2019;

Carraro et al. 2020), suggesting that BHAR is enhanced relative to SFR in the most massive galaxies.

Testing whether AGN accretion behaves differently between galaxies on and above the MS relation requires us to dissect the AGN X-ray luminosity function (XLF) into those two galaxy classes and study how they evolve through cosmic time.

This work aims to constrain the relative contribution of MS and SB galaxies to the XLF since $z \sim 3$. In order to avoid selection biases that might arise from collecting AGNs at a particular wavelength or from flux-limited samples, we model the XLF as the convolution between the galaxy $M_*$ function and a large set of Eddington ratio ($\lambda_{EDD}$) distributions that mimics the stochastic nature of AGN activity (e.g., Aird et al. 2013; Conroy & White 2013; Caplar et al. 2015; Jones et al. 2017, 2019; Weigel et al. 2017). Previous works attempting to achieve this goal (Bernhard et al. 2018) successfully reproduced the observed XLF (Aird et al. 2015) out to $z \sim 1.75$ by assuming a $M_*$-dependent shape of the $\lambda_{EDD}$ distribution for SF and quiescent galaxies with relatively complex shapes.

In this work we tackle a simpler approach, showing that a $M_*$-independent shape of the $\lambda_{EDD}$ distribution, scaled with a $M_*$-dependent normalization, is fully able to reproduce the observed XLF out to $z \sim 3$. This method strongly reduces the number of free parameters, while being fully motivated by recent observational grounds (Section 2.2). Moreover, we are able to predict the relative incidence of AGNs of a given $L_X$ and redshift, separately within MS and SB galaxies, putting constraints on the typical SMBH duty cycle on and above the MS. This analysis serves as an important test case for making predictions on the expected SMBH growth rate at different redshifts, $M_*$, and MS offsets.

The manuscript is structured as follows. Section 2 illustrates our initial assumptions and the statistical approach adopted in this work. The best XLF prediction for MS and SB galaxies is presented in Section 3, quantifying its uncertainties and dissecting its evolution with redshift and $L_X$. We further infer the relative contribution of MS and SB galaxies to the global SMBH accretion history since $z \sim 3$. In Section 4 we test our modeling, interpret our findings, and discuss the implications of this study in the framework of SMBH–galaxy evolution since $z \sim 3$. Finally, we list our concluding remarks in Section 5.

Throughout this paper, we adopt a Chabrier (2003) initial mass function (IMF), and we assume a flat cosmology with $\Omega_m = 0.30$, $\Omega_{\Lambda} = 0.70$, and $H_0 = 70 \text{ km s}^{-1} \text{ Mpc}^{-1}$.

## 2. Methodology

The main goal of the present work is to infer how the average BHAR evolves relative to the host-galaxy mass and star formation activity, while matching the observed evolution of the X-ray emission. This analysis further enables us to constrain the occurrence of AGN activity in galaxies since $z \sim 3$ and to dissect the relative contributions of MS and SB populations to the global XLF.

### 2.1. Prior Assumptions

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 \sim 0$ (e.g., Pellegrini 2010) and at $z \sim 2$ (Olsen et al. 2013; Civano et al. 2014). These works generally found low-level X-ray AGN activity, with $L_X < 10^{43} \text{ erg s}^{-1}$, predominantly attributed to free-free emission from hot ($T \sim 10^{6-7} \text{ 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 \sim 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 $\lambda_{\text{EDD}}$ ($<10^{-3}$) than X-ray and MIR-selected AGNs ($>10^{-2}$). This is also supported by studies on the intrinsic $\lambda_{\text{EDD}}$ distribution in quiescent versus star-forming galaxies, reporting systematically lower mean $\lambda_{\text{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 $\lambda_{\text{EDD}}$ distribution of MS galaxies at each redshift. Then we rescaled the $\lambda_{\text{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 $\lambda_{\text{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 $\lambda_{\text{EDD}}$ distributions do not differ between MS and SB galaxies (Section 2.6), with only the corresponding break values ($\lambda_{\text{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_{\text{EDD}}^*$. The main reason is that a simple shift in $\lambda_{\text{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 $\lambda_{\text{EDD}}$ function. We briefly discuss the effect of relaxing this condition in Section 2.7. More details on the $\lambda_{\text{EDD}}$ profiles for the MS and SB populations are given in Section 2.6.

## 2.2. Our Approach

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 \text{BHAR} \rangle$ (or mean $\langle L_X \rangle$) by multiplying the SFR by various $M_*$-dependent BHAR/SFR relations from the literature.

4. We assume a large set of $\lambda_{\text{EDD}}$ distributions, each normalized to have a mean value equal to the corresponding $\langle L_X \rangle$ expected from (3), at a given $M_*$ and redshift.

5. Each simulated $\lambda_{\text{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.

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.

## 2.3. Galaxy Stellar Mass Function

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$^2$, see Laigle et al. 2016). They provide the $M_*$ function separately between SF and passive galaxies, based on the $[\text{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_{\text{SB}}$) appears to evolve linearly, from $f_{\text{SB}} = 0.015$ ($z = 0$) to $f_{\text{SB}} = 0.03$ ($z = 1$), while it stays flat at higher redshifts. By simply scaling the SF galaxy $M_*$ function down by $f_{\text{SB}}$ (or $1 - f_{\text{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
<figure><img src="image_3.png" /><figcaption>**Figure 1.** Sketch of the convolution model used in this work to derive the XLF. The five steps are summarized as follows. (1) We parameterize the galaxy $M_*$ function of SF galaxies at each redshift ($0.5 < z < 3$). (2) At each $M_*$ and redshift, we read the corresponding SFR from a lognormal SFR kernel centered at the mean MS or SB relation. (3) We derive the expected average (BHAR) (or $\langle L_X \rangle$) from $M_*$-dependent BHAR/SFR relations. (4) We assume a large set of $\lambda_{EDD}$ distributions, each normalized to match the corresponding mean $\langle L_X \rangle$ based on (3) at a given $M_*$ and redshift. (5) Each simulated $\lambda_{EDD}$ distribution is convolved with the $M_*$ function (as highlighted in the blue box), yielding the predicted XLF. Each step is described in the corresponding section and run separately for MS and SB galaxies. Our predicted XLF (combining MS and SB galaxies) will then be compared with the observed XLF of Aird et al. (2015) in Section 3.1. See text for details.</figcaption></figure>

corresponding $M_*$ function at redshift $z = 0.5$, $z = 1$, $z = 2$, and $z = 3$ across a $M_*$ range of $10^8 < M_* < 10^{12} M_\odot$.

## 2.4. The MS and the SB Loci

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 \sim 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 \times \text{SFR\_MS}$, while SB galaxies are centered at $5.3 \times \text{SFR\_MS}$ (Schreiber et al. 2015). Both relations were rescaled to a Chabrier (2003) IMF. The $1\sigma$ 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 SFR. 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.

## 2.5. BHAR/SFR Relation with $M_*$

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 $K_s$-selected galaxy sample to investigate the BHAR/SFR relationship in the GOODS-South field. Interestingly, they found a roughly constant BHAR/SFR $\sim 10^{-3}$ with redshift (at $0.5 < z < 2.5$), which nicely reproduced the local $M_{BH}-M_{\text{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 \sim 2$, split between MS, SB, and passive systems (Daddi et al. 2004), in the COSMOS field (Scoville et al. 2007).
<figure><img src="image_4.png" /><figcaption>**Figure 2.** Galaxy $M_*$ function at various redshifts for SF galaxies, taken from Davidzon et al. (2017). Solid and dashed lines mark MS and SB galaxies, respectively.</figcaption></figure>

<figure><img src="image_5.png" /><figcaption>**Figure 3.** Evolving MS relation in the SFR–$M_*$ plane, taken from Schreiber et al. (2015) and scaled to a Chabrier (2003) IMF. Solid and dashed lines highlight the central loci of MS and SB galaxies, respectively. These were defined as $0.87 \times \text{SFR}_{\text{MS}}$ for MS galaxies and $5.3 \times \text{SFR}_{\text{MS}}$ for the SB population (Schreiber et al. 2015).</figcaption></figure>

The authors found that MS galaxies display a $M_*$-dependent BHAR/SFR relation, as BHAR/SFR $\propto M_*^{0.44}$. In addition, they argued that SB galaxies at $z \sim 2$ show $2 \times$ 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 $\lambda_{\text{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

<figure><img src="image_6.png" /><figcaption>**Figure 4.** Relationship between BHAR and SFR as a function of $M_*$. Dashed lines represent a least squares fitting (in the log–log space) of the data presented in Mullaney et al. (2012, red triangles), Rodighiero et al. (2015, blue squares), and Aird et al. (2019, black circles). The best-fit slope and normalization of each fit are reported in the legend. The green dotted–dashed lines indicate the 18 different BHAR/SFR relations explored in this work: 15 of them are taken within $2\sigma$ around the prescription of Aird et al. (2019), while the remaining three are taken to match the Rodighiero et al., Mullaney et al., and flat BHAR/SFR $= 10^{-3}$ trends. See Section 2.5 for more details.</figcaption></figure>

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\sigma$ 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^{-3}$ 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.

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 \times$ 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, \text{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 (⟨BHAR⟩ hereafter). This is
<table><caption>**Table 1**

List of Free Parameters, Ranges, and Relative Assumptions Made in This Work (See Also Section 2.7)</caption><thead><tr><th>Parameter<br>(1)</th><th>Range<br>(2)</th><th>Assumptions<br>(3)</th><th>Effects<br>(4)</th></tr></thead><tbody><tr><td>α</td><td>[0.01; 0.55]</td><td>α[z = 0.5] = 0.55<br>and evolves as (1 + z)γ<br>with γ = [-4.22; 0]<br>(Section 2.6.1 and Figure 5)<br>α(MS) = α(SB)<br>(Section 2.1 and Figure 6)<br>independent of M<sub>⋆</sub><br>(Section 2.6.1)</td><td>Reduce the parameter space<br>in line with empirical studies<br>(Section 2.7)<br>Reduce the number of free parameters<br>that could not be constrained<br>Link sSFR and sBHAR variations (Section 2.6.1)<br>Simplify the shape of p(log λ<sub>EDD</sub>)</td></tr><tr><td>β</td><td>[2, 3, 4, 5]</td><td>β(MS) = β(SB)<br>(Section 2.1 and Figure 6)<br>independent of M<sub>⋆</sub><br>(Section 2.6.1)</td><td>Reduce the number of free parameters<br>that could not be constrained<br>Link sSFR and sBHAR variations (Section 2.6.1)<br>Simplify the shape of p(log λ<sub>EDD</sub>)</td></tr><tr><td>log λ<sup>*</sup><sub>MS</sub></td><td>[-1; +0.5]</td><td>full range explored at each z<br>(Section 2.6.1)<br>independent of M<sub>⋆</sub><br>(Section 2.6.1)</td><td>The positive shift with redshift<br>is genuine (Section 3.4)<br>The z-evolution of λ<sup>*</sup> is mirrored in L<sup>*</sup><sub>X</sub><br>at each M<sub>⋆</sub> (Section 3.2)</td></tr><tr><td>log λ<sup>*</sup><sub>SB</sub></td><td>[-1; +0.5]</td><td>λ<sup>*</sup><sub>SB</sub> > λ<sup>*</sup><sub>MS</sub><br>(Section 2.6.1)<br>independent of M<sub>⋆</sub><br>(Section 2.6.1)</td><td>The positive shift with redshift<br>is induced by λ<sup>*</sup><sub>MS</sub> (Section 3.4)<br>The SB-MS offset in λ<sup>*</sup> is mirrored in L<sup>*</sup><sub>X</sub><br>at each M<sub>⋆</sub> (Section 3.2)</td></tr><tr><td>BHAR/SFR<br>slope(**) with log M<sub>⋆</sub></td><td>[0; 1.05]</td><td>18 values: 15 around A19<br>+3 to match M12, R15<br>and a flat trend at 10<sup>-3</sup><br>(Section 2.5 and Figure 4)<br>The mean ⟨BHAR⟩ anchors<br>the mean p(log λ<sub>EDD</sub>) value at each M<sub>⋆</sub><br>(Section 2.6.2)<br>same relation for MS and SB<br>(Section 2.5)<br>full range explored at each z<br>(Section 2.5)</td><td>M<sub>⋆</sub>-dependent BHAR/SFR ratios<br>are favored, but a flat trend is rejected<br>at ~3σ (Section 4.2)<br>The minimum λ<sub>MIN</sub> changes with M<sub>⋆</sub><br>to accommodate a M<sub>⋆</sub>-independent<br>p(log λ<sub>EDD</sub>) shape (Section 2.6.2)<br>The constant SFR<sub>SB</sub>/SFR<sub>MS</sub> induces a constant<br>mean ⟨BHAR<sub>SB</sub>/BHAR<sub>MS</sub>⟩ ≈ 0.8 dex, at each M<sub>⋆</sub> and z<br>(Section 3.3 and Figure 9)<br>The nonevolution of this relation<br>with redshift is genuine (Section 3.4)</td></tr></tbody></table>

**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.

connected to the mean X-ray luminosity $\langle L_X \rangle$ as follows:

$$ \langle L_X \rangle = \frac{\epsilon c^2}{1 - \epsilon} \cdot \frac{\langle \text{BHAR} \rangle}{k_{\text{BOL}}} \quad (1) $$

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

We acknowledge that the $k_{\text{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_{\text{BOL}}$. More simply, in order to scale the average BHAR back to $L_X$, we need to adopt the same $k_{\text{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_{\text{BOL}}$ as previous studies to get rid of the $k_{\text{BOL}}$ dependence when computing the XLF.

## 2.6. Eddington Ratio Distribution of AGN

In this work, we express the Eddington ratio $\lambda_{\text{EDD}}$ as a proxy for $L_X/M_\star$ (or BHAR/$M_\star$), 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_\star$ 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 $\lambda_{\text{EDD}}$ would appear more luminous than a less massive galaxy at the same $\lambda_{\text{EDD}}$. In particular, assuming a fixed $k_{\text{BOL}} (=22.4$, see Section 2.5) and
a fixed $M_{BH}/M_{\star}$ ratio of 1/500 (Häring & Rix 2004), the $\lambda_{EDD}$ can be linked to $L_X/M_{\star}$ via

$$ \lambda_{\text{EDD}} = \frac{k_{\text{BOL}} L_X}{1.3 \times 10^{38} \times 0.002 M_{\star}/M_{\odot}}. \qquad (2) $$

We will briefly discuss in Section 4.5 how a $M_{\star}$-dependent $M_{BH}/M_{\star}$ ratio (see e.g., Delvecchio et al. 2019) would affect the derived $\lambda_{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_{\star}$-invariant distribution, they observed a steep redshift increase of its normalization, as $\propto(1+z)^{3.8}$ (see also Bongiorno et al. 2012).

In order to reproduce the XLF at $z \gtrsim 1$, a $M_{\star}$-dependent $\lambda_{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 $\lambda_{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 $\lambda_{EDD}$ distribution. As mentioned in Section 2, we assume a broken power law with faint- and bright-end slopes $\alpha$ and $\beta$, respectively, that meet at the break $\lambda^{*}$. In order to mitigate the parameter degeneracy, we assume that MS and SB galaxies share the same slopes ($\alpha, \beta$), while the corresponding breaks $\lambda_{MS}^{*}$ and $\lambda_{SB}^{*}$ are allowed to vary independently within the range $[-1; +0.5]$ in log space (with steps of 0.1) at each redshift. This $\lambda^{*}$ range was chosen to be consistent with the typical position of the knee found in recent determinations of the $\lambda_{EDD}$ distribution at $z \lesssim 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_{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_{\star}$ (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_{SB}^{*} > \lambda_{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 $\lambda^{*}$ values are assumed to be $M_{\star}$-invariant. Although we acknowledge that the intrinsic $\lambda_{EDD}$ shape might be more complex ($M_{\star}$-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 $\lambda_{EDD}$ to a rigid shift in $\lambda^{*}$. The only foreseen $M_{\star}$ 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 $\alpha$ 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_{\star}$ function with some $M_{\star}$-dependent $p(\log \lambda_{EDD})$ distribution of AGNs. Similarly for our XLF, if the faint end of the $M_{\star}$ function steepens with redshift (Figure 2) while the faint-end XLF flattens with redshift, this latter feature can be reproduced if the $p(\log \lambda_{EDD})$ at low $\lambda_{EDD}$ intrinsically flattens with redshift (Bongiorno et al. 2016; Bernhard et al. 2018; Caplar et al. 2018).

We parameterize the redshift flattening of $\alpha$ as follows:

$$ \alpha(z) = \alpha(z_0) \cdot \left( \frac{1+z}{1+z_0} \right)^{\gamma}. \qquad (3) $$

For a given $\lambda_{EDD}$ distribution, the bright-end slope $\beta$ 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_{\star}$ function at the high-$M_{\star}$ end. Previous studies found that $\beta$ 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 $\beta$, we assume $\beta$ takes the values [2, 3, 4, 5]. As pointed out in Caplar et al. (2018), we stress that changing $\beta$ 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 $\lambda_{EDD}$, we measure the probability density function $p(\log \lambda_{EDD}|M_{\star}, z)$ as a function of ($M_{\star}, z$), which is defined as follows (see Aird et al. 2019):

$$ \int_{-\infty}^{\infty} p(\log \lambda_{\text{EDD}}|M_{\star}, z) d \log \lambda_{\text{EDD}} = 1. \qquad (4) $$

This approach assumes that all SMBHs are accreting, however weak their accretion rate might be. Therefore, $p(\log \lambda_{EDD})$ reflects the entire distribution of specific $L_X/M_{\star}$ encompassed by SMBHs during their life cycle. According to this formalism, the mean $\lambda_{EDD}$ of the model ($\langle \lambda_{\text{mod}} \rangle$) defines the “typical” $\langle L_X/M_{\star} \rangle$ averaged over the entire SMBH life cycle. This quantity can be written as

$$ \langle \lambda_{\text{mod}} \rangle = \int_{\log \lambda_{\text{MIN}}}^{\log \lambda_{\text{MAX}}} \lambda_{\text{EDD}} \cdot p(\log \lambda_{\text{EDD}}) d \log \lambda_{\text{EDD}}. \qquad (5) $$

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

First, at fixed ($M_{\star}, 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_X \rangle$), namely $\langle \lambda_{\text{exp}} \rangle$.

Second, in order to match $\langle \lambda_{\text{mod}} \rangle$ and $\langle \lambda_{\text{exp}} \rangle$, we scan each simulated $\lambda_{EDD}$ distribution backward from the maximum (log
$\lambda_{MAX} = 2$), assuming a logarithmic step $\Delta(\log \lambda_{EDD}) = 0.02$. At each iteration, we calculate the corresponding $\langle \lambda_{mod} \rangle$ (Equation (5)) and compare it with the expected $\langle \lambda_{exp} \rangle$ taken from Section 2.5. We stop when $\langle \lambda_{mod} \rangle$ equals $\langle \lambda_{exp} \rangle$ within 0.02 dex, which sets $\lambda_{MIN}$. Below this value, we impose $p(\log \lambda_{EDD}) = 0$. When $\log \lambda_{MIN} < -6$ (i.e., $L_X \lesssim 10^{40}$ erg s⁻¹), we truncate the $\lambda_{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 $\lambda_{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 \lambda_{EDD})$ parameters: $\alpha$ (or equivalently $\gamma$), $\beta$, $\lambda_{MS}^*$, and $\lambda_{SB}^*$.

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

## 2.7. Free Parameters

In this section, we summarize the five free parameters introduced in our analysis: ($\alpha, \beta, \lambda_{MS}^*, \lambda_{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 ($\alpha, \beta$) of the $\lambda_{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 $\alpha$(SB) and $\beta$(MS) at low and high $L_X$, respectively. Specifically, our prior assumption on $\alpha$ consists in a progressive flattening of $p(\lambda_{EDD})$ with redshift, in order to reproduce the faint-end flattening of the XLF toward higher redshifts (Section 2.6.1). We start from $\alpha = 0.55$ at $z = 0.5$, which is consistent with the faint-end $\lambda_{EDD}$ slope presented in previous studies at $z < 1$ ($\alpha = 0.65 \pm 0.05$, Aird et al. 2012; $\alpha = 0.45$, Caplar et al. 2018).

The bright-end slope $\beta$ 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_{MS}^*$, $\lambda_{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_{SB}^* > \lambda_{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 $\chi^2$.

## 3. Results

In this section, we present the results of our modeling to reproduce the XLF since $z \sim 3$. The galaxy $M_*$ function (Section 2.3) was convolved with a set of $M_*$-independent Eddington ratio parameters (slopes and break, see Sections 2.6.1 and 2.6.2), but with a $M_*$-dependent normalization that matches the mean BHAR from several BHAR/SFR relations found in the literature (Section 2.5). This analysis was run separately among MS and SB galaxies, which allowed us to infer the relative contribution of each class to the total XLF. With this formalism, the XLF $\Phi(L_X, z)$ was derived as follows:

$$ \Phi(L_X, z) = \int_{M_*} \Phi_*(M'_*, z) \otimes p(L_X|M'_*, z) dM'_* \quad (6) $$

where $\Phi_*(M'_*, z)$ is the galaxy $M_*$ function of the corresponding population, and $p(L_X|M'_*, z)$ is the likelihood distribution of log $L_X$ as a function of ($M_*, z$). The total XLF split between MS and SB galaxies is shown in Section 3.1. The best parameters, along with their uncertainties and confidence ranges, are listed in Table 2. The degeneracy and the evolution of each free parameter are discussed in Section 3.4, where we also present the SMBH accretion rate density dissected for the first time among those two populations.

### 3.1. Total XLF of MS and SB Galaxies since $z \sim 3$

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 ($\alpha$ values), 4 ($\beta$ values), 18 (BHAR/SFR trends with $M_*$), and $(16 \times 15)/2$ combinations of $\lambda^*$ (accounting for the condition $\lambda_{SB}^* > \lambda_{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^{42}$ erg s⁻¹ in the other bins.

We select the best XLF model prediction via a simple $\chi^2$ minimization. Starting from $\alpha(z=0.5) = 0.55$, we first identify the $\gamma$ value that best describes the observed XLF across all redshift bins (i.e., minimizing the global $\chi^2$ at $0.5 \le z \le 3$). This led us to $\gamma = -3.16_{-0.00}^{+0.79} \quad (14$ (at $1\sigma$ level), which defines the flattening trend with redshift. Second, among the pool of models within $1\sigma$ from the best $\gamma$, we searched for the best XLF at each redshift, based on $\chi^2$ 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

<footnote n="14">¹⁴ Zero errors are due to our discrete grid and should be interpreted as smaller than the closest value (see Table 2).</footnote>
<figure><img src="image_7.png" /><figcaption>**Figure 5.** Colored lines show the 15 $\alpha$ vs. redshift trends assumed in this work. Each line marks a different $\gamma$ value, spanning from steep to flat trends with redshift, reaching to $\alpha \sim 0$ at $z=3$ in the most extreme case. We set $\alpha = 0.55$ at $z = 0.5$ to be consistent with previous studies ($\alpha = 0.65 \pm 0.05$ at $0.2 < z < 1.0$, Aird et al. 2012; $\alpha = 0.45$, Caplar et al. 2018).</figcaption></figure>

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 $\pm 1\sigma$ confidence intervals is enclosed within the corresponding dotted-dashed lines. Such a range is delimited by all the predicted XLFs within a certain $\Delta\chi^2$ threshold with $N_{\text{dof}}$ degrees of freedom from the best XLF$^{15}$. The confidence interval around the best XLF also incorporates the propagation of the uncertainties on $\gamma$.

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 \sim 3$ without invoking complex $\lambda_{\text{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 $\lambda_{\text{EDD}}$, independent of $M_*$ (Sections 2.1 and 2.6.1). The only difference in $\lambda_{\text{EDD}}$ is driven by the corresponding break $\lambda^*$. With this simple formalism, the flattening of the XLF with redshift is reproduced through a significant flattening of the $\alpha$ 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_{\text{EDD,SB}} > \lambda_{\text{EDD,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

<figure><img src="image_8.png" /><figcaption>**Figure 6.** Set of $\lambda_{\text{EDD}}$ probability distributions that best reproduce the observed XLF (see Section 3.1), here shown only at $M_* = 10^{10.5} M_{\odot}$ for illustrative purposes. We indicate MS and SB galaxies with solid and dashed lines, respectively. Colors mark our four redshift bins. All functions are calculated down to their minimum $\lambda_{\text{MIN}}$ and normalized to unity, as detailed in Section 2.6.2.</figcaption></figure>

$\lambda^* \propto 0.048(1+z)^{2.5}$ at $z < 2$ and constant at $z \gtrsim 2$ (see Figure 10).

Figure 7 clearly shows that MS galaxies dominate $\Phi(L_X)$ at $L_X \lesssim 10^{44.5}$ erg s$^{-1}$, while SB galaxies tend to take over at higher luminosities, with a crossover $L_X^{\text{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

As mentioned in Section 2.1, here we quantitatively address the contribution of quiescent galaxies to the global XLF at various redshifts. For consistency, we adopt an approach similar to that presented for star-forming galaxies (Section 2.2). We convolve the galaxy $M_*$ function of quiescent galaxies (Davidzon et al. 2017) with a set of $\lambda_{\text{EDD}}$ distributions. Instead of plugging into our modeling five additional free parameters for the quiescent galaxy population, we conservatively assume that they share the same best $\lambda_{\text{EDD}}$ parameters of MS galaxies ($\alpha, \beta, \lambda_{\text{MS}}^*$) inferred from Section 3.1. The normalization of the corresponding $\lambda_{\text{EDD}}$ distribution is, however, set differently to match empirical studies of quiescent galaxies. We base our validation on recent observational grounds by Carraro et al. (2020), who stacked deep Chandra images of the COSMOS field, including a $M_*$ complete sample of quiescent galaxies. They inferred mean $L_X$ as a function of $M_*$ and redshift out to $z \sim 3$, which we use to anchor the mean $\lambda_{\text{EDD}}$ of the corresponding distributions (as presented in Section 2.6.2). Our convolution yields the XLF of quiescent galaxies out to $z \sim 3$ (orange dashed-dotted lines in Figure 7). While we acknowledge that our approach is not formally the same as that adopted for MS and SB galaxies, we stress that our reasoning is largely supported by observational studies finding lower mean BHAR (or $\lambda_{\text{EDD}}$) in quiescent galaxies compared to star-forming analogs (e.g., Rodighiero et al. 2015; Bernhard et al. 2018; Yang et al. 2018; Aird et al. 2019). We therefore

<footnote>$^{15}$ $N_{\text{dof}}$ is the difference between the observed data points $N_d$ and the number of free parameters of each redshift bin (i.e., $N_d-4$ at $z = 0.5$, $N_d-5$ in the other bins, corresponding to $\Delta\chi^2 = 4.71$ and 5.88, respectively).</footnote>
<figure><img src="image_9.png" /><figcaption>**Figure 7.** The best 2–10 keV AGN X-ray luminosity function predicted by our modeling at various redshifts (red solid lines). The XLF is dissected between MS (green dashed lines) and SB (blue dashed lines) galaxies at each redshift, with their ±1σ confidence intervals enclosed by the corresponding dotted–dashed lines. The upper limit XLF made by quiescent galaxies (orange dotted–dashed line) is detailed in Section 3.1.1. Data points are from the compilation of Aird et al. (2015), containing data in both the soft (0.5–2 keV, magenta points) and hard (2–10 keV, blue points) X-ray bands.</figcaption></figure>

interpret the resulting quiescent XLF as an upper limit. Nevertheless, if the intrinsic $\lambda_{EDD}$ distribution differs from that of MS galaxies, we stress that the normalization and break $L_X$ of the quiescent XLF would change in opposite directions. Given the generally negligible contribution of quiescent galaxies displayed at all $L_X$ (Figure 7), we might expect them to become potentially comparable to SB galaxies only at low $L_X$, where SBs are already subdominant and poorly constrained. We stress that the mean stacked $L_X$ reported by Carraro et al. display a positive dependence on both $M_*$ and redshift, in accordance with previous studies (see, e.g., Wang et al. 2017; Aird et al. 2018; Bernhard et al. 2018). The subdominant role highlighted by this test might be caused by the steep drop of the galaxy quiescent $M_*$ function toward higher redshift, which counterbalances the increasing X-ray AGN fraction (but see Georgakakis et al. 2014).

It is worth noticing that neglecting the quiescent galaxy population does not contradict the observed prevalence of X-ray AGNs with SFR $2\times$ below the MS (e.g., Mullaney et al. 2015). Indeed, the definition of MS adopted in this work accounts for a $1\sigma$ scatter of a factor of two (Section 2.4). This implies that X-ray AGNs lying $2\times$ below the MS are within the MS locus, and therefore would not nominally contribute to quiescent galaxies. In addition, X-ray AGNs found below the MS display Seyfert-like luminosities ($L_X < 10^{44} \text{ erg s}^{-1}$), thus consistent with MS galaxies hosting moderately luminous AGNs. In summary, our check demonstrates that quiescent galaxies make a very minor ($<10\%$) contribution to the space

density of X-ray AGNs at all $L_X$ and redshifts analyzed in this work. We believe this justifies not plugging them into our modeling.

### 3.2. XLF Split in $M_*$ Bins

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}$, $10^{10} < M_* < 10^{11}$, and $M_* > 10^{11} M_\odot$) 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 $\lambda_{EDD}$ distribution: because we assumed a $M_*$-independent shape of the $\lambda_{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 $\lambda_{EDD}$ with $M_*$, which translates to a different break $L_X$ with $M_*$.

Given these considerations, unsurprisingly we see that galaxies of different $M_*$ dominate at different $L_X$. Particularly, $M_* < 10^{10}$ galaxies make a negligible contribution to the XLF, in both MS and SB populations, out to $z \sim 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^{10}$ galaxies, especially at faint $L_X$. As for galaxies at $10^{10} < M_* < 10^{11} M_\odot$, they tend to dominate the XLF up to the knee $L_X$ of the corresponding population, thus contributing to the bulk X-ray
<figure><img src="image_10.png" /><figcaption>Figure 8. Best 2–10 keV AGN X-ray luminosity function dissected among three $M_*$ bins: $M_* < 10^{10}$ (dashed lines), $10^{10} < M_* < 10^{11}$ (dotted–dashed lines), and $M_* > 10^{11}$ M$_\odot$ (solid lines), separately for MS (green) and SB (blue) galaxies.</figcaption></figure>

emission at all redshifts. Lastly, galaxies with $M_* > 10^{11}$ M$_\odot$ 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 $\lambda_{EDD}$ distribution is shifted toward higher values.

### 3.3. Dissecting the SMBH Accretion History among MS and SB Galaxies

Given the derived XLF, we are able to derive the black hole accretion rate density (BHARD or $\Psi_{BHAR}$) since $z \sim 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:

$$ \Psi_{BHAR}(z) = \int_{0}^{\infty} \frac{1 - \epsilon_{rad}}{\epsilon_{rad} c^2} L_{AGN} \phi(L_{AGN}) d \log L_{AGN} \quad (7) $$

where $\epsilon_{rad}$ is the matter-to-radiation conversion efficiency, which is assumed for simplicity to take the constant value $\epsilon_{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 $\Psi_{BHAR}(z)$ estimates. Nevertheless, we will discuss below the effect that a constant $k_{BOL} = 22.4$ would have on the estimated $\Psi_{BHAR}(z)$.

Figure 9 shows the $\Psi_{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 $\pm 1\sigma$ uncertainties by simply propagating the $\pm 1\sigma$ 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 $\Psi_{BHAR}(z)$ (orange downward arrows). As displayed in Figure 9, the MS population makes the bulk $\Psi_{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 $\Psi_{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
<figure><img src="image_11.png" /><figcaption>**Figure 9.** SMBH accretion rate density $\Psi_{BHAR}(z)$ since $z \sim 3$ (red) split between MS (green) and SB (blue) galaxies. The MS population makes the bulk $\Psi_{BHAR}$ at all cosmic epochs, while SBs are subdominant but evolve in a similar fashion. Downward arrows mark the upper limits on $\Psi_{BHAR}$ obtained from quiescent galaxies (Section 3.1.1). The overall BHARD agrees by design with the derivation by Aird et al. (2015, purple stars) and displays a shape similar to the SFRD (Madau & Dickinson 2014), here scaled down by 3300 to ease the comparison (gray dashed line). As displayed in the bottom panel, we find $f_{SB} \sim 20\%–30\%$ of the BHARD. Similarly, $f_{SB} \sim 10\%–15\%$ of the full SFRD at $z \sim 2$ and marginally consistent the local value (Sargent et al. 2012, magenta squares). Error bars are provided at the $1\sigma$ level.</figcaption></figure>

cosmic epochs (see also Gruppioni et al. 2013; Magnelli et al. 2013).

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 \times$ 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 $\approx 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 $\Psi_{BHAR}$, being $<4.3\%$ ($z = 0.5$), $<6.5\%$ ($z = 1$), $<1.8\%$ ($z = 2$), and $<1.6\%$ ($z = 3$). If the total $\Psi_{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).

### 3.4. Uncertainties and Evolution of Free Parameters

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_{EDD}$ distribution ($\alpha, \beta$), the break Eddington ratios of MS and SB galaxies ($\lambda_{MS}^*, \lambda_{SB}^*$), and the BHAR/SFR relation with $M_*$.

The combination of these parameters identifies the $\pm 1\sigma$ 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 $\alpha$ 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_{EDD}$ distribution at $z = 3$ (top panel of Figure 10). The bright-end slope $\beta$ is unconstrained from our analysis. Among the $\beta$ values initially allowed by our grid (2, 3, 4, and 5), the resulting $\pm 1\sigma$ confidence intervals range between 3 and 5 (see Table 2). Conversely, the flatter value of $\beta = 2$ is marginally (at $>1\sigma$) 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_{EDD}$ distributions are here assumed to contribute to the bright-end XLF, with their relative breaks being placed at different $L_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_{EDD}$ of MS galaxies, $\lambda_{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_{MS}^*$ with redshift, from $\lambda_{MS}^* = -0.7$ at $z = 0.5$ to $\lambda_{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):

$$ \lambda_{MS}^* = 10^{-0.9 \pm 0.4} \cdot (1 + z)^{1.12 \pm 0.19}. \quad (8) $$

Not only is $\lambda_{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_*$ density (from Aird et al. 2015 and Ilbert et al. 2013, respectively), from which they identified the corresponding knee $L_{AGN}$ and $M_*$. For the whole star-forming galaxy population (i.e., MS and SB), they obtained $\lambda^* \propto (1+z)^{2.5}$ out to $z \sim 2$, and constant at $z > 2$.

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

$$ \lambda_{SB}^* = 10^{0.0 \pm 1.0} \cdot (1 + z)^{0.80 \pm 0.57}. \quad (9) $$

As previously mentioned, our prior assumption that MS and SB galaxies share the same BHAR/SFR ratio, at each $M_*$ and redshift, produces a constant gap between $\lambda_{SB}^*$ and $\lambda_{MS}^*$ of the order of 0.8 dex (i.e., a factor of 6, Section 3.3). This gap mimics the $\times 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_{SB}^*$,
<table><caption>**Table 2**
Uncertainties and Confidence Ranges (at 1σ Level) of the Free Input Parameters Assumed in This Work (α, β, λ<sup>*</sup><sub>MS</sub>, λ<sup>*</sup><sub>SB</sub>, and the BHAR/SFR Relation with M<sub>*</sub>)</caption><thead><tr><th>Redshift<br>(1)</th><th>χ<sup>2</sup><sub>red</sub><br>(2)</th><th>α<br>(3)</th><th>β<br>(4)</th><th>log λ<sup>*</sup><sub>MS</sub><br>(5)</th><th>log λ<sup>*</sup><sub>SB</sub><br>(6)</th><th>Slope<br>BH/SF<br>(7)</th><th>Norm<br>BH/SF<br>(8)</th><th>L<sup>cross</sup><sub>X</sub><br>(9)</th><th>f(λ<sub>EDD,MS</sub> > 0.1)<br>(10)</th><th>f(λ<sub>EDD,SB</sub> > 0.1)<br>(11)</th></tr></thead><tbody><tr><td>z = 0.5</td><td>0.60</td><td>0.55</td><td>4<sup>+1</sup><sub>-0</sub><sup>**</sup></td><td>-0.7<sup>+0.1</sup><sub>-0.2</sub></td><td>0.0<sup>+0.5</sup><sub>-0.5</sub><sup>**</sup></td><td>0.95<sup>+0.00</sup><sub>-0.30</sub></td><td>-13.4<sup>+3.1</sup><sub>-0.0</sub></td><td>44.59<sup>+0.09</sup><sub>-0.23</sub></td><td>0.004<sup>+0.001</sup><sub>-0.002</sub></td><td>0.030<sup>+0.006</sup><sub>-0.011</sub></td></tr><tr><td>z = 1</td><td>0.29</td><td>0.22<sup>+0.06</sup><sub>-0.00</sub></td><td>5<sup>+0</sup><sub>-2</sub><sup>**</sup></td><td>-0.5<sup>+0.1</sup><sub>-0.3</sub></td><td>0.5<sup>+0.0</sup><sub>-1.0</sub><sup>**</sup></td><td>0.73<sup>+0.22</sup><sub>-0.24</sub></td><td>-11.2<sup>+2.5</sup><sub>-2.2</sub></td><td>44.74<sup>+0.09</sup><sub>-0.26</sub></td><td>0.016<sup>+0.007</sup><sub>-0.004</sub></td><td>0.038<sup>+0.019</sup><sub>-0.021</sub></td></tr><tr><td>z = 2</td><td>0.85</td><td>0.06<sup>+0.05</sup><sub>-0.00</sub></td><td>4<sup>+0</sup><sub>-1</sub></td><td>-0.4<sup>+0.1</sup><sub>-0.1</sub></td><td>0.4<sup>+0</sup><sub>-0.8</sub></td><td>0.73<sup>+0.22</sup><sub>-0.29</sub></td><td>-11.2<sup>+3.0</sup><sub>-2.2</sub></td><td>44.96<sup>+0.09</sup><sub>-0.12</sub></td><td>0.044<sup>+0.006</sup><sub>-0.015</sub></td><td>0.102<sup>+0.001</sup><sub>-0.056</sub></td></tr><tr><td>z = 3</td><td>1.32</td><td>0.03<sup>+0.02</sup><sub>-0.00</sub></td><td>5<sup>+0</sup><sub>-0</sub><sup>**</sup></td><td>-0.2<sup>+0.2</sup><sub>-0.0</sub></td><td>0.5<sup>+0.0</sup><sub>-0.3</sub><sup>**</sup></td><td>0.91<sup>+0.00</sup><sub>-0.01</sub></td><td>-13.0<sup>+0.2</sup><sub>-0.0</sub></td><td>45.15<sup>+0.31</sup><sub>-0.15</sub></td><td>0.065<sup>+0.002</sup><sub>-0.001</sub></td><td>0.153<sup>+0.001</sup><sub>-0.041</sub></td></tr></tbody></table>

**Note.** In addition, we report the crossover $L_X$ ($L_X^{cross}$) and the fraction of the $\lambda_{EDD}$ distribution spent above 10% Eddington ($f[\lambda_{EDD} > 0.1]$), for both MS and SB galaxies. The double vertical line separates free parameters (left) from other byproduct quantities (right). Zero lower or upper errors are due to our discrete parameter grid, and they should be interpreted as smaller than the closest value. (**) This symbol denotes those values touching the upper edge of our input grid, and should be taken as lower limits.

<figure><img src="image_12.png" /><figcaption>**Figure 10.** Top panel: redshift evolution of the faint-end slope of the $\lambda_{EDD}$ (α) in the form $\alpha \propto (1+z)^\gamma$. The best solution is given by $\gamma = -3.16^{+0.79}_{-0.00}$. Bottom panel: redshift evolution of the break $\lambda_{EDD}$ for MS and SB galaxies (green and blue points, respectively). We fit a power-law trend, finding $\lambda_{MS}^* \propto (1+z)^{1.12}$ (MS) and $\lambda_{SB}^* \propto (1+z)^{0.80}$ (SB). Uncertainties are given at the 1σ level. For comparison, the curve from Caplar et al. (2018) is shown (red hatched area), valid for AGN host galaxies close to the knee of the galaxy M* function and AGN bolometric LF, respectively. See Section 3.4 for details.</figcaption></figure>

empirical studies do support the condition $\lambda_{SB}^* > \lambda_{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.22}_{-0.29}$, 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^{-3}$, would lead to a higher than observed faint-end XLF, due to excessive mean BHAR values arising from low-mass galaxies.

## 4. Discussion

In this section, we try to interpret our results in a broader context of AGN and galaxy evolution. First, we validate our model predictions against known observed trends that have been reported in the literature (Section 4.1). Next, we discuss the broad implications of our findings within a twofold framework: (1) the role of cold gas content in driving SMBH accretion within MS and SB galaxies (Section 4.2), and (2) the AGN duty cycle as a function of MS offset and redshift (Section 4.3).

### 4.1. Testing Our Modeling against the Observed $L_X$-SFR Relation

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<sup>44</sup> erg s<sup>-1</sup>) 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<sup>44</sup> erg s<sup>-1</sup>) 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
<figure><img src="image_13.png" /><figcaption>**Figure 11.** Linear mean SFR 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$. Solid lines represent our predicted trends, while data points are from observational samples of X-ray-selected AGNs. The color bar indicates the redshift of each data set. We note that the SFR of each sample was properly scaled 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_\odot$). Our modeling shows a good agreement with the observed $L_X$–SFR trends at all redshifts. The “bump” predicted by our curves at high $L_X$ is driven by the gradual predominance of SB galaxies, with the position of the bump scaling with redshift (see Figure 12).</figcaption></figure>

$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 \sim 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 \sim 3.5$. The color bar indicates the redshift of each data set.

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_\odot$). 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 $\lambda_{\text{EDD}}$ distribution, which incorporates the stochasticity of SMBH accretion. However, we are able to circumvent this issue by constraining

<figure><img src="image_14.png" /><figcaption>**Figure 12.** Ratio between SB-related XLF and total (MS and SB) XLF, namely $f_{\text{SB}}(L_X)$, as a function of $L_X$ (solid lines) and redshift (different colors). Dashed lines indicate the corresponding $\pm 1\sigma$ confidence XLF propagated to $f_{\text{SB}}$. While the SB fraction accounts for just a few percent at $L_X < 10^{14.5} \text{ erg s}^{-1}$, it increases to nearly 100% at $L_X > 10^{45} \text{ erg s}^{-1}$, due to the gradual predominance of SB galaxies at the highest $L_X$. The black horizontal line marks $f_{\text{SB}} = 0.5$, which identifies the crossover $L_X^{\text{cross}}$.</figcaption></figure>

the main $\lambda_{\text{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.

## 4.2. What Drives the Evolving SMBH Growth in MS and SB Galaxies?

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^{10} < M_* < 10^{11} M_\odot$) MS galaxies, the SB population appears to take over at relatively high $L_X$, enabling us to reproduce the bright-end XLF since $z \sim 3$. Specifically, Figure 12 displays the fractional contribution of SB galaxies to the XLF ($f_{\text{SB}}$) as a function of $L_X$ and redshift. Solid lines mark the best $f_{\text{SB}}(L_X)$, with colors used to indicate our four redshift values. Dashed lines indicate the corresponding $\pm 1\sigma$ interval propagated from the XLF of SBs. We observe that $f_{\text{SB}}$ increases with $L_X$ from a few percent ($L_X < 10^{43.5} \text{ erg s}^{-1}$) to nearly 100% ($L_X > 10^{45} \text{ erg s}^{-1}$). As expected, the scatter becomes narrower toward higher $L_X$, where SBs increasingly dominate. The black horizontal line marks $f_{\text{SB}} = 0.5$, which identifies the crossover $L_X^{\text{cross}}$ (see top panel of Figure 13). This value scales as

$$L_X^{\text{cross}} (\text{erg s}^{-1}) = 10^{44.36 \pm 0.20} \cdot (1+z)^{1.28 \pm 0.33}. \quad (10)$$

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 $f_{\text{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.
<figure><img src="image_15.png" /><figcaption>**Figure 13.** Top panel: crossover $L_X$ ($L_X^{cross}$) as a function of redshift. The dashed line corresponds to a power-law fit as $L_X^{cross} \propto (1+z)^{1.28}$, suggesting a roughly linear increase. Bottom panel: fraction of the $\lambda_{EDD}$ distribution above 10% Eddington ($f[\lambda_{EDD} > 0.1]$), for both MS (green points) and SB (blue points) galaxies. Error bars are given at the $1\sigma$ level. Dashed lines indicate the power-law fit of both populations. For comparison, data taken from Grimmett et al. (2019) and Aird et al. (2019) are reported, based on more complex $\lambda_{EDD}$ functions. Though our $f[\lambda_{EDD} > 0.1]$ estimates stand slightly higher than their data, we consistently observe a monotonic increase with redshift, with SB galaxies displaying higher mean $\lambda_{EDD}$ (see Section 4.3 for details).</figcaption></figure>

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$\sigma$ 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 $\lambda_{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 $\lambda^*$ 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 $\lambda_{EDD}$ distribution support a progressive flattening at low $\lambda_{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 $\lambda_{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$\sigma$ significance.

Our results reveal a significant evolution of $\lambda_{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 $\lambda_{EDD}$ is higher toward earlier times, for both populations. A qualitatively similar cosmic evolution of the active AGN fraction and $\lambda_{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 \sim 2-3$ (as $f_{gas} \propto (1+z)^2$, 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 $\lambda_{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 ($\lambda_{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 $\lambda_{EDD}$ distribution.

### 4.3. AGN Duty Cycle in MS and SB Galaxies

The evolution of the characteristic $\lambda^*$ 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 $\lambda_{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[\lambda_{EDD} > 0.1]$) and how this evolves with redshift across MS and SB galaxies. We consider the $\lambda_{EDD}$ distribution corresponding to the best XLF at each redshift and $M_*$, separately for MS and SB, and integrate each function from $\lambda_{MAX} = 100$ (our maximum value) down to
$\lambda_{\text{EDD}} = 0.1$. Because the $\lambda_{\text{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 $\lambda_{\text{EDD}}$ as being proportional to the time spent above that $\lambda_{\text{EDD}}$. At each redshift, we weight all the $f[\lambda_{\text{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[\lambda_{\text{EDD}} > 0.1]$ at each redshift, separately for MS and SB galaxies. The bottom panel of Figure 13 shows $f[\lambda_{\text{EDD}} > 0.1]$ for both MS (green points) and SB (blue points) galaxies. Error bars are given at the $1\sigma$ 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:

$$f[\lambda_{\text{EDD,MS}} > 0.1] = 10^{-2.54 \pm 0.19} \cdot (1+z)^{2.26 \pm 0.40} \quad (11)$$

$$f[\lambda_{\text{EDD,SB}} > 0.1] = 10^{-1.82 \pm 0.57} \cdot (1+z)^{1.68 \pm 1.02}, \quad (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 $\lambda_{\text{EDD}}$ profiles. Though our $f[\lambda_{\text{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_{\text{EDD}}$.

Interestingly, our trend broadly follows the evolution of the molecular gas fraction $f_{\text{gas}}$ with redshift, namely $f_{\text{gas}} \propto (1+z)^2$ (e.g., Saintonge et al. 2013; Tacconi et al. 2013). We speculate that a link between $f_{\text{gas}}$ and typical AGN Eddington ratio $(\propto L_X/M_{\text{BH}} \sim L_X/M_*$ if assuming a fixed $M_{\text{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_{\text{dep}} = M_{\text{gas}}/\text{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; Cibelin 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.

### 4.4. Extrapolating the XLF out to $z \sim 5$

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

<figure><img src="image_16.png" /><figcaption>Figure 14. Best predicted XLFs at $0.5 < z < 3$ (solid lines) and their extrapolation out to $z \sim 5$ (dashed lines). Observational data from Vito et al. (2018) are shown for comparison at $z \sim 3.3$ (yellow circles) and $z \sim 4.1$ (red circles). See text for details.</figcaption></figure>

$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 \sim 3.3$ ($3 < z < 3.6$, yellow circles) and $z \sim 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 \sim 3.3$ and $z \sim 4.1$, with the corresponding $\lambda_{\text{EDD}}$ distributions. These latter are derived by extrapolating the redshift evolution of the slopes $(\alpha, \beta)$ and break $\lambda^*$ 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 \sim 3.3$ and $z \sim 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.

As discussed in Section 2.6, the behavior of the XLF out to $z \sim 5$ is driven by a progressive flattening of the faint-end and a positive shift of the bright-end $\lambda_{\text{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 $\pm 1\sigma$ confidence intervals) out to $z \sim 5$ in the online supplementary material, for both MS and SB galaxies.

### 4.5. From sBHAR to Intrinsic Eddington Ratio: Super-Eddington Growth?

We briefly explore the implications of relaxing the assumption of a constant $M_{\text{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,
<figure><img src="image_17.png" /><figcaption>**Figure 15.** Comparison between $\lambda_{\text{EDD}}$ distributions of MS galaxies at $M_{*} = 10^{10.5} M_{\odot}$. Solid lines mark the $p(\log \lambda_{\text{EDD}})$ used in this work, based on the assumption that $M_{\text{BH}}/M_{*} = 1/500$ (e.g., Häring & Rix 2004, see Figure 6). Dashed lines show the equivalent distribution shifted in the intrinsic $\lambda_{\text{EDD}}$ space by assuming a $M_{*}$-dependent $M_{\text{BH}}/M_{*}$ ratio for typical MS galaxies (Delvecchio et al. 2019), incorporating the scatter of the $M_{\text{BH}}-M_{*}$ conversion (shaded area). Colors mark our four redshift bins.</figcaption></figure>

recent studies support a $M_{*}$-increasing $M_{\text{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_{\text{BH}}-M_{*}$ scaling relation. In agreement with the above-mentioned literature, we found that in galaxies with $M_{*} \gtrsim 10^{10} M_{\odot}$ the $M_{\text{BH}}/M_{*}$ ratio increases with $M_{*}$ as $\log(M_{\text{BH}}/M_{*}) = -11.14 + 0.70 \times \log M_{*}$. In less massive galaxies, instead, the relation may flatten out depending on the assumed BH seed mass, which typically ranges between $10^2$ and $10^6 M_{\odot}$ (e.g., Begelman & Rees 1978).

In this work, we have shown that the bulk of SMBH growth occurs in MS galaxies with $10^{10} < M_{*} < 10^{11} M_{\odot}$ (Section 3.2 and Figure 8). In that range, the relation obtained in Delvecchio et al. (2019) yields $M_{\text{BH}}/M_{*} \approx 1/5000$. From Equation (2), this implies that the intrinsic Eddington ratio ($\lambda_{\text{EDD,INTR}}$) would shift up by a factor of 10. Figure 15 shows the comparison between the standard $\lambda_{\text{EDD}}$ distributions ($M_{\text{BH}}/M_{*} = 1/500$, solid lines) obtained in this work against the $\lambda_{\text{EDD,INTR}}$ distributions (dashed lines) for MS galaxies at different redshifts. The effect of adopting empirical $M_{\text{BH}}/M_{*}$ ratios is that the break $\lambda_{\text{EDD,INTR}}$ shifts above the Eddington limit. This trend is significant, given the relatively small uncertainties on our break $\lambda_{\text{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_{\text{EDD,INTR}}$ rises from 0.03 ($z=0.5$) to 0.19 ($z=3$), thus well below the Eddington limit.

This effect would be further amplified in $M_{*} < 10^{10} M_{\odot}$ galaxies, albeit it could be partly mitigated by assuming quite massive BH seeds ($M_{\text{BH,SEED}} \gtrsim 10^{5} M_{\odot}$) 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_{\text{BH}}/M_{*} = 1/500$, a lower radiative efficiency $\epsilon < 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_{\text{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 $\lambda_{\text{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_{\text{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_{\text{BH}}$ measurements) would be useful for observationally tying down the average break $\lambda_{\text{EDD,INTR}}$ in a statistical manner.

## 5. Summary and Conclusions

In this paper, we decipher the evolution of the AGN duty cycle in galaxies from the XLF, separating the contribution of MS and SB galaxies since $z \sim 3$. While these populations are known to display profound differences in structure and gas content, still open are the questions of whether the rate and incidence of SMBH accretion depend on MS offset and how they evolve over cosmic time. In order to account for the stochasticity of AGN activity and mitigate possible selection effects, we modeled the XLF as the convolution between the galaxy $M_{*}$ function and a large set of simulated $\lambda_{\text{EDD}}$ distributions, as done in a number of previous works (e.g., Aird et al. 2012; Bongiorno et al. 2012; Caplar et al. 2015; Jones et al. 2017; Weigel et al. 2017; Bernhard et al. 2018). In contrast to most studies, we assumed a very simple modeling, characterized by $M_{*}$-independent $\lambda_{\text{EDD}}$ parameters (slopes and break), normalized with $M_{*}$-dependent BHAR/SFR relations reported in the literature (Section 2). This allows us to derive a large set of predicted XLFs, separately between MS and SB galaxies, with a simple statistical approach and well anchored to empirical grounds.

Our analysis relies on three prior assumptions (Section 2.1): (1) the XLF is predominantly made by MS and SB galaxies, while passive systems have a negligible contribution (as later confirmed in Section 3.1.1); (2) we parameterize the $\lambda_{\text{EDD}}$ distribution as a broken power law with slopes ($\alpha, \beta$) that meet at the break $\lambda^{*}$; (3) the values of $\alpha$ and $\beta$ are assumed not to differ between MS and SB galaxies at fixed redshift.

The comparison between our model predictions and the observed XLF (Aird et al. 2015) reveals a very good agreement at all redshifts (Section 3.1), which leads us to the following main results.

(1) We reproduce the observed XLF through a continuous flattening of the faint-end $\lambda_{\text{EDD}}$ distribution, as well as a positive shift of the break $\lambda^{*}$ with redshift, consistent with previous studies (Caplar et al. 2015, 2018). Driven by our empirically motivated assumptions, SB galaxies stand above by a constant offset of 0.8 dex, reaching break $\lambda^{*}$
close to or slightly above the Eddington limit (Section 3.4 and Figure 10).

(2) By splitting the XLF into $M_*$ bins, we find that the bulk XLF is made by massive galaxies ($10^{10} < M_* < 10^{11} M_\odot$) on the MS, while merging-driven BH accretion in SB galaxies becomes dominant only in bright quasars with $L_X > 10^{44.36}(1+z)^{1.28}$ erg s⁻¹ (Figure 13). The inferred BHARD traced by the MS population shows a peak at $z \sim 2$ and declines at lower redshifts in a fashion similar to the SFRD (Madau & Dickinson 2014). Quiescent galaxies are estimated to contribute <6% of the integrated BHARD at each redshift (Section 3.3).

(3) We underline that a $M_*$-dependent relation between BHAR and SFR is strongly favored by our modeling and in line with recent studies (Aird et al. 2019). The best solution corresponds to BHAR/SFR $\propto M_*^{0.73[+0.22, -0.29]}$, while a constant BHAR/SFR trend is rejected at $\sim 3\sigma$ significance because it would overpredict the XLF at low $L_X$ arising from low-$M_*$ galaxies (see Bernhard et al. 2018). This finding implies that the cosmic buildup of SMBH and galaxy mass does not occur in lockstep at all epochs, but it evolves nonlinearly as the galaxy grows in $M_*$ (Section 4.2).

(4) Our modeling successfully reproduces the relatively flat $L_X$-SFR relation observed in X-ray-selected AGNs (Stanley et al. 2015) since $z \sim 3$. This bolsters the reliability of our approach in predicting realistic SFR estimates for X-ray AGNs across a wide $L_X$ and redshift range (Section 4.1).

(5) Finally, we argue that the probability of finding highly accreting ($\lambda_{EDD} > 10\%$) AGNs notably increases with redshift, from 0.4% (3.0%) at $z = 0.5\% - 6.5\%$ (15.3%) at $z = 3$ for MS (SB) galaxies (Figure 13), which supports a longer AGN duty cycle in the early universe, especially in dusty starbursting galaxies (Section 4.3). This is expectable if the level of SMBH accretion is tightly linked to the amount of usable cold gas in the host.

Our proposed framework serves as an important toy model for predicting the incidence of AGN activity in star-forming galaxies on and above the MS, the typical SFR and $M_*$ of X-ray AGNs, and the fraction of AGNs lying within MS and SB galaxies, at different luminosities and cosmic epochs. This modeling also opens potential questions about super-Eddington BH growth and different $\lambda_{EDD}$ prescriptions for explaining the assembly of the $M_{BH}-M_*$ relation.

Our key results broadly support a long-lasting interplay between SMBH accretion and star formation in galaxies, both showing enhanced activity at earlier epochs. This scenario is plausible if the evolution of cold gas content drives the triggering and maintenance of both phenomena over cosmic time. We speculate that merger-driven (or massive cold gas inflow-driven) SMBH accretion might be widespread in high-redshift ($z > 1$) SB galaxies, explaining the onset of Eddington-limited activity and the longer AGN duty cycle relative to MS analogs.

We thank the referee for his/her careful reading of the manuscript. I.D. is grateful to Iary Davidzon for useful input on the galaxy stellar mass function. I.D. is supported by the European Union's Horizon 2020 research and innovation program under the Marie Skłodowska-Curie grant agreement No. 788679. A.C. acknowledges the support from grant PRIN MIUR 2017. R.C. acknowledges financial support from CONICYT Doctorado

Nacional No. 21161487 and CONICYT PIA ACT172033. F.V. acknowledges financial support from CONICYT and CASSACA through the fourth call for tenders of the CAS-CONICYT Fund, and CONICYT grant Basal-CATA AFB-170002.

## ORCID iDs

I. Delvecchio https://orcid.org/0000-0001-8706-2252

E. Daddi https://orcid.org/0000-0002-3331-9590

J. Aird https://orcid.org/0000-0003-1908-8463

J. R. Mullaney https://orcid.org/0000-0002-3126-6712

A. Cimatti https://orcid.org/0000-0002-4409-5633

G. Zamorani https://orcid.org/0000-0002-2318-301X

N. Caplar https://orcid.org/0000-0003-3287-5250

F. Vito https://orcid.org/0000-0003-0680-9305

G. Rodighiero https://orcid.org/0000-0002-9415-2296

## References

Aird, J., Alexander, D. M., Ballantyne, D. R., et al. 2015, ApJ, 815, 66

Aird, J., Coil, A. L., & Georgakakis, A. 2017, MNRAS, 465, 3390

Aird, J., Coil, A. L., & Georgakakis, A. 2018, MNRAS, 474, 1225

Aird, J., Coil, A. L., & Georgakakis, A. 2019, MNRAS, 484, 4360

Aird, J., Coil, A. L., Moustakas, J., et al. 2012, ApJ, 746, 90

Aird, J., Coil, A. L., Moustakas, J., et al. 2013, ApJ, 775, 41

Alexander, D. M., & Hickox, R. C. 2012, NewAR, 56, 93

Aversa, R., Lapi, A., de Zotti, G., Shankar, F., & Danese, L. 2015, ApJ, 810, 74

Begelman, M. C., & Rees, M. J. 1978, MNRAS, 185, 847

Behroozi, P. S., Wechsler, R. H., & Conroy, C. 2013, ApJ, 770, 57

Benson, A. J., Bower, R. G., Frenk, C. S., et al. 2003, ApJ, 599, 38

Bernhard, E., Grimmett, L. P., Mullaney, J. R., et al. 2019, MNRAS, 483, L52

Bernhard, E., Mullaney, J. R., Aird, J., et al. 2018, MNRAS, 476, 436

Bernhard, E., Mullaney, J. R., Daddi, E., Ciesla, L., & Schreiber, C. 2016, MNRAS, 460, 902

Béthermin, M., Daddi, E., Magdis, G., et al. 2012, ApJL, 757, L23

Béthermin, M., Wu, H.-Y., Lagache, G., et al. 2017, A&A, 607, A89

Bongiorno, A., Merloni, A., Brusa, M., et al. 2012, MNRAS, 427, 3103

Bongiorno, A., Schulze, A., Merloni, A., et al. 2016, A&A, 588, A78

Bournaud, F., Dekel, A., Teyssier, R., et al. 2011, ApJL, 741, L33

Bower, R. G., Schaye, J., Frenk, C. S., et al. 2017, MNRAS, 465, 32

Calabrò, A., Daddi, E., Puglisi, A., et al. 2019, A&A, 623, A64

Caplar, N., Lilly, S. J., & Trakhtenbrot, B. 2015, ApJ, 811, 148

Caplar, N., Lilly, S. J., & Trakhtenbrot, B. 2018, ApJ, 867, 148

Carraro, R., et al. 2020, A&A, submitted

Chabrier, G. 2003, ApJL, 586, L133

Chen, C.-T. J., Hickox, R. C., Alberts, S., et al. 2013, ApJ, 773, 3

Cibinel, A., Daddi, E., Sargent, M. T., et al. 2019, MNRAS, 485, 5631

Civano, F., Fabbiano, G., Pellegrini, S., et al. 2014, ApJ, 790, 16

Conroy, C., & White, M. 2013, ApJ, 762, 70

Cowie, L. L., Barger, A. J., Bauer, F. E., & González-López, J. 2020, ApJ, 891, 69

Cresci, G., Mainieri, V., Brusa, M., et al. 2015, ApJ, 799, 82

Croton, D. J., Springel, V., White, S. D. M., et al. 2006, MNRAS, 365, 11

Daddi, E., Bournaud, F., Walter, F., et al. 2010, ApJ, 713, 686

Daddi, E., Cimatti, A., Renzini, A., et al. 2004, ApJ, 617, 746

Daddi, E., Elbaz, D., Walter, F., et al. 2010, ApJL, 714, L118

Davidzon, I., Ilbert, O., Laigle, C., et al. 2017, A&A, 605, A70

Delvecchio, I., Daddi, E., Shankar, F., et al. 2019, ApJL, 885, L36

Delvecchio, I., Gruppioni, C., Pozzi, F., et al. 2014, MNRAS, 439, 2736

Delvecchio, I., Lutz, D., Berta, S., et al. 2015, MNRAS, 449, 373

Di Matteo, T., Khandai, N., DeGraf, C., et al. 2012, ApJL, 745, L29

Dubois, Y., Peirani, S., Pichon, C., et al. 2016, MNRAS, 463, 3948

Dubois, Y., Volonteri, M., & Silk, J. 2014, MNRAS, 440, 1590

Duras, F., Bongiorno, A., Piconcelli, E., et al. 2017, A&A, 604, A67

Elbaz, D., Dickinson, M., Hwang, H. S., et al. 2011, A&A, 533, A119

Elmegreen, D. M., Elmegreen, B. G., Ravindranath, S., & Coe, D. A. 2007, ApJ, 658, 763

Fabian, A. C. 2012, ARA&A, 50, 455

Ferrarese, L., & Merritt, D. 2000, ApJL, 539, L9

Förster Schreiber, N. M., Genzel, R., Bouché, N., et al. 2009, ApJ, 706, 1364

Geach, J. E., Smail, I., Moran, S. M., et al. 2011, ApJL, 730, L19

Genzel, R., Tacconi, L. J., Gracia-Carpio, J., et al. 2010, MNRAS, 407, 2091
Georgakakis, A., Aird, J., Schulze, A., et al. 2017, MNRAS, 471, 1976

Georgakakis, A., Pérez-González, P. G., Fanidakis, N., et al. 2014, MNRAS, 440, 339

Glikman, E., Simmons, B., Mailly, M., et al. 2015, ApJ, 806, 218

Gobat, R., Daddi, E., Magdis, G., et al. 2018, NatAs, 2, 239

Goulding, A. D., Forman, W. R., Hickox, R. C., et al. 2014, ApJ, 783, 40

Grimmett, L. P., Mullaney, J. R., Jin, S., et al. 2019, MNRAS, 487, 4071

Gruppioni, C., Pozzi, F., Rodighiero, G., et al. 2013, MNRAS, 432, 23

Gültekin, K., Richstone, D. O., Gebhardt, K., et al. 2009, ApJ, 698, 198

Habouzit, M., Volonteri, M., & Dubois, Y. 2017, MNRAS, 468, 3935

Häring, N., & Rix, H.-W. 2004, ApJL, 604, L89

Heckman, T. M., & Best, P. N. 2014, ARA&A, 52, 589

Hickox, R. C., Jones, C., Forman, W. R., et al. 2009, ApJ, 696, 891

Hickox, R. C., Mullaney, J. R., Alexander, D. M., et al. 2014, ApJ, 782, 9

Hopkins, P. F., Hernquist, L., Cox, T. J., & Kereš, D. 2008, ApJS, 175, 356

Hopkins, P. F., Richards, G. T., & Hernquist, L. 2007, ApJ, 654, 731

Ilbert, O., McCracken, H. J., Le Fèvre, O., et al. 2013, A&A, 556, A55

Jones, M. L., Hickox, R. C., Mutch, S. J., et al. 2017, ApJ, 843, 125

Jones, M. L., Hickox, R. C., Mutch, S. J., et al. 2019, ApJ, 881, 110

Kawaguchi, T., Aoki, K., Ohta, K., & Collin, S. 2004, A&A, 420, L23

Kennicutt, R. C., Jr. 1998, ApJ, 498, 541

Kim, D.-W., & Fabbiano, G. 2013, ApJ, 776, 116

Kocevski, D. D., Faber, S. M., Mozena, M., et al. 2012, ApJ, 744, 148

Kormendy, J., & Ho, L. C. 2013, ARA&A, 51, 511

Laigle, C., McCracken, H. J., Ilbert, O., et al. 2016, ApJS, 224, 24

Le Fèvre, O., Cassata, P., Cucciati, O., et al. 2013, A&A, 559, A14

Lee, N., Sanders, D. B., Casey, C. M., et al. 2015, ApJ, 801, 80

Leroy, A. K., Walter, F., Brinks, E., et al. 2008, AJ, 136, 2782

Luo, B., Brandt, W. N., Xue, Y. Q., et al. 2017, ApJS, 228, 2

Lupi, A., Volonteri, M., Decarli, R., et al. 2019, MNRAS, 488, 4004

Lusso, E., Comastri, A., Simmons, B. D., et al. 2012, MNRAS, 425, 623

Lynden-Bell, D. 1969, Natur, 223, 690

Madau, P., & Dickinson, M. 2014, ARA&A, 52, 415

Magdis, G. E., Rigopoulou, D., Helou, G., et al. 2013, A&A, 558, A136

Magnelli, B., Popesso, P., Berta, S., et al. 2013, A&A, 553, A132

Magorrian, J., Tremaine, S., Richstone, D., et al. 1998, AJ, 115, 2285

Marconi, A., Risaliti, G., Gilli, R., et al. 2004, MNRAS, 351, 169

Merloni, A., & Heinz, S. 2008, MNRAS, 388, 1011

Morganti, R., Oosterloo, T. A., Emonts, B. H. C., et al. 2003, ApJL, 593, L69

Morganti, R., Tadhunter, C. N., & Oosterloo, T. A. 2005, A&A, 444, L9

Mullaney, J. R., Alexander, D. M., Aird, J., et al. 2015, MNRAS, 453, L83

Mullaney, J. R., Daddi, E., Béthermin, M., et al. 2012, ApJL, 753, L30

Netzer, H., Lani, C., Nordon, R., et al. 2016, ApJ, 819, 123

Noeske, K. G., Weiner, B. J., Faber, S. M., et al. 2007, ApJL, 660, L43

Olsen, K. P., Rasmussen, J., Toft, S., & Zirm, A. W. 2013, ApJ, 764, 4

Page, M. J., Symeonidis, M., Vieira, J. D., et al. 2012, Natur, 485, 213

Pellegrini, S. 2010, ApJ, 717, 640

Reines, A. E., & Volonteri, M. 2015, ApJ, 813, 82

Ricci, C., Trakhtenbrot, B., Koss, M. J., et al. 2017, Natur, 549, 488

Rodighiero, G., Brusa, M., Daddi, E., et al. 2015, ApJL, 800, L10

Rodighiero, G., Daddi, E., Baronchelli, I., et al. 2011, ApJL, 739, L40

Rodighiero, G., Enia, A., Delvecchio, I., et al. 2019, ApJL, 877, L38

Rosario, D. J., Santini, P., Lutz, D., et al. 2012, A&A, 545, A45

Rosario, D. J., Trakhtenbrot, B., Lutz, D., et al. 2013, A&A, 560, A72

Saintonge, A., Lutz, D., Genzel, R., et al. 2013, ApJ, 778, 2

Sanders, D. B., Soifer, B. T., Elias, J. H., et al. 1988, ApJ, 325, 74

Santini, P., Rosario, D. J., Shao, L., et al. 2012, A&A, 540, A109

Sargent, M. T., Béthermin, M., Daddi, E., & Elbaz, D. 2012, ApJL, 747, L31

Schawinski, K., Koss, M., Berney, S., & Sartori, L. F. 2015, MNRAS, 451, 2517

Schmidt, M. 1959, ApJ, 129, 243

Schmidt, M. 1963, Natur, 197, 1040

Scholtz, J., Alexander, D. M., Harrison, C. M., et al. 2018, MNRAS, 475, 1288

Schreiber, C., Pannella, M., Elbaz, D., et al. 2015, A&A, 575, A74

Schulze, A., Bongiorno, A., Gavignaud, I., et al. 2015, MNRAS, 447, 2085

Schulze, A., Silverman, J. D., Daddi, E., et al. 2019, MNRAS, 488, 1180

Scoville, N., Aussel, H., Brusa, M., et al. 2007, ApJS, 172, 1

Scoville, N., Lee, N., Vanden Bout, P., et al. 2017, ApJ, 837, 150

Shankar, F., Bernardi, M., Sheth, R. K., et al. 2016, MNRAS, 460, 3119

Shankar, F., Bernardi, M., Sheth, R. K., et al. 2019, MNRAS, 485, 1278

Shankar, F., Weinberg, D. H., & Miralda-Escudé, J. 2009, ApJ, 690, 20

Shankar, F., Weinberg, D. H., & Miralda-Escudé, J. 2013, MNRAS, 428, 421

Shao, L., Lutz, D., Nordon, R., et al. 2010, A&A, 518, L26

Shirakata, H., Kawaguchi, T., Oogi, T., Okamoto, T., & Nagashima, M. 2019, MNRAS, 487, 409

Soltan, A. 1982, MNRAS, 200, 115

Speagle, J. S., Steinhardt, C. L., Capak, P. L., & Silverman, J. D. 2014, ApJS, 214, 15

Stanley, F., Alexander, D. M., Harrison, C. M., et al. 2017, MNRAS, 472, 2221

Stanley, F., Harrison, C. M., Alexander, D. M., et al. 2015, MNRAS, 453, 591

Tacconi, L. J., Genzel, R., Neri, R., et al. 2010, Natur, 463, 781

Tacconi, L. J., Neri, R., Genzel, R., et al. 2013, ApJ, 768, 74

Ueda, Y., Akiyama, M., Hasinger, G., Miyaji, T., & Watson, M. G. 2014, ApJ, 786, 104

Vasudevan, R. V., & Fabian, A. C. 2007, MNRAS, 381, 1235

Veilleux, S., Kim, D.-C., & Sanders, D. B. 2002, ApJS, 143, 315

Vito, F., Brandt, W. N., Yang, G., et al. 2018, MNRAS, 473, 2378

Vito, F., Maiolino, R., Santini, P., et al. 2014, MNRAS, 441, 1059

Wang, T., Elbaz, D., Alexander, D. M., et al. 2017, A&A, 601, A63

Weigel, A. K., Schawinski, K., Caplar, N., et al. 2017, ApJ, 845, 134

Williams, R. J., Quadri, R. F., Franx, M., van Dokkum, P., & Labbé, I. 2009, ApJ, 691, 1879

Yang, G., Brandt, W. N., Vito, F., et al. 2018, MNRAS, 475, 1887