AGN host galaxy mass function in COSMOS:
Key Words.:Galaxies:active, Galaxies: fundamental parameters, Galaxies: evolution
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 i.e., , distribution function (SARDF) up to z2.5 with 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 slope flattens with increasing redshift while the normalization increases. This implies that, for a given stellar mass, higher objects have a peak in their space density at earlier epoch compared to the lower ones, 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 a almost constant and a low mass slope that flattens with redshift. Compared to the stellar mass function, we find that the HGMF has a similar shape and that, up to 11.5 the ratio of AGN host galaxies to star forming galaxies is basically constant (10%). Finally, the comparison of the AGN HGMF for different luminosity and specific accretion rate sub-classes with the phenomenological model prediction by Peng et al., 2010 for the “transient” population, i.e. 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, i.e. , feedback from luminous AGN ( [erg/s]) may be responsible for the quenching of star formation in the host galaxy.
Super-massive black hole (SMBH) growth, nuclear activity, and galaxy evolution, have been found to be closely related. In fact, over the last 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 in which Active Galactic Nuclei (AGN) are key players in the process of galaxy formation and evolution. Several theoretical models (e.g. Somerville et al. 2001; Granato et al. 2004; Monaco & Fontanot 2005; Springelet 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 to find the mechanism responsible for the simultaneous fuelling of the central BH and the formation of new stars in the host galaxy, as well as the quasi-simultaneous shut-off of both 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), they usually propose the presence of an energetic AGN-driven feedback, i.e. a strong wind originated 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 shut off both 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 galaxies heating and/or pushing away the gas which is forming stars.
Star formation quenching via some mechanism is required also to prevent the overgrowth of massive galaxies, hosted in the most massive dark-matter haloes (e.g. Read & Trentham 2005). Such a “mass quenching” mechanism, irrespective 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 such a role for AGN, 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 been often contradictory, i.e. 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 by the models including 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 interpret this as a physical link between these two phenomena which however do not directly imply a causal link.
In fact, irrespective of AGN feedback, an essential pre-requisite to understand the role of black hole activity in galaxy evolution is to have a accurate and unbiased census of the AGN population an 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, i.e. 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 the study of their stellar mass function and/or the active fraction or duty cycle of AGN occurrence in galaxies of 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 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). In fact, Aird et al. (2012) showed that the intrinsic distribution of specific accretion rates at follows a power law, whose shape does not evolve with redshift, independent of stellar mass. This result has been confirmed and extended out to by Bongiorno et al. (2012).
In this work, we build upon the aforementioned 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 . 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 and/or be responsible for mass quenching using the model prediction from Peng et al. (2010) for the mass function of the ’transient population” (i.e. galaxies in the process of being mass-quenched).
The paper is organized as follows: In Sec. 2 we present the X-ray selected sample we are using. Sec 3 presents the method used to derive the specific accretion rate distribution function and the AGN host galaxy mass function (Sec. 3.2 and 3.3) and their results (Sec. 3.4). In Sec 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 subsamples, with the model prediction for quenching galaxies by Peng et al. (2010).
Throughout this paper, a standard cosmology (=0.3, =0.7 and H=70 km s Mpc) has been assumed. The stellar masses are given in units of solar masses for a Chabrier IMF (Chabrier 2003).
2 The Sample
The AGN sample considered here has been extracted from the XMM-COSMOS point-like source catalogue (Hasinger et al. 2007; Cappelluti et al. 2009) whose optical identifications and multiwavelength properties have been presented by Brusa et al. (2010). The catalog contains 1800 X-ray sources detected above flux limits of 510, 310 and 710 erg cms 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 is chosen since the soft band can be affected by obscuration that can lead to a redshift-dependent incompleteness (i.e. flux limited surveys pick up more obscured objects at higher redshift, see e.g. Gilli et al. 2010). However, this band may still suffers from incompleteness due to heavily obscured and Compton Thick (CT, ) AGN, whose detection probability is strongly reduced because the intrinsic emission can be significantly suppressed due to repeated Compton scattering and photoelectric absorption.
Out of the full 1800 sources, we identify a final sample of 927 hard X-ray selected AGN in the redshift range 0.3z2.5. All hard X-ray sources have accurate photometric redshifts (Salvato et al. 2011) while half (581/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 have been 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, by applying a K-correction computed by assuming an intrinsic X-ray spectrum with a photon index =1.8, an exponential cut-off at E=200 keV and, a photoelectric absorption corresponding to the observed N column density. The [2-10] keV luminosity is given by:
where D is the luminosity distance and is the term which accounts for the K-correction and absorption correction. The absorbing column density N for our sample has been derived as in Merloni et al. (2014). For the brightest sources (above 200 pn counts in the 0.5-10 keV band of XMM-Netwton) N is obtained from the full spectral analysis of Mainieri et al. (2011), which is available for 195/927 of the AGN. For the remaining sources, N is estimated in a statistical fashion, by assessing the value of the ‘observed’ spectral slope from the hardness ratio and assessing the value of the ’observed’ spectral slope drawn from a normal distribution with mean and dispersion of . 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 as 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 the reader to this paper for a detailed description of the method.
Following Aird et al. (2012) and Bongiorno et al. (2012) we define “specific accretion rate” (see also Brusa et al. 2009; Georgakakis et al. 2014) as a directly measurable quantity which can be regarded as a proxy for the black hole growth rate relative to the stellar mass of the host galaxy, , 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 SMBH’s Eddington ratio, , 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 =25 (Marconi et al. 2004; Lusso et al. 2012) and a constant host stellar to black hole mass ratio of 500 (Marconi & Hunt 2003; Häring & Rix 2004), approximately corresponds to the Eddington limit, while would give 1% of Eddington. The bivariate distribution M - for the analyzed sample is shown in Fig. 1 where different colors correspond to different redshift ranges as labeled.
For the determination of the mass function, we further restrict our sample in stellar mass and specific accretion rate , applying the following cuts: and . The latter criterion is motivated by the requirement of having a clear cut in above which we define the AGN as active (see below). The chosen minimum value corresponds to the lowest observed value in our intermediate redshift bin and furthermore corresponds approximately to 1% of Eddington which we chose in the following as our minimum 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 .
3 The AGN Host Galaxy Mass Function and specific accretion rate distribution function
In order 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.
In fact, completeness in L does not directly ensure completeness in . As previously reported, AGN show a wide range of Eddington ratios (Kauffmann & Heckman 2009; Schulze & Wisotzki 2010), and thus also a wide range of (), 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 towards high mass BHs and high galaxy mass i.e. since an AGN with low Eddington ratio will be included in the sample only if its M is high enough to be above the given luminosity (L) limit, given the relation between M-M, a bias towards high-mass black holes induces a bias toward high-mass galaxies. 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 ; and (3) the stellar mass completeness down to our threshold in units of specific accretion rate .
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 which have been missed in the sample due to their high column density . For this correction we use the distribution as a function of and published by Ueda et al. (2014) based on several X-ray AGN surveys (see their eq. (5) and (6)). We integrate over the distribution between , i.e. we do not include Compton thick AGN in our HGMF determination. The fraction of CT AGN is still uncertain and the distribution above 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 % (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, we additionally suffer from significant incompleteness due to the fact that a broad range of can be associated to a given luminosity . To account for this effect in the HGMF, we need to include an additional term to the incompleteness function based on the distribution of . Using this distribution function, we correct for incompleteness down to a fixed threshold in , which we set at . The HGMF is therefore defined as the mass function of all AGN above this threshold. The most rigorous and self consistent approach to do this is by determining the HGMF and the SARDF simultaneously, e.g. via the maximum likelihood method described in the next section.
3.2 Maximum Likelihood method
We here present the methodology of determining the SARDF and the HGMF simultaneously as a bivariate distribution function of stellar mass and specific accretion rate, i.e. , where gives the space density of AGN with stellar mass host galaxies between and and a specific accretion rate between and at the redshift . The HGMF, SARDF and the X-ray AGN LF (XLF) can be derived as different marginalizations over this bivariate distribution function. We use the maximum likelihood method developed by Schulze & Wisotzki (2010) and extended by Schulze et al. (2015) to compute . 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 is implemented here for the joint determination of the HGMF and SARDF.
The technique minimizes the likelihood function , where the probability distribution for each object is given by:
where is the bivariate distribution function of stellar mass and specific accretion rate that we want to derive, is the X-ray selection function given by the sensitivity map in the 2-10 keV band and is the absorption distribution function, taken from Ueda et al. (2014). We use the estimates presented in Sec. 2.1 to compute (and therefore ) and for our sample. The factor corresponds to the total number of objects in the sample predicted by the model and it is given by integrating over , , N, and , i.e.
where we integrate over the distribution between , while our integration ranges in , and are , and , as discussed in Sec. 2.1.
Our sample also contains 12 AGN without measurements, due to poor quality photometry. However, we account for these sources using their luminosity and redshift information integrated over the entire mass range, i.e.
Our XMM-COSMOS based sample covers only a limited dynamical range in , 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 include as additional observational data the XLF. In this way, we ensure consistency with the XLF observations over its full observationally determined luminosity range. In particular, we use the binned XLF from Miyaji et al. (2015) and compute the value for the comparison with the XLF implied by the HGMF and SARDF. We then add this likelihood to that of the XMM-COSMOS sample. The study by Miyaji et al. (2015) uses the same distribution as Ueda et al. (2014) for the determination of the XLF which we also employ 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 and , which are not accounted for in our bivariate distribution function fit. This may lead to an overestimate of the space density at [erg/s].
The total likelihood to minimize is given by:
where is the number of AGN with measurements in our sample and is the number of AGN with only 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 assume that the bivariate distribution function is separable, i.e. the specific accretion rate distribution is mass independent and vice versa. Under this assumption, the bivariate distribution function is given by:
where is the normalization of the bivariate distribution function, is the SAR-term, is the M-term and 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:
where contains now also a dependence on the mass. We use this more general parametrization as our default model. We want to point out that the SAR-term and the M-term are not equal to the SARDF and HGMF.
The HGMF and the SARDF are calculated by integrating over and over , respectively. To be specific:
In case of separable SAR- and M-terms, as in Eq. 6, the SARDF (HGMF) has the same shape as () and only the absolute normalization is determined by the marginalisation. However, in the more general case of Eq. 7, this is not necessarily the case, which is why the HGMF and SARDF cannot then be explicitly expressed as analytic functions.
We here consider the following parametric models for the individual terms: the M-term is modeled using a Schechter function:
While a model with a low mass slope evolving with redshift has been included, we find that the best fit parameters are indeed consistent with no z-evolution in .
The SAR-term is instead described by a double power law:
where the low slope , with set at , and the break with .
The assumption of a double power law for , allows to recover 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 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 parameterize the redshift evolution of the normalization of the space density as:
where we fixed , 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 is shown in Fig. 2 while 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 initilised the MCMC ”walkers” around the best fit maximum likelihood solution. The quoted uncertainties represent the 16 and 84%-tile of the parameter distribution, marginalized over all other parameters apart from . The latter is not determined by the Maximum likelihood fit and their error is given by .
As mentioned above, our best fit HGMF and SARDF given by Eq. 8 and Eq. 9 cannot be expressed as simple analytic functions, due to the entanglement of and 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, at each redshift, we performed a least-squares fit to the HGMF (computed via Eq. 8) with a standard Schechter function with normalisation , break and low mass slope , and the SARDF (computed via Eq. 9) with a double power law with normalisation , break and slopes , . We provide the best fit parameters in Tab. 2 and 3.
3.3 V method
An additional consistency check can be obtained by computing the AGN host galaxy mass function using the V method. The V for each individual object is given by:
where is the effective area as a function of and given by the total survey area times the incompleteness function. We emphasize here that the values used are not identical to the values that would be used for the computation of the AGN luminosity function. This is because, as discussed above, we have to account in the incompleteness function also 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:
where is the X-ray selection function given by the sensitivity map in the [2-10] keV band, is the absorption distribution function from Ueda et al. (2014) and is the SARDF term in . The latter term is required for the mass dependent incompleteness function in addition to the ones needed for the computation of the luminosity dependent incompleteness function .
While the V 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, due to 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 . On the contrary, the maximum likelihood provides a parametric estimate of the mass function, and determines the HGMF and SARDF simultaneously and self-consistently. Therefore we only use the V method as a consistency check. For the function we assume the best fit -dependent SAR-term determined above (Eq. 11), normalized within , which again defines our lower integration limit.
The AGN Host Galaxy Mass Function is thus computed in three redshift bins as:
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 values.
As shown in the figure, we overall find a good agreement between the V 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.
In the upper panles of Fig. 4, we show the SAR-term f (Eq. 11), described by a double power-law with mass dependent, but redshift independent break . 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 characteristic slope flattens from -1.35 to -0.54 from the lowest to the highest redshift bin. The overall normalization on the contrary increases for increasing redshift (see Tab. 3). The increasing normalization with redshift was already noted in Aird et al. (2012) and Bongiorno et al. (2012). In those works, the specific accretion rate distribution was parametrized with a single power law over the full redshift range, but already Bongiorno et al. (2012) noticed the presence of a break above log. 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 do not report a change in the shape of the specific accretion rate distribution with redshift, we find a SARDF clearly flattening towards higher redshift. It is important to note that, compared to the aforementioned works, there are some differences. First, here we determine the SARDF, i.e. the absolute space density as a function of , while the previous studies present , i.e. the AGN fraction in the galaxy population. Furthermore, we account for obscuration by integration over the distribution, which generally steepens our low slope.
The work by Aird et al. (2012) is refers to , and thus it did not cover a sufficiently large redshift range to constrain this shape evolution. 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), due to the simpler single power-law parametrisation. Finally, we model the bivariate distribution function of and and not only the SARDF and include 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 and a low mass slope flattening with redshift (i.e. in the first redshift bin, in the second and in the third one; see Eq. 2). We compare 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 11.5, the HGMF but also the total and SF galaxy mass functions are not well constrained by the data (see Fig. 1) due to the limited volume sampled in both cases. This region is indicated by the dashed lines in Fig. 3. Furthermore, in the highest -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 to constrain 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.
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. Such ratio indicates the active fraction or duty cycle of AGN activity in the galaxy population, if we consider AGN with (% Eddington) 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 , the active fraction is approximately constant at , while at it increases over our three redshift bins from to to . This trend is in qualitative agreement 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 could 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 , the mass range tracing the bulk of the population, in all three redshift bins. At the high mass end for the AGN/SF galaxy ratio and for also the active fraction appear to increase with stellar mass. Future studies will be 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 behaviour, i.e. the luminosity-dependent evolution, seen in the XLF out to . 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 (left panel), (central panel) and (right panel). The dependence, based on the XLF from Miyaji et al. (2015) shows the well known AGN downsizing behaviour (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 dependence, we see that higher objects () have a peak in their space density at an earlier cosmic epoch compared to the lower objects (), i.e. also showing a clear downsizing trend. The 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 towards 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 Fig. 6 upper panels, we show the AGN HGMF for different luminosity sub-classes, i.e. 43 [erg/s] (magenta), 43 [erg/s] (cyan), 43.5 [erg/s] (green), 44 [erg/s] (red), and 44.5 [erg/s] (blue). As expected the high mass end is dominated by luminous AGN ( [erg/s]), while the low mass bins are mainly populated by low luminosity objects ( [erg/s]) whose contribution above 11 is negligible. Our definition threshold of directly excludes any AGN with [erg/s] above . This also implies that when applying an AGN definition by a luminosity threshold, as usually done, you will tend to find an active fraction increasing with mass, 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 bins i.e., log (magenta), log (cyan), log (green), and log (red). Overall, the mass distributions of AGN of different specific accretion rate have a similar shape, only mildly affected by the dependence in our SARDF model.
4 The mass function of galaxies in the process of being mass-quenched
According to the model described in Peng et al. (2010), the quenching process, i.e. the process which 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/red galaxies are separable.
In Peng et al. (2010) paper, it is speculated that the environment quenching occur 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) confirm the expectation on the environment quenching as due to satellite galaxies, studying the mass function of central and satellite galaxies. Here we want to test whether the mass quenching process can be linked to AGN feedback.
The strength of the Peng et al. (2010) approach is that this phenomenological model is based on simple observational inputs, which allow one to successfully reproduce many of the features of the galaxy population. Moreover, the model is able to give a clear prediction for 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):
where , and are the parameters of the Schechter function which describes the star-forming galaxy mass function and is the exponent in the power law relation that links the specific star formation rate (sSFR) and the stellar mass (see Eq. LABEL:eq:sSFR). Here we use the data for star-forming galaxies from Ilbert et al. (2013), and force 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 slope of the high-mass end, which is the most difficult part of the stellar MF to be constrained, as we will point out later in this section. The value is the period of time the “transient” signature is visible, and is not constrained by the Peng et al. (2010) model. Here, we assume that the transient phase corresponds to the active feedback/blow-out phase i.e. the gas depletion time-scale associated with the outflow. Current observations suggest this time-scale to be of the order of (Maiolino et al. 2007; Feruglio et al. 2010; Cicone et al. 2014). Finally, is the evolving specific star formation rate. Here we consider the recent measurement of the from Lilly et al. (2013, eq. (2)):
Starting from the star-forming galaxy mass function (green line in Fig. 7), we then derive, using 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 we show in Fig. 7 the predictions for a range of ; the blue solid line is for while the blue dashed lines correspond to and (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 (Menci et al. 2008) and that galaxy-scale outflows are energy-driven, i.e., 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 , i.e. the more luminous the AGN is, the more powerful outflows are produced. This means that 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, i.e. , and of different sub-samples. We test the agreement using sub-samples applying in addition different cuts on either or , as shown in Fig. 6. We do not consider more complicated cuts or for example a luminosity dependent transition time-scale, which could improve the agreement between the two mass functions, in order to keep the comparison as simple as possible. We find that the class of objects that best reproduces, in terms of both shape and normalization, the expected mass function are: [erg/s] (red solid line and yellow shaded area) at , and [erg/s] at . Reducing the threshold in leads to a space density in the HGMF higher than the expected for the ”transient” objects at the low mass end. On the contrary, specific accretion rate based sub-samples do not seem to reproduce the expected mass function particularly well. This is because, within the Peng et al. (2010) model, for a constant , the fractional density of the ”transition” population strongly decreases at low masses: only very few low-mass galaxies experience quenching at any redshift. On the other hand, the poluation of AGN above any given threshold increases towards low stellar masses (see the bottom panel of Fig. 6): rapidly growing high Eddington ratio objects can be found in galaxies of any mass, at all redshifts. Thus, any model that invokes a fixed threshold in 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 Peng et al. (2010) finding.
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 is due to the fact that, to apply the recipe from Peng et al. (2010), we used a single Schechter function to fit the star-forming galaxy mass function 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 ( [erg/s]) at stellar masses 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 to 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 via another mechanism. Lower luminosity AGN might contribute here, if their AGN feedback mechanism would operate on a different transition time-scale . 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 less than 10M. 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
In this paper, we have studied the host galaxy stellar mass function of a sample of 1000 AGN detected in the XMM-COSMOS field in the 2-10keV band at 0.32.5. We derived the SARDF and the HGMF simultaneously as a bivariate distribution function of stellar mass and specific accretion rate , using the maximum likelihood method developed by Schulze & Wisotzki (2010) and extended by Schulze et al. (2015).
Our results can be summarized as follows:
The SARDF is best described by a double power-law with a mass dependent but redshift independent break and a low characteristic slope which flattens from to with increasing redshift. The overall normalization on the contrary increases for increasing redshift.
The AGN HGMF is described by a Schechter function with constant and a low mass slope flattening with redshift from =-0.41 at to =-0.03 at . We derived the active fraction of AGN activity by comparison with the stellar mass function by Ilbert et al. (2013) and we find a redshift evolution in its mass dependence at the high mass end, where the the fraction of AGN in massive galaxies increases from 3% at to 20% at .
The redshift evolution of the SARDF and AGN HGMF allows us to gain a deeper understanding into the physical drivers of the AGN downsizing behaviour, seen in the XLF out to . We find that 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 objects have a peak in their space density at earlier epoch compared to the lower AGN.
We compare 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 sub-samples obtained with different luminosity and -cuts. We find at the high masses (i.e. ) that the population that agrees with the model prediction is that of luminous AGN having [erg/s] (i.e. [erg/s]). Both 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 be in agreement and complementary with the recent findings by Peng et al. (2015) who suggested “strangulation” as the primary quenching mechanism at lower masses (i.e. ).
Acknowledgements.We thank the referee for the careful reading and precious suggestions which 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).
Appendix A Model comparison
|Model||fixed parameters||AIC||relative likelihood|
|default (A)||/ , free||0||1.0|
|(A) + -evolution in||/ , , free||10|
|-evolution in (C)||, , free||14|
|(C) + -evolution in||, , , free||28|
|no dependence in (B)||, / free||311|
|(B) + -evolution in||/ , free||311|
|no -evolution in||/ , free||531|
|no -evolution in and||, / free||731|
As mentioned in Sec. 3.2, we derived uncertainties via MCMC computation of the posterior distribution function (PDF). In Fig. 8 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, e.g. between the break and the slope of the -term or between several of the redshift evolution parameters (, , ).
Besides the default parametric model we presented above, we also explored additional 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 use the Akaike information criterion (AIC; Akaike 1974). It is given by AIC, where is the likelihood as defined above, is the number of parameters in the model and is the size of the sample. While the AIC penalizes against overfitting, it is known to be less penalizing than e.g. the Bayesian information criterion (BIC). Nevertheless, in general the difference in AIC between the models tested below is significant enough to draw firm conclusions.
Our default model (hereafter model A) allows for a mass dependence in the SARDF, while previous studies assumed a distribution, independent of mass (e.g. Aird et al. 2012; Bongiorno et al. 2012). In fact, we first started our analysis with a model without such a mass dependent term, i.e. (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. 10) and fits well over most of the plane, as shown in the middle panels of Fig. 9. We provide the best fit parameters in Table 4. However, it is not able to recover the observed number of objects at low masses and high- (upper left corner). In the upper panels of Fig. 9 we show the plane for our default model A. This mass dependent SARDF model is able to match the observations also in this region (since the SARDF has a higher break at low mass). The AIC ratio between the two models prefers our model A with a relative likelihood of , providing strong evidence that this model provides a better description for 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. In fact, if we neglect the XLF and only use the data from our XMM-COSMOS sample, such a mass dependence is not strongly required and also the flattening in the SARDF is less strong, while still present. We give the best fit solution without XLF data as model ”no XLF” in Tab. 4. It also matches the XMM-COSMOS sample in the plane well (see lower panels in Fig. 9). However, this is achieved by proposing a high space density of objects with low and low , 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. 10, 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 , where XMM-COSMOS still probes a relative 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 which is consistent with XLF measurements outside of the luminosity range directly covered and thus constrained by our sample. On the other hand, the XLF alone is degenerate between and , 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 largely driven by the number of few low mass and high 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 -term , with fixed to 1.1 (hereafter model C). However, we found the best fit redshift evolution parameter to be consistent with zero , thus the inclusion of this additional parameter is not statistically justified. This is also confirmed by the AIC ratio between the two models. Similarly, we found that the addition of another parameter which allows for redshift evolution in the break of the -term, , is not statistically justified (see Table 5).
We also performed the same test of goodness to justify the redshift dependence in the slope of the SAR-term and found this term to be statistically justified. We provide the relative AIC values and the corresponding relative likelihood for different models, compared to our default model (A) in Table 5.
- Aird, J., Coil, A. L., Georgakakis, A., et al. 2015, ArXiv e-prints
- 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
- Akaike, H. 1974, IEEE Transactions on Automatic Control, 19, 716
- Bongiorno, A., Merloni, A., Brusa, M., et al. 2012, MNRAS, 427, 3103
- Bongiorno, A., Zamorani, G., Gavignaud, I., et al. 2007, A&A, 472, 443
- Boyle, B. J. & Terlevich, R. J. 1998, MNRAS, 293, L49
- Brusa, M., Civano, F., Comastri, A., et al. 2010, ApJ, 716, 348
- Brusa, M., Fiore, F., Santini, P., et al. 2009, A&A, 507, 1277
- 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, ArXiv e-prints
- 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., Fiore, F., Carniani, S., et al. 2015, ArXiv e-prints
- Feruglio, C., Maiolino, R., Piconcelli, E., et al. 2010, A&A, 518, L155
- 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., Cappelluti, N., Brunner, H., et al. 2007, ApJS, 172, 29
- Hasinger, G., Miyaji, T., & Schmidt, M. 2005, A&A, 441, 417
- Hopkins, P. F., Hernquist, L., Cox, T. J., et al. 2006, ApJS, 163, 1