CONTINUITY EQUATION AND ABUNDANCE MATCHING

# Black Hole and Galaxy Coevolution from Continuity Equation and Abundance Matching

R. Aversa11affiliation: SISSA, Via Bonomea 265, 34136 Trieste, Italy 44affiliation: INFN-Sezione di Trieste, via Valerio 2, 34127, Trieste, Italy , A. Lapi11affiliation: SISSA, Via Bonomea 265, 34136 Trieste, Italy 22affiliation: Dip. Fisica, Univ. ‘Tor Vergata’,Via Ricerca Scientifica 1, 00133 Roma, Italy 33affiliation: INAF-Osservatorio Astronomico di Trieste, Via Tiepolo 11, 34131, Trieste, Italy 44affiliation: INFN-Sezione di Trieste, via Valerio 2, 34127, Trieste, Italy , G. de Zotti11affiliation: SISSA, Via Bonomea 265, 34136 Trieste, Italy 55affiliation: INAF-Osservatorio Astronomico di Padova, Vicolo Osservatorio 5, 35122 Padova, Italy , F. Shankar66affiliation: School of Physics and Astronomy, University of Southampton, Southampton SO17 1BJ, UK , L. Danese11affiliation: SISSA, Via Bonomea 265, 34136 Trieste, Italy 33affiliation: INAF-Osservatorio Astronomico di Trieste, Via Tiepolo 11, 34131, Trieste, Italy 44affiliation: INFN-Sezione di Trieste, via Valerio 2, 34127, Trieste, Italy
###### Abstract

We investigate the coevolution of galaxies and hosted supermassive black holes throughout the history of the Universe by a statistical approach based on the continuity equation and the abundance matching technique. Specifically, we present analytical solutions of the continuity equation without source term to reconstruct the supermassive black hole (BH) mass function from the AGN luminosity functions. Such an approach includes physically-motivated AGN lightcurves tested on independent datasets, which describe the evolution of the Eddington ratio and radiative efficiency from slim- to thin-disc conditions. We nicely reproduce the local estimates of the BH mass function, the AGN duty cycle as a function of mass and redshift, along with the Eddington ratio function and the fraction of galaxies with given stellar mass hosting an AGN with given Eddington ratio. We exploit the same approach to reconstruct the observed stellar mass function at different redshift from the UV and far-IR luminosity functions associated to star formation in galaxies. These results imply that the buildup of stars and BHs in galaxies occurs via in-situ processes, with dry mergers playing a marginal role at least for stellar masses and BH masses , where the statistical data are more secure and less biased by systematic errors. In addition, we develop an improved abundance matching technique to link the stellar and BH content of galaxies to the gravitationally dominant dark matter component. The resulting relationships constitute a testbed for galaxy evolution models, highlighting the complementary role of stellar and AGN feedback in the star formation process. In addition, they may be operationally implemented in numerical simulations to populate dark matter halos or to gauge subgrid physics. Moreover, they may be exploited to investigate the galaxy/AGN clustering as a function of redshift, mass and/or luminosity. In fact, the clustering properties of BHs and galaxies are found to be in full agreement with current observations, so further validating our results from the continuity equation. Finally, our analysis highlights that (i) the fraction of AGNs observed in slim-disc regime, where anyway most of the BH mass is accreted, increases with redshift; (ii) already at substantial dust amount must have formed over timescales yr in strongly starforming galaxies, making these sources well within the reach of ALMA surveys in (sub-)millimeter bands.

black hole physics - galaxies: formation - galaxies: evolution - methods: analytical - quasars: supermassive black holes

## 1. Introduction

Kinematic and photometric observations of the very central regions in local, massive early-type galaxies strongly support the almost ubiquitous presence of black holes (BHs) with masses (Dressler 1989; Kormendy & Richstone 1995; Magorrian et al. 1998; for a recent review Kormendy & Ho 2013). Their formation and evolution is a major problem in astrophysics and physical cosmology.

The correlations between the central BH mass and galaxy properties such as the mass in old stars (Kormendy & Richstone 1995; Magorrian et al. 1998; Marconi & Hunt 2003; Häring & Rix 2004; McLure & Dunlop 2004; Ferrarese & Ford 2005; Graham 2007; Sani et al. 2011; Beifiori et al. 2012; McConnell & Ma 2013; Kormendy & Ho 2013), the velocity dispersion (Ferrarese & Merritt 2000; Gebhardt et al. 2000; Tremaine et al. 2002; Gültekin et al. 2009; McConnell & Ma 2013; Kormendy & Ho 2013; Ho & Kim 2014), and the inner light distribution (Graham et al. 2001; Lauer et al. 2007; Graham & Driver 2007; Kormendy & Bender 2009) impose strong ties between the formation and evolution of the BH and that of the old stellar population in the host galaxy (Silk & Rees 1998; Fabian 1999; King 2005; for a recent review see King 2014).

A central role in this evolution is played by the way dark matter (DM) halos and associated baryons assemble. So far it has been quite popular, e.g. in most semi-analytic models, to elicit merging as the leading process; as to the baryons, ‘wet’ and ‘dry’ mergers or a mixture of the two kinds have been often implemented (for a recent review see Somerville & Davé 2015). On the other hand, detailed analyses of DM halo assembly indicate a two-stage process: an early fast collapse during which the central regions reach rapidly a dynamical quasi-equilibrium, followed by a slow accretion that mainly affects the halo outskirts (e.g., Zhao et al. 2003; Wang et al. 2011; Lapi & Cavaliere 2011). Thus one is led to consider the rapid starformation episodes in the central regions during the fast collapse as the leading processes in galaxy formation (e.g., Lapi et al. 2011, 2014; Cai et al. 2013). Plainly, the main difference between merging and fast collapse models relates to the amount of stars formed in-situ (e.g., Moster et al. 2013).

While body simulations of DM halo formation and evolution are nowadays quite robust (though details of their results are not yet fully understood), the outcomes of hydrodynamical simulations including star formation and central BH accretion are found to feature large variance (Scannapieco et al. 2012; Frenk & White 2012). This is expected since most of the relevant processes involving baryons such as cooling, gravitational instabilities, angular momentum dissipation, star formation and supermassive BH accretion occur on spatial and temporal scales well below the current resolution.

On the other hand, observations of active galactic nuclei (AGNs) and galaxies at different stages of their evolution have spectacularly increased in the last decade at many wavelengths. In particular, the AGN luminosity function is rather well assessed up to though with different uncertainties in the X-ray (Ueda et al. 2014; Fiore et al. 2012; Buchner et al. 2015; Aird et al. 2010, 2015), UV/optical (Richards et al. 2006; Croom et al. 2009; Masters et al. 2012; Ross et al. 2012; Fan et al. 2006; Jiang et al. 2009; Willott et al. 2010a), and IR bands (Richards et al. 2006; Fu et al. 2010; Assef et al. 2011; Ross et al. 2012); these allow to infer the BH accretion rate functions at various redshifts. In addition, luminosity functions of galaxies are now available up to in the ultraviolet (UV; Bouwens et al. 2015; Finkelstein et al. 2014, Weisz et al. 2014; Cucciati et al. 2012; Oesch et al. 2010; Reddy & Steidel 2009; Wyder et al. 2005) and up to in the far-infrared (FIR) band (Lapi et al. 2011; Gruppioni et al. 2013; Magnelli et al. 2013); these allow to infer the star formation rate (SFR) function at various redshifts.

As for galaxies selected by their mid- and far-IR emission, the distribution function of the luminosity associated to the formation of massive stars shows that at the number density of galaxies endowed with star formation rates yr is Mpc. The density is still significant, Mpc, for yr. On the other hand, the UV selection elicits galaxies forming stars at much lower rates yr up to . The complementarity between the two selections is ascribed to the increasing amount of dust in galaxies with larger SFRs (Steidel et al. 2001; Mao et al. 2007; Bouwens et al. 2013, 2015; Fan et al. 2014; Cai et al. 2014; Heinis et al. 2014). From deep, high resolution surveys with ALMA at (sub-)mm wavelengths there have been hints of possible source blending at fluxes mJy (Karim et al. 2013; Ono et al. 2014; Simpson et al. 2015b). On the other hand, observations at high spatial resolution of sub-mm selected, high redshift galaxies with the SMA and follow-ups at radio wavelengths with the VLA show that galaxies exhibiting a few yr have a number density Mpc (Barger et al. 2012, 2014), fully in agreement with the results of Lapi et al. (2011) and Gruppioni et al. (2013) based on Herschel (single dish) surveys.

Studies on individual galaxies show that several sub-mm galaxies at high redshift exhibit yr concentrated on scales kpc (e.g., Finkelstein et al. 2014; Neri et al. 2014; Rawle et al. 2014; Riechers et al. 2014; Ikarashi et al. 2014; Simpson et al. 2015a; Scoville et al. 2014). Size ranging from a few to several kpc of typical high strongly star forming galaxies has been confirmed by observations of many gravitational lensed objects (e.g., Negrello et al. 2014). In addition, high spatial resolution observations around optically selected quasars put in evidence that a non negligible fraction of host galaxies exhibits yr (Omont et al. 2001, 2003; Carilli et al. 2001; Priddey et al. 2003; Wang et al. 2008; Bonfield et al. 2011; Mor et al. 2012).

The clustering properties of luminous sub-mm selected galaxies (Webb et al. 2003; Blain et al. 2004; Weiss et al. 2009; Hickox et al. 2012; Bianchini et al. 2015) indicate that they are hosted by large halos with masses several and that the star formation timescale is around Gyr.

The statistics on the presence of AGNs along the various stages of galaxy assembling casts light on the possible reciprocal influence between star formation and BH accretion (for a recent review, see Heckman & Best 2014 and references therein), although the fine interpretation of the data is still debated. On one side, some authors suggest that star formation and BH accretion are strongly coupled via feedback processes, while others support the view that the two processes are only loosely related and that the final relationships among BH mass and galaxy properties are built up along the entire Hubble time with a relevant role of dry merging processes.

Most recently, Lapi et al. (2014) have shown that the wealth of data at strongly support the view that galaxies with final stellar mass proceed with their star formation at an almost constant rate over Gyr, within a dusty interstellar medium (ISM). At the same time several physical mechanisms related to the star formation, such as gravitational instabilities in bars or dynamical friction among clouds of starforming gas or radiation drag (Norman & Scoville 1988; Shlosman et al. 1989, 1990; Shlosman & Noguchi 1993; Hernquist & Mihos 1995; Noguchi 1999; Umemura 2001; Kawakatu & Umemura 2002; Kawakatu et al. 2003; Thompson et al. 2005; Bournaud et al. 2007, 2011; Hopkins & Quataert 2010, 2011), can make a fraction of the ISM loose angular momentum and flow into a reservoir around the seed BH. The accretion from the reservoir to the BH can be as large as times the Eddington rate, leading to slim-disc conditions (Abramowicz et al. 1988; Watarai et al. 2001; Blandford & Begelman 2004; Li 2012; Begelman 2012; Madau et al. 2014; Volonteri & Silk 2014), with an Eddington ratio and an average radiative efficiency . This results in an exponential increase of the BH mass and of the AGN luminosity, with an folding timescale ranging from a few to several years. Eventually, the AGNs at its maximum power can effectively transfer energy and momentum to the ISM, removing a large portion of it from the central regions and so quenching the star formation in the host. The reservoir around the BH is no more fed by additional gas, so that even the accretion and the nuclear activity come to an end.

More in general, we can implement lightcurves for the luminosity associated to the star formation and to the BH accretion in a continuity equation approach. In the context of quasar statistics, the continuity equation has been introduced by Cavaliere et al. (1971) to explore the optical quasar luminosity evolution and its possible relation with the radiosource evolution. Soltan (1982) and Choksi & Turner (1992) exploited the mass-energy conservation to derive an estimate of the present mass density in inactive BHs. The extension and the derivation of the BH mass function has been pioneered by Small & Blandford (1992), who first attempted to connect the present day BH mass function to the AGNs luminosity evolution. A simplified version in terms of mass-energy conservation has been used by Salucci et al. (1999), who have shown that the distribution of the mass accreted onto central BHs during the AGNs activity well matches the mass function of local inactive BHs. A detailed discussion of the continuity equation in the quasar and central supermassive BH context has been presented by Yu & Lu (2004, 2008). In the last decade the continuity equation has been widely used, though the lightcurve of the AGN, one of the fundamental ingredients, was largely based on assumptions (e.g., Marconi et al. 2004; Shankar et al. 2009, 2013; Merloni & Heinz 2008; Wang et al. 2009; Cao 2010). Results on the BH mass function through the continuity equation have been reviewed by Kelly & Merloni (2012) and Shankar (2013).

We will also implement the continuity equation for the stellar content of galaxies. This has become possible because the UV surveys for Lyman Break Galaxies, the wide surveys HerMES and H-ATLAS obtained with the Herschel space observatory made possible to reconstruct the star formation rate function in the Universe out to for SFRs yr. Therefore, we can exploit the continuity equation, in an analogous manner as routinely done for the AGNs; the BH mass is replaced by the mass in stars and the bolometric luminosity due to accretion is replaced by the luminosity generated by the formation of young, massive stars.

As for the stellar mass function, it is inferred by exploiting the observed luminosity function in the wavelength range of the SED dominated by the emission from older, less massive stars. The passage from stellar luminosity to mass is plagued by several problems, which result in uncertainties of order of a factor of , increasing for young, dusty galaxies (e.g., Cappellari et al. 2013; Conroy et al. 2014). Therefore the mass estimate is more robust for galaxies with quite low star formation and/or in passive evolution. All in all, the stellar mass function of galaxies is much easier to estimate, and hence much better known, than the BH mass function, particularly at high redshift. Reliable stellar mass functions are available both for local and redshift up to galaxy samples (e.g., Bernardi et al. 2013; Maraston et al. 2013; Ilbert et al. 2013; Santini et al. 2012; Stark et al. 2009; González et al. 2011; Duncan et al. 2014). The comparison between the observed stellar mass function and the results from the continuity equation sheds light on the relative contribution of dry merging and of in-situ star formation. In the present paper we will solve the continuity equation for AGN and for the stellar component after inserting the corresponding lightcurves derived from the data analysis of Lapi et al. (2011, 2014, see above).

Once the stellar and the BH mass functions at different redshifts are known, these can also be compared with the abundance of DM halos to obtain interesting relationships between halo mass and galaxy/BH properties. Such a technique, dubbed abundance matching, has been exploited by several authors (e.g., Vale & Ostriker 2004; Shankar et al. 2006; Moster et al. 2010, 2013; Behroozi et al. 2013; Behroozi & Silk 2015). In this paper, the technique is refined and used in connection with the outcomes of the continuity equation to tackle the following open issues in galaxy formation and evolution:

• Is the BH mass function reflecting the past AGN activity? What was the role of merging in shaping it? (see Sect. 2.1 and Appendix B)

• How does the BH duty cycle evolve? What can we infer on the radiative efficiency and on the Eddington ratio of active BHs? (see Sect. 2.1.4)

• Is there any correlation between the central BH mass and the halo mass and how does it evolve with time? (see Sect. 3.1.1)

• Which is the relationship between the AGN bolometric luminosity and the host halo mass? Can we use this relationship with the duty cycle to produce large simulated AGN catalogs? (see Sect. 3.1.1)

• Which are the bias properties of AGNs? Do they strongly depend on luminosity and redshift? (see Sect. 3.1.1)

• Can the evolution of the stellar mass function be derived through the continuity equation as in the case of the BH mass function, by replacing the accretion rate with the SFR? Does dry merging play a major role in shaping the stellar mass function of galaxies? What is the role of the dust in the star formation process within galaxies? (see Sect. 2.2, App. B and C)

• Which is the relationship between the SFR, the stellar mass of the galaxies, and the mass of the host halo? Does the star formation efficiency (i.e., the fraction of baryons going into stars) evolve with cosmic time? (see Sect. 3.1.2)

• Which is the relationship between the bolometric luminosity of galaxies due to star formation and the host halo mass? Can we use this relationship with the stellar duty cycle to produce large simulated catalogs of star forming galaxies? (see Sect. 3.1.2)

• Which are the bias properties of star forming and passively evolving galaxies? (see Sect. 3.1.2)

• How does the specific star formation rate evolve with redshift and stellar mass? (see Sect. 3.1.3)

• Which is the relationship between the BH mass and the stellar mass at the end of the star formation and BH mass accretion epoch? Does it evolve with time? (see Sect. 3.1.4)

• How and to what extent can we extrapolate the relationships for both galaxies and hosted AGNs to higher, yet unexplored, redshift? (see Sects. 3 and 4)

To answer these questions we have organized the paper as follows. In Sect. 2 we present the statistical data concerning AGN and starforming galaxies, introduce and motivate the corresponding lightcurves, and solve the continuity equation to derive the BH and stellar mass functions at different redshifts. In Sect. 3 we exploit the abundance matching technique to infer relationships among the properties of the BH, stellar, and dark matter component in galaxies. In Sect. 4 we discuss and summarize our findings.

Throughout this work we adopt the standard flat concordance cosmology (Planck Collaboration 2014, 2015) with round parameter values: matter density , baryon density , Hubble constant km s Mpc with , and mass variance on a scale of Mpc. Stellar masses and luminosities (or SFRs) of galaxies are evaluated assuming the Chabrier’s (2003) IMF.

## 2. Continuity equation

Given an evolving population of astrophysical sources, we aim at linking the luminosity function tracing a generic form of baryonic accretion (like that leading to the growth of the central BH or the stellar content in the host galaxies) to the corresponding final mass function . To this purpose, we exploit the standard continuity equation approach (e.g., Small & Blandford 1992; Yu & Lu 2004), in the integral formulation

 N(L,t)=∫∞0dM[∂tN(M,t)−S(M,t)]∑idτidL(L|M,t) ; (1)

here is the time elapsed since the triggering of the activity (the internal clock, different from the cosmological time ), and is the time spent by the object with final mass in the luminosity range given a lightcurve ; the sum allows for multiple solution of the equation . In addition, is a source term due to ‘dry’ merging (i.e., adding the whole mass content in stars or BHs of merging objects without contributing significantly to star formation or BH accretion). In solving Eq. (1) we shall set the latter to zero, and investigate the impact of dry merging in the dedicated App. B. Note that by integrating Eq. (1) in from to the present time , one recovers Eq. (18) of Yu & Lu (2004).

If the timescales (that encase the mass-to-energy conversion efficiency) are constant in redshift and luminosity, then a generalized Soltan (1982) argument concerning the equivalence between the integrated luminosity density and the local, final mass density can be straightforwardly recovered from Eq. (1) without source term by multiplying by and integrating it over and

 ∫t00dt∫∞0dLLN(L,t)=∫∞0dMN(M,t)∫∞0dLL∑idτidL=const×∫∞0dMMN(M,t) , (2)

where the last equivalence holds since const . Specifically, for the BH population the constant is equal to in terms of an average radiative efficiency . We shall see that an analogous expression holds for the stellar component in galaxies.

More in general, Eq. (1) constitutes an integro-differential equation in the unknown function , that can be solved once the input luminosity function and the lightcurve have been specified. Specifically, we shall use it to derive the mass function of the supermassive BH and stellar component in galaxies throughout the history of the Universe. Remarkably, we shall see that Eq. (1) can be solved in closed analytic form under quite general assumptions on the lightcurve.

#### 2.0.1 Connection with standard approaches for BHs

It is useful to show the connection of Eq. (1) with the standard, differential form of the continuity equation for the evolution of the BH mass function as pioneered by Small & Blandford (1992) and then used in diverse contexts by many authors (e.g., Salucci et al. 1999; Yu & Tremaine 2002; Marconi et al. 2004; Shankar et al. 2004, 2009; Merloni & Heinz 2008; Wang et al. 2009; Cao 2010). Following Small & Blandford (1992), BHs are assumed to grow in a single accretion episode, emitting at a constant fraction of their Eddington luminosity erg s in terms of the Eddington time yr. The resulting lightcurve can be written as

 LAGN(τ|MBH,t)=λMBHc2tEdde(τ−τlife)/τef    τ≤τlife; (3)

here is the final BH mass, is the total duration of the luminous accretion phase, and is the folding time in terms of the mass-energy conversion efficiency . Then one has

 dτdLAGN=τefLAGNΘH[LAGN≤LAGN(MBH)] , (4)

where the Heaviside step function specifies that a BH with final mass cannot have shone at luminosity exceeding . Equivalently, only BHs with final masses exceeding can have attained a luminosity and so can contribute to the integral on the right hand side of Eq. (1). Hence such an equation can be written as

 LAGNN(LAGN,t)=∫∞MBH(LAGN)dMBH[∂tN(MBH,t)−S(MBH,t)]τef . (5)

Differentiating both sides with respect to and rearranging terms yields

 ∂tN(MBH,t)+1τef∂MBH[LAGNN(LAGN,t)]|LAGN(MBH)=S(MBH,t) . (6)

Now one can formally write that

 N(LAGN,t)=N(MBH,t)dMBHdLAGNδAGN(MBH,t)          ⟨LAGN⟩=δAGN(MBH,t)LAGN (7)

in terms of the BH duty cycle . Since by definition , one finally obtains the continuity equation in the form

 ∂tN(MBH,t)+∂MBH[⟨˙MBH⟩N(MBH,t)]=S(MBH,t) ; (8)

the underlying rationale is that, although individual BHs turn on and off, the evolution of the BH population depends only on the mean accretion rate .

### 2.1. The BH mass function

We now solve Eq. (1) to compute the BH mass function at different redshifts.

#### 2.1.1 BHs: luminosity function

Our basic input is constituted by the bolometric AGN luminosity functions, that we build up as follows. We start from the AGN luminosity functions at different redshifts observed in the optical band by Richards et al. (2006), Croom et al. (2009), Masters et al. (2012), Ross et al. (2012), Fan et al. (2006), Jiang et al. (2009), Willott et al. (2010a), and in the hard X-ray band by Ueda et al. (2014), Fiore et al. (2012), Buchner et al. (2015), Aird et al. (2015).

Then we convert the optical and X-ray luminosities to bolometric ones by using the Hopkins et al. (2007) corrections111Most of the optical data are given in terms of magnitude at Å. First, we convert them to band ( Å) using the relation (Fan et al. 2001), then we pass to band luminosities in solar units , and finally we go to bolometric luminosities in solar units after . For this last step we recall that the band luminosity of the Sun erg s is about half its bolometric one erg s. In some other instances the original data are expressed in terms of a -corrected band magnitude . We adopt the relation with the magnitude (Richards et al. 2006) and then convert to bolometric as above.. Note that in the literature several optical and X-ray bolometric corrections have been proposed (see Marconi et al. 2004; Hopkins et al. 2007; Shen et al. 2011; Lusso et al. 2012; Runnoe et al. 2012). Those by Marconi et al. (2004) and Lusso et al. (2012) are somewhat smaller by in the optical and by in the hard X-ray band with respect to Hopkins et al. (2007). In fact, since bolometric corrections are intrinsically uncertain by a factor (e.g., Vasudevan & Fabian 2007; Lusso et al. 2012; Hao et al. 2014), these systematic differences between various determinations are not relevant. We shall show in Sect. 2.1.4 that our results on the BH mass function are marginally affected by bolometric corrections. In addition, we correct the number density for the fraction of obscured (including Compton thick) objects as prescribed by Hopkins et al. (2007) for the optical data and according to Ueda et al. (2014; see also Ueda et al. 2003) for the hard X-ray data. We stress that both the bolometric and the obscuration correction are rather uncertain, with the former affecting the luminosity function mostly at the bright end, and the latter mostly at the faint end.

Given the non-homogeneous nature and the diverse systematics affecting the datasets exploited to build up the bolometric luminosity functions, a formal minimum fit is not warranted. We have instead worked out an analytic expression providing a sensible rendition of the data in the relevant range of luminosity and redshift. For this purpose, we use a modified Schechter function with evolving characteristic luminosity and slopes. The luminosity function in logarithmic bins writes

 N(logLAGN,z)=Φ(z)[LAGNLc(z)]1−α(z)exp⎧⎨⎩−[LAGNLc(z)]ω(z)⎫⎬⎭ (9)

The normalization , the characteristic luminosity , and the characteristic slopes and evolve with redshift according to the same parametrization

 p(z)=p0+kp1χ+kp2χ2+kp3χ3 (10)

with and . The parameter values are reported in Table 1. The functional form adopted here is similar to the widely-used double powerlaw shape (e.g., Ueda et al. 2014; Aird et al. 2015), but with a smoother transition between the faint and bright end slopes; all in all, it provides a data representation of comparable quality. In fact, we stress that the results of the continuity equation approach are insensitive to the specific parameterization adopted for the luminosity function and its evolution, provided that the quality in the rendition of the data be similar to ours. For example, in Sect. 2.1.4 we shall show explicitly that our results on the BH mass function are marginally affected when using a double powerlaw shape in place of a modified Schechter to represent the AGN luminosity functions.

In Fig. 1 we illustrate the bolometric AGN luminosity function at various redshifts, including both our data collection and our analytic parameterization of Eq. (9), with an estimate of the associated uncertainty; the extrapolation is also shown for illustration. In the inset we plot the evolution with redshift of the AGN luminosity density, computed as

 ρLAGN(z)=∫dlogLAGNN(logLAGN,z)LAGN (11)

and the contribution to the total by specific luminosity ranges.

#### 2.1.2 BHs: lightcurve

As a further input to the continuity equation, we adopt the following lightcurve (Yu & Lu 2004)

 LAGN(τ|MBH,t) = λ0MBH,Pc2/tEdde(τ−τP)/τef    0≤τ≤τP = λ0MBH,Pc2/tEdde−(τ−τP)/τD    τP≤τ≤τP+ζτD = 0                                   τ>τP+ζτD

This includes two phases: an early one up to a peak time when the BH grows exponentially with a timescale to a mass and emits with an Eddington ratio until it reaches a peak luminosity ; a late phase when the luminosity declines exponentially on a timescale up to a time when it shuts off. With , , we denote the average Eddington ratio and radiative efficiency during the early, ascending phase. The folding time associated to them is .

The lightcurve in Eq. (2.1.2) has been set in Lapi et al. (2014) in order to comply with the constraints imposed by large observational datasets concerning:

• the fraction of X-ray detected AGNs in FIR/K-band selected host galaxies (e.g., Alexander et al. 2005; Mullaney et al. 2012a; Wang et al. 2013a; Johnson et al. 2013);

• the fraction of FIR detected galaxies in X-ray AGNs (e.g., Page et al. 2012; Mullaney et al. 2012b; Rosario et al. 2012) and optically selected quasars (e.g., Mor et al. 2012; Wang et al. 2013b; Willott et al. 2015);

• related statistics via stacking of undetected sources (e.g., Basu-Zych et al. 2013).

The same authors also physically interpreted the lightcurve according to a specific BH-galaxy coevolution scenario.

In a nutshell, the scenario envisages that the early growth of the BH occurs in an interstellar medium rich in gas and strongly dust-enshrouded (Lapi et al. 2014; also Chen et al. 2015). The BH accretes in a demand-limited fashion with values of Eddington ratios appreciably greater than unity, though the radiative efficiency may keep to low values because slim-disc conditions develop. Since the BH mass is still small, the nuclear luminosity, though appreciably super-Eddington, is much lower than that of the starforming host galaxy, and the whole system behaves as a sub-mm bright galaxy with an X-ray nucleus. On the other hand, close to the peak of the AGN lightcurve, the BH mass has grown to large values, and the nuclear emission becomes comparable or even overwhelms that of the surrounding galaxy. Strong winds from the nucleus remove gas and dust from the ambient medium stopping the star formation in the host, while the whole system shines as an optical quasar. If residual gas mass is still present in the central regions, it can be accreted in a supply-driven fashion so originating the declining part of the lightcurve; this phase corresponds to the onset of the standard thin disk accretion, which yields the observed SEDs of UV/optically-selected type-1 AGNs (Elvis et al. 1994; Hao et al. 2014). Actually, the data concerning the fraction of starforming galaxies in optically-selected quasar samples suggest such a descending phase to be present only for luminous objects, while in low-luminosity ones tiny residual mass is present and the AGN fades more drastically after the peak. When the accreting gas mass ends, the BH becomes silent, while the stellar populations in the galaxy evolve passively. For the most massive objects, the outcome will be a local elliptical-type galaxy with a central supermassive BH relic.

All in all, we set the fiducial values of the parameters describing the BH lightcurve on the basis of the Lapi et al. (2014) analysis. We shall discuss the effects of varying them in Sect. 2.1.4. Specifically, we fiducially adopt and for luminous AGNs with peak luminosity , while , i.e., the declining phase is almost absent for low-luminosity objects. To interpolate continuously between these behaviors, we use a standard -function smoothing

 τDτef=3[1−12erfc(12logLAGN1013L⊙)] (13)

which is illustrated in Fig. 3 (bottom panel). We note that our results will be insensitive to the detailed shape of the smoothing function. The value of is fiducially adopted, since after a time after the peak the BH mass has almost saturated to its final value. Results are unaffected by modest variation of this parameter.

We also fiducially assume that the Eddington ratio of the ascending phase depends on the cosmic time (or redshift ) after

 λ0(z)=4[1−12erfc(z−23)] , (14)

as illustrated in Fig. 3 (top panel). As shown by Lapi et al. (2006, 2014), such moderately super-Eddington values at high redshift are necessary to explain the bright end of the quasar luminosity function (see also Shankar et al. 2009, 2013). During the demand-limited, ascending phase of the lightcurve, exceeds the characteristic value for the onset of a slim accretion disc (Laor & Netzer 1989). On the other hand, during the declining phase of the lightcurve, the Eddington ratio declines almost exponentially, so that after the characteristic time the transition to a thin accretion disc takes place. At high redshift where , the thin-disc regime sets in only after a time after the peak, while at low redshift where it sets in about after the peak. We notice that statistically the fraction of slim discs should increase toward high-, as suggested by the data analysis of Netzer & Trakhtenbrot (2014), paving the way for their use as standard candles for cosmological studies (Wang et al. 2013c). The time-averaged value of during the declining phase is, to a good approximation, given by , while the time-averaged value during the thin disk regime ranges from at to at . We will see that such values averaged over the Eddington distribution associated to the adopted lightcurve reproduce well the observational determinations (Vestergaard & Osmer 2009; Kelly & Shen 2013).

As to the radiative efficiency, we take into account the results of several numerical simulations and analytic works (Abramowicz et al. 1988; Mineshige et al. 2000; Watarai et al. 2001; Blandford & Begelman 2004; Li 2012; Begelman 2012; Madau et al. 2014), that indicate a simple prescription to relate the efficiency and the Eddington ratio in both slim and thin disc conditions

 ϵ=ϵthin2λeλ/2−1 ; (15)

here is the efficiency during the thin disc phase, which may range from for a non-rotating to for a maximally rotating Kerr BH (Thorne 1974). We will adopt as our fiducial value (see Davis & Laor 2011). In conditions of mildly super-Eddington accretion with a few the radiative efficiency applies, while in sub-Eddington accretion regime with it quickly approaches the thin disc value . We also take into account that along the declining portion of the lightcurve increases given the almost exponential decrease of . The time averaged values of the efficiency during the declining phase and during the thin disc regime are illustrated in Fig. 3. We expect that the redshift dependence of the average efficiency is negligible during the thin disc regime; this is in qualitative agreement with the findings by Wu et al. (2013) based on spectral fitting in individual type-1 quasars (see also Davis & Laor 2011 for a low- determination), and by Cao (2010) based on continuity equation analysis. However, we caution the reader that the determination of the radiative efficiency is plagued by several systematic uncertainties and selection effects (see discussion by Raimundo et al. 2012). Large samples of AGNs with multiwavelength SEDs and BH masses are crucial in fully addressing the issue.

The final BH mass is easily linked to the mass at the peak appearing in Eq. (2.1.2). One has

 MBH=∫τP+ζτD0dτ′1−ϵϵc2LAGN(τ′)=MBH,P[1+fϵτDτef(1−e−ζ)] . (16)

The correction factor takes into account the modest change of the quantity along the declining phase. We have checked that for any reasonable value of . Notice that at high redshift where , the fraction of BH mass accumulated in thin-disc conditions is only of the total, while it can be as large as at low redshift where . This is relevant since most of the BH mass estimates at high- are based on single-epoch method, that probes the UV/optical bright phase.

The evolution with the internal time of the AGN luminosity, mass, and Eddington ratio are sketched in Fig. 2. We also schematically indicate with colors the stages (according to the framework described below Eq. 2.1.2) when the galaxy is detectable as a FIR-bright source, and the nucleus is detectable as an X-ray AGN and as an optical quasar.

#### 2.1.3 BHs: solution

Given the lightcurve in Eq. (2.1.2), the fraction of the time spent by the BH per luminosity bin reads

 ∑idτidLAGN=τef+τDLAGNΘH[LAGN≤LAGN(MBH)] , (17)

where is the maximum luminosity corresponding to a final BH mass , that can be written as

 LAGN(MBH)=λMBHc2tEdd[1+fϵτDτef(1−e−ζ)]−1 ; (18)

the expression stresses the relevance of the mass accretion during the AGN descending phase. This implies that the time spent in a luminosity bin is longer by a factor than on assuming a simple growing exponential curve, and that Eq. (18) is implicit since is itself a function of the luminosity.

Using Eq. (17) in the continuity equation (neglecting dry merging, i.e., no source term) yields

 LAGNN(LAGN,t)=∫∞MBH(LAGN)dMBH∂tN(MBH,t)[τef+τD] , (19)

where the minimum final mass that have shone at is given by the inverse of Eq. (18). We proceed by differentiating both sides with respect to and rearranging terms to find

 LAGNfBH,LAGN∂LAGN[LAGNN(LAGN,t)]τef+τD=−[∂tN(MBH,t)MBH]|MBH(LAGN) ; (20)

in deriving this equation we have defined , which is not equal to unity since in Eq. (18) depends on .

Finally, we integrate over cosmic time and pass to logarithmic bins. The outcome reads

 N(logMBH,t)=−∫t0dt′fBH,LAGN∂lnLAGN[N(logLAGN)]τef(LAGN,t′)+τD(LAGN,t′)|LAGN(MBH,t′) . (21)

Note that in practice we have started the integration at assuming that the BH mass function at that time was negligibly small. This solution constitutes a novel result. In the case when , and when and are constant with redshift and luminosity, the above equation reduces to the form considered by Marconi et al. (2004).

#### 2.1.4 BHs: results

In Fig. 4 we illustrate our results on the supermassive BH mass function at different representative redshifts. The outcomes of the continuity equation can be fitted by the functional shape of Eq. (9) with replaced by , and with the parameter values reported in Table 1; the resulting fits are accurate to within in the redshift range from to and over the BH mass range from a few to a few .

We also illustrate two determinations of the local mass function. One is from the collection of estimates by Shankar et al. (2009), that have been built by combining the stellar mass or velocity dispersion functions with the corresponding (Haring & Rix 2004) or (Tremaine et al. 2002) relations of elliptical galaxies and classical bulges. The other is the determination by Shankar et al. (2012) corrected to take into account the different relations followed by pseudobulges. In addition, we present the determination at by Vika et al. (2009) based on an object-by-object analysis and on the (McLure & Dunlop 2003) relationship.

The BH mass function at from the continuity equation provides and almost perfect rendition of the local estimates by Shankar et al. (2009) and Vika et al. (2009) when is adopted. At we find a BH mass function which is very similar to the local determination. Our result is in good shape with, though on the high side of, the determination by Li et al. (2011), based on luminosity (or stellar) mass functions and mild evolution of the (or ) relationship. The same also holds at , which is not plotted for clarity.

At we find a BH mass function which at the knee is a factor about below the local data. We are in good agreement with the determination by Ueda et al. (2014) based on continuity equation models. This is expected since we adopt similar bolometric luminosity functions, and around we have similar values of the Eddington ratio and radiative efficiency. At we find a BH mass function which is about orders of magnitude smaller than the local data. We compare our result with the estimate by Willott et al. (2010b) in the range . This has been derived by combining the Eddington ratio distribution from single-epoch BH mass estimates to the optical quasar luminosity functions corrected for obscured objects. At the knee of the mass function we find a good agreement with our result based on the continuity equation, while at lower masses we predict a slightly higher number of objects.

The reasonable agreement with previous determinations in the redshift range validates our prescriptions for the lightcurves, the redshift evolution of , and the relation of Eq. (15). Besides, we recall that these were already independently tested against the observed fractions of AGNs hosted in sub-mm galaxies and related statistics (Lapi et al. 2014; see Sect. 1).

Note that during the slim-accretion regime, where most of the BH mass is accumulated, the effective efficiency amounts to given our assumed value in Eq. (15), see also Fig. 3. This requires a bit of discussion. In principle, during a coherent disk accretion, the BH is expected to spin up very rapidly, and correspondingly the efficiency is expected to attain values (Madau et al. 2014), corresponding to after Eq. (15). However, such a high value of the efficiency would produce a local BH mass function in strong disagreement with the data. This can be understood just basing on the standard Soltan argument. In fact, the BH mass density inferred from the AGN luminosity density would amount to Mpc Mpc. Plainly, the result would fall short of the local observational determinations, that yields a fiducial mass density of Mpc (using the Shankar et al. 2009 local mass function). The discrepancy is even worse if one considers the local mass function obtained by combining the velocity dispersion or stellar mass function with the recently revised or relations by McConnell & Ma (2013) and Kormendy & Ho (2013), which feature a higher overall normalization.

In App. B we have also tested the relevance of dry merging processes (contributing via the source term in the continuity equation) in shaping the BH mass function. At BH merging effects are found to be statistically negligible (see also Shankar et al. 2009), although smaller mass BHs may undergo substantial merging activity with possible impact on the seed distribution (for a review, see Volonteri 2010). At our tests indicate that the mass function is mildly affected only for .

Thus an average slim-disc efficiency is required. During the slim-disc accretion, such a low efficiency can be maintained by, e.g., chaotic accretion, efficient extraction of angular momentum by jets, or similar processes keeping the BH spin to low levels (King & Pringle 2006; Wang et al. 2009; Cao 2010; Li et al. 2012; Barausse 2012; Sesana et al. 2014). We also remark that an efficiency eases the formation of supermassive BHs at very high redshift , so alleviating any requirement on initial massive seeds (Volonteri 2010). On the other hand, the supermassive BH mass function only poorly constrains the values of the BH spin during the final thin-disc phase, which the current estimates suggest to be rather high (Reynolds et al. 2013).

Bolometric corrections and obscured accretion can concur to alleviate the requirement of a low slim-disc efficiency. Bolometric corrections are based on studies of SEDs for large samples of AGNs (e.g., Elvis et al. 1994; Richards et al. 2006; Hopkins et al. 2007; Lusso et al. 2010, 2012; Vasudevan et al. 2010; Hao et al. 2014). In fact, the SEDs depend on the main selection of the objects (e.g., X-ray, UV, optical, IR), possibly on the Eddington ratio (Vasudevan et al. 2010; Lusso et al. 2012), and on bolometric luminosity (Hopkins et al. 2007). The recent analysis of Hao et al. (2014) finds no significant dependencies on redshift, bolometric luminosity, BH mass and Eddington ratio of the mean SEDs for a sample of about X-ray selected type-1 AGNs, although a large dispersion is signalled. A large fraction of objects with accretion obscured at wavelengths ranging from X-ray to optical bands has been often claimed, also in connection with their contribution to the X-ray background (Comastri et al. 1995). The fraction compatible with it at substantial X-ray energies has been recently discussed by Ueda et al. (2014) and properly inserted in our AGN bolometric luminosity functions.

Concerning the overall evolution of the BH mass function, we find that most of the BH mass growth occurs at higher redshifts for the more massive objects (see the inset of Fig. 4). The overall BH mass density at amounts to Mpc, in excellent agreement with observational determinations.

In Fig. 5 we show how our results on the mass function depend on various assumptions. The top and middle panels illustrate the effect of changing the parameters of the lightcurve: radiative efficiency , Eddington ratio , timescale of the descending phase , and duration of the descending phase . For clarity we plot results only at and . We illustrate our fiducial model, and compare it with the outcome for values of the parameters decreased or increased relative to the reference ones.

To understand the various dependencies, it is useful to assume a simple, piecewise powerlaw shape of the luminosity function in the form , with at the faint and at the bright end. Then it is easily seen from Eq. (21) that the resulting mass function behaves as

 N(logMBH)∝1−ϵϵλ1−η[1+(τD/τef)(1−e−ζ)]η1+τD/τefηM−ηBH . (22)

Thus the BH mass function features an almost inverse dependence on at given BH mass. The dependence on is inverse at the high-mass end, which is mostly contributed by high luminosities where . On the other hand, it is direct at the low-mass end, mainly associated to faint sources with . The opposite applies to the dependence on , since roughly . Finally, the dependence on is mild, and practically irrelevant for since the exponential in Eq. (16) tends rapidly to zero. Differences in the results are more evident in the than in the mass function, since this is an integrated quantity, as expressed by Eq. (21).

In the bottom left panel we illustrate the effect of changing the optical/X-ray bolometric corrections: the black lines refer to our reference one by Hopkins et al. (2007), while the blue and red lines refer to the ones proposed by Marconi et al. (2004) and by Lusso et al. (2012), respectively. It is easily seen that the impact on the BH mass function is limited, actually well within the uncertainties associated to the input luminosity functions, and to the observational determinations of the local BH mass function.

In the bottom right panel, we illustrate the effect of changing the functional form used to analytically render the AGN luminosity functions: the black lines refer to our fiducial rendition via a modified Schechter function (cf. Eq. 9), while the green lines refer to a standard double powerlaw representation (e.g., Ueda et al. 2014; Aird et al. 2015). It is seen that our results on the BH mass function are marginally affected; this is because both shapes render comparably well the input AGN luminosity functions.

In Fig. 6 we illustrate the Eddington ratio distribution associated to the overall adopted lightcurve at different redshift and for different final BH masses. Typically at given redshift and BH mass, the distribution features a Gaussian peak at high values of , and then a powerlaw increase toward lower values of before an abrupt cutoff. The peak reflects the value of in the ascending part of the lightcurve. Actually since is constant there, the peak should be a Dirac function. However, small variations around the central value and observational errors will broaden the peak to a narrow Gaussian as plotted here; a dispersion of dex has been safely adopted. The powerlaw behavior reflects the decrease of during the declining part of the lightcurve at late times, and the cutoff in the distribution mirrors that of the lightcurve. The relative contribution of the Gaussian peak at high and of the powerlaw increase at low depends on the relative duration of the declining and ascending phases. Thus at given redshift, small mass BHs feature a much more prominent peak and a less prominent powerlaw increase relative to high mass ones. This is because in small mass objects the descending phase is shorter. At given BH mass, the distributions shift to the left, i.e., toward smaller values of , as the redshift decreases. This is because the initial value decreases with redshift, as prescribed by Eq. (14).

Such a distribution has been computed under the assumption that the overall lightcurve can be sampled. However, from an observational perspective, the Eddington ratio distribution is usually determined via single-epoch BH mass estimates of type-1 AGNs. This implies that only a portion of the descending phase can be sampled. To ease the comparison with observations, we present in the middle and bottom panels of Fig. 6 the expected distribution considering only the descending phase (including both the final portion of the slim-disc phase and the while thin-disc phase, with ), and only the thin-disc phase (i.e., the portion with ). The resulting distributions feature a powerlaw shape, whose slope depends on the portion of the declining phase that can be effectively sampled: the shorter this portion, the steeper the powerlaw. The result is roughly consistent with the observational determinations by, e.g., Kelly & Shen (2013), although a direct comparison is difficult due to observational selection effects. In fact, different observations are likely to sample diversely the initial part of the declining phase, and this will possibly make the expected and the observed distributions even more similar. Note that especially at , BH reactivations, which are not included in our treatment (both in terms of lightcurve descriptions and of stochasticity of the events), can contribute to broaden the Eddington ratio distribution toward very low values of as estimated in the local Universe (e.g., Kauffmann & Heckman 2009; Brandt & Alexander 2015).

In Fig. 7 we present the AGN duty cycle averaged over the Eddington ratio distribution associated to the adopted lightcurve. Specifically, this has been computed as

 ⟨δAGN⟩(MBH,t)≡NAGN(logMBH,t)N(logMBH,t)=1N(logMBH,t)∫dlogλ P(logλ|MBH,z)× (23) ×N(logLAGN,t)|LAGN(MBH,λ)

where is given by Eq. (18). In our approach based on the continuity equation, the duty cycle is a quantity derived from the luminosity and mass functions. It provides an estimate for the fraction of active BHs relative to the total. At given redshift, the average duty cycle increases with the BH mass, since more massive BHs are typically produced by more luminous objects, that feature the descending phase of the lightcurve. On the contrary, small mass BHs are originated mainly by low-luminosity objects for which the descending phase is absent. At given BH mass, the duty cycle increases with the redshift, essentially because to attain the same final mass, BHs stay active for a larger fraction of the shorter cosmic time. This is especially true for BHs with high masses, up to the point that they are always active () for . This agrees with the inferences from the strong clustering observed for high-redshift quasars (Shen et al. 2009; Shankar et al. 2010; Willott et al. 2010b; Allevato et al. 2014); we will further discuss the issue in Sect. 3.1.1. The increase of the duty cycle with BH mass is consistent with the active fraction measured by Bundy et al. (2008), Xue et al. (2010), although the issue is still controversial and strongly dependent on obscuration-corrections (see Schulze et al. 2015). On the other hand, we again remark that our approach does not include AGN reactivations, which may strongly enhance the duty cycle for low-luminosity objects especially at , accounting for the estimates by, e.g., Ho et al. (1997), Green & Ho (2007), Goulding & Alexander (2009), Schulze & Wisotzki (2010).

In Fig. 8 (top panel) we present the AGN Eddington ratio averaged over the lightcurve, computed as

 ⟨λ⟩(MBH,t)≡1N(logMBH)∫dlogλ λP(logλ|MBH,z)N(logLAGN)|LAGN(MBH,λ) . (24)

At given final BH mass, the Eddington ratio decreases with the redshift, as a consequence of the dependence adopted in Eq. (14). The average values are consistent with those observed for a sample of quasars by Vestergaard & Osmer (2009). Note that to take into account the observational selection criteria, we have used the Eddington ratio distribution associated to the descending phase, presented in the middle panel of Fig. 6.

In Fig. 8 (bottom left panel) we show the Eddington ratio function, that has been computed as

 N(logλ,z)≡∫dlogMBH P(logλ|MBH,z)NAGN(logMBH,z) ; (25)

the outcome is mildly sensitive to the lower integration limit, and a value has been adopted to compare with data (see Schulze et al. 2015). For the sake of completeness, we present the results when using the Eddington ratio distribution associated to the thin disk phase (cf. bottom panel of Fig. 6) or to the whole descending phase (cf. middle panel of Fig. 6), being the outcome for the overall lightcurve very similar to this latter case. Our results from the continuity equation are confronted with the estimates from Schulze & Wisotzki (2010) at , and from Schulze et al. (2015) and Nobuta et al. (2012) at , finding a nice agreement within the observational uncertainties and the clear systematic differences between datasets.

In Fig. 8 (bottom right panel) we present a related statistics, i.e., the fraction of galaxies with given stellar mass hosting an AGN (active fraction) with a given Eddington ratio. This has been computed simply on dividing the quantity by the number of galaxies with given stellar mass (the stellar mass function, cf. Sect. 2.2). Plainly, must be set to the BH mass corresponding to , the result being rather insensitive to the ratio adopted; we further take into account a scatter of dex in this relationship, whose effect is to make the active fraction to depend on more weakly than the Eddington ratio distribution depends on . We illustrate the outcome for a range of stellar masses from to ; it turns out to be only mildly dependent on , and especially so at low redshift , as also indicated by current observations.

In fact, our results can be compared with the observational estimates at by Aird et al. (2012) and Bongiorno et al. (2012). The latter authors actually provide the active fraction as a function of the observable quantity ; this can be converted into an Eddington ratio by assuming an X-ray bolometric correction and a value for the ratio. Bongiorno et al. (2012) suggest an overall conversion factor (here cgs units are used for the quantities on the r.h.s.). We also plot their data points when using a conversion (corresponding, e.g., to a larger ratio or a lower bolometric correction), giving more consistency with the determination by Aird et al. (2012).

All in all, our results from the continuity equation are found to be in good agreement with the observational estimates, reproducing their mild dependence on stellar mass and their shape for down to an Eddington ratio a few . On the other hand, AGNs at with tiny accretion rates corresponding to are likely triggered by reactivations, which are not included in our lightcurve, and can contribute to maintain a powerlaw shape of the Eddington ratio function and of the active fraction down to as observed by Aird et al. (2012).

### 2.2. The stellar mass function

Now we turn to the evolution of the stellar mass function from the SFR-luminosity function.

#### 2.2.1 Stars: SFR function from UV and FIR luminosity

SFR can be inferred by lines (mainly Ly and H) and by continuum emission (mainly UV, FIR, radio and X-ray) of star forming galaxies (for a review, see Kennicutt & Evans 2012). The SFR is directly proportional to the UV (chiefly far-UV, FUV) luminosity, which traces the integrated emission by young, massive stars. On the other hand, even a small amount of dust causes large extinction of the UV emission. The absorbed luminosity is re-emitted at longer wavelengths, mostly in the m interval. Therefore ideal estimates would be based on both the UV () and the FIR () observed luminosities for large galaxy samples at relevant redshifts. This would allow to derive the total luminosity proportional to the SFR

 LSFR=LcorrUV=LUV+fLIR ; (26)

here the fraction is meant to subtract from the budget the FIR luminosity contributed by diffuse dust (cirrus) absorbing the light from less massive older stars.

Actually, the SFR functions are inferred from UV or FIR-selected samples. In both cases calibrations and corrections come in (Calzetti et al. 2000; Hao et al. 2011; Murphy et al. 2011; Kennicutt & Evans 2012). The calibration constants between SFR and luminosity in UV and FIR are practically the same, as expected on energy conservation arguments (Kennicutt 1998; Kennicutt & Evans 2012); for FIR luminosity we have

 logSFRM⊙ yr−1≈−9.81+logfLIRL⊙ , (27)

while for extinction-corrected UV luminosity we have

 logSFRM⊙ yr−1≈−7.42−0.4McorrνUV≈−9.76+logνUVLcorrνUVL⊙ , (28)

being the frequency corresponding to Å. 222Some UV data are given at a restframe wavelength different from Å; Eq. (28) still holds provided that on the right hand side the correction Å is added. E.g., for Å the correction amounts to and the zero point calibration becomes ..

The FIR luminosity ascribable to the diffuse dust emission (cirrus) depends on several aspects such as stellar content (mass, age and chemical composition), dust content and spatial distribution. The cirrus emission is characterized by dust temperature lower than the emission associated with star formation in molecular clouds (Silva et al. 1998; Rowlands et al. 2014). There are local galaxies with quite low SFR, whose FIR luminosity is dominated by cirrus emission. E.g., Hao et al. (2011) found for a sample of nearby starforming galaxies with SFR yr. However, the fraction strongly reduces with increasing star formation (e.g., Clemens et al. 2013). For strong local starbursting galaxies with yr and , such as Arp 220, we get a few percent (Silva et al. 1998; Rowlands et al. 2014). Hereafter we will assume that for and that at such large luminosities the SFR can be soundly inferred from the FIR luminosity.

At low luminosity, the SFR is better estimated from UV emission. For this purpose it is essential to allow for dust absorption, that may drastically reduce the UV luminosity to a few percent or even less of its intrinsic value. When only UV data are available, the correlation between the UV slope and the IRX ratio is largely used to infer the dust attenuation (Meurer et al. 1999). While initially proposed only for low redshift galaxies, the method has been tested and applied also to high redshifts (Reddy et al. 2010; Overzier et al. 2011; Hao et al. 2011; Bouwens et al. 2013, 2015). However, for large values of the slope and of the attenuation, the spread around the correlation becomes huge (Overzier et al. 2011; Reddy et al. 2010) and the estimate of attenuation becomes quite uncertain even in local samples (e.g., Howell et al. 2010). On the other hand, the estimate of attenuation for UV-selected samples is less dispersed for galaxies with SFRs yr. In such instances the correction to UV luminosity is more secure and relatively small on average (Bouwens et al. 2013). In fact, this is also suggested by the UV attenuation inferred by combining the H attenuation and the Calzetti extinction curve (Hopkins et al. 2001; Mancuso et al. 2015).

Given all these considerations, we build up the overall SFR-luminosity function as follows. We start from the luminosity function at different redshifts observed in the FIR band by Magnelli et al. (2013), Gruppioni et al. (2013), Lapi et al. (2011), and in the UV band by Bouwens et al. (2015), Oesch et al. (2010), Reddy & Steidel (2009), and Wyder et al. (2005). The data are reported in Fig. 9. In passing note that the SFR and the SFR-luminosity scales on the upper and lower axis have been related assuming the approximate relation , and so the number density per unit SFR or per unit luminosity is the same.

For the FIR samples we assume , while for the UV samples at redshift we have exploited the dust correction suggested by Meurer et al. (1999) and Bouwens et al. (2013, 2015). At the attenuation has been kept to the values found by Bouwens et al. (2013) for galaxies. This assumption at produces an UV attenuation somewhat in between those proposed by Wyder et al. (2005) and Cucciati et al. (2012), and the one proposed by Hopkins et al. (2001). However, we stress that for galaxies with the correction is anyway smaller than a factor .

Fig. 9 shows that at any redshift we lack a robust determination of the SFR-luminosity function at intermediate luminosities; this occurs for two reasons: first, UV data almost disappear above (see also Reddy et al. 2010) because of dust extinction, while FIR data progressively disappear below because of current observational limits. Second, the UV correction for or intrinsic SFR yr becomes progressively uncertain, as discussed above.

To fill in the gap, we render the overall SFR distribution with a continuous function, whose shape is basically determined at the bright end by the FIR data and at the faint end by the UV data. Specifically, we exploit the same modified-Schechter functional shape of Eq. (9), with replaced by and the parameter values reported in Table 1. The UV data at the faint and FIR data at the bright are smoothly connected by our analytic renditions at various redshifts. We also illustrate an estimate of the associated uncertainty. In the inset we plot the evolution with redshift of the SFR-luminosity density and the contribution to the total by specific luminosity ranges.

It happens that our rendition of the data closely follows the model proposed by Mao et al. (2007) and Cai et al. (2014), wherein the extinction is strongly differential with increasing SFR (and gas metallicity). In such models, the faint end of the UV luminosity function at high redshift is dictated by the rate of halo formation, while the bright end is modeled by the dust content in rapidly starforming galaxies. At reliable statistics concern only UV-selected galaxies endowed with low SFR. At high luminosity we have extrapolated the behavior from lower redshift , finding a good agreement with the model proposed by Cai et al. (2014). This extrapolation implies, at , a significant fraction of dusty galaxies with SFR yr, which are missed by UV selection. Clues of such a population are scanty, but not totally missing. Riechers et al. (2014) detected a dust obscured galaxy at with SFR yr, and Finkelstein et al. (2014) a second one at with SFR yr. The large SFR end at will be probed in the near future by ALMA and JWST observations.

In passing, we have also reported the extrapolation of the SFR-luminosity function to (cyan line in Fig. 9). It is interesting to compare this with the recent estimates from UV observations by Bouwens et al. (2015). At the extrapolation matches the observed number density around Mpc. For smaller we remark that the slope of the luminosity function is highly uncertain; data extrapolation suggests a slope in the range from to , as illustrated by the cyan shaded area. At the other end, for , the extrapolated number density is around Mpc, a factor around times larger than that observed in the UV. This possibly suggests that dust already at affects the UV data toward the bright end, as it happens at lower redshift.

#### 2.2.2 Stars: lightcurve

There are three time-honored assumptions regarding the behavior of the SFR as a function of galactic age: exponentially increasing (up to a ceiling value), exponential decreasing, constant.

Here we specialize to a very simple, constant lightcurve, motivated by the recent FIR data from the Herschel satellite concerning high-redshift, luminous starbursting galaxies, and their physical interpretation on the basis of the BH-galaxy coevolution model by Lapi et al. (2014). Hence we adopt

 LSFR(τ|M⋆,t) = κ⋆˙M⋆            τ≤τburst = 0                   τ>τburst

where is a dimensional constant converting SFR into bolometric luminosity. For a Chabrier IMF we have yr erg s yr (see Sect. 2.2.1). The constant SFR represents an average over the fiducial period of the burst , with the total mass of formed stars amounting to .

Since the more massive stars restitute most of their mass to the ISM, the total amount of surviving mass is , where is the restituted fraction that depends on the IMF and on the time elapsed from the burst. For the Chabrier IMF the mass in old, less massive stars approaches to when the time elapsed is larger than a few Gyrs. Since we shall exploit the continuity equation also at relatively short cosmic times at we adopt the value corresponding to Gyr (see below).

The most recent observations by ALMA have undoubtedly confirmed that the SFR in massive high-redshift galaxies must have proceeded over a timescale around Gyr at very high rates a few yr under heavily dust-enshrouded conditions (e.g., Scoville et al. 2014, 2015, their Table 1). The observed fraction of FIR detected host galaxies in X-ray (e.g., Page et al. 2012; Mullaney et al. 2012b; Rosario et al. 2012) and optically selected (e.g., Mor et al. 2012; Wang et al. 2013b; Willott et al. 2015) AGNs points toward a SFR abruptly shutting off after this period of time. In the analysis by Lapi et al. (2014) this rapid quenching is interpreted as due to the energy feedback from the supermassive BH growing at the center of the starbursting galaxy. In the first stages of galaxy evolution the BH is still rather small and the nuclear luminosity is much less than that associated to the star formation in the host. The SFR is then regulated by feedback from SN explosions, and stays roughly constant with time, while the AGN luminosity increases exponentially. After a period Gyr in massive galaxies the nuclear luminosity becomes dominant, blowing away most of the gas and dust from the ambient medium and hence quenching abruptly the star formation in the host. On the other hand, longstanding data on stellar population and chemical abundances of galaxies with final stellar masses indicate that star formation have proceeded for longer times regulated by supernova feedback (see reviews by Renzini et al. 2006; Conroy 2013; Courteau et al. 2014; and references therein).

On this basis, we adopt a timescale for the duration of the starburst given by

 τburst(t)=1Gyr(1+z3.5)−3/2[1+2erfc(43logLSFR1010.5L⊙)] ; (30)

the dependence on the cosmic time mirrors that of the dynamical/condensation time, in turn reflecting the increase of the average density in the ambient medium. In addition, the -smoothing connects continuously the behavior for bright and faint objects expected from the discussion above. We tested that our results are insensitive to the detailed shape of the smoothing function. At high redshift, as noted by Lapi et al. (2014), such a timescale is around folding time of the hosted BH (i.e., Gyr). The luminosity scale corresponds to SFR yr.

#### 2.2.3 Stars: solution

Given the lightcurve in Eq. (2.2.2), the fraction of the time spent per luminosity bin reads

 ∑idτidLSFR=τburstδD[LSFR−LSFR(M⋆)] (31)

with the SFR-luminosity associated to the final stellar mass ; the Dirac delta-function specifies that, since the lightcurve is just a constant, the luminosity associated to a stellar mass must be in the luminosity bin .

Using this expression in the continuity equation Eq. (1) without source term yields

 LSFRN(LSFR,t)τburstf⋆,LSFR=[∂tN(M⋆,t)M⋆]|M⋆(LSFR,t) , (32)

where the final stellar mass that have shone at is given by

 M⋆(LSFR,t)=(1−R)LSFRτburstκ⋆ . (33)

In deriving these equation, we have used . On the same line of Sect. 2.1.3, we integrate over cosmic time, and turn to logarithmic bins. The outcome reads

 N(logM⋆,t)=∫t0dt′[N(logLSFR)f⋆,LSFRτburst]|LSFR(M⋆,t′) . (34)

This solution constitutes a novel result. Note that our approach exploits in the continuity equation the full SFR-luminosity function, and is almost insensitive to initial conditions; in these respects it differs from the technique developed by Leja et al. (2015; see also Peng et al. 2010) to evolve the stellar mass function backwards from basing on the observed SFR