Stellar Mass Function of Active and Quiescent Galaxies via the Continuity Equation
A. Lapi Claudia Mancuso A. Bressan L. Danese




1SISSA, Via Bonomea 265, 34136 Trieste, Italy2INFN-Sezione di Trieste, via Valerio 2, 34127 Trieste, Italy3INAF-Osservatorio Astronomico di Trieste, via Tiepolo 11, 34131 Trieste, Italy
4INAF-IRA, Via P. Gobetti 101, 40129 Bologna, Italy

The continuity equation is developed for the stellar mass content of galaxies, and exploited to derive the stellar mass function of active and quiescent galaxies over the redshift range z∼ 0-8. The continuity equation requires two specific inputs gauged on observations: (i) the star formation rate functions determined on the basis of the latest UV+far-IR/sub-mm/radio measurements; (ii) average star-formation histories for individual galaxies, with different prescriptions for discs and spheroids. The continuity equation also includes a source term taking into account (dry) mergers, based on recent numerical simulations and consistent with observations. The stellar mass function derived from the continuity equation is coupled with the halo mass function and with the SFR functions to derive the star formation efficiency and the main sequence of star-forming galaxies via the abundance matching technique. A remarkable agreement of the resulting stellar mass function for active and quiescent galaxies, of the galaxy main sequence and of the star-formation efficiency with current observations is found; the comparison with data also allows to robustly constrain the characteristic timescales for star formation and quiescence of massive galaxies, the star formation history of their progenitors, and the amount of stellar mass added by in-situ star formation vs. that contributed by external merger events. The continuity equation is shown to yield quantitative outcomes that must be complied by detailed physical models, that can provide a basis to improve the (sub-grid) physical recipes implemented in theoretical approaches and numerical simulations, and that can offer a benchmark for forecasts on future observations with multi-band coverage, as it will become routinely achievable in the era of JWST.








§ INTRODUCTION


Some recent findings have significantly rekindled interest in the field of galaxy formation and evolution. The first concerns the discovery of an abundant population of dusty star-forming galaxies at redshifts z 1, that has been shown to be responsible for the bulk of the cosmic star formation history, in particular around the crucial redshifts z≈ 2-3 where it peaks (e.g., Gruppioni et al. 2013; Rowan-Robinson et al. 2016; Lapi et al. 2017; Bourne et al. 2017; Dunlop et al. 2017; Novak et al. 2017), and to be present even out to z 6 (e.g., Cooray et al. 2014; Riechers et al. 2017; Zavala et al. 2017). Such achievement has become feasible only recently thanks to wide-area far-IR/sub-mm surveys conducted by Herschel, ASTE/AzTEC, APEX/LABOCA, JCMT/SCUBA2, and ALMA-SPT (e.g., Gruppioni et al. 2013, 2015; Lapi et al. 2011; Weiss et al. 2013; Strandet et al. 2016; Koprowski et al. 2014, 2016), in many instances eased by gravitational lensing from foreground objects (e.g., Negrello et al. 2014, 2017; Nayyeri et al. 2016). In fact, galaxies endowed with star formation rates Ṁ_⋆ a few tens M_⊙ yr^-1 at redshift z 2 were largely missed by rest-frame optical/UV surveys because of heavy dust obscuration, difficult to correct for with standard techniques based only on UV spectral data (e.g., Bouwens et al. 2016, 2017; Mancuso et al. 2016a; Pope et al. 2017; Ikarashi et al. 2017; Simpson et al. 2017).

High-resolution, follow-up observations of these galaxies in the far-IR/sub-mm/radio band via ground-based interferometers, such as SMA, VLA, PdBI, and recently ALMA, have revealed star formation to occur in a few collapsing clumps distributed over spatial scales smaller than a few kpcs (see Simpson et al. 2015; Ikarashi et al. 2015; Straatman et al. 2015; Spilker et al. 2016; Barro et al. 2016; Tadaki et al. 2017). A strongly baryon-dominated stellar core with high ongoing SFR is often surrounded out to 15 kpc by a clumpy, unstable gaseous disk in nearly keplerian rotation (e.g., Genzel et al. 2017; Swinbank et al. 2017).

Observations of dusty star-forming galaxies in the optical and near/mid-IR band from Spitzer, WISE, and HST have allowed to characterize their stellar mass content. The vast majority feature stellar masses strongly correlated to the SFR, in the way of an almost linear relationship dubbed 'Main Sequence', with a normalization
steadily increasing as a function of redshift and a limited scatter around 0.25 dex (see Daddi et al. 2007; Rodighiero et al. 2011, 2015; Speagle et al. 2014; Whitaker et al. 2014; Renzini & Peng 2015; Salmon et al. 2015; Tasca et al. 2015; Kurczynski et al. 2016; Tomczak et al. 2016; Bourne et al. 2017; Dunlop et al. 2017; Schreiber et al. 2017).

Another relevant piece of news concerns the discovery by deep near-IR surveys of an increasing number of massive galaxies M_⋆ several 10^10  M_⊙ at high redshift z 2 (see Bernardi et al. 2013, 2017; Ilbert et al. 2013; Duncan et al. 2014; Tomczak et al. 2014; Caputi et al. 2015; Grazian et al. 2015; Thanjavur et al. 2016; Song et al. 2016; Davidzon et al. 2017). Even more interestingly, some of them are found to be already in passive evolution at z 2-3, and to feature chemical properties similar to local early-type galaxies, including a (super)solar metallicity and a pronounced α-enhancement. There is the intriguing yet still debated possibility that the dusty star-forming objects seen in the far-IR/sub-mm band constitute the progenitors of the massive (quiescent) galaxies increasingly detected at high redshifts via deep near-IR surveys (Straatman et al. 2014, 2016; Lonoce et al. 2015; Kriek et al. 2016; Mawatari et al. 2016; Michalowski et al. 2016; Davidzon et al. 2017; Glazebrook et al. 2017).

Relevant model-independent information on the cosmic star formation and mass growth history can be inferred by comparing the observed SFR function, stellar mass function, and main sequence for active and quiescent galaxies (e.g., Leja et al. 2015; Contini et al. 2016; Tomczak et al. 2016; Mancuso et al. 2016a,b; Steinhardt et al. 2017). This procedure can provide stringent constraints, e.g., on the typical timescales for star formation and quiescence, on the overall star formation efficiency, on the initial mass function (IMF), and on the amount of stellar mass added by in-situ star formation vs. that contributed by external merger events. Such outcomes can also be helpful to improve the (sub-grid) physical recipes implemented in theoretical models and numerical simulations, that currently face some difficulties in reproducing the observed abundances of strongly star-forming and massive quiescent galaxies at z 2-3 (e.g., Wellons et al. 2015; Behroozi & Silk 2017; Dave et al. 2017; Furlong et al. 2017; Rong et al. 2017; Hopkins et al. 2017).

In the present paper we pursue the above strategy, for the first time in a quantitative way, by exploiting the specific tool constituted by the `continuity equation'. Being originally devised to connect quasar statistics to the demographics of supermassive black hole relics (Cavaliere et al. 1971; Soltan 1982; Small & Blandford 1992; Salucci et al. 1999; Yu & Lu 2004, 2008; Marconi et al. 2004; Merloni & Heinz 2008; Shankar et al. 2009, 2013; Aversa et al. 2015), here we develop it for the stellar mass content of galaxies, in order to derive the stellar mass function of active and passive galaxies at different redshifts from the SFR functions and average star formation histories for individual objects. Our approach includes in the continuity equation a source term taking into account dry mergers and tidal stripping effects, gauged on observations and on state-of-the-art numerical simulations. With the term dry mergers we refer to events adding the whole mass content in stars of merging objects without contributing significantly to in-situ star formation; starbursts triggered by wet mergers, although included as star-forming objects populating the SFR functions, are expected to contribute little to the final stellar mass, and especially so for massive galaxies.

Moreover, we will exploit the abundance matching technique to derive the star formation efficiency and the main sequence of star-forming galaxies, and compare the outcome to recent observational determinations. Specifically, we will demonstrate via the continuity equation that the dusty, strongly star-forming galaxies at z 2 are indeed the progenitors of massive quiescent galaxies, and that the latter's mass growth is dominated by in-situ star formation with an overall efficiency of less than 20%.

The plan of the paper is straightforward. In  <ref> we describe the basic ingredients of our analysis: the SFR functions and the adopted star formation histories for individual galaxies; in  <ref> we solve the continuity equation for the stellar mass function of active and passive galaxies, and describe how to derive from those the star formation efficiency and the main sequence of star-forming galaxies; in  <ref> we present our results and compare them to observations, discussing the relevant implications for galaxy formation and evolution; in  <ref> we summarize our findings.

Throughout the work we adopt the standard flat cosmology (Planck Collaboration XIII 2016) with round parameter values: matter density Ω_M = 0.32, baryon density Ω_b = 0.05, Hubble constant H_0 = 100  h km s^-1 Mpc^-1 with h = 0.67, and mass variance σ_8 = 0.83 on a scale of 8  h^-1 Mpc. Stellar masses and SFRs (or luminosities) of galaxies are evaluated assuming the Chabrier's (2003) IMF.



§ BASIC INGREDIENTS


Our analysis relies on two basic ingredients: (i) an observational determination of the SFR function at different redshifts; (ii) deterministic evolutionary tracks describing the average star formation history of individual galaxies. In this section we recall the notions relevant for the investigation of the stellar mass function, deferring the reader to the papers by Mancuso et al. (2016a,b) and Lapi et al. (2017) for more details.



 §.§ SFR functions and cosmic SFR density


The first ingredient is constituted by the intrinsic SFR function dN/ dlogṀ_⋆, namely the number density of galaxies per logarithmic bin of SFR [logṀ_⋆,logṀ_⋆+ dlogṀ_⋆] at given redshift z. This has been accurately determined by Mancuso et al. (2016a,b) and Lapi et al. (2017) by exploiting the most recent determinations of the evolving galaxy luminosity functions from (dust-corrected) UV, far-IR, sub-mm, and radio data.



The SFR function can be described as a smooth Schechter function

    dN dlogṀ_⋆(Ṁ_⋆,z) = 𝒩(z) [Ṁ_⋆Ṁ_⋆, c(z)]^1-α(z) e^-Ṁ_⋆/Ṁ_⋆, c(z) ,

with three parameters: the overall normalization 𝒩, the characteristic SFR Ṁ_⋆, c and the faint end slope α. The redshift evolution of each parameter has been measured via an educated fit to the observed data in unitary redshift bins by Mancuso et al. (2016a,b).
As extensively discussed by the latter authors, the SFR function is mainly determined by (dust-corrected) UV data for SFR Ṁ_⋆ 30  M_⊙ yr^-1 since in this range dust emission is mainly due to the diffuse (cirrus) dust component and standard UV dust-corrections based on the UV slope are reliable (see Meurer et al. 1999; Calzetti 2000; Bouwens et al. 2015, 2016, 2017); here we use the Meurer/Calzetti extinction law, but note that switching to a Small Magellanic Cloud (SMC) extinction law affects mildly the SFR function at the faint end (see also Sect. <ref> and Fig. <ref>). On the other hand, the SFR function is mainly determined by far-IR/sub-mm/radio data for SFRs Ṁ_⋆ 10^2  M_⊙ yr^-1 since in this range dust emission is largely dominated by molecular clouds, and UV corrections are wildly dispersed and statistically fail (see Silva et al. 1998; Efstathiou et al. 2000; Coppin et al. 2015; Reddy et al. 2015; Fudamoto et al. 2017).

The resulting SFR functions at representative redshifts are illustrated along with the relevant data collection in Fig. <ref>. In Mancuso et al. (2016a,b; 2017) and Lapi et al. (2017) we have validated them against independent datasets, including integrated galaxy number counts at significative far-IR/sub-mm/radio wavelengths, counts/redshift distributions of strongly gravitationally-lensed galaxies, main sequence of star-forming galaxies and AGNs, redshift evolution of the cosmic SFR, and high-redshift observables including the history of cosmic reionization.

All in all, our determination of the SFR functions implies a significant number density of dusty star-forming galaxies with SFR Ṁ_⋆ 10^2  M_⊙ yr^-1, currently missed by (dust-corrected) UV data. To highlight more clearly this point, in Fig. <ref> we also report at z 1 the SFR function that would have been inferred basing solely on UV data, dust corrected via the UV slope. The UV data considerably underestimate the SFR function for SFRs Ṁ_⋆ 30  M_⊙ yr^-1, because of strong dust extinction. Interestingly, the shape of the SFR function for Ṁ_⋆ 10^2  M_⊙ yr^-1, which so far has been probed only indirectly at z 4 due to sensitivity limits in current wide-areas far-IR surveys, is found to agree out to z 6 with the constraints from the recent VLA-COSMOS radio survey (Novak et al. 2017) and from the few individual galaxies detected at z 5 with ALMA and SMA (e.g., Riechers et al. 2017; Zavala et al. 2017). We shall demonstrate via the continuity equation that a robust probe on the bright end of the SFR function at high-redshift z 4 is provided by the galaxy stellar mass function.



For the analysis in the present paper, we shall assume that at z 1 active galaxies populate the total (UV+far-IR/sub-mm/radio) SFR function, and feature a spheroid-like star-formation history; the latter envisages a nearly constant beahvior of the SFR as a function of galaxy age, with a timescale of 1 Gyr at high SFR Ṁ_⋆ 30  M_⊙ yr^-1, increasing to a few Gyrs for lower SFRs (see Sect. <ref> for details). At high SFRs such a population comprises dusty starforming, far-IR/sub-mm selected galaxies, which will turn out to be the progenitors of local massive dead spheroids with masses M_⋆ a few 10^10  M_⊙; at low SFRs, it comprises mildly obscured, UV selected galaxies (e.g., Lyman break galaxies), that will end up in objects with stellar masses M_⋆ 10^10  M_⊙.

At z 1 we will show that a bimodal star-formation history is required. On the one hand, the bright end of the total (UV+far-IR/sub-mm/radio) SFR function is assumed to be populated by galaxies with the same spheroid-like star-formation history sketched above; since the knee of the SFR function toward z∼ 0 recedes a lot, the star-formation timescales are on average appreciably longer than at z∼ 1, attaining up to a few Gyrs (possibly splitted in many recurrent, shorter bursts); this population comprises a mixed bag of objects, including low-mass spheroids (e.g., bulges), irregulars, and reactivations of massive galaxies. On the other hand, the UV-inferred SFR function is assumed to be populated by galaxies with a disc-like star-formation history, i.e., exponentially declining SFR as a function of galaxy age with long timescales of several Gyrs (see also Cai et al. 2013, 2014); these objects will end up in disc-dominated galaxies with stellar masses M_⋆ several 10^10  M_⊙. In Sect. <ref> we will describe in detail the adopted spheroid-like or disc-like star-formation histories for individual galaxies.

From the SFR function, we can straightforwardly compute the cosmic SFR density as

    ρ_ SFR(z) = ∫ dlogṀ_⋆  dN dlogṀ_⋆ Ṁ_⋆ ,

integrated down to a limit Ṁ_⋆ 10^-1  M_⊙ yr^-1 for fair comparison to observational data, in particular with current blank-field UV surveys at high z 4; the outcome is illustrated in Fig. <ref>. The result from the (dust-corrected) UV-inferred SFR functions is in good agreement with the UV data by Schiminovich et al. (2005) at z≲ 4 and by Bouwens et al. (2015, 2016, 2017) at z≳ 4. It also agrees with the estimate by ALMA observations of UV-selected galaxies in the HUDF (see Dunlop et al. 2017); this is because the rather small area of the HUDF survey allows to pick up only moderately star forming galaxies with mild dust obscuration, on which the UV slope-based corrections still work pretty well.

However, the cosmic SFR density from (dust-corrected) UV data is inconsistent with other datasets both at low and high redshift. Specifically, at redshifts z 4 it falls short with respect to the multiwavelength determination by Hopkins & Beacom (2006) based on UV/optical, radio, Hα and mid-IR 24 μm data, to the far-IR measurements from Herschel by Magnelli et al. (2013) and Gruppioni et al. (2013), and to the recent estimate from deblended data from Herschel, JCMT/AzTEC and JCMT/SCUBA-2 in the GOODS field by Liu et al. (2017). At redshifts z 4 it underestimates the determinations based on stacking of far-IR data from Herschel by Rowan-Robinson et al. (2016), the measurements from radio data by Novak et al. (2017), and the estimates based on long GRB rates from Swift by Kistler et al. (2009; 2013). This mostly reflects the fact, already mentioned above, that the UV-inferred SFR functions (even corrected for dust extinction via the UV slope) appreciably underestimate the number density of dusty galaxies with Ṁ_⋆ 30  M_⊙ yr^-1.

The agreement with all these datasets is substantially improved when basing on the cosmic SFR density computed from the UV+far-IR/sub-mm/radio SFR functions. We also illustrate the contribution to the total density from objects with spheroid-like and disc-like star-formation histories. We remark that at z 1 most of the SFR density is contributed by dusty starforming progenitors of local massive quiescent spheroids, while at z 1 it is contributed both by disc-dominated galaxies and by low-mass spheroids, irregulars, and reactivated massive galaxies.



 §.§ Star-formation history of individual galaxies


The second ingredient of our analysis is constituted by deterministic evolutionary tracks for the history of star formation in individual galaxies. The relevant quantity Ṁ_⋆(τ|M_⋆,t) is the behavior of the SFR as a function of the internal galactic age τ (i.e., the time since the beginning of significant star formation activity) for a galaxy with relic stellar mass M_⋆ at cosmological time t (corresponding to redshift z).

For high z 1 strongly star-forming galaxies (that will turn out to be the progenitors of local dead massive spheroids), we base on the indications emerging from many SED-modeling studies (e.g., Papovich et al. 2011; Smit et al. 2012; Moustakas et al. 2013; Steinhardt et al. 2014; Citro et al. 2016; Cassará et al. 2016); these suggest a slow power-law increase of the SFR Ṁ_⋆∝τ^κ with κ 1 over a timescale τ_ sphe Gyr, then followed by a rapid quenching, at least for massive objects. For the sake of simplicity (cf. Sect. <ref>), here we adopt the law

    Ṁ_⋆(τ|M_⋆,t) = κ+1 1-ℛ M_⋆τ_ sphe (τ/τ_ sphe)^κ Θ_ H[τ≤τ_ sphe] ,

where Θ(·) is the Heaviside step function. The quantity ℛ is the fraction of mass restituted to the interstellar medium by massive stars, computed in the instantaneous recycling approximation; for a Chabrier IMF and star-formation timescales  Gyr, ℛ≈ 0.4 applies.

As to the parameters involved in the above expressions, recent observations by ALMA have shown that in high-redshift galaxies the star formation occurred within a compact region  a few kpcs over timescales τ_ sphe 0.5-1 Gyr at violent rates Ṁ_⋆ a few 10^2  M_⊙ yr^-1 under heavily dust-enshrouded conditions (e.g., Scoville et al. 2014, 2016; Simpson et al. 2015; Ikarashi et al. 2015; Straatman et al. 2015; Spilker et al. 2016; Tadaki et al. 2017). A duration of the main star formation episode τ_ sphe 0.5-1 Gyr in high-redshift dusty star-forming galaxies, which are the candidate progenitors of massive spheroids, is also confirmed by local observations of the α-enhancement, i.e., iron underabundance compared to α elements. This occurs because star formation is stopped, presumably by some form of energetic feedback (e.g., due to the central supermassive black hole), before type Ia surpernova explosions can pollute the interstellar medium with substantial iron amounts (e.g., Romano et al. 2002; Thomas et al. 2005; Gallazzi et al. 2006; for a review see Renzini 2006). Contrariwise, in low-mass spheroidal galaxies with M_⋆ 10^10  M_⊙ data on the age of stellar population and on chemical abundances indicate that star formation has proceeded for longer times, mainly regulated by supernova feedback and stellar winds (see review by Conroy 2013).

On this basis, we parameterize the timescale for the duration of the main SFR episode in objects with spheroid-like star-formation history
as a function of the peak SFR value Ṁ_⋆(τ_ sphe|M_⋆,t)=(κ+1)  M_⋆/(1-ℛ) τ_ sphe via the implicit equation

    τ_ sphe   =   τ_ sphe^+ +τ_ sphe^- 2+τ_ sphe^+ -τ_ sphe^- 2×
    
       ×   tanh[Ṁ_⋆(τ_ sphe|M_⋆,t) 5  M_⊙  yr^-1] ,
    
    τ_ sphe^+    =    0.6  Gyr (1+z 3)^-3/2 ,      τ_ sphe^- = t_z ;

this has to be solved on a grid of M_⋆ and t (or z; see details in Sect. <ref>). The tanh(·) function interpolates smoothly between the short timescale τ_ sphe^+ 1 Gyr for high star-forming galaxies, and a long timescale τ_ sphe^-∼ t_z of the order of the cosmic time for galaxies with low SFRs. Note that in τ_ sphe^+ the dependence on redshift matches that of the dynamical time ∝ 1/√(G ρ)∝ (1+z)^-3/2, in turn following the increase in average density ρ∝ (1+z)^3 of the ambient medium. Our results will be insensitive to the specific shape of the smoothing function. With a similar parameterization, Mancuso et al. (2016b) have been able to reproduce the main sequence of star-forming galaxies at z≈ 2. We recall that at z 1, since the knee of the SFR functions recedes a lot, most of the objects with spheroid-like star-formation histories are characterized by moderate SFR Ṁ_⋆ 10  M_⊙ yr^-1, hence rather long star-formation timescales up to a few Gyrs (cf. Fig. <ref>, bottom panel).

As to the quenching timescale, the observed fraction of far-IR detected host galaxies in X-ray (e.g., Mullaney et al. 2012; Page et al. 2012; Rosario et al. 2012; Barger et al. 2015; Stanley et al. 2015; Harrison et al. 2016) and optically selected AGNs (e.g., Mor et al. 2012; Wang et al. 2013; Willott et al. 2015; Xu et al. 2015; Netzer et al. 2016; Harris et al. 2016) points toward a SFR abruptly stopping, at least in massive galaxies, after τ_ sphe over a short timescale 10^8 yr due to the action of feedbacks (see Lapi et al. 2014). To avoid introducing an additional parameter, in Eq. (<ref>) we truncate the SFR abruptly after τ_ sphe^+; we checked that an exponential quenching over a short timescale τ_ sphe^+/ζ with ζ a few will add an additional, small delay and produce very similar outcomes on the stellar mass function.



On the other hand, in low redshift z 1 disc-dominated galaxies, it is well known that on average star formation declines exponentially as a function of the galactic age, with a long characteristic timescale of several Gyrs; for example, for our Milky Way it amounts to ≈ 6-7 Gyrs (see Chiappini et al. 1997; Courteau et al. 2014; Pezzulli & Fraternali 2016). In view of these classic evidences, we adopt

    Ṁ_⋆(τ|M_⋆,t)    =   1 1-ℛ M_⋆τ_ disc  e^-τ/τ_ disc ,
    
    τ_ disc   =    6  Gyr (1+z 2)^-3/2 .

In Fig. <ref> (top panel) we illustrate an example of the resulting spheroid-like and disk-like star-formation histories; the relative star-formation timescales as a function of SFR and cosmic time are also shown (bottom panel). For both disc- and spheroid-like histories, we assume a dispersion of 0.25 dex around the average star-formation timescales; this value is inspired by the scatter observed in the specific SFRs Ṁ_⋆/M_⋆ (the inverse of a mass doubling time) of active galaxies at different redshift (e.g., Madau & Dickinson 2014), and it will turn out to produce the observed dispersion in the resulting star formation main sequence (cf. Sect. <ref> and Fig. <ref>).

We caveat the reader that in the literature there have been attempts to parameterize with a unique shape the average star-formation history of galaxies. A classic way involves the so-called 'delayed exponential' model Ṁ_⋆(τ)∝τ^κ  e^-τ/τ_⋆, with two parameters κ and τ_⋆ controlling the early powerlaw rise and the late exponential decline. More recently, analogy with the behavior of the cosmic SFR density (see Gladders et al. 2013) and indications from numerical simulations (see Diemer et al. 2017) have suggested a lognormal shape Ṁ_⋆(τ)∝ e^-(lnτ/τ_⋆)^2/2σ_⋆^2/√(2πσ_⋆^2) τ with  the parameters τ_⋆ and σ_⋆ controlling peak time and width. Other descriptions with more complex parametric form have also been proposed based on observations (e.g., Leitner & Kravtsov 2011) or empirical models (e.g., Behroozi et al. 2013; Moster et al. 2013). All these shapes can be useful to describe the star-formation history averaged over the entire population of a galaxy survey; however, chemical and photometric data
require to differentiate between disc-like and spheroid-like star-formation histories, making the parametric models for each class (e.g., see Fig. 2 in Diemer et al. 2017) essentially indistinguishable from our simple adopted shapes. For example, to describe the history of a starforming disc the timescale of the early rise has to be much faster than that of the late decline, to mirror the exponential model of Eq. (<ref>); contrariwise, in a massive spheroid progenitor the SFR must be nearly constant and then abruptly quenched, to mirror the power-law truncated model of Eq. (<ref>).



§ THE CONTINUITY EQUATION


The continuity equation has been originally devised for connecting the AGN statistics to the demographics of both active and dormant supermassive black holes (Cavaliere et al. 1971; Soltan 1982; Small & Blandford 1992; Salucci et al. 1999; Yu & Lu 2004, 2008; Marconi et al. 2004; Merloni & Heinz 2008; Shankar et al. 2009, 2013). Aversa et al. (2015) have been the first to show that it can be also applied to the stellar component in galaxies, to link the evolution across cosmic times of the SFR function to the stellar mass functions. The continuity equation in integral formulation is written


    dN dlogṀ_⋆(Ṁ_⋆,t) = ∫ dlog M_⋆ ∂_t [ dN dlog M_⋆(M_⋆,t)-S(M_⋆,t)]  dτ dlogṀ_⋆(Ṁ_⋆|M_⋆,t) ;


here the term on the l.h.s. is the (known) SFR function, while under the integral on the r.h.s. the first factor is the cosmic time derivative of the (unknown) stellar mass function minus a source term due to dry mergers (i.e., adding the whole mass content in stars of merging objects without contributing significantly to in-situ star formation), and the second factor is the overall time spent by a galaxy in a bin of SFR obtained from the star formation history. The interested reader can find in Aversa et al. (2015) an extended discussion of how and under which hypothesis the standard differential form of the continuity equation is recovered.

In general, the continuity equation above is integro-differential and has to be solved numerically. If the source term due to dry merging is negligible (as it turns out to be indeed for z 1 according to simulations, see Sect. <ref> for details) and the star formation histories have simple shapes like in Eqs. (<ref>) and (<ref>), the continuity equation can be solved analytically along the following lines (see Aversa et al. 2015). First, the time lapses spent by the galaxy in a logarithmic bin of SFR read


    dτ dlogṀ_⋆(Ṁ_⋆|M_⋆,t)    =   1κ (κ+1)^1/κ Ṁ_⋆^1/κ M_⋆^1/κ τ_ sphe^1+1/κ ln(10) Θ_ H[Ṁ_⋆≤(κ+1) M_⋆ (1-ℛ) τ_ sphe] ,
          
       =   τ_ disc ln(10) Θ_ H[Ṁ_⋆≤M_⋆ (1-ℛ) τ_ disc] ,


for galaxies with spheroid- and disc-like star-formation histories, respectively. In both expressions the Heaviside step function Θ_ H(·) specifies the maximum SFR contributing to a given final stellar mass.

Inserting these expressions in the continuity equation Eq. (<ref>), differentiating with respect to Ṁ_⋆ and then integrating over cosmic time yield the closed form solutions


    dN(log M_⋆,t) dlog M_⋆   =    -κ (1+κ)^1/κ  M_⋆^1/κ ∫_0^t dt' ∂_lnṀ_⋆
    f_τ_ sphe [Ṁ_⋆^-1/κ τ_ sphe^-1-1/κ  dN dlogṀ_⋆(Ṁ_⋆,t')]_|Ṁ_⋆ = (1+κ) M_⋆ (1-ℛ) τ_ sphe ;
          
       =    -∫_0^t dt' ∂_lnṀ_⋆
    f_τ_ disc [τ_ disc^-1  dN dlogṀ_⋆(Ṁ_⋆,t')]_|Ṁ_⋆ = M_⋆ (1-ℛ) τ_ disc ,


again for galaxies with spheroid- and disc-like star-formation histories, respectively; in both expression we have used the shorthand f_τ≡ 1+∂_logṀ_⋆logτ, which is not trivially equal to one when τ depends explicitly on the SFR (as in Eq. <ref>). The above equation is numerically solved on a grid in M_⋆ and z (or cosmic time t_z). We use a grid of 100 equally-spaced points in log M_⋆ [M_⊙] within the range [8,13] and a grid of 1000 equally-spaced points in redshift z within the range [0,20]; for optimal interpolation, the SFR functions and the star-formation timescales have been defined on the same grid of redshift and on a grid of 100 equally-spaced points in logṀ_⋆ [M_⊙ yr^-1] within the range [-2,4].



 §.§ Dry merging


In presence of mergers, the source S(M_⋆,t)=S_+-S_- actually includes the difference between a creation S_+ and a destruction S_- term. The former depends on the merger rate of objects with smaller masses into the descendant mass M_⋆, while the latter depends on the merger rates of the mass M_⋆ into more massive objects. Given the merger rate dN_ merg/ dlog M_⋆  dμ  dt for the production of a descendant mass M_⋆ by the merging of two progenitors with (smaller to higher) mass ratio μ, the creation term reads

    S_+(M_⋆,t) = 1 2 ∫_μ_ min^1 dμ  dN_ merg dlog M_⋆  dμ  dt(M_⋆,μ,t)

while the destruction term is written

    S_-(M_⋆,t)    =   1 2∫_μ_ min^1 dμ [ dN_ merg dlog M_⋆  dμ  dt(M_⋆(1+μ)/μ,μ,t)+.
    
       +   .  dN_ merg dlog M_⋆  dμ  dt(M_⋆ (1+μ),μ,t)] .

In the above μ_ min is the minimum progenitors' mass ratio; typically μ_ min=0.3 includes only 'major mergers', 0.1 includes major and minor mergers, 0.1 practically includes all mergers. We take μ_ min=0.01 in the following.

We base on the outcomes of the Illustris simulations by Rodriguez-Gomez et al. (2015, 2016), who provide a handy fitting function for the merger rate per descendant galaxy

    dn_ merg dμ  dt(M_⋆,μ,t) = A(z) (M_⋆ 10^10  M_⊙)^ω(z)×
    
          ×[1+(M_⋆ 2× 10^11  M_⊙)^δ(z)] μ^β(z)+γ log(M_⋆/10^10  M_⊙)

where

    A(z)   =   A_0 (1+z)^η ,    ω(z)=ω_0 (1+z)^ω_1 ,
    
    β(z)   =   β_0 (1+z)^β_1,    δ(z)=δ_0 (1+z)^δ_1

with A_0≈ 10^-2.2287 Gyr^-1, η≈ 2.4644, ω_0≈ 0.2241, ω_1≈-1.1759, β_0≈-1.2595, β_1≈ 0.0611, γ≈ -0.0477, δ_0≈ 0.7668, δ_1≈ -0.4695. The above authors have validated this expression against various datasets, including observations of galaxy pairs.
Two remarks are in order here. First, we caveat that at z 1 the Illustris simulation does not perfectly reproduce the observed galaxy stellar mass function; however, this does not concern much the merger rates, since as we shall demonstrate the growth in stellar mass at high redshift is mainly dominated by in situ star formation. Second, our results will turn out to be robust against other choices of the merger rate; we checked this by exploiting the galaxy merger rates extracted from the hydrodynamic simulations by Stewart et al. (2009), and the halo merger rates based on the N-body simulation by Fakhouri et al. (2010) coupled with empirical relationships connecting halo and stellar mass (e.g., Moster et al. 2013).

Multiplying the merger rates above by the stellar mass function yields the quantity

    dN_ merg dlog M_⋆  dμ  dt(M_⋆,μ,t)   =    dn_ merg dμ  dt(M_⋆,μ,t)×
    
       ×    dN dlog M_⋆  dt(M_⋆,t)

entering the source terms in Eqs. (<ref>) and (<ref>). As discussed by Rodriguez-Gomez et al. (2015, 2016) the above expression also takes into account stellar mass stripping from satellites prior to dry mergers. Plainly, when merging is introduced, the continuity equation becomes fully integro-differential and must be solved numerically. We have computed the full solution and found that, to a good approximation, one can solve the problem iteratively, using Eq. (<ref>) as the zero-th order solution and then updating it with the correction due to the merging terms.



 §.§ Cosmic stellar mass density and Soltan argument


Once the redshift-dependent stellar mass function is known, the cosmic stellar mass density is obtained as

    ρ_M_⋆(t) ≡∫ dlog M_⋆  M_⋆  dN dlog M_⋆(M_⋆,t) ,

where the integration is typically performed over stellar masses above 10^8  M_⊙ for fair comparison with observational determinations (see discussion by Madau & Dickinson 2014).

Interestingly, if τ_ sphe,disc is independent of, or only weakly dependent on Ṁ_⋆, a Soltan (1982) argument holds for the stellar content of galaxies (see Aversa et al. 2015); classically, this connects the cosmic luminosity density to the relic mass density of a population via an average conversion efficiency. In the present context, the Soltan argument can be easily found by multiplying both sides of Eq. (<ref>) by Ṁ_⋆ and integrating over it and over cosmic time, to obtain

    ρ_M_⋆ = (1-ℛ) ∫_0^t  dt' ∫ dlogṀ_⋆ Ṁ_⋆  dN dlogṀ_⋆(Ṁ_⋆,t') .

We highlight that here the IMF-dependent factor 1-ℛ plays the role of the radiative efficiency in the classic Soltan argument for black holes. Note that for conventional IMFs most of the stellar mass in galaxies resides in stars with mass 1  M_⊙; since these stars emit most of their luminosity in the near-IR, the galaxy stellar mass M_⋆ can be inferred by the near-IR luminosity functions. On the other hand, the SFR function is determined from UV and far-IR/sub-mm/radio observations, as discussed in Sect. <ref>. Thus in principle accurate determinations of the SFR and stellar mass functions in these independent manners could be exploited via the Soltan argument above to constrain the average galaxy IMF at different redshifts. In practice, however, the dependence on the IMF is weak, and current observational uncertainties do not allow to fulfill this program now.



 §.§ Star formation efficiency


We now connect the stellar mass function to the underlying, gravitationally dominant DM component, with the aim of deriving the star formation efficiency f_⋆≡ M_⋆/f_ b  M_ H. This represents the fraction of the baryonic mass f_ b  M_ H≈ 0.16  M_ H initially associated to a DM halo of mass M_ H that has been eventually converted into stars. To this purpose, we exploit the abundance matching technique, a standard way of deriving a monotonic relationship between galaxy and halo properties by matching the corresponding integrated number densities (e.g., Vale & Ostriker 2004; Shankar et al. 2006; Moster et al. 2013; Behroozi et al. 2013).



For fair comparison with the determination of the star formation efficiency f_⋆≡ M_⋆/f_ b ⟨ M_ H⟩ via weak gravitational lensing (e.g., Velander et al. 2014; Hudson et al. 2015; Mandelbaum et al. 2016) and galaxy kinematics (e.g., More et al. 2011; Wojtak & Mamon 2013), that are based on galaxy samples selected by stellar mass, we aim at deriving the average halo mass ⟨ M_ H⟩ (M_⋆,z) associated to a given M_⋆. In the abundance matching formalism, this relationship is obtained via the equation (see Aversa et al. 2015 for details)

    ∫_log⟨ M_ H⟩ (M_⋆,z)^∞ dlog M_ H'  dN dlog M_ H(M_ H',z) = 
    
       =   ∫_-∞^+∞ dlog M_⋆'  dN dlog M_⋆(M_⋆',z) 1 2  erfc{log[M_⋆/M_⋆']√(2) σ_log M_ H} ,

holding when a lognormal distribution of M_ H at given M_⋆
with dispersion σ_log M_ H is assumed.

We follow previous studies based on various semi-empirical methods of galaxy and halo connection (see Rodriguez-Puebla et al. 2015, their Fig. 10) and adopt σ_log M_ H≈max[0.05,0.05+0.15 (log M_⋆ [M_⊙]-10)] for log M_⋆ [M_⊙] within the range [8.5,12.5]. In Eq. (<ref>) the quantity dN/ dlog M_ H is usually taken as the halo mass function from N-body simulations (e.g., Tinker et al. 2008; Watson et al. 2013; Bocquet et al. 2016; Comparat et al. 2017), that includes galaxy groups and clusters. This is particularly suitable when comparing with observational determinations of the star formation efficiency based on weak gravitational lensing (see references above), that integrate all the DM mass along the line of sight, including that associated to the surrounding galaxy environment.

However, in order to infer the star formation efficiency of individual galaxies, and not of a galaxy system like a group or a cluster, it would be more appropriate to use the galaxy halo mass function, i.e., the mass function of halos hosting one individual galaxy. This can be built up from the overall halo mass function by adding to it the contribution of subhalos,
and by probabilistically removing from it the contribution of halos corresponding to galaxy systems via halo occupation distribution modeling. We defer the reader to Appendix A of Aversa et al. (2015) for details on such a procedure.





 §.§ Galaxy main sequence


The vast majority of galaxies is endowed with stellar masses strongly correlated to the ongoing SFR, in the way of an almost linear relationship dubbed 'Main Sequence', with a normalization steadily increasing as a function of redshift, and with a limited scatter around 0.25 dex (see Daddi et al. 2007; Rodighiero et al. 2011, 2015; Speagle et al. 2014; Whitaker et al. 2014; Renzini & Peng 2015; Salmon et al. 2015; Tasca et al. 2015; Kurczynski et al. 2016; Tomczak et al. 2016; Bourne et al. 2017; Dunlop et al. 2017; Schreiber et al. 2017).

We exploit the abundance matching between the SFR functions (Fig. <ref>) and the stellar mass functions (Fig. <ref>) self-consistently derived from the continuity equation (cf. Eq. <ref>) to compute the average SFR ⟨Ṁ_⋆⟩(M_⋆,z) associated to a given stellar mass M_⋆. This reads

    ∫_log⟨Ṁ_⋆⟩ (M_⋆,z)^∞ dlogṀ_⋆'  dN dlogṀ_⋆(Ṁ_⋆',z) =
    
       =   ∫_-∞^+∞ dlog M_⋆'  dN dlog M_⋆(M_⋆',z) 1 2  erfc{log[M_⋆/M_⋆']√(2) σ_logṀ_⋆} ,

holding when a lognormal distribution of Ṁ_⋆ at given M_⋆ with dispersion σ_logṀ_⋆≈ 0.15 dex is adopted (see Aversa et al. 2015).

The comparison of the resulting main sequence with the observational data will actually constitute an additional constraint on the assumed star formation histories for individual galaxies (see Eqs. <ref> and <ref>), on the star formation timescales and the associated scatter, and on the robustness of our results to other aside assumptions discussed in previous sections.



§ RESULTS


In Fig. <ref> we present the stellar mass function at different redshifts obtained via the continuity equation, including both in situ star formation and (dry) mergers. We highlight the average result as solid lines, and the 1σ dispersion expected from the scatter in the star formation timescales and merging histories as shaded areas (see Sect. <ref> for details). We compare our results to recent observational data (Moffett et al. 2016; Thanjavur et al. 2016; Bernardi et al. 2017; Davidzon et al. 2017; Tomczak et al. 2014; Grazian et al. 2015; Song et al. 2016), finding an excellent agreement.

We stress that in-situ star formation within galaxies dominates over dry mergers in building up the stellar mass function at high redshifts, all the way down to z∼ 1, while at lower redshifts z 1 dry mergers can contribute appreciably to the stellar mass growth.  This is highlighted on comparing the solid and dashed lines in Fig. <ref>, which illustrate the mass function at different redshifts when including or not dry mergers, respectively (actually the two sets of curves are superimposed for z 1.5). The effect of dry mergers on the stellar mass function is twofold: the number of low mass galaxies is decreased appreciably because of the merging into larger units, whereas the high-mass end of the stellar mass function is boosted toward larger masses because of mass additions from smaller objects; dry mergers mainly affect the most massive galaxies that are typically dominated by the spheroidal component. Such a picture is in agreement with what is recently emerging from state-of-the-art numerical simulations (see Schaye et al. 2015; Rodriguez-Gomez et al. 2016), semiempirical models (see Behroozi et al. 2013) and analysis of observations based on density-matching arguments (see Hill et al. 2017).



The dotted lines show the contribution to the stellar mass function from galaxies with spheroid-like star-formation history that featured SFRs Ṁ_⋆ 100  M_⊙ yr^-1; this corresponds to the limiting value currently sampled in wide-area far-IR surveys out to z 4 (e.g., the Herschel-ATLAS, see Lapi et al. 2011). It is seen that the descendants of these galaxies populate the high-mass end of the local stellar mass function, and thus are mainly present-day massive dead spheroids (e.g., Moffett et al. 2016). This demonstrates on a statistical basis that strongly starforming galaxies observed in the far-IR/(sub-)mm band constitute the progenitors of massive spheroids. By the same token, we stress that to test at z 4 the outcomes of the continuity equation, and better constrain the input SFR functions and the parameters of the star-formation history for spheroid progenitors, it will be extremely relevant to improve the accuracy in the determination of the stellar mass function at the high-mass end for M_⋆ a few 10^10  M_⊙ out to z 6 via wide-area near-IR surveys.



In Fig. <ref> we focus on the stellar mass function of galaxies with disc-like star-formation histories at z≈ 0. Our result is in excellent agreement with the observed stellar mass function of disc-dominated galaxies from decomposed data (Moffett et al. 2016; Bernardi et al. 2017); thus we find a good correspondence between objects populating the UV-inferred SFR function, to which we assigned disc-like star formation histories, and galaxies with observed disc-dominated morphology in the stellar mass function. We highlight that discs contribute considerably to the total stellar mass function for stellar masses M_⋆ a few 10^10  M_⊙, and that the effects of mergers on their stellar mass function are negligible. It is seen from Eq. (<ref>) that the observed steepness for M_⋆ a few 10^10  M_⊙ in the local stellar mass function of disc-dominated galaxies mirrors that in the input UV-inferred SFR functions at z 1. We caveat that assigning a disc-like star-formation history (with long star-formation timescales) even to objects populating the UV+far-IR/sub-mm/radio SFR functions at z 1 would considerably overproduce the number of massive discs; this is because the UV+far-IR/sub-mm/radio functions are much higher than the UV-inferred ones at given SFR. In fact, the UV+far-IR/sub-mm/radio SFR functions at z 1  must be populated by objects with spheroid-like star-formation history; the continuity equation shows these star formation events to change little the total stellar mass function at z∼ 0 with respect to that at z∼ 1, mildly affecting the number density of galaxies with M_⋆ 10^10  M_⊙.

In Fig. <ref> we show how our resulting stellar mass
function depend on the input SFR function and on the parameters of the star formation history. To highlight such dependencies in simple terms it is convenient to assume a piecewise powerlaw shape of the SFR function dN/ dlogṀ_⋆∝Ṁ_⋆^-χ, with χ 1 at the faint and χ>1 at the bright end. Then it is easily seen from Eq. (<ref>) that the resulting stellar mass function (in absence of mergers) behaves as

    dN dlog M_⋆   ∝   1+κ χ (1+κ)^χ  (1-ℛ)^χ+1/κ M_⋆^-χ τ_ sphe^χ-1
          
       ∝   χ  (1-ℛ)^χ  M_⋆^-χ τ_ disc^χ-1 ,

for galaxies with spheroid-like and disc-like star formation histories, respectively. Thus, the stellar mass function features an almost direct dependence on the star-formation timescales τ_ sphe,disc at the high-mass end, which is mostly contributed by high SFRs where χ> 1; on the other hand, the dependence is inverse but mild at the low-mass end, mainly contributed by low-SFR galaxies with χ 1; note, however, that the value of the SFR where χ appreciably exceeds unity is much lower for the UV-inferred than for the UV+far-IR/sub-mm/radio SFR functions. The dependence on the parameter κ entering the star-formation history Ṁ_⋆(τ)∝τ^κ is mild, direct at the low-mass and inverse at the high-mass end. The dependence on the IMF is encapsulated in the restituted fraction 1 - ℛ, and in the factor used to convert the observed UV+far-IR/sub-mm/radio luminosity function into the SFR function; e.g., passing from the Chabrier to the Salpeter (1955) IMF, the high mass end of the stellar mass function is increased somewhat, while a strong suppression is originated when basing on a top-heavy IMF (e.g., Lacey et al. 2010). Finally, adopting a SMC extinction law in place of the Calzetti for the determination of the input SFR function amounts to alter somewhat the exponent χ and thus changes little the final outcome on the stellar mass function; on the other hand, adopting the UV-inferred SFR function at any redshift (i.e., neglecting far-IR/sub-mm/radio data) would imply to strongly underestimate the stellar mass function for large stellar masses (see discussion by Mancuso et al. 2016a,b) that are indeed built up in dusty star-forming galaxies with violent SFRs.



In Fig. <ref> we focus on the stellar mass function of quiescent (passively evolving) galaxies; these systems have been increasingly observed with appreciable number density out to high redshift z 4 after selection via color-color diagrams in deep near-IR surveys (see Tomczak et al. 2014; Davidzon et al. 2017; Lonoce et al 2017; Glazebrook et al. 2017). Typically, these selections tend to pick up galaxies that have been quenched since, and then passively evolving over, a quiescence time interval Δ t_ qui∼ 250-500 Myr. Thus we compute the associated stellar mass function from the continuity equation by replacing the upper limit of integration in Eq. (<ref>) with τ_ sphe-Δ t_ qui. The result for two different values of Δ t_ qui≈ 250 and 500 Myr encompasses very well with the observational determinations out to z 4. Plainly, higher quiescence time Δ t_ qui imply a lower mass function, especially toward higher redshift where the cosmic time is smaller and progressively closer to τ_ sphe. The decrease of the mass function at the low mass end is due to the fact that small galaxies are still actively forming stars, since they feature longer star formation timescales; as a consequence, the fraction of galaxies in passive evolution decreases rapidly with stellar mass. The downturn shifts toward larger masses toward higher z, passing from 10^10 to a few 10^10  M_⊙ from z≈ 0 to z 3.



In Fig. <ref> we show the cosmic stellar mass density, obtained according to Eq. (<ref>). Our result from integrating the overall stellar mass function from the continuity equation is compared with the data collection by Madau & Dickinson (2014) and with the recent estimates by Song et al. (2016) and Davidzon et al. (2017) at high redshift. The agreement between our results and the data is remarkably good. We also highlight the contribution to the total stellar mass density from galaxies with disc-like and spheroid-like star-formation histories; the latter dominates the overall mass density at any redshift, though at z 1 discs brings an appreciable contribution around 40%. Note that the fraction of galaxies that are quiescent (here we use Δ t_ qui≈ 250 Myr, see above discussion), which are basically massive spheroids, constitute only a fraction 50% of the total mass density (contributed also by small spheroids/irregulars and discs that are still active) in the local Universe and rapidly declines to values 10% at higher redshift z 2. The overall shape agrees well with the estimates by Muzzin et al. (2013), Straatman et al. (2014) and Davidzon et al. (2017).

For reference, in the figure we also report the mass density (scaled down by a factor 10^-2) of galaxy halos with mass M_ H 10^8.5  M_⊙, the minimal threshold for efficient star formation required to solve the missing  satellite problem (see Boylan-Kolchin et al. 2014; Wetzel et al. 2016; Lapi et al. 2017). The evolution in halos and in the stellar mass content of galaxies differs both in shape and in normalization; these differences stem from: the inefficiency of galaxy formation due to feedback processes (e.g., supernovae, stellar winds, active galactic nuclei); the decrease at high z in the number density of massive halos, that are the hosts of the most massive galaxies; the inability to grow massive galaxies at high redshift since the growth timescales become comparable to the age of the Universe.



In Fig. <ref> we show the star-formation efficiency, computed according to Eq. (<ref>). The green solid line and shaded area illustrate the outcome when matching the total stellar mass function from the continuity equation to the overall halo mass function at z≈ 0; note that the shaded area takes into account the uncertainty in the determination of the local stellar mass function from the continuity equation and that arising from the rather flat shape of the average ⟨ M_ H⟩ (M_⋆) correlation at the high mass end. Our result is compared with the local data for early and late type galaxies by various authors, determined via weak lensing (see Mandelbaum et al. 2016; Velander et al. 2014; Rodriguez-Puebla et al. 2015; Hudson et al. 2015) and satellite kinematics (see More et al. 2011; Wojtak & Mamon 2013). We stress that the abundance matching results must be confronted with the data of spheroidal galaxies for stellar masses above, and with the data of disc-dominated galaxies below, a few 10^10  M_⊙; this is because spheroids and discs mostly contribute to the local stellar mass function in such stellar mass ranges (see Moffett et al. 2016). Provided that, we find a very good agreement, within the uncertainties of the respective datasets.

The dotted green line is instead the outcome when matching the local stellar mass function with the galactic halo mass function. This highlights that the decrease in star-formation efficiency at large stellar masses is somewhat spurious, being related to the fact that the most massive galaxies tend to live at the center of group/cluster halos, which contain a lot of DM. Considering instead only the DM mass belonging to individual galactic halos would imply the efficiency to stay almost constant or increase somewhat at large masses out to M_⋆∼ 10^12  M_⋆.

The resulting values and shape of the star formation efficiency as a function of stellar mass is easily understood in terms of feedback processes. It is apparent that, because of feedbacks, galaxy formation is a very inefficient process: at most 20-30% of the original baryonic content of halos is converted into stars; this occurs for galaxies with final stellar mass around a few 10^10  M_⊙ (corresponding to halos with mass M_ H≈ 10^12  M_⊙). At small stellar masses, the action of supernova feedback is predominant, while for large stellar masses AGN feedback is likely more relevant; the mass of maximum efficiency corresponds approximately to the transition between supernova and AGN feedbacks (see Shankar et al. 2006; Moster et al. 2013; Aversa et al. 2015).

The green dashed line is the outcome of matching the stellar mass function and the overall halo mass function at z≈ 2, and compares well with the efficiencies measured at the same redshift from Hα observation by Burkert et al. (2016). The outcome is also similar, within a factor of 2, to the determination via abundance matching by Moster et al. (2013), Behroozi et al. (2013), and Aversa et al. (2015). The similarity of the efficiency at z≈ 2 to the local value is indicative that star formation is mainly an in-situ process (see Lilly et al. 2013; Moster et al. 2013; Aversa et al. 2015; Mancuso et al. 2016a).

We also present as a red solid line the outcome of matching the overall halo mass function with the stellar mass function of passively evolving galaxies, again finding a pleasingly agreement with the local data for spheroidal galaxies. It is extremely interesting to notice that the global efficiency at z≈ 2 can be brought on the efficiency of quiescent galaxies at z≈ 0 by allowing:  (i) an evolution of the stellar mass by a factor 50% due to late star formation or dry mergers (see Rodriguez-Puebla et al. 2017; also Sect. <ref>); (ii) an halo mass evolution by a factor 4 (M_ H/10^14  M_⊙)^0.12 due to late smooth accretion or tidal stripping (see McBride et al. 2009; Fakhouri et al. 2010; Lapi et al. 2013). We stress that such result is again indicative of the in-situ nature of the star formation in spheroid progenitors, and is also extremely relevant for understanding the evolution of the specific angular momentum in galaxies (see Shi et al. 2017).

In Fig. <ref> we show the main sequence of star-forming galaxies at different redshifts, as obtained by matching the SFR function and the stellar mass functions from the continuity equation after Eq. (<ref>). The outcome at z≈ 2 is in pleasing agreement with the observational determination from large statistics of mass-selected galaxy samples by Rodighiero et al. (2015).
This further substantiate our assumed star formation histories for individual galaxies, that are illustrated on three representative cases by red dotted lines; their shape is dictated by the slowly increasing SFR Ṁ_⋆∝τ^1/2 and appreciably rising stellar mass M_⋆∝τ^3/2, which imply Ṁ_⋆∝ M_⋆^1/3. Then the main sequence corresponds to the portions of such tracks where galaxies spend most of their lifetime in
logarithmic bins of M_⋆.

To highlight the relevance of observational selections different from that based on stellar mass, in Fig. <ref> we also report data points for individual, far-IR selected galaxies by Koprowski et al. (2016), Ma et al. (2015b), Negrello et al. (2014), along with Dye et al. (2015), da Cunha et al. (2015), and Dunlop et al. (2017) mainly at redshifts z∼ 1-4.
An appreciable fraction of the individual, far-IR selected galaxies around z≈ 2 (highlighted in red) lie above the main sequence, i.e., at SFR values higher than expected on the basis of the average relationship at given M_⋆. These off-main-sequence objects can be simply interpreted (see Mancuso et al. 2016b) as galaxies caught in an early evolutionary stage, and still accumulating their stellar mass. Thus young star-forming galaxies are found to be preferentially located above the main sequence or, better, to the left of it. As time goes by and stellar mass increases, the galaxy moves toward the average main sequence relationship, around which it will spend most of its lifetime. Afterwards, the SFR is quenched by feedbacks and the galaxy will then evolve passively to become a local early-type; then it will populate a region of the SFR versus stellar mass diagram substantially below the main sequence. These loci of 'red and dead' galaxies are indeed observed locally (see Renzini & Peng 2015), and start to be pinpointed even at high redshift (see Man et al. 2016).



§ SUMMARY


We have developed the continuity equation for the stellar mass content of galaxies, and have exploited it to derive the stellar mass function of active and quiescent galaxies at redshifts z∼ 0-8 from the observed SFR functions and disc-like or spheroid-like star-formation histories for individual galaxies. Our approach based on the continuity equation includes a source term due to dry merging gauged on state-of-the-art numerical simulations and consistent with observations. We have then used the abundance matching technique to investigate the star formation efficiency and the main sequence of star-forming and quiescent galaxies. By comparing these outcomes to current observational estimates, we have inferred constraints on the characteristic timescales for star formation and quiescence, on the overall star formation efficiency, and on the amount of stellar mass added by in-situ star formation vs. that contributed by external (dry) merger events.

Our main findings are the following:




  * We have found that the stellar mass function computed from the continuity equation is in excellent agreement with current observational constraints in the extended redshift range z∼ 0-8.  At high redshift z 1 the mass function is produced by galaxies with spheroid-like star-formation histories, featuring an approximately constant (or slowly increasing) behavior of the SFR as a function of galactic age; the SFR must last for a time τ_ sphe≈ fraction of Gyr in strongly star-forming galaxies, while it can proceed over a longer time interval up to a few Gyrs for less massive objects: this reflects the differential action of supernova and AGN feedbacks in systems with different mass. We stressed the relevance of using as input of the continuity equation the SFR function estimated from far-IR/sub-mm/radio, in addition to UV, observations. This is because strongly star-forming galaxies are heavily dust-enshrouded, and as such their intrinsic SFR is considerably underestimated by UV observations, even when corrected for dust extinction according to standard prescriptions based on the UV slope. We have highlighted that the mass growth of spheroids is dominated by in-situ star formation for z 1, while at lower redshift dry mergers contribute a mass budget 50% especially in the most massive objects.


  * At low redshift z 1, we have shown that the stellar mass function of disc-dominated galaxies is well reproduced in our approach when using as input the UV-inferred SFR functions and an exponentially declining SFR history with a long characteristic timescales τ_ disc≈ several Gyrs. On the other hand, we have noted that assigning such a disc-like star-formation history to the UV+far-IR/sub-mm/radio SFR functions would considerably overproduce the number of massive discs; this is because obscuration is mild in starforming discs, so the SFR function from dust-corrected UV data must be effectively exploited as input of the continuity equation. The effects of mergers on the stellar mass function of discs are negligible.


  * We have found that the stellar mass function of quiescent galaxies from the continuity equation is in excellent agreement with current observational constraints for z 4. We thus have demonstrated quantitatively via the continuity equation that the dusty, strongly star-forming galaxies recently discovered thanks to wide area far-IR/sub-mm surveys at z 1 are indeed the progenitors of the massive quiescent galaxies increasingly detected out to high redshift z 4 via deep near-IR surveys. We have estimated that the typical time of quiescence (i.e., with absent or negligibly small SFR) for these galaxies is around Δ t_ qui≈ 250-500 Myr.
    To further test the outcomes of the continuity equation, and better constrain the input SFR functions and the parameters of the star-formation history for individual galaxies, it will be crucial to improve the accuracy in the determination of the stellar mass function at the high-mass end M_⋆ 10^11  M_⊙ out to z 6 via wide areas near-IR surveys.


  * We have determined the cosmic mass density, finding it in excellent agreement with observational determination out to z∼ 0-8, both for active and quiescent galaxies. The continuity equation implies an analogue of the Soltan argument for the stellar component, in such a way that the cosmic stellar mass density is by construction consistent with the cosmic time-integrated star formation history, besides a factor depending on the IMF.


  * We have determined the star-formation efficiency of galaxies as a function of the stellar mass in the local Universe, finding it in good agreement with diverse observations. We have found, in line with previous studies, that the efficiency of star-formation is lower than f_⋆≈ 20-30%, with the maximum value being attained around a characteristic stellar mass of a few 10^10  M_⊙. The behavior as a function of stellar mass can  be ascribed to different form of feedbacks regulating star formation in galaxies, with supernovae and stellar winds dominating for stellar masses below the characteristic one, and AGN feedback dominating above.
    We have also pointed out that the decline of the efficiency for large
    masses is somewhat spurious, being related to the fact that the most
    massive galaxies tend to live at the center of group/cluster halos, which contain a lot of DM; considering instead only the DM mass belonging to individual galactic halos would imply the efficiency to stay almost constant or increase somewhat at large stellar masses. Finally, we have stressed that the similarity of the efficiency at z 2 to the local one is indicative of the early, in-situ nature of the star formation process, at least for massive spheroidal galaxies; in fact, we have noted that the global efficiency at z≈ 2 can be brought on that observed locally for quiescent galaxies by letting the stellar mass to evolve of a modest factor 50% due to late star formation or dry mergers and the halo mass to evolve by a factor of a few due to late smooth accretion and/or tidal stripping.


  * We have computed the main sequence of star forming galaxies via abundance matching of the input SFR function and of the stellar mass function self-consistently derived from the continuity equation. We have found a remarkable agreement with the observational determinations at different redshifts, so further constraining our input star-formation histories and timescales. We have highlighted how off-main sequence galaxies (located above the average relation) can be simply interpreted in the light of our star formation histories as
    young objects, caught when their stellar mass is still to be accumulated; they will then progressively move onto the main sequence, where they will spend most of their lifetime as active galaxies, before being quenched.



Finally, we conclude by stressing that the added value of the continuity equation, developed here on the stellar component of galaxies, is to provide quantitative, yet largely model-independent outcomes which must be complied by detailed physical models. In particular, the continuity equation allows a full exploitation of the redshift-dependent SFR functions, stellar mass functions, and galaxy main sequence, in order to determine the average star-formation histories and timescales of individual galaxies. Our analysis highlights that a bimodal star-formation history is required for spheroids and discs: the former must be characterized by a nearly constant SFR over short timescales  Gyr (increasing somewhat for less star-forming objects), and the latter must feature a SFR exponentially declining over long timescales of several Gyrs. Such outcomes of the continuity equation can provide inspiring hints on ways to improve the (sub-grid) physical recipes implemented in theoretical models and numerical simulations; moreover, they can offer a benchmark for forecasts on future observations at very high redshift with multi-band coverage on medium and wide areas, as it will become routinely achievable with the advent of the JWST.


We are grateful to J. Beacom, K. Glazebrook, M. Massardi, F. Pozzi, and P. Salucci for stimulating discussions. We thank the referee for a constructive report. Work partially supported by PRIN MIUR 2015 `Cosmology and Fundamental Physics: illuminating the Dark Universe with Euclid' and PRIN INAF 2014 `Probing the AGN/galaxy co-evolution through ultra-deep and ultra-high-resolution radio surveys'. AL acknowledges the RADIOFOREGROUNDS grant (COMPET-05-2015, agreement number 687312) of the European Union Horizon 2020 research and innovation programme.




Alavi, A., Siana, B., Richard, J., et al. 2016, ApJ, 832, 56

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

Barger, A. J., Cowie, L. L., Owen, F. N., et al. 2015, ApJ, 801, 87

Barro, G., Kriek, M., Pérez-González, P. G., et al. 2016, ApJ, 827, L32

Behroozi, P., & Silk, J. 2017, MNRAS, in press [arXiv:1609.04402]

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

Bernardi, M., Meert, A., Sheth, R. K., et al. 2017, MNRAS, 467, 2217

Bernardi, M., Meert, A., Sheth, R. K., et al. 2013, MNRAS, 436, 697

Bocquet, S., Saro, A., Dolag, K., & Mohr, J.J. 2016, MNRAS, 456, 2361

Bourne, N., Dunlop, J. S., Merlin, E., et al. 2017, MNRAS, 467, 1360

Bouwens, R.J., Oesch, P.A., Illingworth, G.D., Ellis, R.S., &
Stefanon, M. 2017, ApJ, 843, 129

Bouwens, R. J., Aravena, M., De Carli, R., et al. 2016, ApJ, 833, 72

Bouwens, R. J., Illingworth, G. D., Oesch, P. A., et al. 2015, ApJ, 803, 34

Burkert, A., Forster Schreiber, N. M., Genzel, R., et al. 2016, ApJ, 826, 214

Cai, Z.-Y., Lapi, A., Bressan, A., et al. 2014, ApJ, 785, 65

Cai, Z.-Y., Lapi, A., Xia, J.-Q., et al. 2013, ApJ, 768, 21

Calzetti, D., Armus, L., Bohlin, R. C., et al. 2000, ApJ, 533, 682

Caputi, K. I., Ilbert, O., Laigle, C., et al. 2015, ApJ, 810, 73

Cassará, L. P., Maccagni, D., Garilli, B., et al. 2016, A&A, 593, A9

Cavaliere, A., Morrison, P., & Wood, K. 1971, ApJ, 170, 223

Chabrier G. 2003, ApJ, 586, L133

Chiappini, C., Matteucci, F., & Gratton, R. 1997, ApJ, 477, 765

Citro, A., Pozzetti, L., Moresco, M., & Cimatti, A. 2016, A&A, 592, A19

Comparat, J., Prada, F., Yepes, G., & Klypin, A. 2017, MNRAS, 469, 4157

Conroy, C. 2013, ARA&A, 51, 393

Contini, E., Kang, Xi, Romeo, A. D., & Xia, Q. 2017, ApJ, 837, 27

Cooray, A., Calanog, J., Wardlow, J. L., et al. 2014, ApJ, 790, 40

Coppin, K. E. K., Geach, J. E., Almaini, O., et al. 2015, MNRAS,
446, 1293

Courteau, S., Cappellari, M., de Jong, R. S., et al. 2014, RvMP, 86, 47

Cucciati, O., Tresse, L., Ilbert, O., et al. 2012, A&A, 539, A31

da Cunha, E., Walter, F., Smail, I. R., et al. 2015, ApJ, 806, 110

Daddi, E., Alexander, D. M., Dickinson, M., et al. 2007, ApJ, 670, 173

Davé, R., Rafieferantsoa, M.H., Thompson, R.J., & Hopkins, P.F. 2017, MNRAS, 467, 115

Davidzon, I., Ilbert, O., Laigle, C., et al. 2017, ApJ, in press [arXiv:1701.02734]

Diemer, B., Sparre, M., & Torrey, P. 2017, ApJ, 839, 26

Duncan, K., Conselice, C. J., Mortlock, A., et al. 2014, MNRAS, 444, 2960

Dunlop, J. S., McLure, R. J., Biggs, A. D., et al. 2017, MNRAS, 466, 861

Dye, S., Furlanetto, C., Swinbank, A. M., et al. 2015, MNRAS, 452, 2258

Efstathiou, A., Rowan-Robinson, M., & Siebenmorgen, R. 2000, MNRAS, 313, 734


Finkelstein, S. L., Ryan, R. E., Jr., Papovich, C., et al. 2015, ApJ, 810, 71

Fudamoto, Y., Oesch, P.A., Schinnerer, E., et al. 2017, MNRAS, in press [arXiv:1705.01559]

Furlong, M., Bower, R. G., Crain, R. A., et al. 2017, MNRAS, 465, 722

Gallazzi, A., Charlot, S., Brinchmann, J., & White, S. D. M. 2006, MNRAS, 370, 1106

Genzel, R., Forster Schreiber, N. M., Ubler, H., et al. 2017, Natur, 543, 397

Gladders, M.D., Oemler, A., Dressler, A., Poggianti, B., Vulcani, B., & Abramson, L. 2013, ApJ, 770, 64

Glazebrook, K., Schreiber, C., Labbé, I., et al. 2017, Natur, 544, 71

Grazian, A., Fontana, A., Santini, P., et al. 2015, A&A, 575, A96

Gruppioni, C., Calura, F., Pozzi, F., et al. 2015, MNRAS, 451, 3419

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

Harris, K., Farrah, D., Schulz, B., et al. 2016, MNRAS, 457, 4179

Harrison, C. M., Simpson, J. M., Stanley, F., et al. 2016, MNRAS, 457, L122

Hill, A.R., Muzzin, A., Franx, M., et al. 2017, ApJ, 837, 147

Hopkins, P.F., Wetzel, A., Keres, D., et al. 2017, MNRAS, in press [arXiv:1702.06148]

Hudson, M. J., Gillis, B. R., Coupon, J., et al. 2015, MNRAS, 447, 298

Ikarashi, S., Ivison, R. J., Caputi, K.I., et al. 2017, ApJ, 835, 286

Ikarashi, S., Ivison, R. J., Caputi, K.I. 2015, ApJ, 810, 133

Ilbert, O., McCracken, H. J., le Fevre, O., et al. 2013, A&A, 556, A55

Kennicutt, R. C. 1983, ApJ, 272, 54

Koprowski, M., Dunlop, J. S., Michalowski, M. J., et al. 2016, MNRAS, 458, 4321

Koprowski, M. P., Dunlop, J. S., Michalowski, M. J., Cirasuolo, M., & Bowler, R. A. A. 2014, MNRAS, 444, 117

Kriek, M., Conroy, C., van Dokkum, P.G., et al. 2016, Natur, 540, 248

Kurczynski, P., Gawiser, E., Acquaviva, V., et al. 2016, ApJL, 820, L1

Lacey, C. G., Baugh, C. M., Frenk, C. S., et al. 2010, MNRAS, 405, 2

Lapi, A., Mancuso, C., Celotti, A., & Danese, L. 2017, ApJ, 835, 37

Lapi, A., Salucci, P., & Danese, L. 2013, ApJ, 772, 85

Lapi, A., Gonzalez-Nuevo, J., Fan, L., et al. 2011, ApJ, 742, 24

Lee, S.-K., Idzi, R., Ferguson, H. C., et al. 2009, ApJS, 184, 100

Leitner, S.N., & Kravtsov, A.V. 2011, ApJ, 734, 48

Leja, J., van Dokkum, P. G., Franx, M., & Whitaker, K. E. 2015, ApJ, 798, 115

Lilly, S. J., Carollo, C. M., Pipino, A., et al. 2013, ApJ, 772, 119

Lonoce, I., Longhetti, M., Maraston, C., et al. 2015, MNRAS, 454, 3912

Ma, J., Gonzalez, A. H., Viera, J. D., et al. 2016, ApJ, 832, 114

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

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

Man, A.W.S., Greve, T.R., Toft, S., et al. 2016, ApJ, 820, 11

Mancuso, C., Lapi, A., Shi, J., et al. 2016a, ApJ, 823, 128

Mancuso, C., Lapi, A., Shi, J., et al. 2016b, ApJ, 833, 152

Mandelbaum, R., Wang, W., Zu, Y., et al. 2016, MNRAS, 457, 3200

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

Mawatari, K., Yamada, T., Fazio, G. G., Huang, J.-S., & Ashby, M. L. N. 2016, PASJ, 68, 46

Mendel, J. T., Simard, L., Palmer, M., Ellison, S.L., & Patton, D.R. 2014, ApJS, 210, 3

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

Meurer, G. R., Heckman, T. M., & Calzetti, D. 1999, ApJ, 521, 64

Michalowski, M.J., Dunlop, J.S., Koprowski, M.P., et al. 2016, MNRAS, 469, 492

Moffett, A.J., Lange, R., Driver, S.P., et al. 2016, MNRAS, 462, 4336

Mor, R., Netzer, H., Trakhtenbrot, B., Shemmer, O., & Lira, P. 2012, ApJL, 749, L25

More, S., van den Bosch, F. C., Cacciato, M., et al. 2011, MNRAS, 410, 210

Moster, B. P., Naab, T., & White, S. D. M. 2013, MNRAS, 428, 3121

Moustakas, J., Coil, A. L., Aird, J., et al. 2013, ApJ, 767, 50

Mullaney, J. R., Daddi, E., Bethermin, M., et al. 2012, ApJL, 753, L30

Nayyeri, H., Keele, M., Cooray, A., et al. 2016, ApJ, 823, 17

Negrello, M., Amber, S., Amvrosiadis, A., et al. 2017, MNRAS, 465, 3558

Negrello, M., Hopwood, R., Dye, S., et al. 2014, MNRAS, 440, 1999

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

Novak, M., Smolcic, V., Delhaize, J., et al. 2017, A&A, 602, A5

Oesch, P. A., Bouwens, R. J., Carollo, C. M., et al. 2010, ApJL, 725, L150

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

Papovich, C., Finkelstein, S. L., Ferguson, H. C., Lotz, J. M., & Giavalisco, M. 2011, MNRAS, 412, 1123

Pezzulli, G., & Fraternali, F. 2016, MNRAS, 455, 2308

Planck Collaboration XIII 2016, A&A, 594, A13

Pope, A., Montana, A., Battisti, A., et al. 2017, ApJ, 838, 137

Reddy, N. A., Kriek, M., Shapley, A. E., et al. 2015, ApJ, 806, 259

Renzini, A. 2006, ARA&A, 44, 141

Renzini, A., & Peng, Y.-J. 2015, ApJL, 801, L29

Riechers, D.A., Daisy Leung, T.K., Ivison, R., et al. 2017, ApJ, submitted [preprint arXiv:1705.09660]

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

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

Rodriguez-Gomez, V., Pillepich, A., Sales, L.V., et al. 2016, MNRAS, 458, 2371

Rodriguez-Gomez, V., Genel, S., Vogelsberger, M., et al. 2015, MNRAS, 449, 49

Rodriguez-Puebla, A., Avila-Reese, V., Yang, X., et al. 2015, ApJ, 799, 130

Romano, D., Silva, L., Matteucci, F., & Danese, L. 2002, MNRAS, 334, 444

Rong, Y., Jing, Y., Gao, L., et al. 2017, MNRAS, 471, L36

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

Rowan-Robinson, M., Oliver, S., Wang, L., et al. 2016, MNRAS, 461, 1100

Salmon, B., Papovich, C., Finkelstein, S. L., et al. 2015, ApJ, 799, 183

Salpeter, E. E. 1955, ApJ, 121, 161

Salucci, P., Szuszkiewicz, E., Monaco, P., & Danese, L. 1999, MNRAS, 307, 637

Schaye, J., Crain, R. A., Bower, R. G., et al. 2015, MNRAS, 446,
521

Schreiber, C., Pannella, M., Leiton, R., et al. 2017, A&A, 599, A134

Scoville, N., Sheth, K., Aussel, H., et al. 2016, ApJ, 820, 83

Scoville, N., Aussel, H., Sheth, K., et al. 2014, ApJ, 783, 84

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

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

Shankar, F., Lapi, A., Salucci, P., de Zotti, G., & Danese, L. 2006, ApJ, 643, 14

Shi, J., Lapi, A., Mancuso, C., Wang, H., & Danese, L. 2017, ApJ, submitted

Silva, L., Granato, G. L., Bressan, A., & Danese, L. 1998, ApJ,
509, 103

Simpson, J. M., Smail, I., Wang, W.-H., et al. 2017, ApJ, 844, L10

Simpson, J. M., Smail, I., Swinbank, A. M., et al. 2015, ApJ, 799, 81

Small, T. A., & Blandford, R. D. 1992, MNRAS, 259, 725

Smit, R., Bouwens, R. J., Franx, M., et al. 2012, ApJ, 756, 14

Soltan, A. 1982, MNRAS, 200, 115

Song, M., Finkelstein, S. L., Ashby, M. L. N., et al. 2016, ApJ, 825, 5

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

Spilker, J. S., Bezanson, R., Marrone, D. P., et al. 2016, ApJ, 832, 19

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

Steinhardt, C.L., Yurk, D., & Capak, P., 2017, MNRAS, 468, 849

Steinhardt, C. L., Speagle, J. S., & Capak, P. 2014, ApJL, 791, L25

Straatman, C. M. S., Spitler, L. R., Quadri, R. F., et al. 2016, ApJ, 830, 51

Straatman, C. M. S., Labbé, I., Spitler, L. R., et al. 2015, ApJ, 808, L29

Straatman, C.M.S., Labbé, I., Spitler, L.R., et al.
2014, ApJ, 783, L14

Strandet, M. L., Weiss, A., Vieira, J. D., et al. 2016, ApJ, 822, 80

Swinbank, M., Harrison, C., Trayford, J., et al. 2017, ApJ, 467, 3140

Tadaki, K.-I., Genzel, R., Kodama, T., et al. 2017, ApJ, 834, 135

Tasca, L. A. M., Le Févre, O., Hathi, N. P., et al. 2015, A&A, 581, A54

Thanjavur, K., Simard, L., Bluck, A. F. L., & Mendel, T. 2016, MNRAS, 459, 44

Thomas, D., Maraston, C., Bender, R., & Mendes de Oliveira, C. 2005, ApJ, 621, 673

Tinker, J. L., Kravtsov, A. V., Klypin, A., et al. 2008, ApJ, 688, 709

Tomczak, A. R., Quadri, R. F., Tran, K. H., et al. 2016, ApJ, 817, 118

Tomczak, A. R., Quadri, R. F., Tran, K.-V. H., et al. 2014, ApJ, 783, 85

Vale, A., & Ostriker, J. P. 2004, MNRAS, 353, 189

van der Burg, R. F. J., Hildebrandt, H., & Erben, T. 2010, A&A, 523, A74

Velander, M., van Uitert, E., Hoekstra, H., et al. 2014, MNRAS, 437, 2111

Wang, R., Wagg, J., & Carilli, C. L. 2013, ApJ, 773, 44

Watson, W.A., Iliev, I.T., D'Aloisio, A., et al. 2013, MNRAS, 433, 1230

Weiss, A., De Breuck, C., Marrone, D. P., et al. 2013, ApJ, 767, 88

Wellons, S., Torrey, P., Ma, C.-P., et al. 2015, MNRAS, 449, 361

Whitaker, K. E., Franx, M., Leja, J., et al. 2014, ApJ, 795, 104


Willott, C. J., Bergeron, J., & Omont, A. 2015, ApJ, 801, 123

Wojtak, R., & Mamon, G. A. 2013, MNRAS, 428, 2407

Wyder, T. K., Treyer, M. A., Milliard, B., et al. 2005, ApJL, 619, L15

Xu, L., Rieke, G. H., Egami, E., et al. 2015, ApJ, 808, 159

Yu, Q., & Lu, Y. 2004, ApJ, 602, 603

Yu, Q., & Lu, Y. 2008, ApJ, 689, 732

Zavala, J.A., Montana, A., Hughes, D.H., et al. 2017, Natur, submitted



