# AGN host galaxy mass function in COSMOS Is AGN feedback responsible for the mass-quenching of galaxies?

A. Bongiorno, A. Schulze, A. Merloni, G. Zamorani, O. Ilbert, F. La Franca, Y. Peng, E. Piconcelli, V. Mainieri, J. D. Silverman, et al.

► To cite this version:

A. Bongiorno, A. Schulze, A. Merloni, G. Zamorani, O. Ilbert, et al. AGN host galaxy mass function in COSMOS Is AGN feedback responsible for the mass-quenching of galaxies?. Astronomy and Astrophysics - A&A, 2016, 588, 10.1051/0004-6361/201527436. [hal-01434353](https://hal.science/hal-01434353)

HAL Id: hal-01434353

https://hal.science/hal-01434353

Submitted on 9 Sep 2021

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
# AGN host galaxy mass function in COSMOS

**Is AGN feedback responsible for the mass-quenching of galaxies?**

A. Bongiorno<fnref n="1" />, A. Schulze<fnref n="2" />, A. Merloni<fnref n="3" />, G. Zamorani<fnref n="4" />, O. Ilbert<fnref n="5" />, F. La Franca<fnref n="6" />, Y. Peng<fnref n="7" />, E. Piconcelli<fnref n="1" />,
V. Mainieri<fnref n="8" />, J. D. Silverman<fnref n="2" />, M. Brusa<fnref n="9" />,<fnref n="4" />, F. Fiore<fnref n="1" />, M. Salvato<fnref n="3" />, and N. Scoville<fnref n="10" />

Received 23 September 2015 / Accepted 5 January 2016

## ABSTRACT

We investigate the role of supermassive black holes in the global context of galaxy evolution by measuring the host galaxy stellar mass function (HGMF) and the specific accretion rate, that is, $\lambda_{SAR}$, the distribution function (SARDF), up to $z \sim 2.5$ with $\sim 1000$ X-ray selected AGN from XMM-COSMOS. Using a maximum likelihood approach, we jointly fit the stellar mass function and specific accretion rate distribution function, with the X-ray luminosity function as an additional constraint. Our best-fit model characterizes the SARDF as a double power-law with mass-dependent but redshift-independent break, whose low $\lambda_{SAR}$ slope flattens with increasing redshift while the normalization increases. This implies that for a given stellar mass, higher $\lambda_{SAR}$ objects have a peak in their space density at earlier epoch than the lower $\lambda_{SAR}$ objects, following and mimicking the well-known AGN cosmic downsizing as observed in the AGN luminosity function. The mass function of active galaxies is described by a Schechter function with an almost constant $M_{*}^{\star}$ and a low-mass slope $\alpha$ that flattens with redshift. Compared to the stellar mass function, we find that the HGMF has a similar shape and that up to $\log(M_{*}/M_{\odot}) \sim 11.5$, the ratio of AGN host galaxies to star-forming galaxies is basically constant ($\sim 10\%$). Finally, the comparison of the AGN HGMF for different luminosity and specific accretion rate subclasses with a previously published phenomenological model prediction for the “transient” population, which are galaxies in the process of being mass-quenched, reveals that low-luminosity AGN do not appear to be able to contribute significantly to the quenching and that at least at high masses, that is, $M_{*} > 10^{10.7} M_{\odot}$, feedback from luminous AGN ($\log L_{bol} \ge 46$ [erg/s]) may be responsible for the quenching of star formation in the host galaxy.

**Key words.** galaxies: active – galaxies: high-redshift – galaxies: evolution

## 1. Introduction

Supermassive black hole (SMBH) growth, nuclear activity, and galaxy evolution have been found to be closely related. In the past 15 years, the discovery of tight correlations between galaxies and their central nuclei properties (see Kormendy & Ho 2013, and references therein) as well as similar evolutionary trends between the growth histories of SMBHs and galaxies (e.g., Boyle & Terlevich 1998; Marconi et al. 2004) have established a new paradigm that shows active galactic nuclei (AGN) to be the key in the process of galaxy formation and evolution. Several theoretical models (e.g., Somerville et al. 2001; Granato et al. 2004; Monaco & Fontanot 2005; Springel et al. 2005; Croton et al. 2006; Hopkins et al. 2006; Schawinski et al. 2006; Cen & Chisari 2011) have been developed to explain this co-evolution and find the mechanism responsible for the simultaneous fuelling of the central BH and the formation of new stars in the host galaxy and the quasi-simultaneous shut-off of these

two processes. While the physical scales of interest (a few pc) cannot be directly resolved in these models and in current numerical simulations (e.g., Sijacki et al. 2015), the models usually propose the presence of an energetic AGN-driven feedback, which is a strong wind originating from the AGN that deposits the energy released by the accretion process within the host galaxy (Faucher-Giguère & Quataert 2012). This mechanism is able to link black hole growth and star formation and the shut-off of these two processes in a self-regulated manner. However, it is still unclear and observationally not proven whether AGN-driven feedback processes do indeed have an effect on the global properties of the galaxy population, in particular in suppressing the star formation (SF) in their host galaxy heating and/or pushing away the gas that forms stars.

Star formation quenching by some mechanism is also required to prevent the overgrowth of massive galaxies that are hosted in the most massive dark-matter haloes (e.g., Read & Trentham 2005). This “mass quenching” mechanism, regardless

<footnote n="1">¹ INAF-Osservatorio Astronomico di Roma, via Frascati 33, 00040 Monteporzio Catone, Italy
e-mail: angela.bongiorno@oa-roma.inaf.it (OAR)</footnote>

<footnote n="2">² Kavli Institute for the Physics and Mathematics of the Universe, Todai Institutes for Advanced Study, the University of Tokyo, 277-8583 Kashiwa, Japan (Kavli IPMU, WPI)</footnote>

<footnote n="3">³ Max-Planck-Institut fuer Extraterrestrische Physik (MPE), Postfach 1312, 85741 Garching, Germany</footnote>

<footnote n="4">⁴ INAF-Osservatorio Astronomico di Bologna, via Ranzani 1, 40127 Bologna, Italy</footnote>

<footnote n="5">⁵ Aix-Marseille Université, CNRS, LAM (Laboratoire d'Astrophysique de Marseille) UMR 7326, 13388 Marseille, France</footnote>

<footnote n="6">⁶ Dipartimento di Matematica e Fisica, Università Roma Tre, via della Vasca Navale 84, 00146 Roma, Italy</footnote>

<footnote n="7">⁷ Cavendish Laboratory, University of Cambridge, 19 J. J. Thomson Ave., Cambridge CB3 0HE, UK</footnote>

<footnote n="8">⁸ European Southern Observatory, Karl-Schwarzschild-Str. 2, 85748 Garching bei Muenchen, Germany</footnote>

<footnote n="9">⁹ Dipartimento di Fisica e Astronomia, Università di Bologna, viale Berti Pichat 6/2, 40127 Bologna, Italy</footnote>

<footnote n="10">¹⁰ California Institute of Technology, MC 249-17, 1200 East California Boulevard, Pasadena, CA 91125, USA</footnote>
of its physical origin, would suppress the growth of massive galaxies and explain the steep decline of the galaxy mass function above a given characteristic mass. While supernova feedback is not energetic enough in this mass regime, a central AGN would be an efficient mechanism.

To investigate whether AGN could play this role, detailed studies on single objects have been performed to search for signatures of AGN feedback. Massive outflows on several kpc scales have been observed in a few cases (Cano-Díaz et al. 2012; Feruglio et al. 2010; Cresci et al. 2015a; Feruglio et al. 2015), but up to now the evidence that such outflows are indeed responsible for suppressing star formation in the region of the outflow is circumstantial (Cano-Díaz et al. 2012; Cresci et al. 2015a,b). Further progress can be made through statistical studies of the properties of active galaxies (e.g., SFR) compared to normal galaxies. However, results have often been contradictory: some authors found that AGN mainly lie above or on the main sequence (MS) of galaxies (Santini et al. 2012; Mullaney et al. 2012), while others (Bongiorno et al. 2012; Mullaney et al. 2015) found the SFR of AGN hosts to be lower than the average MS galaxies, as expected from the models that include AGN feedback. Bundy et al. (2008) compared the star formation quenching rate with the rate at which AGN activity is triggered in galaxies and showed that these two quantities agree over a range of masses. They interpreted this as a physical link between these two phenomena, which does not directly imply a causal link, however.

Regardless of AGN feedback, an essential pre-requisite for understanding the role of black hole activity in galaxy evolution is to have an accurate and unbiased census of the AGN population and its relation to the properties of their host galaxies. The former is basically provided by the AGN luminosity function, which is now well established over a wide range of redshift and luminosity (Ueda et al. 2014; Buchner et al. 2015; Aird et al. 2015; Miyaji et al. 2015; Silverman et al. 2008). Deep X-ray surveys established a trend of AGN downsizing, which means that the most luminous AGN have the peak in their space density at earlier times than lower luminosity AGN (Ueda et al. 2003; Hasinger et al. 2005), which is also seen in optical surveys (Bongiorno et al. 2007; Croom et al. 2009). This trend is similar to the downsizing in the galaxy population (Cowie et al. 1996), where the most massive galaxies build their mass at earlier times than lower mass galaxies.

Linking black hole growth to their host galaxies requires studying their stellar mass function and/or the active fraction or duty cycle of AGN occurrence in galaxies of a given stellar mass (e.g., Bundy et al. 2008; Xue et al. 2010; Georgakakis et al. 2011; Aird et al. 2012; Bongiorno et al. 2012; Lusso et al. 2012). Most of these studies define AGN activity above a certain X-ray luminosity threshold and found the fraction of AGN at given $L_X$ to increase with stellar mass. However, this may lead to a biased view, since AGN at different masses cover different ranges of Eddington ratios for a given luminosity range, and AGN have been found to show a wide distribution of Eddington ratios (e.g., Kauffmann & Heckman 2009; Schulze & Wisotzki 2010). Aird et al. (2012) showed that the intrinsic distribution of specific accretion rates at $z < 1$ follows a power law whose shape does not evolve with redshift, independent of stellar mass. This result has been confirmed and extended out to $z < 2.5$ by Bongiorno et al. (2012).

We here build upon these previous studies of AGN hosts by establishing the bivariate distribution function of stellar mass and specific accretion rate for a hard X-ray selected AGN sample over the redshift range $0.3 < z < 2.5$. We use the derived

AGN host galaxy stellar mass function to test the hypothesis of AGN feedback as driver of star formation quenching. In particular, we test whether the AGN population can be associated with and/or be responsible for mass quenching using the model prediction from Peng et al. (2010) for the mass function of the transient population (these are galaxies in the process of being mass-quenched).

The paper is organized as follows: in Sect. 2 we present the X-ray selected sample. Section 3 presents the method used to derive the specific accretion rate distribution function and the AGN host galaxy mass function (Sects. 3.2 and 3.3) and their results (Sect. 3.4). In Sect. 4 we address the question of the link between AGN and star formation quenching by comparing the AGN host galaxy mass function, computed for different samples, with the model prediction for quenching galaxies by Peng et al. (2010).

Throughout this paper, a standard cosmology ($\Omega_m = 0.3$, $\Omega_\lambda = 0.7$ and $H_0 = 70 \text{ km s}^{-1} \text{ Mpc}^{-1}$) has been assumed. The stellar masses are given in units of solar masses for a Chabrier IMF (Chabrier 2003).

## 2. Sample

The AGN sample considered here was extracted from the XMM-COSMOS point-like source catalog (Hasinger et al. 2007; Cappelluti et al. 2009), whose optical identifications and multi-wavelength properties were presented by Brusa et al. (2010). The catalog contains ~1800 X-ray sources detected above flux limits of ~$5 \times 10^{-16}$, ~$3 \times 10^{-15}$ and ~$7 \times 10^{-15}$ erg cm$^{-2}$ s$^{-1}$ in the [0.5–2] keV, [2–10] keV and [5–10] keV bands, respectively.

Our analysis is based on objects that have been detected in the hard [2–10] keV band. The restriction to a hard X-ray selected sample was chosen since the soft band can be affected by obscuration that can lead to a redshift-dependent incompleteness (flux-limited surveys pick up more obscured objects at higher redshift, see, e.g., Gilli et al. 2010). However, this band may still suffer from incompleteness that is due to heavily obscured and Compton-thick (CT, $\log(N_H) > 24$ [cm$^{-2}$]) AGN, whose detection probability is strongly reduced because the intrinsic emission can be significantly suppressed by repeated Compton scattering and photoelectric absorption.

Out of the full ~1800 sources, we identified a final sample of 927 hard X-ray selected AGN in the redshift range $0.3 < z < 2.5$. All hard X-ray sources have accurate photometric redshifts (Salvato et al. 2011), while about half (581 of 927) have secure spectroscopic redshifts.

### 2.1. X-ray luminosities, host galaxy stellar masses, and specific accretion rates

Rest-frame, intrinsic X-ray [2–10] keV luminosities for the final sample were derived from the observed hard X-ray flux. Following La Franca et al. (2005), we converted the observed [2–10] keV fluxes to the intrinsic [2–10] keV luminosities for each AGN with a given measured $N_H$ by applying a K-correction computed by assuming an intrinsic X-ray spectrum with a photon index $\Gamma = 1.8$, an exponential cut-off at $E = 200$ keV and, a photoelectric absorption corresponding to the observed $N_H$ column density. The [2–10] keV luminosity is given by

$$L_{[2-10] \text{ keV}}^{\text{rf}} = F_{[2-10] \text{ keV}} 4\pi D_L^2 K(z, N_H), \quad (1)$$

where $D_L$ is the luminosity distance and $K(z, N_H)$ is the term that accounts for the K-correction and absorption correction. The absorbing column density $N_H$ for our sample was derived as in
<figure><img src="image_2.png" /><figcaption>**Fig. 1.** Bivariate distribution for the analyzed hard X-ray selected sample in the $M_{\star} - \lambda_{SAR}$ plane color-coded according to the redshift bins. The horizontal and vertical dashed lines correspond to the lower limit cuts applied in $M_{\star}$ and $\lambda_{SAR}$.</figcaption></figure>

Merloni et al. (2014). For the brightest sources (above 200 pc counts in the 0.5–10 keV band of XMM-Newton) $N_H$ was obtained from the full spectral analysis of Mainieri et al. (2011), which is available for 195 of 927 of the AGN. For the remaining sources, $N_H$ was estimated in a statistical fashion by assessing the value of the observed spectral slope from the hardness ratio and assessing the value of the intrinsic spectral slope drawn from a normal distribution with mean and dispersion of $\Gamma_{int} = 1.8 \pm 0.2$. While this estimate shows a significant scatter, there are no apparent systematic biases, as demonstrated in Merloni et al. (2014). Therefore these estimates can be robustly used for the statistical studies performed in this paper.

Host galaxy stellar masses have been derived in Bongiorno et al. (2012) using a two-component (AGN and galaxy) SED-fitting technique. We refer to this paper for a detailed description of the method.

Following Aird et al. (2012) and Bongiorno et al. (2012), we define the specific accretion rate $\lambda_{SAR} \equiv L_X/M_\star$ (see also Brusa et al. 2009; Georgakakis et al. 2014) as a directly measurable quantity that can be regarded as a proxy for the black hole growth rate relative to the stellar mass of the host galaxy, $M_{BH}/M_\star$, after taking into account the (luminosity-dependent) bolometric correction (e.g., Marconi et al. 2004; Lusso et al. 2012) and a radiative efficiency factor. It is also related to the Eddington ratio of the SMBH, $\lambda_{Edd} = L_{bol}/M_{BH}$, by applying the bolometric correction factor and the scaling relationship between black hole mass and host stellar mass. Assuming as an approximation a mean bolometric correction $k_{bol} = 25$ (Marconi et al. 2004; Lusso et al. 2012) and a constant mass ratio of the stellar host to the black hole of 500 (Marconi & Hunt 2003; Häring & Rix 2004), $\log \lambda_{SAR} = 34$ [erg/s/M$_\odot$] approximately corresponds to the Eddington limit, while $\log \lambda_{SAR} = 32$ [erg/s/M$_\odot$] would give 1% of the Eddington limit. The bivariate distribution $M_\star - \lambda_{SAR}$ for the analyzed sample is shown in Fig. 1, where different colors correspond to different redshift ranges as labeled.

To determine the mass function, we furthermore restricted our sample in stellar mass $M_\star$ and specific accretion rate $\lambda_{SAR}$ by applying the following cuts: $M_\star > 10^{9.5} M_\odot$ and $\lambda_{SAR} > 10^{32} \text{ erg/s/M}_\odot$. The latter criterion is motivated by the

requirement of having a clear cut in $\lambda_{SAR}$ above which we define the AGN as active (see below). The chosen minimum $\lambda_{SAR}$ value corresponds to the lowest observed value in our intermediate-redshift bin and furthermore corresponds approximately to 1% of the Eddington limit, which we chose in the following as our lowest threshold to define an active black hole, consistent with studies of type 1 AGN (Schulze et al. 2015). After applying these limits, our sample is reduced to 877 AGN with $0.3 < z < 2.5$.

### 3. AGN host galaxy mass function and specific accretion rate distribution function

To derive the AGN host galaxy mass function (HGMF) and the specific accretion rate distribution function (SARDF), we have to account for various selection effects in our flux-limited AGN sample. This requires a careful assessment of the incompleteness function.

Completeness in $L_X$ does not directly ensure completeness in $M_\star$. As previously reported, AGN show a wide range of Eddington ratios (Kauffmann & Heckman 2009; Schulze & Wisotzki 2010), and thus also a wide range of $L_X/M_\star$ ($\lambda_{SAR}$), with a distribution falling below the corresponding Eddington limit approximately following a power-law distribution (Aird et al. 2012; Bongiorno et al. 2012).

A luminosity-complete AGN sample will be biased toward high-mass BHs and high galaxy mass. This means that because an AGN with a low Eddington ratio will be included in the sample only if its $M_{BH}$ is high enough to be above the given luminosity ($L_X$) limit, a bias toward high-mass black holes induces a bias toward high-mass galaxies because of the relation between $M_{BH} - M_\star$. This effect has to be carefully taken into account when building a galaxy-mass-complete sample starting from an X-ray flux-limited AGN sample.

#### 3.1. Incompleteness function

Our corrections for incompleteness account for three effects: (1) the X-ray sensitivity function; (2) the absorption correction $f(N_H | L_X, z)$; and (3) the stellar mass completeness down to our threshold in units of specific accretion rate $\log \lambda_{SAR} = 32$ [erg/s/M$_\odot$].

The first selection effect to consider is the position-dependent X-ray flux limit based on the sensitivity map computed by Cappelluti et al. (2009). The absorption correction accounts for the sources that were missed in the sample because of their high column density $N_H$. For this correction we used the $N_H$ distribution as a function of $z$ and $L_X$ that was published by Ueda et al. (2014) based on several X-ray AGN surveys (see their Eqs. (5) and (6)). We integrated over the $N_H$ distribution between $20 < \log N_H < 24$, i.e., we did not include Compton-thick AGN in our HGMF determination. The fraction of CT AGN is still uncertain and the $N_H$ distribution above $\log N_H = 24$ is poorly known (Ueda et al. 2014; Buchner et al. 2015; Aird et al. 2015). The contribution of CT AGN to the AGN space density is expected to lie between ~10–40% (Gilli et al. 2007; Treister et al. 2009; Vignali et al. 2014; Buchner et al. 2015; Lansbury et al. 2015). These two corrections applied to the flux-limited sample result in a luminosity-complete sample.

As described above, the sample is also significantly incomplete because a broad range of $M_\star$ can be associated with a given luminosity $L_X$. To account for this effect in the HGMF, we need to include an additional term to the incompleteness function based on the distribution of $\lambda_{SAR}$. Using this distribution
function, we corrected for incompleteness down to a fixed threshold in $\lambda_{SAR}$, which we set at $\log \lambda_{SAR} = 32$ [erg/s/$M_{\odot}$]. The HGMF is therefore defined as the mass function of all AGN above this $\lambda_{SAR}$ threshold. The most rigorous and self-consistent approach to do this is by determining the HGMF and the SARDF simultaneously, for instance, with the maximum likelihood method described in the next section.

## 3.2. Maximum likelihood method

We here present the method of determining the SARDF and the HGMF simultaneously as a bivariate distribution function of stellar mass and specific accretion rate, $\Psi(M_{\star}, \lambda_{SAR}, z)$, where $\Psi(M_{\star}, \lambda_{SAR}, z)d\log M_{\star}d\log \lambda_{SAR}$ gives the space density of AGN with stellar mass host galaxies between $\log M_{\star}$ and $\log M_{\star} + d\log M_{\star}$ and a specific accretion rate between $\log \lambda_{SAR}$ and $\log \lambda_{SAR} + d\log \lambda_{SAR}$ at the redshift $z$. The HGMF, SARDF, and the X-ray AGN LF (XLF) can be derived as different marginalizations over this bivariate distribution function. We used the maximum likelihood method developed by Schulze & Wisotzki (2010) and extended by Schulze et al. (2015) to compute $\Psi(M_{\star}, \lambda_{SAR}, z)$. While these works focused on the joint determination of the active black hole mass function and the Eddington ratio distribution function (using type 1 AGN), the method was implemented here for the joint determination of the HGMF and SARDF.

The technique minimizes the likelihood function $S = -2 \sum \ln p_i$, where the probability distribution $p_i$ for each object is given by

$$p_i(M_{\star}, \lambda_{SAR}, N_{\mathrm{H}}, z) = \frac{1}{N} \Psi(M_{\star}, \lambda_{SAR}, z) \mathcal{I}(M_{\star}, \lambda_{SAR}, z, N_{\mathrm{H}}) \times f(N_{\mathrm{H}} | L_{\mathrm{X}}, z) \frac{dV}{dz}, \quad (2)$$

where $\Psi(M_{\star}, \lambda_{SAR}, z)$ is the bivariate distribution function of stellar mass and specific accretion rate that we wish to derive, $\mathcal{I}(M_{\star}, \lambda_{SAR}, z, N_{\mathrm{H}}) = \mathcal{I}(L_{\mathrm{X}}, z, N_{\mathrm{H}})$ is the X-ray selection function given by the sensitivity map in the 2–10 keV band, and $f(N_{\mathrm{H}} | L_{\mathrm{X}}, z)$ is the absorption distribution function, taken from Ueda et al. (2014). We used the $N_{\mathrm{H}}$ estimates presented in Sect. 2.1 to compute $L_{\mathrm{X}}$ (and therefore $\lambda_{SAR}$) and $f(N_{\mathrm{H}} | L_{\mathrm{X}}, z)$ for our sample. The factor $N$ corresponds to the total number of objects in the sample predicted by the model, and it is given by integrating over $M_{\star}$, $\lambda_{SAR}$, $N_{\mathrm{H}}$, and $z$,

$$N = \iiint\int \Psi(M_{\star}, \lambda_{SAR}, z) \mathcal{I}(M_{\star}, \lambda_{SAR}, z, N_{\mathrm{H}}) \times f(N_{\mathrm{H}} | L_{\mathrm{X}}, z) \frac{dV}{dz} d\log N_{\mathrm{H}} d\log \lambda_{SAR} d\log M_{\star} dz, \quad (3)$$

where we integrated over the $N_{\mathrm{H}}$ distribution between $20 < \log N_{\mathrm{H}} < 24$, while our integration ranges in $M_{\star}$, $\lambda_{SAR}$ and $z$ are $9.5 < \log M_{\star} < \infty$, $32 < \log \lambda_{SAR} < \infty$ and $0.3 < z < 2.5$, as discussed in Sect. 2.1.

Our sample also contains 12 AGN without $M_{\star}$ measurements because of poor-quality photometry. However, we accounted for these sources using their luminosity and redshift information integrated over the entire mass range,

$$p_j(L_{\mathrm{X}}, N_{\mathrm{H}}, z) = \int p_j(M_{\star}, \lambda_{\mathrm{SAR}}, N_{\mathrm{H}}, z) d\log M_{\star}. \quad (4)$$

Our XMM-COSMOS based sample covers only a limited dynamical range in $L_{\mathrm{X}}$, narrower than the full range over which

the XLF is currently determined. This might lead to degenerate solutions for the bivariate distribution function, some of which may be inconsistent with the XLF. Ideally, we would like to construct the HGMF and SARDF including deeper and larger area surveys, but this is beyond the scope of the present work. To reduce this effect, we included the XLF as additional observational data. In this way, we ensured consistency with the XLF observations over its full observationally determined luminosity range. In particular, we used the binned XLF from Miyaji et al. (2015) and computed the $\chi^2$ value for the comparison with the XLF implied by the HGMF and SARDF. We then added this likelihood to that of the XMM-COSMOS sample. The study by Miyaji et al. (2015) used the same $N_{\mathrm{H}}$ distribution as Ueda et al. (2014) to determine the XLF as we employed here. Over our range in redshift and luminosity, the XLF by Miyaji et al. (2015) is consistent with other recent studies (Ueda et al. 2014; Buchner et al. 2015; Aird et al. 2015), thus our results are robust against the specific choice of XLF.

We caution that the faint end of the XLF is not directly constrained by our sample: the XLF will also include AGN below our threshold in $M_{\star}$ and $\lambda_{SAR}$, which are not accounted for in our bivariate distribution function fit. This may lead to an overestimate of the space density at $\log L_{\mathrm{X}} < 43$ [erg/s].

The total likelihood to minimize is given by

$$S_{\mathrm{tot}} = -2 \sum_{i=1}^{N_{M_{\star}}} \ln p_i(M_{\star}, \lambda_{\mathrm{SAR}}, N_{\mathrm{H}}, z) - 2 \sum_{j=1}^{N_{L_{\mathrm{X}}}} \ln p_j(L_{\mathrm{X}}, N_{\mathrm{H}}, z) + \chi^2(\mathrm{XLF}), \quad (5)$$

where $N_{M_{\star}}$ is the number of AGN with $M_{\star}$ measurements in our sample and $N_{L_{\mathrm{X}}}$ is the number of AGN with only $L_{\mathrm{X}}$ known. The absolute normalization of the bivariate distribution function is then determined by scaling to the total observed number of objects in the sample.

Following Aird et al. (2012) and Bongiorno et al. (2012), we first assumed that the bivariate distribution function $\Psi(M_{\star}, \lambda_{SAR}, z)$ is separable, meaning that the specific accretion rate distribution is mass independent and vice versa. Under this assumption, the bivariate distribution function is given by

$$\Psi(M_{\star}, \lambda_{\mathrm{SAR}}, z) = \Psi^* f_{\lambda_{\mathrm{SAR}}} (\lambda_{\mathrm{SAR}}, z) f_{\star}(M_{\star}, z) f_z(z), \quad (6)$$

where $\Psi^*$ is the normalization of the bivariate distribution function, $f_{\lambda_{\mathrm{SAR}}}(\lambda_{\mathrm{SAR}}, z)$ is the SAR term, $f_{\star}(M_{\star}, z)$ is the $M_{\star}$-term and $f_z(z)$ is a redshift evolution term.

However, for the SAR term, we also tested a mass-dependent model and found this model to provide a better description of our data (see Appendix A for more details). The bivariate distribution function is therefore written as

$$\Psi(M_{\star}, \lambda_{\mathrm{SAR}}, z) = \Psi^* f_{\lambda_{\mathrm{SAR}}} (\lambda_{\mathrm{SAR}}, M_{\star}, z) f_{\star}(M_{\star}, z) f_z(z), \quad (7)$$

where $f_{\lambda_{\mathrm{SAR}}}(\lambda_{\mathrm{SAR}}, M_{\star}, z)$ now also contains a dependence on the mass. We used this more general parametrization as our default model. We point out that the SAR term $f_{\lambda_{\mathrm{SAR}}}$ and the $M_{\star}$ term $f_{\star}$ are not equal to the SARDF and HGMF.

The HGMF and the SARDF are calculated by integrating $\Psi(M_{\star}, \lambda_{SAR}, z)$ over $\lambda_{SAR}$ and over $M_{\star}$, respectively. To be specific,

$$\begin{aligned} \Phi_{\star}(M_{\star}, z) &= \frac{dN}{dV d\log M} = \int_{32}^{\infty} \Psi(M_{\star}, \lambda_{\mathrm{SAR}}, z) d\log \lambda_{\mathrm{SAR}} \\ &= \int_{32}^{\infty} \Psi^* f_{\lambda_{\mathrm{SAR}}} (\lambda_{\mathrm{SAR}}, M_{\star}, z) f_{\star}(M_{\star}, z) f_z(z) d\log \lambda_{\mathrm{SAR}} \end{aligned} \quad (8)$$
<figure><img src="image_3.png" /><figcaption>**Fig. 2.** Bivariate distribution function $\Psi(M_*, \lambda_{SAR}, z)$, for our best-fitting parametric model, derived through the maximum likelihood method, in three redshift bins. The orange contours indicate lines of constant space density (from $10^{-10}$ to $10^{-3}$ Mpc$^{-3}$), separated by a factor of 10 each.</figcaption></figure>

<table><caption>**Table 1.** Best-fit model parameters and their errors for the bivariate distribution function of stellar mass and SAR (Eq. (7)).</caption><thead><tr><th rowspan="2">log(Ψ*)</th><th colspan="2">f<sub>*</sub>(M<sub>*</sub>, z) from Eq. (10)</th><th colspan="5">f<sub>λSAR</sub>(λ<sub>SAR</sub>) from Eq. (11)</th><th colspan="3">f<sub>z</sub>(z) from Eq. (12)</th></tr><tr><th>log M<sub>*</sub><sup>*</sup> [M<sub>☉</sub>]</th><th>α</th><th>log λ<sub>SAR,0</sub><sup>*</sup></th><th>k<sub>1</sub></th><th>log M<sub>*,0</sub></th><th>γ<sub>1,0</sub></th><th>k<sub>γ</sub></th><th>γ<sub>2</sub></th><th>p<sub>1</sub></th><th>p<sub>2</sub></th><th>z<sub>0</sub></th></tr></thead><tbody><tr><td>-6.86</td><td>10.99</td><td>0.24</td><td>33.8<sup>*</sup></td><td>-0.48</td><td>11.0<sup>*</sup></td><td>-1.01</td><td>0.58</td><td>-3.72</td><td>5.82</td><td>2.36</td><td>1.1<sup>*</sup></td></tr><tr><td>+/-0.01</td><td>+/-0.03</td><td>+0.07<br/>-0.06</td><td>—</td><td>+0.03<br/>-0.03</td><td>—</td><td>+0.02<br/>-0.02</td><td>+0.02<br/>-0.02</td><td>+0.09<br/>-0.09</td><td>+0.12<br/>-0.13</td><td>+0.08<br/>-0.08</td><td>—</td></tr></tbody></table>

**Notes.** The parameters denoted with an asterisk are kept fixed during the fit.

and

$$ \Phi_{\lambda_{SAR}}(\lambda_{SAR}, z) = \frac{dN}{dV \, d\log \lambda} = \int_{9.5}^{\infty} \Psi(M_{*}, \lambda_{SAR}, z) \, d\log M_{*} \\ = \int_{9.5}^{\infty} \Psi^{*} f_{\lambda_{SAR}}(\lambda_{SAR}, M_{*}, z) f_{*}(M_{*}, z) f_{z}(z) \, d\log M_{*}. \quad (9) $$

In case of separable SAR- and $M_*$-terms, as in Eq. (6), the SARDF (HGMF) has the same shape as $f_{\lambda_{SAR}}(f_*)$ and only the absolute normalization is determined by the marginalization. However, in the more general case of Eq. (7), this is not necessarily the case, which is why the HGMF and SARDF then cannot be explicitly expressed as analytic functions.

We here considered the following parametric models for the individual terms: the $M_*$ term is modeled using a Schechter function:

$$ f_{*}(M_{*}, z) = \left(\frac{M_{*}}{M_{*}^{*}}\right)^{\alpha} e^{-\left(\frac{M_{*}}{M_{*}^{*}}\right)}. \quad (10) $$

While a model with a low-mass slope $\alpha$ evolving with redshift was included, we found that the best-fit parameters are indeed consistent with no z-evolution in $\alpha$.

The SAR term is instead described by a double power-law:

$$ f_{\lambda_{SAR}}(\lambda_{SAR}, M_{*}, z) = \frac{1}{\left(\frac{\lambda_{SAR}}{\lambda_{SAR}^{*}(M_{*})}\right)^{-\gamma_{1}(z)} + \left(\frac{\lambda_{SAR}}{\lambda_{SAR}^{*}(M_{*})}\right)^{-\gamma_{2}}}, \quad (11) $$

where the low $\lambda_{SAR}$ slope $\gamma_1(z) = \gamma_{1,0} + k_\gamma(z - z_0)$, with $z_0$ set at 1.1, and the break log $\lambda_{SAR}^*(M_*) = \log \lambda_{SAR,0}^* + k_\lambda(\log M_* - \log M_{*,0})$ with $\log M_{*,0} = 11$.

The assumption of a double power-law for $f_{\lambda_{SAR}}$, allows recovering the double power-law shape of the XLF with a Schechter function HGMF, as demonstrated by Aird et al. (2013). We fixed the break value to $\log \lambda_{SAR,0}^* = 33.8$ [erg/s/M$_{\odot}$] to limit the number of free parameters. This value is close to the implied Eddington limit, consistent with the approach in the study of Aird et al. (2013), and with the tentative evidence for such a break first reported in Bongiorno et al. (2012).

Finally, we parameterized the redshift evolution of the normalization of the space density as

$$ f_z(z) = \begin{cases} (1+z)^{p_1} & z \le z_0 \\ (1+z_0)^{p_1} \left(\frac{1+z}{1+z_0}\right)^{p_2} & z > z_0, \end{cases} \quad (12) $$

where we fixed $z_0 = 1.1$, motivated by the break redshift used in the LDDE model in the XLF from Miyaji et al. (2015) and approximately corresponding to the central redshift in our sample.

The best-fit bivariate distribution function $\Psi(M_*, \lambda_{SAR}, z)$ is shown in Fig. 2, the best-fitting parameters and their errors are given in Table 1. We computed the uncertainties of each parameter using a Markov chain Monte Carlo (MCMC) sampling of the likelihood function space, using emcee (Foreman-Mackey et al. 2013), a Python implementation of an affine invariant MCMC ensemble sampler as presented by Goodman et al. (2010). We used uniform priors for our free parameters and initialized the MCMC “walkers” around the best-fit maximum likelihood solution. The quoted uncertainties represent the 16 and 84 percentile of the parameter distribution, marginalized over all other parameters except for $\Psi^*$. The latter is not determined by the maximum likelihood fit, and their error is given by $1/\sqrt{N_{\text{tot}}}$.

As mentioned above, our best-fit HGMF and SARDF given by Eqs. (8) and (9) cannot be expressed as simple analytic functions because $M_*$ and $\lambda_{SAR}$ are too entangled in the SARDF term. For a better quantitative representation of the redshift evolution of HGMF and SARDF and for illustrative purposes, we provide an analytic approximation of the two distribution functions, evaluated at the center of our three redshift bins. For this, we performed a least-squares fit to the HGMF (computed via Eq. (8)) at each redshift with a standard Schechter function with normalization $\Phi_M^*$, break $M_*$, and low-mass slope $\alpha$, and the SARDF (computed via Eq. (9)) with a double power-law with normalization $\Phi_\lambda^*$, break $\lambda_{SAR}^*$, and slopes $\gamma_1, \gamma_2$. We provide the best-fit parameters in Tables 2 and 3.

### 3.3. $V_{\text{max}}$ method

An additional consistency check can be obtained by computing the AGN host galaxy mass function using the $V_{\text{max}}$ method. The
<figure><img src="image_4.png" /><figcaption>**Fig. 3.** Upper panels: total AGN HGMF with the associated errors in three redshift bins derived through the maximum likelihood (black line and gray shaded area) and the $V_{\text{max}}$ (data points) methods. In each bin the lowest z fit is reported for reference with a dashed line. The red and blue lines are the total and the star-forming galaxy stellar mass functions with the associated errors (shown in magenta and cyan shaded areas) from Ilbert et al. (2013). Lower panels: ratio of AGN host galaxies to the total (red line) and the star-forming (blue line) galaxy population as a function of stellar mass in the same redshift bins.</figcaption></figure>

<table><caption>**Table 2.** Best-fit model parameters for the AGN host galaxy mass Schechter function, computed in our three redshift bins.</caption><thead><tr><th>⟨z⟩</th><th>log(Φ<sub>M</sub><sup>*</sup>)</th><th>log M<sub>*</sub><sup>*</sup></th><th>α</th></tr></thead><tbody><tr><td>0.55</td><td>-3.83<sup>+0.04</sup><sub>-0.05</sub></td><td>10.99<sup>+0.03</sup><sub>-0.03</sub></td><td>-0.41<sup>+0.04</sup><sub>-0.04</sub></td></tr><tr><td>1.15</td><td>-3.54<sup>+0.04</sup><sub>-0.05</sub></td><td>10.99<sup>+0.03</sup><sub>-0.03</sub></td><td>-0.24<sup>+0.04</sup><sub>-0.04</sub></td></tr><tr><td>2.00</td><td>-3.84<sup>+0.04</sup><sub>-0.04</sub></td><td>10.99<sup>+0.03</sup><sub>-0.03</sub></td><td>-0.03<sup>+0.05</sup><sub>-0.05</sub></td></tr></tbody></table>

<table><caption>**Table 3.** Best-fit model parameters of the AGN specific accretion rate double power-law function, computed in our three redshift bins.</caption><thead><tr><th>⟨z⟩</th><th>log Φ<sub>λ</sub><sup>*</sup></th><th>γ<sub>1</sub></th><th>γ<sub>2</sub></th><th>log λ<sub>SAR</sub><sup>*</sup></th></tr></thead><tbody><tr><td>0.55</td><td>-6.04<sup>+0.08</sup><sub>-0.08</sub></td><td>-1.35<sup>+0.02</sup><sub>-0.02</sub></td><td>-3.64<sup>+0.10</sup><sub>-0.11</sub></td><td>34.33<sup>+0.04</sup><sub>-0.04</sub></td></tr><tr><td>1.15</td><td>-5.22<sup>+0.08</sup><sub>-0.09</sub></td><td>-1.02<sup>+0.02</sup><sub>-0.02</sub></td><td>-3.61<sup>+0.10</sup><sub>-0.10</sub></td><td>34.32<sup>+0.03</sup><sub>-0.03</sub></td></tr><tr><td>2.00</td><td>-4.85<sup>+0.08</sup><sub>-0.09</sub></td><td>-0.54<sup>+0.03</sup><sub>-0.03</sub></td><td>-3.58<sup>+0.10</sup><sub>-0.10</sub></td><td>34.30<sup>+0.03</sup><sub>-0.03</sub></td></tr></tbody></table>

$V_{\text{max}}$ for each individual object is given by

$$V_{\text{max}}(M_{\star}) = \int_{z_{\text{min}}}^{z_{\text{max}}} \mathcal{A}(M_{\star}, z) \frac{dV}{dz} dz, \qquad (13)$$

where $\mathcal{A}(M_{\star}, z)$ is the effective area as a function of $M_{\star}$ and z given by the total survey area $\Omega$ times the incompleteness function. We emphasize here that the $V_{\text{max}}(M_{\star})$ values used are not identical to the $V_{\text{max}}(L_{\text{X}})$ values that would be used to compute the AGN luminosity function. This is because, as discussed above, we also have to account in the incompleteness function for the SARDF in addition to the sensitivity function and the absorption correction. The incompleteness function thus includes three terms and can be written as

$$\bar{I}(M_{\star}, z) = \int_{20}^{24} \int_{\lambda_{\text{SARmin}}}^{\infty} \bar{I}(M_{\star}, \lambda_{\text{SAR}}, z, N_{\text{H}}) f(N_{\text{H}} | L_{\text{X}}, z) \\ \times f_{\lambda_{\text{SAR}}}(\lambda_{\text{SAR}}, M_{\star}, z) \text{dlog } \lambda_{\text{SAR}} \text{dlog } N_{\text{H}}, \qquad (14)$$

where $\bar{I}(M_{\star}, \lambda_{\text{SAR}}, z, N_{\text{H}}) = \bar{I}(L_{\text{X}}, z, N_{\text{H}})$ is the X-ray selection function given by the sensitivity map in the [2–10] keV band, $f(N_{\text{H}}|L_{\text{X}}, z)$ is the absorption distribution function from Ueda et al. (2014), and $f_{\lambda_{\text{SAR}}}(\lambda_{\text{SAR}}, M_{\star}, z)$ is the SARDF term in $\Psi_{\lambda_{\text{SAR}}}(M_{\star}, \lambda_{\text{SAR}}, z)$. The latter term is required for the mass-dependent incompleteness function $\bar{I}(M_{\star}, \lambda_{\text{SAR}}, z, N_{\text{H}})$ in addition to the terms needed to compute the luminosity-dependent incompleteness function $\bar{I}(L_{\text{X}}, z)$.

While the $V_{\text{max}}$ method has the advantage of providing a non-parametric estimate of the AGN host galaxy mass function, it has the disadvantages that it requires a specific assumption for the SARDF term and, furthermore, it does not include the additional constraints from the AGN XLF, which, as a result of the limited luminosity range probed by our sample, makes the results less robust in particular at the low-mass end, where we only probe a limited range in $\lambda_{\text{SAR}}$. In contrast, the maximum likelihood provides a parametric estimate of the mass function and determines the HGMF and SARDF simultaneously and self-consistently. Therefore we only used the $V_{\text{max}}$ method as a consistency check. For the function $f_{\lambda_{\text{SAR}}}(\lambda_{\text{SAR}}, M_{\star}, z)$ we assumed the best-fit $M_{\star}$-dependent SAR-term determined above (Eq. (11)), normalized within $\log \lambda_{\text{SAR}} > 32$ [erg/s/M$_{\odot}$], which again defines our lower integration limit.

The AGN host galaxy mass function is thus computed in three redshift bins as

$$\Phi(M_{\star}) = \frac{1}{\Delta \log M_{\star}} \sum_{i=1}^{N_{\text{obj}}} \frac{1}{V_{\text{max}}}, \qquad (15)$$

and the binned values are shown in Fig. 3 together with the maximum likelihood result. The error bars are determined by bootstrapping of the sample with their $V_{\text{max}}(M_{\star})$ values.

As shown in the figure, we find a good overall agreement between the $V_{\text{max}}$ binned AGN HGMF and the AGN HGMF based on the maximum likelihood method. This confirms the adopted parametric model in the maximum likelihood approach and verifies the robustness of our results.
<figure><img src="image_5.png" /><figcaption>**Fig. 4.** *Upper panels:* the SAR term split into three redshift bins and described as a double power-law with a mass dependent $\lambda_{\text{SAR}}^*$ as in Eq. (11). *Lower panels:* SARDF derived through the maximum likelihood method by integrating the bivariate distribution function over $M_*$ (Eq. (8)). The shaded area includes the errors on the parameters.</figcaption></figure>

### 3.4. Results

In the upper panles of Fig. 4 we show the SAR term $f_{\lambda_{\text{SAR}}}$ (Eq. (11)), described by a double power-law with mass-dependent, but redshift-independent break $\lambda_{\text{SAR}}^*$. The SARDF, shown in the lower panels of the same figure, is obtained by integrating the bivariate distribution function (including the above function) over $M_*$.

The SARDF can be described by a double power-law whose low $\lambda_{\text{SAR}}$ characteristic slope flattens from -1.35 to -0.54 from the lowest to the highest redshift bin. Conversely, the overall normalization $\phi_{\lambda}^*$ increases for increasing redshift (see Table 3). The increasing normalization with redshift was noted in Aird et al. (2012) and Bongiorno et al. (2012). In these works, the specific accretion rate distribution was parametrized with a single power-law over the full redshift range, but Bongiorno et al. (2012) noted the presence of a break above $\log \lambda_{\text{SAR}} > 34$ [erg/s/M$_{\odot}$]. Furthermore, Aird et al. (2013) argued for a break in the specific accretion rate distribution to be consistent with the XLF.

While these previous studies did not report a change in the shape of the specific accretion rate distribution with redshift, we find a SARDF clearly flattening toward higher redshift. It is important to note that compared to the aforementioned works, there are some differences. First, here we determined the SARDF, that is, the absolute space density as a function of $\lambda_{\text{SAR}}$, while the previous studies present $p_{\text{AGN}}(\lambda_{\text{SAR}}|M_*)$, which is the AGN fraction in the galaxy population. Furthermore, we accounted for obscuration by integration over the $N_{\text{H}}$ distribution, which generally steepens our low $\lambda_{\text{SAR}}$ slope.

The sample used in Bongiorno et al. (2012) is instead similar and largely overlaps the one used in this study. A more accurate analysis of the sample used in Bongiorno et al. (2012) could indeed reveal the redshift dependence of the specific accretion rate, which was not included in the parametric model presented in Bongiorno et al. (2012) because of the simpler single power-law parametrization. Finally, we modeled the bivariate distribution function of $\lambda_{\text{SAR}}$ and $M_*$ and not only the SARDF and included additional information on the XLF. We discuss the effect of the latter in more detail in the Appendix.

The best-fit HGMF (black line in Fig. 3) is well described by a Schechter function with constant $M_*^*$ and a low-mass slope $\alpha$ flattening with redshift (i.e., $\alpha = -0.41$ in the first redshift bin, -0.24 in the second, and -0.03 in the third; see Eq. (2)). We compared the AGN HGMF with the total galaxy stellar mass function (red curve and shaded magenta region) and the star-forming galaxy mass function (blue curve and shaded cyan region) by Ilbert et al. (2013). We note that at $\log(M_*/M_\odot) > 11.5$, the HGMF but also the total and SF galaxy mass functions are only poorly constrained by the data (see Fig. 1) because of the limited volume sampled in both cases. This region is indicated by the dashed lines in Fig. 3. Furthermore, in the highest z-bin the galaxy mass function of Ilbert et al. (2013) shows an upturn at low masses, captured in their double Schechter function model, which is not captured in our more restricted single Schechter function model for the HGMF. Our data do not allow constraining such an upturn for our AGN sample, which would require a larger sample and probably a deeper flux limit for the galaxies including lower luminosity AGN.
<figure><img src="image_6.png" /><figcaption>**Fig. 5.** Left panel: redshift evolution of the SARDF space density for different $\lambda_{\text{SAR}}$. Magenta: $\log(\lambda_{\text{SAR}}) = 32$ [erg/s/$M_{\odot}$] (~1% Edd); green: $\log(\lambda_{\text{SAR}}) = 32.5$ [erg/s/$M_{\odot}$] (~3% Edd); yellow: $\log(\lambda_{\text{SAR}}) = 33$ [erg/s/$M_{\odot}$] (~10% Edd); red: $\log(\lambda_{\text{SAR}}) = 33.5$ [erg/s/$M_{\odot}$] (~30% Edd); blue: $\log(\lambda_{\text{SAR}}) = 34.0$ [erg/s/$M_{\odot}$] (~Edd); and black: $\log(\lambda_{\text{SAR}}) = 34.5$ [erg/s/$M_{\odot}$] (>Edd). Central panel: redshift evolution of the HGMF space density for different $M_{\star}$. Right panel: redshift evolution of the XLF space density for different $L_{\text{X}}$ by Miyaji et al. (2015).</figcaption></figure>

The ratio of AGN HGMF over total galaxy mass function is shown by the red line and the shaded magenta area in the lower panels of Fig. 3. This ratio indicates the active fraction or duty cycle of AGN activity in the galaxy population, if we consider AGN with $\log \lambda_{\text{SAR}} > 32$ (~1% of the Eddington limit), which corresponds to the definition of an AGN assumed in this paper.

We find a redshift evolution in the mass dependence of the active fraction. At $M_{\star} = 10^{10} M_{\odot}$, the active fraction is approximately constant at ~10%, while at $M_{\star} = 10^{11.5} M_{\odot}$ it increases over our three redshift bins from ~3% to ~8% to ~20%. This trend qualitatively agrees with the results for the SMBH mass dependence of the active fraction of the black hole mass function presented in Schulze et al. (2015). This might be related to the redshift evolution of the gas reservoir available to fuel the AGN, since in high-redshift galaxies a greater amount of gas can be responsible for triggering AGN activity (Tacconi et al. 2010).

The ratio of AGN HGMF to the star-forming mass function (shown by the blue line in the lower panels of Fig. 3) traces the average relation between star-forming and AGN activity as a function of stellar mass. It extends the well-known average agreement between star formation rate density and black hole accretion density (e.g., Marconi et al. 2004) to its stellar mass dependence. Overall, we find a weaker redshift evolution in the shape of this ratio than for the active fraction, where the ratio stays almost constant over $10^{10} < M_{\star} < 10^{11} M_{\odot}$, the mass range tracing the bulk of the population, in all three redshift bins. At the high-mass end for $z > 0.8$ the AGN/SF galaxy ratio and for $z > 1.5$ also the active fraction appear to increase with stellar mass. Future studies are required to confirm or disprove the reality of this trend.

The redshift evolution of the SARDF and HGMF allows a more detailed look at the AGN downsizing behavior, i.e., the luminosity-dependent evolution, seen in the XLF out to $z \sim 2.5$. They probe the more physically meaningful quantities stellar mass and specific accretion rate distribution, and by inference relate to black hole mass and Eddington ratio. In Fig. 5 we show the global trend of the redshift evolution of the space density in bins of $\lambda_{\text{SAR}}$ (left panel), $M_{\star}$ (central panel) and $L_{\text{X}}$ (right panel). The $L_{\text{X}}$ dependence, based on the XLF from Miyaji et al. (2015) shows the well-known AGN downsizing behavior (e.g., Ueda et al. 2003; Hasinger et al. 2005; La Franca et al. 2005; Bongiorno et al. 2007; Silverman et al. 2008). For the $\lambda_{\text{SAR}}$ dependence, we see that higher $\lambda_{\text{SAR}}$ objects

($\log \lambda_{\text{SAR}} > 33.5$ [erg/s/$M_{\odot}$]) peak in their space density at an earlier cosmic epoch than the lower $\lambda_{\text{SAR}}$ objects ($\log \lambda_{\text{SAR}} < 33.5$ [erg/s/$M_{\odot}$]), which means that they also show a clear downsizing trend. The $M_{\star}$ dependence, based on the HGMF, also indicates a downsizing trend, with AGN in lower stellar mass galaxies showing a steeper decline in their space density toward high redshift than higher stellar mass galaxies, but less pronounced than what is seen in the SARDF. This suggests that the downsizing in the AGN luminosity function is due to the combination of a (weak) mass-dependent evolution of the HGMF and the stronger evolution of the SARDF.

In the upper panels of Fig. 6 we show the AGN HGMF for different luminosity subclasses: $\log(L_{\text{X}}) < 43$ [erg/s] (magenta), $\log(L_{\text{X}}) > 43$ [erg/s] (cyan), $\log(L_{\text{X}}) > 43.5$ [erg/s] (green), $\log(L_{\text{X}}) > 44$ [erg/s] (red), and $\log(L_{\text{X}}) > 44.5$ [erg/s] (blue). As expected, the high-mass end is dominated by luminous AGN ($\log L_{\text{X}} > 43$ [erg/s]), while the low-mass bins are mainly populated by low-luminosity objects ($\log L_{\text{X}} < 43$ [erg/s]) whose contribution above $\log(M_{\star}/M_{\odot}) \sim 11$ is negligible. Our definition threshold of $\log L_{\text{X}}/M_{\star} > 32$ directly excludes any AGN with $\log L_{\text{X}} < 43$ [erg/s] above $M_{\star} > 10^{11} M_{\odot}$. This also implies that when applying an AGN definition by a luminosity threshold, as is usually done, an active fraction increasing with mass is usually found, consistent with previous work (e.g., Bundy et al. 2008; Xue et al. 2010; Aird et al. 2012; Silverman et al. 2009).

In the lower panels, we instead show the total AGN HGMF in $\lambda_{\text{SAR}}$ bins: $\log \lambda_{\text{SAR}} > 32.5$ [erg/s/$M_{\odot}$] (magenta), $\log \lambda_{\text{SAR}} > 33$ [erg/s/$M_{\odot}$] (cyan), $\log \lambda_{\text{SAR}} > 33.5$ [erg/s/$M_{\odot}$] (green), and $\log \lambda_{\text{SAR}} > 34$ [erg/s/$M_{\odot}$] (red). Overall, the mass distributions of AGN of different specific accretion rate have a similar shape, only mildly affected by the $M_{\star}$ dependence in our SARDF model.

### 4. Mass function of galaxies in the process of being mass-quenched

According to the model described in Peng et al. (2010), the quenching process, that is, the process that leads to the transition from star-forming to passive galaxies, independent of its physical origin, can be described by two different modes: mass and environment quenching, whose differential effects on the fraction of passive or red galaxies are separable.
<figure><img src="image_7.png" /><figcaption>Fig. 6. Total AGN HGMF in three redshift bins derived through the maximum likelihood (black line) and compared with the AGN HGMF for different AGN sub-samples in (upper panels) luminosities. log(L<sub>X</sub>) < 43 [erg/s] (magenta); log(L<sub>X</sub>) > 43 [erg/s] (cyan), log(L<sub>X</sub>) > 43.5 [erg/s] (green), log(L<sub>X</sub>) > 44 [erg/s] (red), and log(L<sub>X</sub>) > 44.5 [erg/s] (blue); and (lower panels) specific accretion rates λ<sub>SAR</sub>, i.e., log λ<sub>SAR</sub> > 32.5 [erg/s/M<sub>☉</sub>] (magenta); log λ<sub>SAR</sub> > 33 [erg/s/M<sub>☉</sub>] (cyan), log λ<sub>SAR</sub> > 33.5 [erg/s/M<sub>☉</sub>] (green), and log λ<sub>SAR</sub> > 34 [erg/s/M<sub>☉</sub>] (red).</figcaption></figure>

Peng et al. (2010) speculated that the environment quenching occurs in satellite galaxies, while the mass quenching could reflect a feedback mechanism related to star-formation or AGN. In a subsequent paper, Peng et al. (2012) confirmed the expectation on the environment quenching as due to satellite galaxies from studying the mass function of central and satellite galaxies. Here we wish to test whether the mass quenching process can be linked to AGN feedback.

The strength of the approach of Peng et al. (2010) is that this phenomenological model is based on simple observational inputs, which allows successfully reproducing many of the features of the galaxy population. Moreover, the model is able to clearly predict the mass function of the galaxies in the process of being mass-quenched and the inter-relationships between the Schechter parameters for star-forming and passive galaxies.

The mass function of the transient population can be described by a single Schechter function with parameters (see Eq. (28) of Peng et al. 2010)

$$ 
\begin{aligned}
M_{\star, \text{trans}}^{*} &= M_{\star, \text{blue}}^{*} \\
\alpha_{\text{s,trans}} &= \alpha_{\text{s,blue}} + (1 + \beta) \\
\Phi_{\text{trans}}^{*} &= \Phi_{\text{blue}}^{*} s S F R(M_{\star}, z) |_{M_{\star}} \tau_{\text{trans}},
\end{aligned} 
\quad (16) $$

where $M_{\star,blue}^*$, $\alpha_{s,blue}$, and $\Phi_{blue}^*$ are the parameters of the Schechter function that describes the star-forming galaxy mass function, and $\beta$ is the exponent in the power-law relation that links the specific star formation rate (sSFR) and the stellar mass (see Eq. (17)). Here we used the data for star-forming galaxies from Ilbert et al. (2013) and forced the fit with a single Schechter

function. This parametric choice is required to use the model fits provided by Peng et al. (2010) with a single Schechter function as starting MF. This introduces some uncertainties especially with respect to the slope of the high-mass end, which is the most difficult part of the stellar MF to be constrained, as we point out below. The value $\tau_{trans}$ is the period of time that the transient signature is visible and is not constrained by the model of Peng et al. (2010). Here, we assumed that the transient phase corresponds to the active feedback or blow-out phase, which is the gas-depletion timescale associated with the outflow. Current observations suggest this timescale to be on the order of $1-10 \times 10^7$ yr (Maiolino et al. 2007; Feruglio et al. 2010; Cicone et al. 2014). Finally, sSFR(M<sub>star</sub>, z) is the evolving specific star formation rate. Here we consider the recent measurement of the sSFR from Lilly et al. (2013, Eq. (2)),

$$ 
\begin{aligned}
sSFR(M, z) &= 0.07 \left( \frac{M_{\star}}{10^{10.5} M_{\odot}} \right)^{\beta} (1+z)^3 \text{ Gyr}^{-1} \text{ for } z < 2 \\
sSFR(M, z) &= 0.30 \left( \frac{M_{\star}}{10^{10.5} M_{\odot}} \right)^{\beta} (1+z)^{5/3} \text{ Gyr}^{-1} \text{ for } z > 2
\end{aligned} 
\quad (17) $$

with $\beta \sim -0.1$.

Starting from the star-forming galaxy mass function (green line in Fig. 7), we then derived with the above equations the predicted mass function of the “transient” (i.e. in the process of being mass-quenched) population. Given the uncertainties on the value of $\tau_{trans}$, we show in Fig. 7 the predictions for a range of $\tau_{trans} = 1-10 \times 10^7$; the blue solid line is for $\tau = 5 \times 10^7$ while
<figure><img src="image_8.png" /><figcaption>**Fig. 7.** Predicted MF of transient galaxies as derived starting from a single Schechter function based on the data from Ilbert et al. (2013, green line), using the Peng recipe (Peng et al. 2010) for $\tau = 5 \times 10^7$ yr (blue solid line), $\tau = 1 \times 10^7$ yr and $\tau = 1 \times 10^8$ yr (blue dashed lines). We compare them with the total AGN HGMF (black line) and with the HGMF computed with different luminosity cuts. In particular in the first redshift bin, cuts of $\log L_X > 43.8$ [erg/s] (red solid line) and $\log L_X > 44.2$ and $\log L_X > 43.5$ [erg/s] for the lower and upper boundary, respectively, have been applied; while in the second and third redshift bins, these cuts are $\log L_X > 44$ (brown solid line) and $\log L_X > 44.5$ [erg/s] and $\log L_X > 43.5$ [erg/s] for the lower and upper boundaries.</figcaption></figure>

the blue dashed lines correspond to $1 \times 10^7$ yr and $1 \times 10^8$ yr (lower and upper boundary, respectively).

To test whether AGN can be responsible for the mass-quenching of galaxies, we chose to restrict our analysis to the most luminous objects. Theory indeed predicts that the capability of AGN outflows of perturbing the ISM depends on AGN luminosity as $L_{\text{Bol}}^{1/2}$ (Menci et al. 2008) and that galaxy-scale outflows are energy-driven, meaning that their mechanical energy is proportional to the AGN luminosity (Zubovas & King 2012). This scenario is supported by observations that find that the momentum rate of kpc-scale outflows (Sturm et al. 2011; Cicone et al. 2014; Feruglio et al. 2015) is $\ge 10-20 L_{\text{Bol}}/c$. This means that the more luminous the AGN, the more powerful the outflows that are produced. Accordingly, the AGN-driven feedback mechanism should become increasingly more efficient in halting the star-formation in the host galaxy for higher AGN luminosities.

In Fig. 7 we compare the prediction for the mass function of mass quenching transient objects with the HGMF of the total population, that is, $\log \lambda_{\text{SAR}} > 32$, and of different subsamples. We tested the agreement using subsamples by also applying different cuts on either $L_X$ or $\lambda_{\text{SAR}}$, as shown in Fig. 6. To keep the comparsion simple, we did not consider more complicated cuts or a luminosity-dependent transition timescale, for instance, which could improve the agreement between the two mass functions. We found that the class of objects that best reproduces the expected mass function in terms of shape and normalization are $\log L_X > 43.8^{+0.4}_{-0.3}$ [erg/s] (red solid line and yellow shaded area) at $0.3 < z < 0.8$, and $\log L_X > 44 \pm 0.5$ [erg/s] at $0.8 < z < 2.5$. Reducing the threshold in $L_X$ leads to a space density in the HGMF higher than expected for the transient objects at the low-mass end. Conversely, specific accretion-rate-based subsamples do not seem to reproduce the expected mass function particularly well. This is because within the model of Peng et al. (2010), the fractional density of the transition population strongly decreases at low masses for a constant $\tau_{\text{trans}}$: only very few low-mass galaxies experience quenching at any redshift. On the other hand, the poluation of AGN above any given $\lambda_{\text{SAR}}$ threshold increases toward low stellar masses (see the bottom panel of Fig. 6): rapidly growing objects with a high Eddington ratio can be found in galaxies of any mass, at all redshifts. Thus, any model that invokes a fixed threshold in $\lambda_{\text{SAR}}$ to explain the quenching population (see, e.g., Zubovas & King 2012) would predict a too high fraction of low-mass galaxies in the transition phase, in strong contrast with the finding of Peng et al. (2010).

We note that the disagreement between the Peng model prediction and the AGN HGMF present at the high mass end of the third redshift bin arises because we used a single Schechter function to fit the star-forming galaxy mass function to apply the recipe from Peng et al. (2010), which, as pointed out by Ilbert et al. (2013), is not a good fit of the data, especially at the high mass end. The fit with a double Schechter function (as performed in Ilbert et al. (2013) and shown in Fig. 3) would indeed be steeper, thus reducing the number of predicted high-mass objects and the discrepancy with the AGN HGMF.

Overall, we find the space density of luminous AGN ($\log L_X > 43.5-44.5$ [erg/s]) at stellar masses $M_* > 10^{10.7} M_{\odot}$ to be consistent with the space density of galaxies in the star formation quenching phase. This non-trivial result is consistent with the notion that feedback from luminous AGN can be associated with the mass-quenching of galaxies. At lower masses the difference in space density between the luminous AGN mass function and the quenching mass function leaves room for a contribution from another mechanism. Lower luminosity AGN might contribute here, if their AGN feedback mechanism would operate on a different transition timescale $\tau_{\text{trans}}$. Furthermore, Peng et al. (2015) recently suggested that strangulation (a mechanism for which the supply of cold gas is halted) is the primary mechanism responsible for quenching star-formation in local galaxies with a stellar mass lower than $10^{11} M_{\odot}$. Our results are complementary to this work, proposing AGN-driven outflows as a plausible mechanism for halting star formation at higher redshift and for more massive galaxies, although a causal connection is not substantiated.

## 5. Summary and conclusions

We have studied the host galaxy stellar mass function of a sample of ~1000 AGN detected in the XMM-COSMOS field in the 2–10 keV band at $0.3 < z < 2.5$. We derived the SARDF and the HGMF simultaneously as a bivariate distribution function of stellar mass and specific accretion rate $\lambda_{\text{SAR}} \equiv L_X/M_*$, using the maximum likelihood method developed by Schulze & Wisotzki (2010) and extended by Schulze et al. (2015).

Our results can be summarized as follows:

(i) The SARDF is best described by a double power-law with a mass-dependent but redshift-independent break $\lambda_{\text{SAR}}^*$ and a low $\lambda_{\text{SAR}}$ characteristic slope that flattens from -1.35 to
-0.54 with increasing redshift. The overall normalization $\phi_{\lambda}^{*}$, in contrast, increases for increasing redshift.

(ii) The AGN HGMF is described by a Schechter function with constant $M_{*}^{\star}$ and a low-mass slope $\alpha$ flattening with redshift from $\alpha = -0.41$ at $z = 0.55$ to $\alpha = -0.03$ at $z = 2.0$. We derived the active fraction of AGN activity by comparison with the stellar mass function by Ilbert et al. (2013), and we found a redshift evolution in its mass dependence at the high-mass end, where the fraction of AGN in massive galaxies increases from ~3% at $z \sim 0.55$ to ~20% at $z \sim 2$.

(iii) The redshift evolution of the SARDF and AGN HGMF allows us to gain a deeper understanding of the physical drivers of the AGN downsizing behavior, which is seen in the XLF out to $z \sim 2.5$. We find that the downsizing in the AGN luminosity function is due to the combination of a (weak) mass-dependent evolution of the HGMF and the stronger evolution of the SARDF. In particular, we see that higher $\lambda_{SAR}$ objects have a peak in their space density at earlier epoch than the lower $\lambda_{SAR}$ AGN.

(iv) We compared the mass function of the population in the process of being mass-quenched, predicted by the phenomenological model by Peng et al. (2010), with the HGMF computed for different subsamples obtained with different luminosity and $\lambda_{SAR}$-cuts. We find at the high masses (i.e., $M_{\star} > 10^{10.7} M_{\odot}$) that the population that agrees with the model prediction is that of luminous AGN with log $L_X > 43.5-44.5$ [erg/s] (i.e., log $L_{bol} \ge 46$ [erg/s]). Their number density and stellar mass distribution are consistent with those of the transition galaxy population, a crucial, and non-trivial, result of our analysis. While this agreement does not establish a causal connection between star formation quenching and AGN activity, it suggests AGN feedback by powerful outflows from luminous AGN as a plausible mechanism for the mass-quenching of star forming galaxies. This scenario would agree and be complementary to the recent findings by Peng et al. (2015), who suggested strangulation as the primary quenching mechanism at lower masses (i.e., $M_{\star} < 10^{11} M_{\odot}$).

**Acknowledgements.** We thank the referee for the careful reading and helpful suggestions that helped to improve the manuscript. This work is based on the COSMOS program. The HST COSMOS Treasury program was supported through NASA grant HST-GO-09822. This work is mainly based on observations obtained with XMM-Newton, an ESA Science Mission with instruments and contributions directly funded by ESA Member States and the USA (NASA), and with the European Southern Observatory under Large Program 175-A-0839, Chile. A.B. and E.P. acknowledge financial support from INAF under the contract PRIN-INAF-2012. A.S. acknowledges support by JSPS KAKENHI Grant Number 26800098. MCMC corner plots make use of *triangle.py* (Foreman-Mackey et al. 2013).

## References

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

Aird, J., Coil, A. L., Georgakakis, A., et al. 2015, MNRAS, 451, 1892

Akaike, H. 1974, IEEE Transactions on Automatic Control, 19, 716

Bongiorno, A., Zamorani, G., Gavignaud, I., et al. 2007, A&A, 472, 443

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

Boyle, B. J., & Terlevich, R. J. 1998, MNRAS, 293, L49

Brusa, M., Fiore, F., Santini, P., et al. 2009, A&A, 507, 1277

Brusa, M., Civano, F., Comastri, A., et al. 2010, ApJ, 716, 348

Buchner, J., Georgakakis, A., Nandra, K., et al. 2015, ApJ, 802, 89

Bundy, K., Georgakakis, A., Nandra, K., et al. 2008, ApJ, 681, 931

Cano-Díaz, M., Maiolino, R., Marconi, A., et al. 2012, A&A, 537, L8

Cappelluti, N., Brusa, M., Hasinger, G., et al. 2009, A&A, 497, 635

Cen, R., & Chisari, N. E. 2011, ApJ, 731, 11

Chabrier, G. 2003, PASP, 115, 763

Cicone, C., Maiolino, R., Sturm, E., et al. 2014, A&A, 562, A21

Cowie, L. L., Songaila, A., Hu, E. M., & Cohen, J. G. 1996, AJ, 112, 839

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

Cresci, G., Marconi, A., Zibetti, S., et al. 2015b, A&A, 582, A63

Croom, S. M., Richards, G. T., Shanks, T., et al. 2009, MNRAS, 1439

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

Faucher-Giguère, C.-A., & Quataert, E. 2012, MNRAS, 425, 605

Feruglio, C., Maiolino, R., Piconcelli, E., et al. 2010, A&A, 518, L155

Feruglio, C., Fiore, F., Carniani, S., et al. 2015, A&A, 583, A99

Foreman-Mackey, D., Hogg, D. W., Lang, D., & Goodman, J. 2013, PASP, 125, 306

Georgakakis, A., Coil, A. L., Willmer, C. N. A., et al. 2011, MNRAS, 418, 2590

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

Gilli, R., Comastri, A., & Hasinger, G. 2007, A&A, 463, 79

Gilli, R., Vignali, C., Mignoli, M., et al. 2010, A&A, 519, A92

Goodman, J., & Weare, J. 2010, Commun. Appl. Math. Comput. Sci., 5, 65

Granato, G. L., De Zotti, G., Silva, L., Bressan, A., & Danese, L. 2004, ApJ, 600, 580

Häring, N., & Rix, H. 2004, ApJ, 604, L89

Hasinger, G., Miyaji, T., & Schmidt, M. 2005, A&A, 441, 417

Hasinger, G., Cappelluti, N., Brunner, H., et al. 2007, ApJS, 172, 29

Hopkins, P. F., Hernquist, L., Cox, T. J., et al. 2006, ApJS, 163, 1

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

Kauffmann, G., & Heckman, T. M. 2009, MNRAS, 397, 135

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

La Franca, F., Fiore, F., Comastri, A., et al. 2005, ApJ, 635, 864

Lansbury, G. B., Gandhi, P., Alexander, D. M., et al. 2015, ApJ, 809, 115

Lilly, S. J., Carollo, C. M., Pipino, A., Renzini, A., & Peng, Y. 2013, ApJ, 772, 119

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

Mainieri, V., Bongiorno, A., Merloni, A., et al. 2011, A&A, 535, A80

Maiolino, R., Neri, R., Beelen, A., et al. 2007, ApJ, 472, L33

Marconi, A., & Hunt, L. K. 2003, ApJ, 589, L21

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

Menci, N., Fiore, F., Puccetti, S., & Cavaliere, A. 2008, ApJ, 686, 219

Merloni, A., Bongiorno, A., Brusa, M., et al. 2014, MNRAS, 437, 3550

Miyaji, T., Hasinger, G., Salvato, M., et al. 2015, ApJ, 804, 104

Monaco, P., & Fontanot, F. 2005, MNRAS, 359, 283

Mullaney, J. R., Pannella, M., Daddi, E., et al. 2012, MNRAS, 419, 95

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

Peng, Y.-J., Lilly, S. J., Kovač, K., et al. 2007, ApJ, 721, 193

Peng, Y.-J., Lilly, S. J., Renzini, A., & Carollo, M. 2012, ApJ, 757, 4

Peng, Y., Maiolino, R., & Cochrane, R. 2015, Nature, 521, 192

Read, J. I., & Trentham, N. 2005, Roy. Soc. London Philosophical Trans. Ser. A, 363, 2693

Salvato, M., Ilbert, O., Hasinger, G., et al. 2011, ApJ, 742, 61

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

Schawinski, K., Khochfar, S., Kaviraj, S., et al. 2006, Nature, 442, 888

Schulze, A., & Wisotzki, L. 2010, A&A, 516, A87

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

Sijacki, D., Vogelsberger, M., Genel, S., et al. 2015, MNRAS, 452, 575

Silverman, J. D., Green, P. J., Barkhouse, W. A., et al. 2008, ApJ, 679, 118

Silverman, J. D., Lamareille, F., Maier, C., et al. 2009, ApJ, 696, 396

Somerville, R. S., Primack, J. R., & Faber, S. M. 2001, MNRAS, 320, 504

Springel, V., Di Matteo, T., & Hernquist, L. 2005, ApJ, 620, L79

Sturm, E., González-Alfonso, E., Veilleux, S., et al. 2011, ApJ, 733, L16

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

Treister, E., Urry, C. M., & Virani, S. 2009, ApJ, 696, 110

Ueda, Y., Akiyama, M., Ohta, K., & Miyaji, T. 2003, ApJ, 598, 886

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

Vignali, C., Mignoli, M., Gilli, R., et al. 2014, A&A, 571, A34

Xue, Y. Q., Brandt, W. N., Luo, B., et al. 2010, ApJ, 720, 368

Zubovas, K., & King, A. 2012, ApJ, 745, L34
# Appendix A: Model comparison

As mentioned in Sect. 3.2, we derived uncertainties through an MCMC computation of the posterior distribution function (PDF). In Fig. A.1 we show the 1D marginalized PDF for the free parameters of our model and the 2D marginalized PDF for parameter pairs. The latter shows the covariance between these pairs. We find covariance between several parameters, for instance, between the break log $M_*^*$ and the slope $\alpha$ of the $M_*$ term $f_*$, or between several of the redshift evolution parameters ($p_1, p_2, k_\gamma$).

In addition to the default parametric model we presented above, we explored parametric models for the HGMF and SARDF. We here discuss the results of this exercise and provide a justification for our chosen parameterization. To compare the relative quality of our respective parametric model given our data set, we used the Akaike information criterion (AIC<sub>c</sub>; Akaike 1974). It is given by $AIC_c = S + 2K + 2K(N/(N-K-1))$, where S

is the likelihood as defined above, K is the number of parameters in the model, and N is the size of the sample. While the AIC<sub>c</sub> penalizes against overfitting, it is known to be less penalizing than the Bayesian information criterion (BIC), for example. Nevertheless, the difference in AIC between the models tested below is in general significant enough to draw firm conclusions.

Our default model (hereafter model A) allows for a mass dependence on the SARDF, while previous studies assumed a $\lambda_{SAR}$ distribution, independent of mass (e.g., Aird et al. 2012; Bongiorno et al. 2012). We first started our analysis with a model without such a mass-dependent term, $k_\lambda = 0$ (hereafter model B). Our best-fit model B provides an almost equivalent fit to the XLF as our default model A (see red dashed line in Fig. A.3) and fits well over most of the $M_* - \lambda_{SAR}$ plane, as shown in the middle panels of Fig. A.2. We provide the best-fit parameters in Table A.1. However, it is unable to recover the observed number of objects at low masses and high-$\lambda_{SAR}$ (upper left corner). In the upper panels of Fig. A.2 we show the $M_* - \lambda_{SAR}$ plane for our

<figure><img src="image_9.png" /><figcaption>**Fig. A.1.** Posterior probability distribution for the free parameters of our best-fit model. The histograms at the top show the marginalized probability distribution function for each free parameter, while the contour plots show the marginalized 2D probability distribution function for each pair of parameters, to illustrate their covariances.</figcaption></figure>
<figure><img src="image_10.png" /><figcaption>**Fig. A.2.** Comparison of the bivariate distribution in the $M_*$ – $\lambda_{SAR}$ plane between the observations (symbols) and the prediction from the best-fit model (black contours) in the three redshift bins as labeled. The *upper panels* show the comparison for our default model **A**), the *middle panels* show the prediction from the best-fit $M_*$-independent SARDF model **B**), and the *lower panels* show the best-fit solution if no XLF data are included in the likelihood function ("no XLF"). Open symbols indicate upper limits in $M_*$.</figcaption></figure>

<figure><img src="image_11.png" /><figcaption>**Fig. A.3.** Comparison of the X-ray luminosity function (XLF) predicted by the bivariate distribution $\Psi(M_*, \lambda_{SAR}, z)$ with direct observations. The blue solid line shows the XLF derived from our default model (A), the red dashed line is for the best-fit $M_*$-independent SARDF model (B), and the green solid line is for the best fit without the XLF into the likelihood function (no XLF). The black circles show the binned XLF for our XMM-COSMOS sample using the $V_{max}$ method, indicating the luminosity range covered by XMM-COSMOS. The shaded cyan and yellow areas show the XLF by Miyaji et al. (2015) and Ueda et al. (2014), respectively, where the shaded area includes the variation of the XLF over the redshift range of the bin and the uncertainty of the XLF determination.</figcaption></figure>
<table><caption>**Table A.1.** Best-fit model parameters for the bivariate distribution function of stellar mass and SAR for other models, which have been tested but were not considered further.</caption><thead><tr><th>Model</th><th>log (Ψ*)</th><th>log M<sub>*</sub>[M<sub>☉</sub>]</th><th>α<sub>0</sub></th><th>k<sub>α</sub></th><th>log λ<sub>SAR</sub><sup>*</sup></th><th>k<sub>λ</sub></th><th>log M<sub>*</sub>,0</th><th>γ<sub>1,0</sub></th><th>k<sub>γ</sub></th><th>γ<sub>2</sub></th><th>p<sub>1</sub></th><th>p<sub>2</sub></th><th>z<sub>0</sub></th></tr></thead><tbody><tr><td>No XLF</td><td>-6.69</td><td>11.02</td><td>-0.503</td><td>0.0*</td><td>33.8*</td><td>0.05</td><td>11.0*</td><td>-1.11</td><td>0.17</td><td>-1.97</td><td>5.49</td><td>1.45</td><td>1.1*</td></tr><tr><td>B</td><td>-6.67</td><td>10.88</td><td>-0.192</td><td>0.0*</td><td>33.8*</td><td>0*</td><td>11.0*</td><td>-0.85</td><td>0.82</td><td>-3.35</td><td>6.32</td><td>2.39</td><td>1.1*</td></tr><tr><td>C</td><td>-6.90</td><td>11.00</td><td>0.225</td><td>-0.017</td><td>33.8*</td><td>-0.48</td><td>11.0*</td><td>-1.01</td><td>0.60</td><td>-3.80</td><td>5.90</td><td>2.37</td><td>1.1*</td></tr></tbody></table>

**Notes.** The parameters denoted with an asterisk are fixed during the fit.

<table><caption>**Table A.2.** Comparison between different parametric models and our default model through their difference in AIC<sub>c</sub> and their corresponding relative likelihood.</caption><thead><tr><th>Model</th><th>Fixed parameters</th><th>ΔAIC<sub>c</sub></th><th>Relative likelihood</th></tr></thead><tbody><tr><td>Default (A)</td><td>k<sub>α</sub> = 0/k<sub>λ</sub>, k<sub>γ</sub> free</td><td>0</td><td>1.0</td></tr><tr><td>(A) + z-evolution in M<sub>*</sub><sup>*</sup></td><td>k<sub>α</sub> = 0/k<sub>λ</sub>, k<sub>γ</sub>, k<sub>M<sub>*</sub><sup>*</sup></sub> free</td><td>10</td><td>5.5 × 10<sup>-3</sup></td></tr><tr><td>z-evolution in α (C)</td><td>k<sub>α</sub>, k<sub>λ</sub>, k<sub>γ</sub> free</td><td>14</td><td>9.5 × 10<sup>-4</sup></td></tr><tr><td>(C) + z-evolution in M<sub>*</sub><sup>*</sup></td><td>k<sub>α</sub>, k<sub>λ</sub>, k<sub>γ</sub>, k<sub>M<sub>*</sub><sup>*</sup></sub> free</td><td>28</td><td>1.0 × 10<sup>-6</sup></td></tr><tr><td>No M<sub>*</sub> dependence on f<sub>λSAR</sub> (B)</td><td>k<sub>α</sub> = 0, k<sub>λ</sub> = 0/k<sub>γ</sub> free</td><td>311</td><td>2.7 × 10<sup>-68</sup></td></tr><tr><td>(B) + z-evolution in α</td><td>k<sub>λ</sub> = 0/k<sub>α</sub>, k<sub>γ</sub> free</td><td>311</td><td>2.5 × 10<sup>-68</sup></td></tr><tr><td>No z-evolution in f<sub>λSAR</sub></td><td>k<sub>γ</sub> = 0/k<sub>α</sub>, k<sub>λ</sub> free</td><td>531</td><td>6.4 × 10<sup>-116</sup></td></tr><tr><td>No z-evolution in f<sub>λSAR</sub> and f<sub>*</sub></td><td>k<sub>α</sub> = 0, k<sub>γ</sub> = 0/k<sub>λ</sub> free</td><td>731</td><td>1.7 × 10<sup>-159</sup></td></tr></tbody></table>

default model A. This mass-dependent SARDF model is able to match the observations in this region as well (since the SARDF has a higher break at low mass). The AIC<sub>c</sub> ratio between the two models prefers our model A with a relative likelihood of $10^{-68}$, providing strong evidence that this model better describes our data set.

An important constraint on this mass dependence in the SARDF is set by the inclusion of the XLF information in our likelihood function, since our sample does not cover a very wide dynamical range in luminosity. If we neglect the XLF and only use the data from our XMM-COSMOS sample, such a mass dependence is not strongly required and the flattening in the SARDF is likewise weaker, while still present. We give the best-fit solution without XLF data as model “no XLF” in Table A.1. It also matches the XMM-COSMOS sample in the $M_{\star} - \lambda_{SAR}$ plane well (see lower panels in Fig. A.2). However, this is achieved by proposing a high space density of objects with low $M_{\star}$ and low $\lambda_{SAR}$, below the luminosity limit of the XMM-COSMOS sample. Such a high space density of low-luminosity AGN violates observations of the XLF from deeper surveys (Miyaji et al. 2015; Aird et al. 2015) and is thus deemed not physical. This is demonstrated in Fig. A.3, where we show the XLF predicted by our bivariate distribution function models. The “no XLF” model (shown by the green solid line) is consistent with our default model (blue line) and the XLF at $z < 0.8$, where XMM-COSMOS still probes a relatively wide luminosity range. At higher redshift it clearly overpredicts the XLF. This motivates the additional inclusion of the XLF information in

this work, to find a solution to our data set that is consistent with XLF measurements beyond the luminosity range directly covered and thus constrained by our sample. On the other hand, the XLF alone is degenerate between $M_{\star}$ and $\lambda_{SAR}$ and thus does not constrain the SARDF and HGMF, as demonstrated by the identical XLF for models A and B. When combining the XMM-COSMOS data with the XLF, the discussed mass and redshift dependencies are required by the data. However, the mass dependence of the SARDF is mainly driven by the number of the few low-mass and high-$\lambda_{SAR}$ objects. Future studies, incorporating both deeper surveys (e.g., CDFS, CDFN) and shallower, larger area surveys (e.g., XMM-XXL, Stripe 82X), will be essential to settle this question.

In addition, we also tested a model where we allowed for redshift evolution in the slope of the $M_{\star}$ term $\alpha(z) = \alpha_{0} + k_{\alpha}(z - z_{0})$, with $z_{0}$ fixed to 1.1 (hereafter model C). However, we found the best-fit redshift evolution parameter to be consistent with zero $k_{\alpha} = -0.017 \pm 0.043$, thus the inclusion of this additional parameter is not statistically justified. This is also confirmed by the AIC<sub>c</sub> ratio between the two models. Similarly, we found that the addition of another parameter that allows for redshift evolution in the break of the $M_{\star}$ term, $M_{\star}^{*}(z) = M_{\star,0}^{*} + k_{M_{\star}}(z - z_{0})$, is not statistically justified (see Table A.2).

We also performed the same test of goodness to justify the redshift dependence on the slope of the SAR term $k_{\gamma}$ and found this term to be statistically justified. We provide the relative AIC<sub>c</sub> values and the corresponding relative likelihood for different models, compared to our default model (A), in Table A.2.