The frequency of very young galaxies in the local Universe: I. A test for galaxy formation and cosmological models
Dylan P. Tweed Gary A. Mamon Trinh X. Thuan Andrea Cattaneo Avishai Dekel Nicola Menci Francesco Calura Jospeh Silk



firstpage–lastpage 2017

Identifying Young Kepler Planet Host Stars from Keck-HIRES Spectra of Lithium1
    Travis A. Berger2, Andrew W. Howard2,3, & Ann Merchant Boesgaard2
    
==============================================================================





In the local Universe, the existence of very young galaxies (VYGs), having
formed at least half their stellar mass in the last 1 Gyr, is debated. We
predict the present-day fraction of VYGs among central galaxies as a function
of galaxy stellar mass. For this, we apply to high mass resolution
Monte-Carlo halo merger trees (MCHMTs) three (one) analytical models of galaxy formation,
where the ratio of stellar to halo mass (mass growth rate) is a function of halo mass and redshift. 
Galaxy merging is delayed until orbital decay by
dynamical friction. With starbursts associated with halo mergers, our models
predict typically one percent of VYGs up to galaxy masses of 10^10,
falling rapidly at higher masses, and VYGs are usually associated with
recent major mergers of their haloes. Without these
starbursts, two of the models have VYG fractions reduced by 1 or 2 dex at low
or intermediate stellar masses, and VYGs are rarely associated with major
halo mergers. In comparison, the state-of-the-art semi-analytical model (SAM)
of Henriques et al. produces only 0.01 per cent of VYGs at intermediate
masses. Finally, the Menci et al.  SAM run on MCMHTs with Warm Dark Matter
cosmology generates 10 times more VYGs at masses below 10^8 than when run with Cold
Dark Matter. The wide range in these VYG fractions illustrates the usefulness
of VYGs to constrain both galaxy formation and cosmological models.





galaxies: formation – galaxies: evolution – galaxies: dwarf – galaxies:
statistics – methods: numerical 




§ INTRODUCTION


In the standard Λ Cold Dark Matter (ΛCDM) paradigm,
galaxies form by dissipative collapse inside dark matter haloes: the smaller
haloes are the first to detach from the Hubble expansion and collapse. The larger
most massive structures, clusters, groups and massive ellipticals, form later
through mergers.

The mass function and growth of dark matter haloes is
 well understood thanks to  the  
<cit.> theory and its extensions  
<cit.>, confirmed by large-scale cosmological N-body
simulations <cit.>.

However,
it has been a considerable challenge to  understand how galaxies form stars within
these haloes, because of the numerous physical processes involved. We know
that stars are formed in cold Giant Molecular Clouds of gas.
One first needs
to allow the gas to enter the haloes, but this process becomes inefficient in
haloes at the extremes of the mass function: 1) its entropy is too high to fall
into very
low-mass haloes <cit.>, 
2) its entropy is significantly raised when it is shock-heated near the
virial radius around high-mass <cit.>.
Moreover, in dense environments, the outer gas can be stripped before it can fall
onto the disc and fuel the molecular clouds, from
a) the tides from the group/cluster
potential 
<cit.>, and
b) the ram pressure it feels from its motion relative to the hot intra-group/cluster gas
<cit.>.
One then needs to retain the gas in the disc, against the feedback from
1) supernovae <cit.> 
and 2) active galactic nuclei <cit.>.

The early realization that elliptical galaxies of increasing luminosity have
redder colours <cit.> is now understood as partly due to the fact
that more massive ellipticals have older stellar populations
<cit.>, even if the colours of more massive ellipticals are also
redder because of their higher  metallicity  <cit.>.
This downsizing trend of older stellar populations for massive
galaxies can be explained by a decrease in the efficiency of star formation
above some halo mass <cit.>.

On the opposite end, the youngest galaxies should be those with the lowest
metallicities, as the neutral gas from which present stars are formed has not
been polluted by many previous generations of stars. Low-metallicity
star-forming objects possess strong emission-line spectra, characteristic of
HII regions and indicating the presence of an intense burst of star formation
(e.g. ).
Emission-line galaxies tend to be of low stellar mass (, who
used the near infrared J band as a proxy for stellar mass). Similarly, the
strong positive correlations of metallicity with both the luminosity of ellipticals 
<cit.> and the stellar mass of irregular and blue compact dwarfs (BCDs,
 and with galaxies in general,
) suggest that the youngest galaxies must be of low stellar
mass.

In particular, a prime candidate for a galaxy with a very young stellar
population is I Zw 18, which has a very
low stellar mass around 10^7  M_⊙ <cit.> and
an extremely low  metallicity (1/50th of solar).
Its spectrum shows strong emission-lines, indicative of active star
formation producing thousands of O stars emitting plenty of ionizing
radiation. Using the Hubble Space Telescope (HST) to resolve its stellar content
and construct its Hertzprung-Russell diagram, <cit.> found
that the bulk of the stellar population of I Zw 18 is younger than 500
Myr. Later, <cit.> and 
<cit.> used  deeper HST 
imaging data of I Zw 18 and detected an older
stellar population with age greater than  1 Gyr. 
However, the mass of the old stellar population is not known, as it depends on
the unknown star formation history of the galaxy. It is thus not clear
whether most of the stellar mass of I Zw 18 was formed within the last
Gyr or earlier.
Another very young galaxy candidate  is J0811+4730, 
recently discovered by  <cit.>.
These authors found it to be even more metal-poor than I Zw 18, and estimate
that
three-quarters of its
stellar mass is younger than (only) 5 Myr.

Motivated by the lack of galaxies whose bulk of stellar mass was
undoubtedly 
formed in
the last Gyr, and by the debate on the epoch when I Zw 18 formed half of its
stellar mass, we are led to the following questions.
Can very young galaxies
(hereafter VYGs) exist?
If yes, how frequent are VYGs in the
local Universe? 
In this article, we define VYGs as
galaxies in which more than half of the stellar mass  was formed in the last Gyr.
This critical age of 1 Gyr is motivated by the debate on the age of I Zw 18,
but is otherwise arbitrary, and we will also investigate how the frequency of
VYGs depends on this choice of critical age. 
Since I Zw 18 <cit.> and J0811+4730 <cit.> are very
isolated galaxies, and since the physics of satellite galaxies is more
debated than that of centrals, we choose to focus on central galaxies, thus excluding
satellite galaxies.

However, it is not our aim, in this article, to model in
detail the properties of candidate VYGs I Zw 18 and J0811+4730, but to
generally explore the fractions of VYGs among central galaxies
as a function of their z=0 stellar mass. 

We estimate the fraction of VYGs in bins of present stellar mass by using 
current models of galaxy
formation, both analytical and semi-analytical, and we also
compare the predictions between
a Warm Dark Matter cosmology and the standard ΛCDM.
In a companion article (Trevisan et al., in prep., hereafter Paper II), we
estimate the fractions of VYGs as a function of stellar mass 
in the local Universe, using the Sloan Digital
Sky Survey (SDSS) spectral database, and compare them with the model
predictions presented here.

The models that we constructed to predict the frequencies of VYGs are  described
in
Sect. <ref> and tested in Sect. <ref>.
In Sect. <ref>, we compute the fractions of VYGs, for the
different models.
These results are discussed in Sect. <ref> and summarized in Sect. <ref>.



§ METHODS





 §.§ Basic considerations


Our choice of methods is guided by our requirement of producing large samples of
galaxies with sufficient mass resolution to form galaxies in a range of
stellar masses extending down to include the best two cases for VYGs, I Zw 18
(=7.38, )
and 
J0811+4730 (=6.3, ).
We must then resolve the much lower mass
progenitors of such z=0 galaxies.
One method is to use semi-analytical models of galaxy formation and evolution
(hereafter, SAMs), run on the halo merger trees derived from the dark matter
haloes extracted from high-resolution dissipationless cosmological N-body
simulations.
We analyse here the z=0 output of the recent state-of-the-art SAM of
<cit.>, run on 
two dissipationless cosmological N-body simulations: 
the Millennium Simulation (MS, ) and the Millennium-II
Simulation (MS-II,
).
Both simulations are re-scaled in time and space to the Planck 2014 cosmology, 
with (Ω_ m,Ω_Λ,h,σ_8)=(0.315,0.685,0.673,0.829)
using the technique of <cit.>, updated by
<cit.>.

However, with  a particle mass of 9.5× 10^6  M_⊙,
the mass resolution of the  MS-II
is barely sufficient to resolve the haloes around our lower mass galaxies.
For example, galaxy I Zw 18, whose halo log mass is 
= 8.9 <cit.> or = 8.5 <cit.>, both based upon
the distance of ),[We denote M the halo mass and
  m the galaxy stellar mass.] would only be resolved at z=0 with 90
particles
(the MS simulation, whose resolution is 125 coarser, clearly misses these
haloes).
Perhaps the halo mass of I Zw 18 is underestimated as figure 3 of
<cit.> suggests a minimum value of
10^10.
While no halo mass is available for galaxy J0811+4730, we infer from figure 3
of  that its halo mass may be as low as =9.2,
i.e. this galaxy's halo would be resolved with 170 particles.
However, the
progenitors of I Zw 18 and J0811+4730 would be only marginally resolved in the MS-II.
This led us to also
consider Monte-Carlo halo merger trees rather than only rely on cosmological N-body
simulations to achieve adequate mass resolution for the haloes. So, in
addition to considering the SAM of , we also run
simple, single-equation, galaxy formation models on high mass-resolution Monte-Carlo halo
merger trees to derive the growth of the stellar masses of galaxies.

We discuss these Monte-Carlo halo merger trees in Sect. <ref>
and present the
analytical modelling in Sect. <ref> and the semi-analytical models in Sect. <ref>.



 §.§ Halo merger trees



 

Monte-Carlo halo merger trees are designed to generate realistic
merger histories of a given halo of mass M at a redshift z (usually 0).
These merger trees are built by
generating progenitor masses at a higher redshift, and iterating over those
progenitors.
For each
halo of mass M_0 at redshift z_0, 
the mass M_1 of the main progenitor at redshift z_1 is drawn according to a
probability distribution function that can be written as
P_1(M_1|M_0,z_0,z_1). A secondary progenitor mass can then be drawn
following a probability distribution P_2(M_2|M_1,M_0,z_0,z_1). For
some codes, multiple secondary progenitors can be drawn following
P_2(M_i|∑_j=1^i-1 M_j, M_0, z_0, z_1). 
Monte-Carlo halo merger tree codes handle mass conservation in different ways,
either neglecting smooth accretion by imposing ∑_j
M_j = M_0, or incorporating diffuse mass growth, i.e. ∑_j M_j ≤
M_0. The difference Δ M = M_0 -∑_j M_j can be interpreted as smooth accretion or unresolved
mergers. This process is iterated to increasingly higher redshifts, thus building the
branches of the halo merger tree down to a predefined mass resolution
or up to a maximal redshift. Since those probability distributions do
not depend on the previous (lower redshift) outcome, the entire process is
Markovian.

The first implementations of Monte-Carlo halo merger trees
<cit.> used the extension of the
<cit.> model for the cosmic halo mass function, while modern
implementations use more accurate probability distribution functions.
<cit.> have recently compared 6 implementations of 4 halo merger
tree codes, for z=0 halo masses M_0 ranging from 10^11 to 10^15 h^-1 
M_⊙, with branches of mass >10^-4 M_0.
They concluded that the code of <cit.>,
which is based on extended Press
Schechter theory <cit.>, 
with an additional term that is designed to achieve better mass
conservation,
reproduced best the mass assembly histories, merger rates and unevolved
subhalo mass functions
measured
in cosmological N-body simulations with 
the same cosmological parameters as previously used in the Millennium
simulations: (Ω_m,Ω_Λ,h,σ_8) =
(0.25,0.75,0.73,0.9),
which is close to the WMAP 1^ st year
cosmology <cit.>.
 calibrated the two free parameters of
their algorithm to match the conditional mass functions, as well as the distribution
of the epochs of most recent major mergers, both measured in the Millennium
simulation. 
We therefore adopted the code of <cit.>.

The  code creates a binary tree
with very fine time resolution,
so 
that only the main and one secondary halo are drawn for a given halo
mass.
This code enables the user to 
adopt a custom, coarser, output time resolution of the merger tree.
Thus, our tree outputs may contain non-binary mergers, but these are built
from binary mergers at the fine internal resolution of the code. We chose 101
output timesteps in equal increments of 
log(1+z) from redshift z=0 to redshift z=13.
Our first non-zero redshift is z=0.0267, corresponding to a lookback time
of 350 Myr.

Figure <ref> illustrates the Monte-Carlo halo merger tree obtained
with the 
halo merger tree code. 
While mergers of branches are clearly seen along the main branch, they also 
 occur within the
secondary branches. The figure shows the diversity of halo formation
times (summed over the progenitors): before being merged into more
massive haloes, massive haloes always form early (brown) while low-mass haloes form
early (brown) or late (green).  

A 
 halo
can have only one
descendant per timestep (halo fragmentation is not allowed). 
However, the number of progenitors
a
given halo can have is zero (a newborn
halo), one (a quiescent halo), or several
(mergers occurring between time outputs).

We ran the  halo merger tree code 
with the cosmological parameters of the Millennium simulations, i.e.,
(Ω_m,Ω_Λ,h,σ_8) =
(0.25,0.75,0.73,0.9). We adopted final halo
log masses of
log (M_0/ M_⊙) = 7 to 14 in steps of 0.025 dex, and a mass
resolution M_ min = 10^-4 M_0.
The minimum halo mass of 10^7 is chosen to ensure that we include
the halo mass of I Zw 18 (possibly as low as 10^8.5
).
By going to such low halo masses, we are assuming that
the  tree code remains valid 4 dex below where it was tested. 
For each value of log M_0, we have run the halo merger tree code 1000 times with different
random seeds. In total, we have generated 1000× [(14-7)/0.025+1] = 281 000
halo merger trees.

These merger trees are qualitatively very similar
to analogous merger trees extracted from (cosmological) N-body simulations. 
Our Monte-Carlo trees have superior mass resolution in comparison with the
halo trees extracted from the MS or even from the MS-II: the final halo
masses extend 2 dex lower than the haloes resolved by MS-II with 100
particles, and the progenitors of our haloes have 4 extra dex of resolution,
so that our Monte Carlo halo merger trees reach a progenitor mass resolution
that is 6 orders of magnitude better than that of the MS-II haloes. 

Note that, while the haloes in N-body simulations can decrease in mass from one step to the
next (because of tidal forces during close interactions),
our haloes cannot lose mass, by construction.
 


 §.§ Analytical galaxy formation models





  §.§.§ Basic formalism


We first 
consider very simple galaxy formation models 
where stellar masses m_ stars are assigned to
haloes with a star formation efficiency (SFE) that depends on halo
mass M and redshift z, according to

    m_ stars≡m(M,z) = f_ b M F_ SFE(M,z)  ,

where
the tilde sign
is to denote that this is a model.
In equation (<ref>), f_ b = Ω_ b/Ω_
  m
is
the universal baryonic
fraction (f_ b = 0.18 for our adopted Ω_ b = 0.045 as in the Millennium
simulations), while
F_ SFE represents the ratio between the stellar mass of a halo and the
total mass in baryons expected within the halo.
Our analytical models are entirely based on our choice for F_
  SFE(M,z).
They predict the stellar mass, but not the gas mass.

We consider a physically-motivated model, as well as two empirical
ones based on abundance matching (see Sect. <ref>).
We  also  use another empirical model, based on an equation
similar to equation (<ref>), but where masses are replaced by
mass variations.
This use
of several different galaxy formation 
models allows us to gauge the dependence of our results on the uncertainties
of galaxy formation.

If one specifies a form for F_ SFE(M,z), one can derive the stellar
mass history of every halo. Following <cit.>, who pioneered the
use of equation (<ref>),  and <cit.>, 
we forbid stellar masses to decrease, i.e. the stellar mass of the 
(central) galaxy has to be greater or equal to the sum of its galaxy
progenitors. 



  §.§.§ Physical analytical model: <cit.> with <cit.>  at the low-mass end


<cit.> (hereafter C11) have presented a quasi-physical model for m(M,z)
that combines 1) supernova feedback, 2) a gentle cutoff at the 
high-mass end caused by the virial shock around high-mass haloes <cit.>,
and 3) a sharp cutoff at the low-mass end  (motivated by
early hydrodynamical simulations of ) due to the
entropy barrier that prevents high entropy gas from collapsing onto
its halo (given that its entropy cannot decrease).
This is written as F_ SFE∝ 1-v_ reion^2/v_ c^2, where
v_ reion is a constant, while
v_ c^2 is the squared circular velocity of the halo

    v_ c^2 =  [Δ(z) 2]^1/3 [ H(z) G M ]^2/3 ,

where Δ(z) is the ratio of the mean density within the virial radius to
the critical density of the Universe, while H(z) is the Hubble constant.

The supernova feedback prescription of  is physically
motivated: it is 
based on the idea that a fraction of the accreted gas that is processed into
stars is rapidly ejected as supernova winds, whose velocity matches the
virial velocity of the halo and whose energy is assumed to be entirely
mechanical and proportional to the remaining
stellar mass. This yields F_ SFE∝ 1/(1+v_ SN^2/v_
  c^2),
where v_ SN is another constant.
At large halo masses,
the virial shock quenches the infall of cold gas filaments by heating them up
near the virial radius. 
Star formation in the disc is then limited by the
much longer cooling time of the shock-heated gas.
This is assumed to yield 
 F_ SFE∝ 1/(1+M/M_ shock),
where M_ shock is a third constant.

This leads to equation (8) of <cit.>:

    F_ SFE^ C11(M,z) = 1-v_ reion^2/v_ c^2  1 +
      M/M_ shock  (1+v_
        SN^2/v_ c^2)^-1 ,

where the constants
 v_ reion,
v_ SN, 
and 
M_ shock  respectively 
represent the minimum circular velocity for  gas to overcome the entropy
barrier and collapse with the dark matter and subsequently form stars,
the impact of supernova feedback, 
and
the characteristic minimum mass
for the occurrence of virial shocks.
The first term of equation (<ref>)  describes the ability to accrete
gas, while the second term describes the ability to retain this accreted
gas. 

At
high halo masses, galaxy mergers produce galaxies with higher
masses than predicted with  F_ SFE(M,z) (see ),
but the present work focuses on intermediate- and low-mass galaxies.

Given the expected sensitivity of the fraction of VYGs to the form
of the threshold of star formation efficiency, F_ SFE,
we improve on the model of <cit.> by introducing, at the low-mass end,  a 
 smoother cutoff, derived from hydrodynamical 
simulations <cit.>. The accreted mass (disregarding for now
the virial shock affecting higher masses) is no longer m_ accr = f_ b M (1- v_
  reion^2/v_ c^2) but instead

    m_ accr = [1+(2^2/3-1)(v_
            reion/v_ c)^6]^-3/2
      f_ b M  .

Thus, v_ reion is no longer the halo circular velocity below which no
star formation can occur, but instead the halo circular velocity where the
accreted mass is reduced by a factor 2 by the entropy barrier. There
is no longer a sharp cutoff of star formation at v_ c = v_ reion
as in (eq. [<ref>]), but instead the stellar mass
increases as v_ c^9 (much steeper than the Tully-Fisher relation). 

Using equation (<ref>), equation (<ref>) becomes

    F_ SFE^  C+G(M,z) = [
    1+(2^2/3-1)  ( v_ reion/v_ c)^6
    ]^-3/2 (1+M/M_ shock) (1+v_ SN^2/v_ c^2) .


While the terms for the effects of supernovae and the entropy barrier are
respectively motivated by physical principles and hydrodynamical simulations, 
the term involving the virial shocks is, admittedly,  empirical
(which is why we dub this model `quasi-physical', but later call it
`physical' to distinguish it from the fully empirical models that we will discuss below).

In this `Cattaneo+Gnedin' model (hereafter, C+G), we refer to the values of
v_ reion before and after 
reionization as v_ pre-reion and v_ post-reion, respectively.
We adopt the following parameters, which produce a good fit to the
present-day mass function of galaxies (see Figure <ref> below):
v_ SN = 300   km   s^-1,
  v_ pre-reion = 18   km   s^-1,
  v_ post-reion = 50   km   s^-1,
  M_ shock = 10^12 M_⊙.[There was no
  pre-reionization velocity in <cit.>, who had assumed that
  Universe was fully reionized from the start.]
According to equation (<ref>), gas accretion is suppressed by a factor
10 if v_ c = 0.74  v_ reion. Since v_ c = σ_v/0.7
for NFW halos (derived from eqs. [22] and [24] of ),
v_ pre-reion = 18   km   s^-1 corresponds to a halo
temperature of 
10 500  K, i.e. the temperature of atomic Hydrogen cooling. Similarly, 
v_ post-reion = 50   km   s^-1 corresponds to T=82 000 
K. 
Our adopted value of M_ shock is double that of
<cit.> and 10 per cent lower than the value employed by
<cit.>.
Finally, our adopted value of v_ SN is chosen to roughly fit the z=0
stellar mass function.

The effects of these parameters on the stellar mass function is shown in
preliminary versions of this work <cit.>. 
The baryonic  Tully-Fisher relation is well reproduced by the  model  <cit.>.

We could have assumed that the Universe has reionized instantaneously at a fixed
redshift z_ reion, somewhere between 6 and 12.
Although reionization fronts are thought to have spread fast throughout
the Universe <cit.>, there is observational evidence that reionization
took a time comparable to the age of the Universe at that epoch. This is
suggested  
by the redshift difference between the epoch z≈ 9 where the optical depth to
neutral Hydrogen was unity
<cit.>
 and the latest
epoch, z≈ 6,
when evidence of a substantially neutral intergalactic medium is seen
<cit.>.
For this reason, we  have assumed instead that the Universe reionizes
stochastically.
Thus,
for each halo merger tree, we have drawn log (1+z_ reion) from a Gaussian
of mean corresponding to a median reionization redshift of
⟨ z_ reion⟩ = 7.5 
and a standard deviation σ[log(1+z_ reion)] = 0.05.
This stochastic reionization redshift is applied to all progenitors
within a merger tree. We suppose that merging haloes lie in the same region of
the Universe and thus are
reionized at the same redshift.
The mean SFE (only used for illustrative figures, see Fig. <ref>
below) is then

    ⟨ F_ SFE^ C+G(M,z)⟩   =   
    F_ reion(z)  F_ SFE-pre^ C+G(M,z) 
            + 
     [1-F_ reion(z)] F_ SFE-post^ C+G (M,z) 
     ,

where
F_ SFE-pre^ C+G and F_ SFE-post^ C+G are the SFEs for
v_ reion = v_ pre-reion and v_ reion = v_
  post-reion, respectively,
while

    F_ reion(z) = 12 {1+ erf[log(1+z)-log(1+⟨ z_
          reion⟩)√(2) σ[log(1+z_
          reion)]]}

is the probability that reionization occurs before redshift z.
Equation (<ref>) implies that 
 90 per cent of the
reionization occurs between z=9.3 and 6.0.  




  §.§.§ Empirical model: <cit.>


The relation between stellar and halo mass can also be determined empirically, by
comparing the cumulative distribution function (CDF) of the observed stellar mass
function with the CDF of the halo mass function predicted from cosmological
N-body simulations, a method known as
abundance matching <cit.>.

 (, hereafter MNW) have proposed the following form for the star formation efficiency:
 
    F_ SFE^ MNW(M,z) = 2 N f_ b [
     (M M_1)^-β+
     (M M_1)^γ]^-1 ,

where F_0, M_1, β and γ are positive constants.
 
computed the halo and subhalo mass functions in  the
MS and MS-II
rescaled to WMAP-7
cosmological parameters (Ω_ m=0.272, h=0.704, σ_8=0.81).
Since subhaloes are tidally stripped,  computed the subhalo masses at the
time they entered their halo, and accounted for orphan subhaloes that
were no longer resolved in the simulations but should have survived given
their long expected orbital decay time from dynamical friction.
Performing abundance matching at several 
redshifts from 0 to 4,  determined

    N    =    0.0351 - 0.0247  z/(1+z)  , 
    
    M_1    =    dex [11.59 + 1.195 z/(1+z) ]  ,
    β   =    1.376 - 0.826  z/(1+z)  ,
    γ   =    0.608 + 0.329 z/(1+z) .




  §.§.§ Empirical model: <cit.>


 (, hereafter BWC)  
have also (independently) performed abundance matching, with however several
differences as compared to
<cit.>.
First, their halo mass functions are derived from a different
dissipationless cosmological N-body simulation, Bolshoi <cit.>, with
cosmological parameters (Ω_ m=0.27, h=0.70, σ_8=0.82)
consistent with WMAP5 <cit.>.
The second difference is the more complex form for F_ SFE(M,z):

    F_ SFE^ BWC(M,z) = 
      (ϵ M_1
    f_ b M)   dex{ 
    f [log (M M_1)]-f(0)
    } ,

where

    f(x) = -log (1+10^α x) + δ {log[1+ e^x]}^γ
      1+exp(10^-x) .


 argue that their more refined shape for F_ SFE(M,z)
  improves its accuracy by a factor 4. 

The third difference with  is that

fitted their model and the halo+subhalo mass function to, not only the
  stellar mass function at different redshifts, but also to the 
  dependence with mass of the specific star formation rate at different redshifts, as
  well as to the evolution of the cosmic star formation rate  with
  redshift. 
The fourth difference is that the  model is fit to 
observational data at sufficiently high redshifts to start probing
the epoch of reionization:  from 0 to 8 (while  stops at z=4).
  derive the halo+subhalo mass function by assuming a model for the
dependence with halo/subhalo mass of the 
fraction of haloes or subhaloes that are subhaloes (which they calibrate at high
masses to measurements from simulations).

Their analysis yields the best-fit set of parameters:

    ϵ   =    dex [-1.777 - 0.006 (a-1) ν -0.119 (a-1)]   , 
    
    M_1    =    dex{11.514 -[1.793 (a-1)+0.251 z ] ν} , 
    α   =    -1.412 + [0.731 (a-1)] ν , 
    δ   =    3.508 + [2.608 (a-1)-0.043 z] ν , 
    γ   =    0.316 + [1.319   (a-1)+0.279 z] ν ,

where

    a    =   1 1+z ,
    ν   =   exp [-4 (1+z)^2]  .




  §.§.§ Empirical model: <cit.>


The model of  (, hereafter MCP)
replaces the stellar and halo masses of 
equation (<ref>) by their time derivatives (see also ):

    dm dt = f_b  dM dt F_ SFE^ MCP(M,z)  .


 assume lognormal functions for F_ SFE(M) at
given z, and we adopt their evolving, “M_ vir”-based model, which in our
notation is

    F_ SFE^ MCP(M,z) =  E(z) 
    exp [
      - (log M-log M_ peak(z)σ(z))^2
      ] ,

where

    E(z)    =    0.90 (1+z)^-0.74 ,
    log M_ peak(z)    =    11.6 (1+z)^0.03 ,  
    σ(z)    =    0.56 (1+z)^0.25 .


The parameters of the   model were fit to constrain both the z=0 SMHM and
the stellar mass functions at redshifts up to z=3, split between red and blue galaxies.
 obtained these parameters by running their model on
halo merger trees from the MS, with the same cosmological parameters (but
they express their halo masses using
h=0.7 instead of h=0.73, a 4 per cent relative difference that we, hereafter, neglect).

The  model of equation (<ref>) has the advantage 
of producing, through  stellar mass growth, a natural scatter
in the SMHM from the stochasticity of the halo mass growth. 
Moreover, by construction, stellar mass cannot decrease in the  model.
Since the  model is based on
time derivatives, it does not predict a minimal stellar mass, but
a minimal mass growth. Therefore, higher stellar mass growth can
only be caused by mergers. 
We
decreased their normalization E by 30 per cent  (as they suggest)
to account for the
loss of stellar mass by supernova explosions.




  §.§.§ Comparison of the analytical galaxy formation models





Table <ref> summarizes the analytical models used here, as well
as the semi-analytical models presented in Sect. <ref>.
The cosmological parameters are those of the
simulations that were used to calibrate to the observations.

 

Figure <ref> shows the comparison of the 4 models at redshifts
z=0 and the peak SFE at redshifts 1, 2, 4, 6, and 8. Note that only the BWC
model was fit by the authors up to z=8. 
The four 
models agree to first order, although 
there are some notable differences.
While the  model, by construction, has a symmetric F_ SFE in log v_ circ and in log M,
all other models are asymmetric, with a relatively faster decrease of
SFE at low circular velocity or halo masses, except for the
shallow low-end tail in the
 model at low redshift.
But recall that the F_ SFE of the  model is defined in a
different fashion (eq. [<ref>] vs. eq. [<ref>]). 

The peak SFE 
decreases with time in the   model,
increases with time in the   and  models, and is roughly independent of redshift
in the   model.
The circular velocity (respectively, halo mass) where this peak SFE is reached 
decreases (increases) slightly with time in the   model,
increases sharply with time in the  and  models,
and increases moderately with time in the  model up to z=4 (z=2) then
reverses and decreases with time increasingly faster.

At the high-mass end,
the decrease of F_ SFE with halo circular velocity at z=0
follows a  similar
slope (–2.5) in the physical  and empirical  
models,[Admittedly, the high-end slope of SFE versus M in the
  C11 and C+G models (e.g., the denominator in eq. [<ref>])  was a guess, but it is
 remarkably close to the 
  empirical relation of .] while the  model shows a somewhat shallower
decrease (slope –2). In contrast, 
the  model shows an increasingly rapid decrease, with a much steeper slope,
close to -9 at 1 percent of peak SFE.

At the low-mass end,
the rise of F_ SFE with halo circular velocity is fastest for the
  model (F_ SFE rises as v_ c^9 according to
  eq. [<ref>]) and the  model, for which the slope is even
  steeper than 9 at levels where SFE is at 1 percent of its peak value.
At z=0, the   model gives a low-end slope of 4.1 for F_ SFE vs. v_
  c, the   model has a much shallower slope
of 1.2, while that of the  model
slowly tends to
infinity (because of the lognormal SFE vs. halo mass), and at 1 per cent of the peak
SFE, the slope is 3.3 at z=0 and 2.1 at z=5.
At higher redshifts, while the   model retains the same steep slopes at
the low-end of F_ SFE vs. v_ c, the   model becomes shallower
with increasing redshift, while the   model becomes steeper with redshift.
Thus, at z=5, the 's  F_ SFE vs. v_ c relation is steeper
(slope 2.9)
than that of  (slope 2.1).



Figure <ref> compares directly the evolution of F_ SFE
between the 4 analytical models.
In particular, the  model has decreasing SFE in the last 8 or more Gyr
(depending on final halo mass), whereas the  model shows a flat SFE in the last
Gyr, while the  and  models appear to have increasing SFE at all times.
In models with F_ SFE that decrease with time such as  and
 at late times, a galaxy whose
halo mass remains constant in time (from an unlucky lack of mergers that would otherwise
make it grow) will see a decreasing F_ SFE at constant M, hence a
decrease in stellar mass.
In these cases, our forcing the stellar mass not to decrease comes into effect.
In contrast, the  model, based on time derivatives, is the only one
where stellar mass is guaranteed to
never decrease in time.



Figure <ref> illustrates how F_ SFE varies with both halo mass and
redshift for each of the 4 analytical models. Naturally, the essential
features of Figures <ref> and <ref> are recovered:
The higher peak F_ SFE at low redshift for the  model shows as a
whiter colour in Figure <ref>.
The sharp drops of F_ SFE at low halo mass of the  model at z=0 seen in
Figure <ref> are visible here as rapid transitions to black
colour at all redshifts (again, recall that F_ SFE is defined in a different fashion
for the  model). This is also the case for the  model. However, at high mass, the
drop of F_ SFE with mass at z=0 is only pronounced for the
 model (again seen at all redshifts).


Given that the intergalactic medium should be fairly cold before reionization
and then is gradually heated by photoionization from galaxies and quasars,
one expects that star formation  at low circular velocities (masses) should
only occur at 
higher redshifts.
The
 model has this low-mass effect built-in, with v_ post-reion > v_
  pre-reion. But, a striking feature of Figure <ref> 
is that the  model is the only one for which the halo mass at peak
F_ SFE has a clear maximum at z=3 (see also
Fig. <ref>).
In other words, the
 empirical model suggests that the effects of photoionization apply
to all masses, as the entire F_ SFE(M) is shifted to lower masses at
higher redshifts. 
This “reionization of the Universe” feature strongly seen in  and also
visible at the low end of the  model, is not seen in the  or
 models, because they were tuned to observations that did not extend far
enough in redshift, so the extrapolation to higher redshifts  (grey shaded
regions of Fig. <ref>) is incorrect.

Figure <ref> also expands on Figure <ref> by showing
that not only do models  and  show halo masses with decreasing
F_ SFE, but
models  and  also show decreasing F_ SFE for high halo
masses (12).

Moreover, Figure <ref>  shows that the  model is the only one that
displays a decrease at late times of the lower envelope of the halo masses
with efficient star formation. By definition, VYGs require a very late
jump in stellar mass, which should usually come together with a sudden jump
in halo mass (we will discuss this in more detail in

Sect. <ref>). 
One would therefore expect that, contrary to the other 3 analytical
models, the  model will lead to fewer VYGs, because, at relatively late
times, the
SFE will increase as the halo mass rises, causing star formation in the last
few Gyr before
z=0.08, which should make it difficult to have a very young stellar
population at z=0. 



  §.§.§ Galaxy merging in analytical models


The galaxies in our analytical models can grow in a quiescent mode, via
equation (<ref>) (or (eq. [<ref>] for the  model),
or by galaxy mergers.
The new stellar mass of the galaxy in halo i of mass M_i
is the maximum value of the masses given by the
quiescent growth and by the
growth by galaxy mergers:

    m(i,s)    ≡    Max[m(M_i,s), ∑_ prog m_
        old ] 
       =    Max[m(M_i,s),
    ∑_jm(j,s-1)]  ,
 
where the js are the indices of the progenitors of halo i from the
previous timestep, while
m(M,s)  
is the model stellar mass at timestep s for a galaxy in
a halo of mass M (eq. [<ref>]).
For the  model, the model mass in equation (<ref>) is understood to be
m = m_ old + f_ b F_ SFE(M,z) Δ M.

While our different analytical galaxy formation models generally predict
different stellar masses associated to both progenitors, they are run on the
same halo merger trees, hence  galaxy
mergers occur at the same time in each model (but involve different stellar masses).

Our analytical models cannot directly handle possible starbursts in
merging gas-rich galaxies. Instead, we
have implemented 2 different schemes to handle halo mergers that roughly
reproduce the situations of gas-rich and gas-poor mergers.

In our bursty galaxy merging scheme,
we delay subhalo (galaxy) merging by a suitable dynamical friction time,
t_ df
(see eq. [<ref>] below), measured from the last time when the two
haloes were distinct in the halo merger tree.
We do not model the subhalo
mass evolution (its decrease by tidal stripping), and thus the
satellite galaxy that originated from  this branch has a constant stellar mass
until it
merges with the central galaxy.
We do not consider satellite galaxies that merge later than z=0
because of disk storage issues and because our main focus is on central
galaxies.

With bursty galaxy merging,  
the central galaxy sees a boost in stellar mass at the
time of the halo merger, when the first term within the brackets of
equation (<ref>) dominates the second one.
It recovers the masses of the satellites only  after the dynamical
friction time has elapsed.

This

scheme thus creates a burst of star formation at
the first pericentre.
Note that hydrodynamical simulations of
merging galaxies indicate that while a starburst occurs at first pericentre,
there is usually another stronger one at the second pericentre, when the 2 galaxies
complete their merger (e.g. ), but this cannot be
modelled in our code.

In our quiet galaxy merging scheme,
we not only delay the galaxy merger by the dynamical friction time, but we also
modify the stellar mass at the time of the halo merger
(eq. [<ref>]).
For this, we subtract from the new halo mass the total mass ∑ M_k of all the branches
that 1) merge at the current timestep or before, and 2) contain galaxies
turning into satellites that do not have time to merge into the central
galaxy (after dynamical friction) by the current timestep.
The new central galaxy mass is then

    m(i,s) =  Max[m(M_i-∑ M_k,s),
    ∑_jm(j,s-1)]  ,

where the  js are again (as in the bursty merger scheme) the indices of the progenitors of halo i
from the
previous timestep.
With equation (<ref>), there is no boost of stellar mass at the time
of the halo merger, hence our quiet scheme represents well the situation of
dry galaxy mergers.

Following <cit.>, we adopt a dynamical friction time of

    t_ df(M_1,M_2,z) = 
      A M_1/M_2/ln(1+M_1/M_2) r_ v v_ v ,

where r_ v/v_ v=√(2/Δ(z))/H(z) is the orbital circular time divided by 2π at
radius r_ v, with Δ(z) being the overdensity relative to critical for
collapse, for which we adopt the approximation of <cit.>.
In equation (<ref>), A is a dimensionless constant that depends on
orbit eccentricity.
There is some debate on the value of A. 
<cit.> calibrated A=1.43 using hydrodynamical
cosmological simulations, while 
other values of A range from 0.58
<cit.> to 2.34 <cit.>, with
the latter value adjusted to better fit the observed galaxy optical
luminosity functions with SAMs.
We adopted A=2.34 to be consistent with the  SAM and the  model.


The dynamical friction time of equation (<ref>) is always longer than a
few Gyr. It is shortest for equal mass mergers, for which it
is A/(ln 2) √(2/Δ)/H > 0.36/H at all times, and is >0.47/H  and 
>6.2  Gyr for lookback times less than
1 Gyr (where we assumed the Millennium cosmology). Therefore, galaxies merge
several Gyr (hence 
many timesteps) after their halos merge. 

Since our focus is on VYGs, which should be gas-rich, we adopt the bursty
merging as
our primary galaxy
merging scheme, but will later compare to the quiet scheme.



  §.§.§ Practical issues for the analytical models run on the Monte-Carlo
  merger trees


We first run the halo merger tree code, which writes the trees to
files. The analytical galaxy formation models are applied
using a second code that is run on
each halo within each tree.  
This code starts at the highest
redshift of z=13 and works forward in time, producing the stellar mass assembly
history within the halo mass assembly history. 
We assign initial stellar masses to each halo according to the model mass of
equation (<ref>) for the first three models, and to  zero for the  model.

We analyse the stellar mass history summed up over all
progenitors of a particular halo at redshift zero. Had we only considered the growth of stellar mass of
the main progenitor, we would have measured the stellar mass assembly instead
of the star formation history. We discard all z=0 stellar masses below 10^3.

Ideally, since the 
 analytical models have been calibrated on halo
mass functions from different cosmologies, 
it would have been best for us to run the halo merger tree code separately  
for each of our 
analytical models, adapting the cosmological parameters on
which each model was built. This proved difficult to do in practice because
of the large disc space required, hence we ran the halo merger tree code only
once, i.e. with the same cosmological parameters (see
end of Sect. <ref>) for all 4 analytical models.

The different models employ slightly different definitions of halo mass.
Our Monte-Carlo halo merger tree is built with the 
<cit.> algorithm, which is based on masses of the Friends of Friends
halo membership in the cosmological simulation (S. Cole, private communication).
However, in the halo to stellar mass relations obtained from the empirical
abundance matching technique, the halo mass is defined at the spherical
overdensity of 200 times critical () or at 's
() virial value
(, ).
The mass used with the  model is simply that of the Monte-Carlo
halo merger tree, on which the model parameters were roughly calibrated, with
no pre-calibration on halo mass functions. 
These inconsistencies in the mass definition should have no
significant 
effect on the resulting analysis (e.g., log (M_ vir/M_200) is only  0.1).




 §.§ Semi-analytical galaxy formation models


We also considered two semi-analytical models (SAMs) of galaxy formation.
In SAMs, galaxies are modelled as discs plus bulges, each with their
characteristic sizes.
These codes include complex physical recipes to incorporate
gas cooling, star formation,
stellar evolution, 
feedback from supernovae and active galactic nuclei (AGN).
These models should be more realistic than the physical analytical models,
in particular because they usually
model the galaxy positions from those of the subhaloes in
the Millennium simulations.
However, nearly all SAMs fail to reproduce many of the
relations observed for galaxies, while the empirical analytical models are
calibrated to these observations. 



  §.§.§ <cit.>



We first analysed the recent, state-of-the-art SAM of
<cit.> that they ran on the halo merger trees
of the MS-II, rescaled to the Planck cosmology.
The  SAM continues the series of SAMs developed by the Munich team
<cit.>.
It also includes ram pressure stripping of satellites in cluster-mass haloes
and
 the formulation of <cit.> (eq. [<ref>]) to prevent gas
accretion on low-mass haloes. 

 tuned the free parameters  of their model to match
as best as possible the z=0 to 3 observations of the stellar mass functions
and the fractions of passive
galaxies. Such a tuning solves several problems in previous implementations of
the Munich model.
In particular, low-mass galaxies no longer form too early (in contrast with
), thanks to longer timescales for the galactic winds to
fall back into the discs they originated from.

We
extracted  from the Henriques2015a table of the Virgo – Millennium database of the 
German Astrophysical Virtual 
Observatory (GAVO)[http://gavo.mpa-garching.mpg.de/Millennium/Help]  
the following parameters:
the halo mass (corresponding to M_200) centralMvir, 
stellar mass stellarMass,
and the mean (mass-weighted) age massweightedAge, all at z=0 with
m>10^6   M_⊙.[The <cit.> model galaxy ages are
  not directly comparable to those of the other models, because
   consider mean stellar age, while for the
  other models, we consider the median age.]
This extraction yielded 1 175 987 galaxies.
We, hereafter, refer to this virtual catalogue and the corresponding SAM model as .
Since the analytical models described above fail to consider satellite
galaxies, we limit the  catalogue to the central galaxies.



  §.§.§ <cit.>


Finally, we considered a second SAM, by <cit.>, where
the baryonic physics is built on top of Monte-Carlo halo merger trees.
We chose this  SAM, because it 
has been run on two sets of merger trees: one <cit.> for the
ΛCDM cosmology (with Ω_ m=0.3, h=0.7, Ω_
  b=0.04, and σ_8=0.9) and another (;
) for  
a Λ Warm Dark Matter (ΛWDM or simply WDM) cosmology.
We will only discuss this SAM in the context of estimating the effects of WDM
on the fraction of VYGs.

WDM
is analogous to ΛCDM, but with the power spectrum truncated at high
wavenumbers, corresponding to 
scales below ∼ 1 Mpc, 
as expected for a m_χ∼ 0.75 keV thermal relic particle. 
Note that this choice for m_χ is fairly extreme, as the observed
abundance of dwarf galaxies around the Milky Way
leads to
m_χ >  2.3  keV <cit.>.
This low value of m_χ  sets an upper limit to the effects of
the suppression of primordial density fluctuations at high wavenumber.

For both
ΛCDM and ΛWDM cosmologies,  
the merger trees were run
using 21 final, equal-spaced  halo
log masses (9.5 < log (M_0/ M_⊙) ≤ 14.5), with a mass resolution
of 5× 10^7 for all final halo masses M_0. Each z=0 halo mass was run
100 times, leading to a total of 2100 trees (see Table <ref>).

In this SAM, gas is converted into stars through two main
channels: a quiescent accretion mode, occurring on long timescales (∼ 1 Gyr), 
and an interaction-driven mode, where gas – destabilized during
major and minor mergers and fly-by events – turns into stars on
shorter timescales (typically ∼ 10^7 yr). 
AGN activity triggered by galaxy
interactions and the related feedback processes are also included.

The  SAM was run on the CDM and WDM trees with the exact
same recipes and parameters for the baryonic physics.  We do not include it
in our standard analyses and figures, so as not to overcrowd our figures,
given that
it uses 100 times
fewer trees than our standard ones and the lower mass ones are poorly
resolved,
and is a much less used SAM 
than that of <cit.>, whose
output is publicly available, and which has the additional advantage of
being calibrated to the stellar mass functions at redshifts up to 3.



§ TESTS OF THE MODELS





 §.§ Stellar versus halo mass




In this section, we test our analytical models, by
analysing the z=0 stellar masses (m), halo masses (M),   and the
redshift z_50 (the corresponding stellar age or lookback time,
is called age_50) when half the mass in stars is formed.
For the  SAM, which does not easily provide the epochs when half the
stellar mass is formed, we consider instead
the mass-weighted mean age of the stellar
population, which we also denote as z_50.

The  stellar mass – halo
mass (SMHM) relation is a fundamental test of galaxy formation models.
Figure <ref> illustrates  the SMHM, according to the 4 analytical models
as well as the
Henriques and Menci SAMs,
and how galaxy ages depend on their location in the SMHM relation.
Note that the Menci SAM uses 0.25 dex steps of , but their output is
shown at z=0.12 instead of z=0, which leads to somewhat less discretization.

The first aspect to note is that only the  model leads to an SMHM
that matches well the observations of 11 nearby, highly-inclined,
isolated dwarf irregulars 
(black line in Fig. <ref>)
as compiled by <cit.>.
The , , and especially  models
predict lower stellar masses at given halo mass than observed by
, as also does the  SAM.
On the other hand, the
Menci SAM predicts higher stellar masses at given halo masses than found by
.

The better match to observations with the  model appears to be
a consequence of its inclusion of
a change in slope of the SMHM at the low end, which was
not incorporated in the other analytical models. Interestingly, 
<cit.> noticed that the SMHM they obtained for their sample of 11
dwarf  galaxies follows the
SMHM related predicted from the extrapolation to low masses of the abundance matching of
the SDSS stellar mass function (SMF) of <cit.>
with the cosmic halo mass function. This explains the agreement between
their SMHM and that found for , who had also used the SMF
of <cit.>. In comparison,  adopted the same SMF, but
adopted a too restrictive single-slope low-end SMHM parameterization.
Finally,  adopted a different set of SMFs.

The second interesting feature is that the most massive galaxies have formed
their stars earlier than the other galaxies (except for the very low-mass
galaxies in the  model),
i.e. one is witnessing downsizing. In contrast the , Menci,
and especially the  model also predict that the
smallest galaxies are old, i.e. upsizing.
We will analyse in more detail in Sect. <ref> how the typical stellar ages of  galaxies
vary with their z=0 stellar mass.

The third interesting aspect of Fig. <ref> is that, at intermediate masses,
the galaxy stellar masses  derived by running the analytical models on the Monte-Carlo
halo merger trees (points) are often close to those directly obtained from the
models for the redshift corresponding to z_50 (curves of similar colour).
However, 
at the high-mass end, the stellar masses predicted from running the models on
our halo merger trees
are greater than 
those predicted by these two models at z=z_50 (or at any z). For the
, , and  models, this extra stellar mass at high halo mass, which
is clearly seen in the 3 corresponding insets of Fig. <ref>,
suggests that we are in a regime where the analytical stellar mass decreases in time,
which our implementation forbids (see end of Sect. <ref>). 
In Figure <ref>, the ridge of highest SFE occurs at decreasing
halo mass at later times before z=0, except in the  model (where it is
roughly constant). However, in this physical  model,
galaxy mergers can boost the stellar mass beyond the analytical prediction
<cit.>. 
Thus, the  model highlights the role of mergers at high masses, while
showing negligible effects of mergers at lower masses, where the analytical
prediction matches well the outcome of the model run on the halo merger
trees.
The importance of mergers at high masses was previously noted by
<cit.> and <cit.>, who studied the merger rates
of galaxies, by <cit.> who studied  (at better mass
resolution)
the fraction of
stellar mass
acquired by mergers, as well as 
by <cit.> who analyzed the observations of galaxy properties.

The fourth interesting feature of Figure <ref> is that the scatter
in the SMHM relation strongly depends on the model.
The   and  models produce virtually no scatter for halo
mass below 10^12.

In contrast, the  model and Menci SAM produce a large scatter in the SMHM relation at
low mass. The  SAM produces scatter at both low mass (where mass
resolution effects become significant) and high mass (perhaps from the
stochasticity of galaxy mergers).
While the scatter of the stellar mass of central galaxies in the , and Menci SAMs are
low at intermediate halo masses, it is 
even lower in the EAGLE hydrodynamical simulations <cit.>. 
The large scatter in the  SMHM  is clearly due to its higher SFE of low-mass
haloes before reionization  (see Fig. <ref>) combined with the
stochasticity in recombination epochs in the  model. 
The    model  produces some scatter at all masses, as expected from its
use of mass evolution rates instead of masses themselves.



 §.§ Stellar mass function



In theory, the SMF, dn/ dm,
can be estimated by integrating over
the cosmic z=0 halo
mass function (known from cosmological simulations), dn/ dM:

    dn dlog m   =   ∫ dlog M  d^2n dlog M  dlog m
       =   ∫ dlog M ( dn dlog M) ( dn dlog m)_M  .


Our choice of equal numbers of haloes (trees) in logarithmic bins of mass
produces an unrealistic halo mass function and would thus lead to an
unrealistic galaxy SMF.
We can nevertheless compute
the SMF (in bins j of log stellar mass)
as a sum over bins (k) of halo mass, by
normalizing our flat halo log mass function by the cosmic halo mass function.
We thus weight the galaxies i as the ratio of the z=0 cosmic halo mass function
to the available z=0 halo mass function returned from the Monte-Carlo
trees, i.e.

    w_i(M_k) = ( dn/ dlog M)_M=M_kΔ n / Δlog M_k .

We can then write the measured SMF in bins of constant Δlog m_j as


    g(log m_j)    ≡    ( dn dlog m)_m=m_j
       =   Δlog M_k ∑_k (Δ nΔlog M_k)w_i(M_k)
    (Δ nΔlog m_j)_k 
       =   
    1
    Δlog m_j ∑_i ∈ j w_i(M_k)  .


The second equality of equation (<ref>) makes use of
equations (<ref>) and (<ref>) and of the constant
value of Δlog M_k.

We adopted the z=0 cosmic halo mass function of <cit.>, 
computed using HMFCalc[http://hmf.icrar.org] <cit.>
with
the cosmological parameters used in the Millennium Simulations, and for an
overdensity of 200 relative to critical.
[The choice of the model for the halo mass function has little
  impact on the relative weights used to compute the fraction of young
    galaxies in Sect. <ref> below; it only affects the SMF in
    Fig. <ref>,
    and only slightly so.]
[<cit.> has shown that the weights of
  equation (<ref>) are slightly incorrect, because backsplash
  (sub)haloes are double counted. The effects are relatively small (at most 10 per
  cent of haloes of a given mass above 10^11 are backsplash
  according to fig. 4 of ).] 

 

Chabrier03
Yang+09

The resultant z=0 galaxy SMFs produced by our Monte-Carlo
runs for the analytical models
are shown in Figure <ref>, as well as that of the  SAM
(restricted to central galaxies) and the observed  central galaxy
stellar mass function of <cit.>.
The 
figure indicates that the  model fails to reproduce
the knee of the SMF at log(m/ M_⊙) ≃ 10.8. 
In the range 7 ≤≤ 8, the logarithmic slopes of the model SMFs
are
≈ -1.28, –1.36, –1.66,  -1.18, and -1.51 for the , , ,
, and  models, respectively.
In comparison, the slope of the observed 
SMF of  
is ≈ -1.65 at 8.7 ≤log (m/) ≤ 9.3.
Hence, the observed low-end slope is best (worst) matched by the  ()
model.
With the quiet galaxy merging scheme, the SMFs from the analytical
models
are similar, but
shifted down by 0.4 dex, thus matching less well the observed one.



 §.§ Age versus mass



We now study how the epoch, z_50,
when half the
z=0 stellar mass is formed in galaxies
varies with the z=0 stellar mass.
This is illustrated in
Figure <ref>
for the four analytical models of F_ SFE(M,z),
as well as for the
 SAM (for which the redshift corresponds to the arithmetic mean age of
the stellar population of the galaxy).



All models show downsizing at the high-mass
end, where z_50 increases with stellar mass. This is
 due to the negative  slope of the SFE - halo mass relations at high halo masses
(Fig. <ref>): the galaxies with the highest stellar masses
 today live in high mass haloes that grow faster than low-mass ones
 ( and black curves of Fig. <ref>) and also
 faster than the mass of the peak SFE (whiter ridges in Fig. <ref>),  
leading to early quenching of star formation
(see ).

At low galaxy stellar mass, the behaviour of z_50 versus mass
relation differs between models.
While the
   and     models show 

a fairly flat typical median age versus
mass for log (m/) < 9.5, 
the  and   models lead to typically earlier star formation at increasingly lower galaxy
masses, i.e. upsizing, for log (m/)<9. 
Finally, the  SAM leads to an intermediate behaviour with weak upsizing for
 log (m/)<10.

In particular, the  model leads to very old low-mass (m < 10^7)
galaxies,
most of which are formed by z=8 (see upper left panel
of Fig. <ref>). Since these low-mass galaxies have final halo
masses log (M/ M_⊙) < 10.2 (upper panel of Fig. <ref>), one can see that
 these old galaxies were quenched by the rising halo mass linked to the
 rising entropy and  temperature of the intergalactic medium during the epoch
 of reionization (see upper left panel of Fig. <ref>).
This is
qualitatively consistent with the star formation histories of low-mass
(ultra-faint) Local Group dwarf spheroidals <cit.>. However,
such low-mass
galaxies with very old stellar populations are currently difficult to observe
beyond
the Local Group, and may be limited to satellites that are ram
pressure-stripped by the fairly hot gas of their more massive centrals.

In the same mass range, Figure <ref> also displays a younger
galaxy population for the  model,
pretty similar in age to intermediate mass galaxies. These galaxies were formed
after the epoch of reionization. 
These two (very old and younger) galaxy populations for the  model are
related to 
our two v_ reion values.  Their overlap in mass is caused by our adopted
stochasticity of z_ reion. This model is thus particularly well-suited 
in exploring the epoch and duration of reionization through the age distribution of dwarf galaxies.

The other source of discrepancy between the models is the median redshift of
star formation of
galaxies at intermediate stellar masses (8 ≤log(m/ M_⊙) ≤
10):
roughly 2.5, 1, 0.7, 0.6, and 0.2 (the corresponding lookback times being roughly 
11, 7.5, 6.5, 5.5, and 2.5 Gyr) for the , , , , and  models,
respectively.

Finally, the scatter in z_50 is high in the  model (except at low
mass) and in the  model (at low mass),
intermediate in the  and  models, and small in the  model.



§ FRACTION OF VERY YOUNG GALAXIES





 §.§ Estimation of (weighted) fractions of very young galaxies and their
  uncertainties
 

Given the data of Figure <ref>, we can derive the fraction of
VYGs  (i.e., z_50<0.08,
corresponding to median stellar ages of less than 1 Gyr).
We first need to normalize the raw fractions of VYGs by correcting
our flat  halo mass  function (in terms of log M) to the cosmic halo mass
function, as  described in Sect. <ref>.
The fraction of VYGs in the jth bin of log stellar mass can
easily be obtained from equation (<ref>) to yield (dropping all
the log terms for clarity)

    f_ VYG(m_j) = 
    ∑_i∈ j, VYG w_i(M_k)∑_i∈ j w_i(M_k)

where w_i(M_k) is, again, the weight for the galaxy i in the kth bin among
the N_k halo log mass bins, required to
normalize the measured SMF (eq. [<ref>]).

The uncertainties on f_ VYG cannot be
 derived 
with the usual binomial formula

    σ^2[f_ VYG(m_j)] = √(f_ VYG(m_j) 
    [1-f_ VYG(m_j)] N_ eff) ,

because it is not clear what value to adopt for N_ eff, the effective
(given the weights) number of galaxies in bin j of stellar mass. Instead, we
rely on bootstraps: for each bin, j, of stellar mass,
we bootstrap the sample 
to obtain B values for 
 f_ VYG using equation (<ref>), and 
then determine the uncertainty from the standard deviation of these B
values of f_ VYG. 
Unfortunately, the bootstrap method cannot provide errors when there are no VYGs.
In that case, we replace the numerator of f_ VYG in
equation (<ref>) by w_K/2, where K is the bin of weight (or halo
mass) that contains the youngest of all galaxies in the stellar mass bin
(by construction, the youngest galaxy will have a median age greater than 1 Gyr).

This weighting scheme is evidently not used for the    model, since it
is applied to the Millennium simulations, which have realistic halo mass
functions.
In this case, we adopt the standard binomial errors (eq. [<ref>],
replacing N_ eff by the number N of galaxies in the bin j of
stellar mass). 
When there are no VYGs  in the bin j of stellar mass (i.e.,
f_ VYG = 0), we adopt
the <cit.> upper limit:[See the Wilson score in the
  Wikipedia entry on the Binomial proportion confidence
    interval.]

    f_ VYG^ upper = z^2/N_ eff 1+z^2/N_ eff

for z =1.65 standard deviations, corresponding to a 95 percent confidence
level for this upper limit.



 §.§ Fractions of very young galaxies for the different models


Wilson1927
Figure <ref> shows our predicted fractions of VYGs
 as a function of galaxy stellar mass, for the analytical 
 galaxy formation models
with bursty galaxy merging
 and for the  SAM.

The ,   and  
 models predict a fraction of VYGs in the range of 0.3 to 4
per cent for a
 stellar mass range of 10^6 to 10^10.5, with maxima of 2 to
 4 per cent at log (m/) ≈ 9.5. The 
 physically motivated  model  predicts the lowest fraction of VYGs, 
with a peak at 0.7 per cent at = 9.0, and no plateau at lower masses.

Interestingly, for virtually all masses, the four models are ordered in the same way,
with increasingly higher fractions of VYGs at given
stellar mass for
the , , , and  models.
We will present a simple model in Sect. <ref>
to explain our results for the first three models.

 

In contrast, the  SAM leads to very low fractions (less than 0.03 per
cent) of VYGs (defined here as mean mass-weighted age below 1 Gyr) at all
masses where it is well resolved (≥ 7.5). The  model thus
appears to be discrepant  
with the 4 analytical models, with 30 to 800 times
fewer VYGs predicted by the  SAM at intermediate stellar masses (e.g.
log m = 10^9). 
While these predictions are for the central galaxies in the  SAM, we
will see in Sect. <ref> below that satellites and centrals have
similar fractions of VYGs, at these intermediate masses.
At masses below 10^7.5,
the  satellites
are typically 5 times less likely to be VYGs than are centrals.
However, we do not trust the  results below 10^7, because
the  SMHM relation, displayed in Fig. <ref>, shows a
suspicious flattening at log (m/ M_⊙) ≈ 6-7, which appears to
be caused by the limited mass resolution of the N body simulation on which
the  SAM was run.



 §.§ Dependence on critical age




One may ask how the fraction of VYGs depends on our choice of 1 Gyr for
the greatest allowed age for very young galaxies.
Figure <ref> shows the cumulative distribution function (CDF) of
galaxy ages for the different models (using bursty galaxy merging for the
analytical ones) in three bins of stellar mass.
This allows us to measure how the VYG fraction scales with critical ages not
necessarily equal to 1 Gyr.
The slopes of the curves near 1 Gyr indicate that the fraction of VYGs
roughly scales
as age_50^3/2 for most models, but is even more sensitive to the
critical age  (slope of
roughly 2) for models  and  at high stellar mass.

Figure <ref> also indicates wide variations between the models for
the ages of the youngest, say, 1 per cent of galaxies (horizontal line).
For the intermediate bin of stellar mass (log (m/ M_⊙)=9.5±0.75,
long dashed lines), the ages of the 1 per cent youngest galaxies are 
0.25, 0.55, 0.65, 1.5, and 3 Gyr for the , , , , and
 models, respectively. The discrepancy is as pronounced for the high-mass
galaxies (log (m/ M_⊙)=11.0±0.75, solid lines):
0.65, 1.5, 2.2, 4, and 7 Gyr, in the same order of models, except that
 switches with  for the highest age.



 §.§ Dependence on the galaxy merging scheme





Figure <ref> compares the VYG fractions versus mass for the two
galaxy merging schemes. It clearly shows the lack of sensitivity of the
VYG fractions predicted by the  and  models to the presence of a burst of star formation 
associated with the halo merger.
However, the merging scheme has an important effect
on the  and  models.
For the  model, the fraction of VYGs is above 0.1 per cent only in a
  narrow range of stellar masses around =9.5
  for the quiet galaxy  merging scheme, with
  negligible fractions at masses lower than =8.5.
  For the  model, the fraction of VYGs is mostly reduced at intermediate
  and high stellar masses for the quiet
  merging scheme.
  With the quiet merging scheme, only the  and  analytical models reproduce the
very low VYG fractions obtained with the 
Henriques SAM for stellar masses 7.5 ≤≤ 9 and 9.5 ≤≤
10, respectively,  while the other
analytical models predict
higher VYG fractions in these mass ranges. 
We will discuss the reasons for these behaviours in Sect. <ref>.  



§ DISCUSSION





 §.§ Limitations of the galaxy formation models



The galaxy formation models that we considered to assess the frequency of
VYGs have several limitations.

The analytical and  SAM were calibrated to observational data
in a limited range of redshifts and masses (highlighted in
Figures <ref> and <ref>, respectively, for the
analytical models). The fractions of VYGs at
  stellar log masses (in solar units) lower than 8.8 (), 7.4 ( and ), 10.0 ()
  are obtained by extrapolating the models.  Similarly, the  SAM is
  calibrated at z=0 to > 7.2, and at lower stellar mass there is a
  suspicious flattening of the SMHM in Figure <ref>, which may be
  due 
to insufficient mass resolution.


 
In fact, the mass and spatial resolution of the cosmological
simulation on which the  SAM is run  has a strong effect on the fraction
of VYGs at low galaxy masses.
Figure <ref> compares the fractions of VYGs versus
stellar mass found by the  SAM, when run on the high resolution MS-II
simulation and
on the 125 times lower mass resolution MS simulation.
While the fractions of VYGs are similar at intermediate mass, the MS produces
increasingly higher VYG fractions at increasingly lower masses, where it is
poorly resolved, than does the
MS-II.
At m=10^7.25, at the limit where the MS-II is still properly resolved,
the  SAM run on the MS leads to VYG fractions up to 3 orders of magnitude
higher than found when it is run on the MS-II.  The  SAM run on the poor mass
resolution MS forms galaxies at later times, which are thus more likely
to become VYGs.

The  SAM
was run on dark matter
simulations that were rescaled in space and time to reproduce the
large-scale statistics for a more realistic cosmology (Planck
vs. 1st-year WMAP). Tests by <cit.> suggest that the small-scale
quantities such as halo masses, bulk velocities, and the luminosities of their brightest galaxies
are changed by only 10, 5 and  25 per cent, respectively.

However, this rescaling may not
reproduce the astrophysical processes at small scales, i.e. in the nonlinear
regime, in particular the galaxy orbits including dynamical friction, as well as
tidal stripping of galaxies by their host groups.

Compared to the  SAM, our 4 analytical models have the advantage of a much
superior halo mass resolution: our Monte-Carlo halo merger trees are built
for z=0 halo masses down to 10^7, corresponding to the mass of a
single particle in the MS-II simulation on which the  SAM was run.  
Furthermore, the trees resolve branches down to 10^3.
Note that  our halo merger trees
were 
only tested by <cit.> down to log (h M_0/) = 11.

On
the other hand, our analytical models have also their own drawbacks.
The careful reader will have noted that some of the empirical models we use have
been calibrated assuming slightly different sets of cosmological parameters
(see Table <ref>). 
While models  and  have values of Ω_ m, h and σ_8
that are within 1 per cent of one another,
model  has values of Ω_ m, h, and σ_8 that are
respectively 8 per cent below, 4 per cent above, and 10 per cent above those
from the other two models.
Since the VYG
fractions of the  model are in much better agreement with those of the
 model than with those of the  model  (Fig. <ref>), it appears that the VYG
fractions are more sensitive to the baryonic physics than to the details of the
cosmological parameters. 

Our analytical models do not follow the satellite galaxies that survive by
z=0 against merging into the central galaxy of their halo.
However,  Fig. <ref> shows that there are no large
differences in the VYG fractions of 
centrals (open circles) versus
satellites (triangles) in the  SAM,
especially for the
SAM run on the higher resolution MS-II, except for
log( m / ) ≤ 7.5, where the satellites are typically 5 times less
likely to be VYGs. 

Moreover,
with the bursty galaxy merging scheme applied to our first three
analytical models, 
the stellar mass of the 
central galaxy is boosted at the time of the halo merger instead of when the
satellite mergers into the central, which is at least a third of a current Hubble time Gyr later (see
Sect <ref>).
This boost is the consequence of the
larger halo mass of the central galaxy after the halo merger together with
the increasing stellar mass of the analytical models with halo mass (given
the shallower than –1 slopes in F_ SFE versus M 
in the bottom panel of
Fig. <ref>). This will be discussed in detail in
Sect <ref>. In the fourth
model (MCP), the stellar mass growth is directly linked to the halo mass
growth. So, all four models naturally have star formation associated with
halo mergers.
This should be realistic for haloes that are rich in cold gas and
that merge in nearly head-on orbits. If the mergers are off-center, one ought
to delay the starburst by the dynamical friction time, and if the mergers are
gas-poor, there should be no burst from the halo merger.
Since our focus is on low and intermediate mass galaxies, which tend to be
gas rich (fig. 11 of  and references therein), and whose progenitors
must have been even more  gas-rich, our galaxy merging scheme appears to be
sufficiently realistic. 

A final worry of the analytical models is that they may not properly consider
halo mass accretion histories (MAHs).
Indeed,
galaxy properties appear to be not only be related to halo mass, but also to  halo
MAH.
In cosmological simulations, the halo MAH influences the
clustering of haloes <cit.> and their concentration
<cit.>, a process known as assembly bias.
Such galaxy assembly
  bias is observed from mass modelling of low redshift galaxies traced by their
satellites, which  leads to red galaxies  (i.e., with older stellar
populations) having higher concentration haloes than blue (young stars)
galaxies of the same  stellar (or halo)
mass <cit.>.
Galaxy assembly bias can be detected in the outputs of SAMs <cit.>. 
Our Monte-Carlo halo merger trees incorporate assembly bias to a large extent,
since different haloes have different MAHs.
But it is not clear that the strong dependence of halo assembly bias with the
large scale environment (as predicted by ) is implicitly
incorporated
in our halo merger trees.




 §.§ Simple modelling of the fraction of very young galaxies in the analytical models


The relative importance of VYGs in the analytical models can be assessed
from first principles.
According to our definition, VYGs are produced when the stellar mass, summed
over all progenitors, increases 
by more than a factor of 2 in the last Gyr. This growth in stellar mass
occurs in two manners: 1) the quiescent growth
from the models in the absence of halo mergers (i.e., at fixed halo mass); 2)
the growth by halo mergers and corresponding galaxy mergers.
Below, we build a toy model of the stellar mass growth of galaxies that is
based on the simplifying assumption that this growth is independent of the
past history (in other words our simple model is Markovian).



  §.§.§ Model of quiescent growth during last Gyr




In the absence of mergers, VYGs
are produced if the quiescent growth in stellar mass from the analytical
model is a factor of 2 since 1 Gyr.
The stellar mass will vary during the last Gyr as

    Δ m_ quiescent   =    m(t_ lookback=1  Gyr) - m(t_ lookback=0)
    
       =   ∫_t_0-1^t_0
      Max(0,  dm_ model dt)  dt  ,

where t_0 is the age of the Universe expressed in Gyr,
m_ model(M_0,z) is given in equation (<ref>), and the halo
mass is fixed to the  z=0 value.
If the quiescent growth is over 50 per cent, then all galaxies should
be VYGs!

The gold band in Figure <ref> highlights the quiescent growth of
the analytical models during the last Gyr.
While model  has
quiescent growth in the last few Gyr, models  and  lead to negative
growth in the last Gyr. 
Solving the integral in equation (<ref>), we display in
Figure <ref> (solid lines) the relative jump in stellar mass during the last
Gyr in the absence of halo mergers.
Only the
 model shows significant quiescent growth at some halo masses, with a peak growth of
25 per cent for =11.2 (thick green curves in
Figure <ref>). 
In comparison, the  has no quiescent growth, 
while that of the  model is small for <9.4 and zero at higher halo mass.
The  model leads to zero stellar mass growth when the halo masses are
fixed.

We can also consider the hybrid evolution of stellar mass combining the
analytical model m(M,z) using the median evolution of halo mass
instead of fixing the halo mass to its z=0 value. Although the median growth of halo
mass is small (typically 0.005 dex, increasing with final halo mass), this
halo evolution can make a difference in the stellar mass growth. 
Re-computing the integral of equation (<ref>) assuming the median
evolution of halo mass (black curves in Fig. <ref>,
instead of fixed halo mass), we find 
 (dashed curves in Fig. <ref>)
that
the relative growth of stellar mass is boosted slightly (to 28 per cent
for  at =11.2) or strongly (for the  and  models).

The extremely low quiescent growth of  stellar mass in the analytical models
implies that mergers are necessary to produce VYGs in these 4 models.
But the  model at =11 requires less of a boost from mergers than it
does at other final halo masses, as well as compared to
the other models.



  §.§.§ Model of final Gyr stellar mass growth for bursty merging


We next consider the growth of stellar mass, summed over the progenitors,
during a halo merger
in the case of bursty merging.
To keep things simple, we consider (only) two haloes of
masses M_1 and M_2, with respective stellar masses m_1 and m_2,
merging
at timestep s to reach a new mass of M=M_1+M_2.
This simple model does not account for multiple halo mergers in a given timestep.

Considering that the stellar mass of the satellites are frozen,
the stellar mass of the pair can be written

    Δ m_ pair   =    Max [ m(M_1+M_2,s),
    m(M_1,s-1)
     ]
    - 
    m(M_1,s-1) 
       =    Max[m(M_1+M_2,s)-m(M_1,s-1),
    0]
     .


Since the
dynamical friction time (eq. [<ref>]) is always greater than 6.7 Gyr
for lookback times less than 1 Gyr (Sect. <ref>),
the satellites that were involved in halo mergers in the last Gyr do not have
time to merge with the main galaxy. 




Figure <ref> shows how the stellar mass of the merging pair
evolves from the previous timestep to the halo merger one using
equation (<ref>), for the , , and
 models[We 
do not apply this test to the
        model, because we do not know how to determine the stellar mass
       of the pair in the normalization to relative variations.].
Since we do
not follow satellites,
we thus normalize the jump in stellar mass by the stellar mass
m_1 of
the first galaxy.
While the z=0 SMHM matches quasi-perfectly the model prediction for the  and
 models, it predicts less stellar mass for the
 model at ≤11.2 (Fig. <ref>), because the stellar mass
is frozen when F_ SFE reaches its maximum over 8 Gyr ago
(Fig. <ref>).  
The excess log stellar mass in
the  model is well fit (to 0.06 dex rms accuracy)  by 

    Δlog (m)    =    1.100 - 1.678 log (M
      10^10 )
    
                 + 0.633 log (M 10^10 )^2  .
 

We thus correct our  model stellar mass  using equation (<ref>). 

Given that 1 Gyr corresponds to 3 timesteps, the galaxy has three
chances to boost its stellar mass by a factor 2 by mergers.

The brown horizontal lines in Figure <ref> indicate relative stellar
mass jumps amounting to a doubling of stellar mass in a single halo merger or
in 2 or 3 mergers of equal halo mass ratio.
The intersection of the curves with these lines provides the minimum halo mass
ratio of a single merger, or 2 or 3 halo mergers of the same mass ratio, 
to allow this doubling of
stellar mass.
Three halo masses are shown that, through the SMHM relation of
Fig. <ref> correspond to stellar masses in the rough range 7.5  10.6 (depending on the model). 

For the  and  models, the jumps in halo mass  for
=10.5 (dotted curves in Fig. <ref>)
are too small to reach the minimum necessary boost from 2 identical mergers, even for 1:1
mass ratios.
This means that it is virtually impossible to double
the  stellar mass in the last Gyr for these halo masses.
For the  model, the plateau of the relative stellar mass jump at low
halo mass ratio, for low and intermediate halo masses,
suggests that this model may be able to produce low and intermediate stellar
mass VYGs through
minor halo mergers.


The expected numbers of halo mergers per 350 Myr timestep can be
obtained from theoretical analyses <cit.> or from
the analysis of cosmological simulations <cit.>, both of which
are very similar.
We adopt the
differential halo merger rate per unit redshift of 

    R(M,μ,z)    ≡    d^2N_ merge dμ  dz
       =    
     A (M M_ ref)^b
     μ^c exp[(μμ_ ref)^d] (1+z)^e  ,

where μ=M_2/M_1 is the halo mass ratio (with M_1+M_2=M) and where
A=0.010, M_ ref=10^12, 
b=0.133,
c=-1.995,
μ_ ref=0.00972,
d=0.263,
and
e=0.010.
Integrating over halo mass ratios, we first infer that in the 350 Myr
timestep, 
a halo of mass  M = 10^10, 10^11, or 10^12 is expected to be
involved in respectively ⟨ N_ m⟩ = 2, 3, or 4 mergers
with  halos  of mass ranging from 10^-4 M (the extreme mass ratio allowed in our
Monte-Carlo halo merger tree, see Table <ref>) to M.
Hence, halo mergers are ubiquitous in our halo merger tree, although there
are non-negligible probabilities (100 exp(⟨-N_ m⟩) = 14,
5, and 2 per cent respectively for 
=10, 11, and 12) that a halo suffers no merger in a given timestep.
We note here that according to the tests of <cit.>, the
<cit.> code reproduces the halo merger rates of  to better than 20
per cent for all mass ratios for low- and intermediate-mass haloes 
(M=10^11 h^-1
and 10^13 h^-1, at z=0).
 
We can go one step further and integrate the halo merger rate. 
Let g_ q(M) and g_ m(M,μ) respectively be the quiescent
relative growth over 1 Gyr and the relative growth through a halo merger.
The predicted fractions of VYGs obtained in two mergers can then be written
(dropping the dependence of halo mass and redshift in the expression for the
halo merger rate R for clarity)

    f_ VYG = (0.08 3)^2 ∫_0.0001^1  R(μ_1)   dμ_1  ∫_μ_2^ min^1 
        R(μ_2)  dμ_2  ,

where the factor in front of the integral in equation (<ref>)  is
for the conversion from Δ z = 1 to 
one-third of a Gyr,
while
the lower limit of the first integral is the resolution of the Monte-Carlo
halo merger tree that we have used (Table <ref>).
In equation (<ref>),  
 μ_2^ min is the solution of the equation 

    (1+η g_ q) {1+ g_ m[(M 1+μ_2),μ_1]} [1  + g_ m(M,μ_2)]=2  ,

for μ_2,
with η=1/3 to account for the fact that only one-third of the Gyr
interval is available for quiescent growth (since we do not bother to perform
the triple integral required for the three timesteps).




The curves in Figure <ref> show the fractions of VYGs predicted by our
simple model (which involves no Monte-Carlo halo merger trees), as given by
equation (<ref>).
The figure captures the basic features of VYG fractions we found with the analytical
models applied to the Monte-Carlo
halo merger trees
(Fig. <ref> and shown as symbols in Fig. <ref>).
In particular, the  and  models are both remarkably well reproduced:
the  model shows a few per cent VYGs at low and intermediate 
stellar masses and a sharp drop at higher masses
(with the simple model overpredicting the Monte Carlo halo merger tree
calculations by a factor of 2).
Similarly, our simple model catches the details of the  model.
However, our model is less successful in predicting  the VYG fractions with the  model,
as it misses the 0.001 fraction of VYGs at < 8.
(Had we not corrected for the excess stellar mass of
the  model galaxies compared to their model predictions,
eq. [<ref>], the predicted  VYG fraction would be a factor 3 above
that of the  model at low mass.)
These discrepancies may be attributable to the simplicity of
our model and to the large scatter in the SMHM of the  model at low mass.
Nevertheless, the agreement for the  and  models of our simple VYG model with the VYG fraction
measured by the full model on the Monte-Carlo halo merger trees suggests that
the galaxy histories before 1 Gyr matter little in predicting the  and
 VYG
fractions, while the  VYG fractions, on the contrary, appear to depend on past history. 




  §.§.§ Model of final Gyr stellar mass growth for quiet merging


For the quiet merging scheme, we can better understand the jump in stellar
mass with a toy model based on equation (<ref>).

Consider a halo of mass M_0 at
lookback time t=t_0 (corresponding to z=0) and assume that its stellar
mass m_0 is given by its model stellar mass m_0=m(M,0).
This corresponds to the first
term within the brackets of 
equation (<ref>) dominating the second one.
Let F_i(M) be the fraction of  halo mass M at lookback
time t_i (corresponding to redshift z_i)
that came from branches that merged at lookback time ≥ t_i and whose host
satellites are not expected to 
merge (after dynamical friction) before lookback time t.
In other words, F measures the fraction of halo mass coming from
branches with surviving satellites. 
The stellar mass at z=0 can then be written

    m_0 ≡m_0 = m ([1- F_0(M_0)] M_0, 0)  .

At lookback time t_1=1  Gyr, corresponding to redshift z_1,
we write the stellar mass summed over the
progenitors of the z=0 galaxy as

    m_1 = m ([1- F_1(M_1)] M_1, z_1)
      + ∑ m_ sats ,

where M_1 is the mass of its main progenitor at lookback time t_1.
The last term in equation (<ref>)  is the total stellar mass in satellites that merge with
the central between lookback times t_1 and 0
(recall that the stellar masses of satellites are frozen). 



Let  F_i(M) be the median value expected for the values of F_i(M)
and η_i(M) = F_i(M)/ F_i(M), for i=1 and 0.
Thus, η_i(M) represents the excess, relative to the median, of the mass fraction  of halo of mass M at
lookback time t_i that originates from branches whose galaxies would survive as
satellites until a lookback time t < t_i. A halo with η > 1 is one
that could be
considered fortunate to have a higher fraction of its mass than usual from branches with surviving satellites.

The condition for obtaining a VYG is m_0 > 2 m_1 using
equations (<ref>) and (<ref>). 
Solving this inequality is difficult, because of the range of redshifts
when haloes merged with the main progenitor of the current halo, yet their
satellites survived. So, for simplicity, we assume that the satellite term in
equation (<ref>) contributes to a fraction f_ sats of the mass
m_1.
Moreover, we adopt for M_1 the median mass M_1 expected for the
z=0 halo mass M_0. The black curves of Figure <ref> indicate a
very slow typical growth of the main halo progenitor in the last Gyr, of order of
0.005 dex (Sect. <ref>). We shall see (Sect. <ref> below) that contrary to the case of
bursty galaxy merging, VYGs are not related to rapid jumps in halo mass with the quiet merging scheme.  

With these assumptions and
equations (<ref>) and (<ref>),
the condition for a VYG is

    .
         m([1-η_0(M_0)  F_0(M_0)] M_0,0)
        
          >  2 1-f_ sats m([1-η_1(M_1)  F_1(M_1)] M_1,z_1)
       .
    .


Of course, since VYGs are predicted to be rare at all masses,
typical values of F(M) (i.e. η_0 = η_1 = 1)
will not lead to VYGs according to equation (<ref>). 
Instead, the requirement for obtaining a VYG is to have
η_0 < 1 and/or η_1 > 1, i.e. a final halo with a low mass fraction
from branches with surviving satellites and/or progenitors at 1 Gyr lookback time
with high mass fraction from  branches with surviving satellites (to 1 Gyr
lookback time).
Here, we assume
η_1=1 and determine the maximum value η_0^ VYG of η_0 for obtaining a
VYG. If  η_0^ VYG is only slightly below unity, VYGs will be easier to
form than if it is much lower than unity.



The median value F(M)
is not straightforward to extract, e.g. from extended Press Schechter theory
(e.g., ).

Instead, we estimate  F(M) by considering the median halo growth:

    F_i(M) = 1 M ∑_t'>t_i 
    τ_ df > t'-t_i [ M(t'-Δ t_i) - M(t')  ]
     ,

where M(t') is the median mass expected for a halo at lookback
time t',  whose mass at
lookback time t_i is M,
Δ t is 
the previous timestep,[In equation (<ref>), Δ t is subtracted
  (not added) to t' because
t is a lookback time, which decreases with time.] 
and where the dynamical friction time is
τ_ df≡τ_ df[M(t'),M(t'-Δ t) - M(t'), t'].
The expression for F given in equation (<ref>) 
corresponds to summing up all the masses of secondary progenitors
merging with the primary ones.
It is an underestimate, because it
assumes that halo merging is binary. If, instead, many branches merge at
once, each branch will correspond to a lower mass progenitor and its
dynamical friction time would be longer. We thus assume that the lower mass
branches bring negligible mass to the final halo.

Figure <ref> shows that F(M) is of
order one-half, weakly increasing with halo mass, 
and with very little
dependence on the final considered lookback time.
The linear fit in log-log space yields

    log F(M) ≃ a + b log (M)
       ,

with a = -0.686 and b=0.0389 for z=0,
and a = -0.674 and b=0.0383 for lookback time of 1 Gyr.
This median surviving fraction F(M)
increases with halo mass, because large halo masses grow faster,
hence the values of M_2 are relatively more important (this offsets the
shorter dynamical friction times, since they are always greater than one-third
of the current Hubble time).




We can combine equations (<ref>) and (<ref>) to estimate the
maximum value of η_0^ VYG to obtain VYGs.
Figure <ref> shows the maximum values of η_0 to obtain a
VYG, for the , , and  models. We assumed η_1 = 1 and converted the final halo masses to stellar
masses using the SMHM predicted by the model.
Figure <ref> can be compared to the hatched shaded regions of 
Figure <ref>, corresponding to the VYG fractions obtained
with the quiet merging scheme applied to the Monte Carlo halo merger tree. For all choices of f_ sats, we recover
the sharp drop of VYG fractions at ≃ 10.7 for all 3 models.
We also recover most of the hierarchy between the 3 galaxy formation models:
 should produce a much higher fraction of VYGs than  at low mass,
 should overtake  at intermediate mass, and
the  model should dominate the others at intermediate masses.
The only weakness of our simple model is  that it fails to recover the higher
fraction of VYGs expected with the  model at low masses in comparison
to the  model. This might be explained by a higher value of f_
  sats for the  model, whose quiescent
evolution is frozen earlier than the other models by our prevention of
decreasing stellar masses, as seen in Figure <ref>.




 §.§ Evolutionary histories




Figure <ref> shows (grey curves) the median evolution of the main halo
progenitor in mass and circular velocity (top panels) and the median
evolution  of
the stellar mass, summed over all progenitors, for the 4 analytical models.
One notices that for intermediate halo masses (=13), the median
evolution of haloes since z=1 is at constant circular
velocity (see ),
while at low masses halo mass growth is so weak that halo mass is
nearly constant.

As seen in the lower 4 panels of Fig. <ref>,
the median stellar mass evolution (whose time derivative corresponds to the
star formation rate) can be growing or nearly frozen depending on the
analytical model and the final halo mass.

The bottom left panel of Figure <ref> displays the
evolution of 
stellar mass for 
representative very young (blue) and old (red) galaxies for the  model
with delayed
galaxy merging
(which produces no VYGs for massive haloes, ≥ 13).

The VYGs can be easily spotted by a sharp growth in the
stellar mass at low redshift, 
and are
caused by the rapid halo mass growth 
(top left panel).
The bottom left panel of Figure <ref> shows examples of VYGs
produced with unusual long-term histories (the example ending at =10.5) and
others with very typical histories (the examples ending at =8.7, 7.3,
and 5.9).
This illustrates the Markovian nature of the  model
(Sect. <ref>): 
the production of VYGs is not determined by the long-term histories of the
galaxies.



 §.§ Statistics of halo growth in last Gyr





One may ask whether, in the analytical models, the final rapid growth of
stellar mass (summed over all progenitors) is indeed linked with a corresponding rapid growth in
the most massive progenitor of a galaxy's halo.

Figure <ref> 
shows the statistics of halo mass growth in the last Gyr.
The abscissae (slightly shifted for clarity) above 3/4 indicate no mergers   with mass ratios
between 1:1 and 3:1, while those below 1/2 indicate more halo growth than a
single 1:1 major merger. 
In the  and  models, the rapid recent growth in total stellar mass
for the VYGs indeed comes
hand-in-hand with rapid mass growth of the most massive progenitor of the halo.
In contrast, the  and, especially,  models  show fairly frequent cases of rapid
final growth in stellar mass that occurs with slow halo mass growth.
To summarize, the median halo relative mass growth from 1 Gyr to the present is of order
0.015 for all galaxies, but for VYGs it is as high as 0.4 (), 0.6 (),
0.8 (), and 0.9 ().
This suggests the important role of major halo mergers. However, 
in nearly one-quarter of the VYGs in the  model, the haloes grow by less than 25 per
cent, indicating that more minor halo mergers suffice. 



Figure <ref> shows a very different halo merging history
for the VYGs obtained with quiet galaxy  merging: most of the VYGs are
associated with haloes that only slowly increase their mass in the last Gyr.
The distribution of the values of halo mass growth in the last Gyr are
indistinguishable between the VYGs and all galaxies.
This is the consequence of the stellar mass no longer being a function of
the current halo mass, but instead of the current halo mass minus the mass in
the  subhaloes that host
surviving satellites (eq. [<ref>]).



 §.§ The effect of the galaxy environment


Our analytical models completely miss the effects
of a galaxy's environment. Analyses of large observational samples of
galaxies by <cit.> and <cit.> indicate that the
fraction of star forming galaxies is reduced in both massive global environments
(high-mass clusters) and in the inner local environments (the inner regions
of groups and clusters).  Since most galaxies belong to
groups or clusters (e.g., ), one may expect that the lack of
environmental effects in our analytical models may cause substantial
modifications to the fraction of VYGs. 

The abundance matching analytical models mix the different environments to
yield a global SMHM relation. However, moving from high to low-mass groups, 
their central galaxies become increasingly likely to be of 
 spiral morphology
<cit.> and with efficient star formation <cit.>. 
One needs to go beyond standard abundance matching and also incorporate the
link between stellar ages and halo formation times (as indirectly observed by
, see Sect. <ref>, above), as this allows 
 reproducing correctly several environmental
trends <cit.>, in particular the fraction of star forming
central galaxies as a 
function of their luminosity (which is known to correlate strongly with
stellar and group mass).

Therefore, the quenching of star formation by the environment may explain why
the  SAM leads to typically 100  times lower fractions of VYGs at
intermediate mass (8 ≤≤ 10) in comparison with the analytical
models with bursty galaxy merging.
  However, one cannot unreservedly believe the results of the  SAM,
because of its limited mass resolution and also because the treatment of the
physics of quenching in SAMs is still very approximate.
Moreover, Figure <ref> indicates that the typical ages of galaxies
in the  SAM are of the same order as those found by , , and
 
models, but their spread is very narrow at intermediate stellar mass, leading
to much lower VYG fractions in this mass range.



 §.§ Effects of the primordial density fluctuation spectrum


Recently, there has been much interest in allowing for a warmer primordial
density fluctuation spectrum than allows the ΛCDM model.
Moving from  ΛCDM to  ΛWDM means that high
wavenumber primordial density fluctuations are suppressed. This decreases
(but does not suppress) the variance of low-mass fluctuations, making them
rarer at given epochs. This in turn means that haloes of given z=0 mass
collapse later. 
This effectively reduces the numbers of low-mass haloes at all epochs, 
potentially solving (; ) 
the problem of the overabundance of
such low mass haloes in cosmological N-body simulations
<cit.>. 
Since present-day stellar mass varies monotonically with halo
mass, one concludes that low-mass galaxy formation is delayed in the
ΛWDM model.
Therefore,
the fraction of VYGs should be enhanced in the ΛWDM
cosmology.



We have extracted the fractions of VYGs from the <cit.>
SAM (first presented in , see Sect. <ref>), 
run  on two forests of Monte-Carlo halo merger trees, one
for ΛCDM, the other for ΛWDM.
The Menci et al. SAM was run with the exact same parameters on both forests.
Figure <ref> shows that, for the CDM model, the fraction of VYGs among
intermediate mass galaxies is typically 50 times greater than predicted by
the  SAM run on the high-resolution MS-II simulation.
The VYG fraction  among
low-mass satellites is significantly greater in the  SAM (run on both
CDM and WDM halo merger trees) than among centrals in the same mass range,
which is opposite to the trend found in the  SAM (Fig. <ref>).

Figure <ref> also illustrates the effects of the primordial density
fluctuation spectrum on the fractions of VYGs.
We see that the power spectrum has little effect on the fraction of VYGs  at
intermediate masses, while the passage from CDM to WDM boosts the fraction of
VYGs among high-mass satellites. On the other hand,  moving from CDM to WDM leads to 10 to 30
times more VYGs at low masses (m < 10^9, except for one of the four
mass bins – at = 8.25), both for centrals and
satellites, although this increase is not statistically significant for the
centrals (except in the bin at =8.75).
This basically confirms the basic trend for younger galaxies with ΛWDM 
found by <cit.>
for this SAM. The boost in the fraction of VYGs among low-mass galaxies but
not higher mass ones is an illustration of the WDM power
spectrum that suppresses high wavenumbers. Star formation in low-mass
galaxies is thus delayed, 
making a higher fraction of low-mass galaxies
appear very young at z=0.
With the WDM particle mass of 0.75 keV, the boost in the VYG fraction occurs at masses m
< 10^8. Had we considered a higher WDM particle mass, the boost in
the VYG fraction would
have occurred at even lower galaxy mass.

These conclusions must be somewhat tempered by the fact that the Menci Monte Carlo
halo merger trees are not well resolved at low masses.
At stellar masses = 7±0.25 and 8±0.25,
the corresponding halo masses have a median of = 9.5 and 10,
respectively (see also the bottom panel of Fig. <ref>), so given the 
minimum branch mass of =7.7 (Table <ref>), the mass resolution of
the trees is only respectively 60 and 200, compared to 10^4 for our four analytical models. 

Note also that the <cit.> SAM, based upon Monte-Carlo halo merger trees,
includes a rough treatment of the environmental effects on galaxies.
It leads to intermediate fractions of VYGs at
intermediate mass (i.e. for their CDM run) in comparison with the high VYG fractions of the analytical
models and the low VYG fractions of the  SAM (compare
Figs. <ref> and <ref>). 




§ CONCLUSIONS






 §.§ Method and models

In this work, we modelled the frequency of VYGs (defined to be those
with over half the stellar
mass having formed within the last 1 Gyr, corresponding to z=0.08), among
central galaxies, as a
function of present-day stellar mass. For this, we first  produced
 over a quarter million Monte-Carlo  halo
merger trees derived from the code of <cit.>. We then
ran  four simple models of galaxy formation on these halo merger trees.
The first three models give the galaxy stellar mass as a function of
halo mass and redshift. These are the physically-motivated model (
with the refined reionization feedback of ), and
the empirical models of <cit.> and <cit.>
obtained by abundance matching of the halo mass function with the observed
SMF (with, for the latter, additional constraints on the cosmic star
formation rate and the specific star formation rate as a function of mass).
We also considered a fourth analytical model by <cit.>, where
the stellar mass
growth rate is proportional to halo mass growth rate and a simple function of
halo mass and redshift.
 
These analytical models
present differences in the  star formation
efficiency as a function of halo mass and redshift
(Figs. <ref> and <ref>), because
i) one () is physically motivated while the other 3 are empirical;
ii) only one
() is calibrated to data extending to the epoch of reionization;
iii)
only one () is based on mass growth rates instead of masses.
We also considered two SAMs, which incorporate much more complex physics: one
(Henriques) with the advantage of being based on realistic galaxy positions (based on
subhaloes from a dark matter cosmological simulation), the other (Menci) that
was run on both CDM and WDM cosmologies.
Our analytical models run on Monte-Carlo halo merger trees have
the advantage of having more statistics at the high halo mass end and higher
mass resolution.

We updated the stellar masses of our analytical models using
two
galaxy
merging
schemes,
one
involving a starburst at the time of the halo merger,
and one without it.
These two schemes should bracket the evolution of stellar masses.

The galaxy formation models applied to the halo merger trees produce
somewhat different halo to stellar mass relations (Fig. <ref>). Yet, they
generally match well (especially the  and  models) 
the analytical stellar to halo mass relations that they predict at the redshift where half
the mass in stars is formed.

These models predict stellar mass functions that generally agree with the
observations, but with noticeable differences in the low-end slopes.

The dependence on stellar mass of the epoch when half the final stellar
mass is formed (measuring the stellar mass evolution by considering all the
progenitors, not just the main one) shows differences between the galaxy
formation models ran with the bursty merging scheme.
All four analytical models and the SAM show downsizing at the high mass end: 
the median stellar ages decrease with decreasing stellar masses, reaching
a minimum around log (m/) ≃ 9.5 (see Fig. <ref>). 
The  model shows strong
 upsizing (stellar age increasing with decreasing stellar mass) at the
 low end, while the , , and  models show weak upsizing.
The typical median
ages at m = 10^8-10 range from 2.5 to 11 Gyr (Fig. <ref>), according to the
model. 
These differences between models can be explained by differences in their
star formation efficiencies as a function of halo mass and redshift
(Fig. <ref>) combined with our preventing stellar masses to
decrease in time.



 §.§ Frequency of very young galaxies


At z=0, the fractions of galaxies that are very young  depend on the galaxy stellar mass, as well
as on the model and the scheme for galaxy mergers.
For bursty galaxy merging,
the predicted  fraction of VYGs is roughly flat up
to 10^10 for the 
analytical models,  with  peaks at typically a few percent
at stellar mass near 10^9
(Fig. <ref>). 
For quiet galaxy merging, the VYG fractions are similar to those with
bursty merging for two models
( and ), but suppressed for two others, at respectively intermediate
() and low () stellar masses.

With bursty galaxy merging,
VYGs are always  ( and ) or often
( and ) associated with rapid late halo
growth
(Fig. <ref>),
indicating that most
VYGs
are associated with recent major halo mergers.
However, the , and
especially  models
also allow for stellar mass buildup with halo mergers of intermediate mass
ratios
(Fig. <ref>).
With quiet galaxy merging, we found, instead, no relation between halo mass
growth and galaxy stellar mass growth
(Fig. <ref>).

The VYG fractions
with bursty galaxy merging
are very well reproduced for the  and  galaxy
formation models 
(Fig. <ref>) 
by a simple, quasi-analytical, 
Markovian model (Sect. <ref>)
computing the growth of stellar mass in the last
Gyr by combining the growth during halo mergers with the quiescent growth.
This simple model reproduces less well the VYG fractions for the  model,
for which the past
history appears to play a role in the stellar mass buildup in the last Gyr,
hence the VYG fraction.
The trends of VYG fractions with stellar mass for the case of quiet galaxy merging are
qualitatively reproduced with another simple model
(Sect. <ref>),
which involves a lower
fraction of z=0 halo mass in subhaloes hosting surviving satellites than haloes 1
Gyr ago (Fig. <ref>).
The lower fractions of VYGs with the  and  models, regardless of the
merging scheme, thus appear to be
linked to the slow recent growth of their star formation efficiency (Fig. <ref>).

The semi-analytical model of galaxy formation
run by <cit.> on the MS-II produces 30 to 800 times lower
fractions of VYGs in comparison to the
analytical models run on Monte-Carlo halo merger trees with bursty galaxy
merging,
predicting less than 0.03 percent of
VYGs for 10^8 ≤ m ≤ 10^10 (lacking mass resolution at
the low end and statistics at the high end). These very low VYG fractions are
consistent with one () or another () analytical model with quiet
galaxy merging, but not both at once (Fig. <ref>).
These discrepancies may highlight the importance of more accurate
modeling such as in the SAM. Conversely, they may indicate that the  SAM
underestimates the effects of starbursts during galaxy mergers or
quenches too much galaxies within groups.

Finally, only at low masses (m < 10^8) is the fraction of VYGs 
significantly boosted with  Warm Dark Matter compared to Cold Dark Matter,
but the mass resolution of the tree code used for the SAM of <cit.> 
comparing CDM and
WDM may not be sufficient for accurate results.



 §.§ Final remarks



The  semi-analytical models used here may not have the necessary
mass resolution to probe the growth of stellar mass to form dwarf galaxies such as
I Zw 18,  while the physical analytical model is probably too simple and
missing some important astrophysical processes that may be present in the
semi-analytical model. The empirical analytical models should be more accurate,
but the 3 such models tested here differ in their predictions, because of
subtle differences in their stellar to halo mass relations, as confirmed with
the simple models we introduced in Sects. <ref> and <ref>.

We conclude that modelling the fraction of VYGs  is a promising sensitive
test of galaxy formation models,
given the wide range of predictions among the analytical models and with the state-of-the-art SAM. 
In Paper II, we confront these models with observations, 
measuring the fractions of VYGs in the local Universe  as a function of
stellar mass, using the spectral database of the SDSS.



§ ACKNOWLEDGMENTS

We thank the 
anonymous
 referee for many useful comments and additional
 references, which considerably improved our manuscript.
We acknowledge Hannah Parkinson for making her halo merger tree
code publicly available and to Shaun Cole for helpful comments on it.
We thank Steve Murray for making his HMFCalc code publicly available.
Thanks also to Gerry Williger for a critical reading.
DT and TXT are grateful to the Institut d'Astrophysique de Paris for
its hospitality, where a large part of
this work was performed, while GAM thanks the University of Virginia for
its hospitality.
This work was supported by the French Ministère de l'Enseignement
Supérieur et de la Recherche, the Ministère de Affaires Européennes
et Etrangères and the Israeli Ministry Of Science and Technology
within a France–Israel research program  for the project
Central
  Issues on Galaxy Formation, awarded to AD and GAM.
FC acknowledges funding from
the INAF PRIN-SKA 2017 program 1.05.01.88.04.
AD acknowledges support from a Lagrange Fellowship.
The Millennium Simulation databases used in this paper and the web
application providing online access to them were constructed as part of the
activities of the German Astrophysical Virtual Observatory (GAVO).







