# The evolution of compact massive quiescent and starforming galaxies derived
from the $R_{e}-R_{h}$ and $M_{\rm star}-M_{h}$ relations

L. Zanisi1,2 , F. Shankar1 , H. Fu1, A. Rodriguez-Puebla3, V. Avila-Reese3, A.
Faisst4, E. Daddi5, L. Boco6, A. Lapi6,7, M. Giavalisco8, P. Saracco9, F.
Buitrago10,11, M. Huertas-Company12,13,A. Puglisi14, A. Dekel15
1 Department of Physics and Astronomy, University of Southampton, Highfield,
SO17 1BJ, UK
2 DISCnet centre for doctoral training, University of Southampton, Highfield,
SO17 1BJ, UK
3 Instituto de Astronomía, Universidad Nacional Autónoma de México, A. P.
70-264, 04510, México, D.F., México
4 PAC, California Institute of Technology 1200 E California Boulevard,
Pasadena, CA 91125, USA
5 CEA, IRFU, DAp, AIM, Universitè Paris-Saclay, Universitè Paris Diderot,
Sorbonne Paris Citè, CNRS, F-91191 Gif-sur-Yvette,France
6 SISSA, Via Bonomea 265, I-34136 Trieste, Italy
7 IFPU - Institute for fundamental physics of the Universe, Via Beirut 2,
34014 Trieste, Italy
8 Department of Astronomy, University of Massachusetts Amherst, 710 N.
Pleasant St., Amherst, MA 01003, USA
9 INAF - Osservatorio Astronomico di Brera, via Brera 28, 20121 Milano, Italy
10 Departamento de Física Teórica, Atómica y Óptica, Universidad de
Valladolid, 47011 Valladolid, Spain
11 Instituto de Astrofísica e Ciências do Espaço, Universidade de Lisboa, OAL,
Tapada da Ajuda, PT1349-018 Lisbon, Portugal
12 Instituto de Astrofísica de Canarias (IAC); Departamento de Astrofísica,
Universidad de La Laguna (ULL), E-38200, La Laguna, Spain
13 LERMA, Observatoire de Paris, CNRS, PSL, Université Paris Diderot, France
14 Center for Extragalactic Astronomy, Durham University, South Road, Durham
DH1 3LE, United Kingdom
15 Racah Institute of Physics, The Hebrew University, Jerusalem 91904, Israel
E-mail<EMAIL_ADDRESS>
(Accepted for publication in MNRAS)

###### Abstract

The mean size ( effective radius $R_{e}$) of Massive Galaxies (MGs, $M_{\rm
star}>10^{11.2}M_{\odot}$) is observed to increase steadily with cosmic time.
It is still unclear whether this trend originates from the size growth of
individual galaxies (via, e.g., mergers and/or AGN feedback) or from the
inclusion of larger galaxies entering the selection at later epochs
(progenitor bias). We here build a data-driven, flexible theoretical framework
to probe the structural evolution of MGs. We assign galaxies to dark matter
haloes via stellar mass-halo mass (SMHM) relations with varying high-mass
slopes and scatters $\sigma_{\rm SMHM}$ in stellar mass at fixed halo mass,
and assign sizes to galaxies using an empirically-motivated, constant and
linear relationship between $R_{e}$ and the host dark matter halo radius
$R_{h}$. We find that: 1) the fast mean size growth of MGs is well reproduced
independently of the shape of the input SMHM relation; 2) the numbers of
compact MGs grow steadily until $z\gtrsim 2$ and fall off at lower redshifts,
suggesting a lesser role of progenitor bias at later epochs; 3) a time-
independent scatter $\sigma_{\rm SMHM}$ is consistent with a scenario in which
compact starforming MGs transition into quiescent MGs in a few $10^{8}$yr with
a negligible structural evolution during the compact phase, while a scatter
increasing at high redshift implies significant size growth during the
starforming phase. A robust measurement of the size function of MGs at high
redshift can set strong constraints on the scatter of the SMHM relation and,
by extension, on models of galaxy evolution.

###### keywords:

galaxies: abundances – galaxies: high redshift – galaxies: star formation –
galaxies: fundamental parameters

††pagerange: The evolution of compact massive quiescent and starforming
galaxies derived from the $R_{e}-R_{h}$ and $M_{\rm star}-M_{h}$
relations–References

##  1 Introduction

There is now substantial evidence that galaxies of a given stellar mass are
smaller at higher redshift than in the local Universe (e.g., Daddi et al.
2005, Trujillo et al. 2007,Buitrago et al. 2008, van Dokkum et al. 2010,
Cassata et al. 2011, Cimatti et al. 2012,Newman et al. 2012, Huertas-Company
et al. 2013, van Dokkum et al. 2015, Kawamata et al. 2015, Shibuya et al.
2015). The size evolution of the galaxy population in a given stellar mass bin
is well fitted by a relation of the type

$R_{e}(z)\propto(1+z)^{-\alpha}$ (1)

where $R_{e}$ is defined as the radius that encloses half of the galaxy light
(see, e.g., van der Wel et al. 2014 for a different fitting function). It is
found that in general starforming galaxies follow shallower trends (lower
values of $\alpha$) than quiescent galaxies (e.g., van der Wel et al. 2014).
Notably, the size growth rate of starforming galaxies increases with stellar
mass, becoming comparable to that of quiescent galaxies (with $\alpha\sim 1$)
for $M_{\rm star}>10^{11.2}M_{\odot}$ (Faisst et al. 2017, Mowla et al. 2018).
Moreover, in this mass regime, as pointed out in several studies (e.g.,
Bernardi et al. 2011a, b; Cappellari 2016), the behaviour of the scaling
relations and the stellar kinematics differ from that of less massive
galaxies. Thus, the mass scale $M_{\rm star}\approx 10^{11.2}M_{\odot}$ is
critical to understanding galaxy evolution. In this paper we focus on the
structural evolution of galaxies in this high mass regime, which we simply
label in what follows as “Massive Galaxies (MGs)”.

There is no consensus yet as to why MGs were a factor of 3 to 5 smaller in the
past. Minor dry mergers have been invoked as an efficient channel to promote
substantial size increase with relatively modest change in stellar mass (e.g.,
Naab et al. 2009,Oser et al. 2010,Shankar et al. 2013, van Dokkum et al. 2015)
to accomodate for the limited evolution in the high-mass end of the stellar
mass function (SMF) since $z\sim 1.5$ (e.g.,Andreon 2013,Muzzin et al. 2013,
McDermid et al. 2015, Kawinwanichakij et al. 2020). However, the rate of minor
dry mergers may not be sufficient by themselves to account for the entire size
evolution of MGs through cosmic time (Newman et al., 2012; Nipoti et al.,
2009, 2012). More generally, the exact contribution of dry mergers to the mass
assembly of massive galaxies is still a matter of intense debate among both
theoretical studies (e.g., De Lucia & Blaizot 2007; Hopkins et al. 2010b;
Wilman et al. 2013; Rodriguez-Gomez et al. 2015; Qu et al. 2017; Tacchella et
al. 2019; O’Leary et al. 2020; Grylls et al. 2020) and observational works
(e.g., Man et al. 2016; Mundy et al. 2017; Mantha et al. 2018; Duncan et al.
2019).

It has often been debated in the literature whether the size evolution of
galaxies of a given stellar mass stems from the size growth of individual
galaxies or it is a consequence of a “population effect” where newly formed,
larger galaxies enter the mass selection at later epochs thus increasing the
mean size distribution (e.g., Carollo et al. 2013, Shankar et al. 2015,
Gargiulo et al. 2017). This “progenitor bias" effect (van Dokkum & Franx,
1996) has been usually invoked to explain the size evolution of passive
galaxies with $M_{\rm star}<10^{11}M_{\odot}$ (e.g. Faisst et al. 2017;
Fagioli et al. 2016). Most studies agree on the lesser role of progenitor bias
in the size evolution of MGs at $z\lesssim 1$, in favour of a more predominant
role of (dry) mergers in increasing the sizes of individual galaxies
(e.g.,Saglia et al. 2010, Carollo et al. 2013, van der Wel et al. 2014,
Fagioli et al. 2016, Faisst et al. 2017, but see also Gargiulo et al. 2017).
In particular, the disappearance of compact (e.g., Cassata et al. 2011; Barro
et al. 2013) galaxies as the Universe ages is interpreted as a sign that they
must have grown in size individually (van der Wel et al., 2014) while a
constant abundance of compact galaxies implies that progenitor bias dominates
the size growth (Saracco et al., 2010; Gargiulo et al., 2016; Gargiulo et al.,
2017). In this respect, the full distribution of galaxy sizes at fixed stellar
mass, i.e., the size function $\phi(R_{e}|M_{star})$, is an invaluable tool to
disentangle galaxy evolution scenarios, providing simultaneous information on
the mean size $R_{e}$ and the number density of compact galaxies (e.g.,
Shankar et al. 2010; Carollo et al. 2013, Z20).

In addition to mergers and progenitor bias, Active Galactic Nuclei (AGN)
feedback during the compact starforming stages of the evolution of MGs (both
in a submm-FIR phase, e.g. Barro et al. 2016, and an optical “blue nugget”
phase, e.g. Martig et al. 2009; Damjanov et al. 2011; Barro et al. 2013; Fang
et al. 2013; Zolotov et al. 2015; Tacchella et al. 2016, which are potentially
linked in an evolutionary sequence, e.g. Gómez-Guijarro et al. 2019; Puglisi
et al. 2021) may also contribute to both size growth and quenching (Fan et
al., 2008; Fan et al., 2010; Kocevski et al., 2017; Lapi et al., 2018a; van
der Vlugt & Costa, 2019), and the relative evolution of compact starforming
and quiescent galaxies can provide tight constraints on these processes, as we
will further discuss below.

Due to their flexibility and relatively lower number of free parameters, semi-
empirical models have become a popular route to study the mass assembly, star
formation, and merger histories of galaxies (e.g., Conroy & Wechsler 2009;
Hopkins et al. 2010a, b; Behroozi et al. 2013; Moster et al. 2013; Shankar et
al. 2013; Gu et al. 2016; Matthee et al. 2017; Tinker 2017; Rodríguez-Puebla
et al. 2017; Lapi et al. 2018a; Chen et al. 2020; Grylls et al. 2019a, 2020;
O’Leary et al. 2021). The main ingredient in semi-empirical models is the
input stellar mass-halo mass (SMHM) relation, which is extracted from the
cumulative equivalence between the number densities of the stellar mass and
(sub)halo mass functions (e.g., Vale & Ostriker, 2006; Shankar et al., 2006;
Guo et al., 2011; Dutton et al., 2010; Firmani & Avila-Reese, 2010; Leauthaud
et al., 2012; Rodríguez-Puebla et al., 2013; Mandelbaum et al., 2016;
Pillepich et al., 2018; Behroozi et al., 2019; Moster et al., 2018; Kravtsov
et al., 2018; Erfanianfar et al., 2019). However, the systematic uncertainties
in the input data, most notably in the stellar mass function111Ultimately, the
reasons for these discrepancies are thought to originate from the way stellar
masses are estimated. The Initial Mass Function, dust attenuation curve,
stellar population synthesis models and assumed star formation histories, the
inclusion of intra-cluster light and even the photometry choice and background
subtraction algorithms all contribute to the various determinations of the
stellar mass functions found in the literature (see, e.g., Bernardi et al.
2010; Bernardi et al. 2013, 2016, 2017; Kravtsov et al. 2018; Leja et al.
2020; Lower et al. 2020; Guarnieri et al. 2019; Kawinwanichakij et al. 2020
amongst many others), which result in different estimates of the SMHM relation
(e.g., Shankar et al. 2017), have yielded discrepant results along the years.
For example, some groups (e.g., Moster et al., 2013; Behroozi et al., 2013;
Rodríguez-Puebla et al., 2017) proposed a shallower high-mass slope in the
SMHM relation, while others have argued in favour of a steeper slope (e.g.,
Shankar et al., 2014; Tinker et al., 2017; Kravtsov et al., 2018; Grylls et
al., 2019b). In turn, the high-mass end of the SMHM relation plays a crucial
role in, e.g., the number of galaxy pairs, and thus on galaxy merger rates
(Grylls et al., 2020; O’Leary et al., 2021).

Some works also highlighted a further correlation between the size of a
galaxy, $R_{e}$, and that of its host dark matter halo, $R_{h}$ (e.g. Fall
1983; Mo et al. 1998; Kravtsov 2013; Huang et al. 2017; Desmond 2017; Hearin
et al. 2019; Desmond et al. 2018; Somerville et al. 2018; Lapi et al. 2018b;
Mowla et al. 2019; Zanisi et al. 2020, but see also Desmond et al. 2017 ).
When coupled together in a semi-empirical model, the SMHM and $R_{e}-R_{h}$
relations become powerful tools to simultaneously probe the mass and
structural evolution of galaxies in a full cosmological context. For instance,
Somerville et al. 2018 found a non-trivial redshift and stellar mass
dependence of the mean $R_{e}-R_{h}$ relation at $M_{\rm star}\lesssim
10^{11.2}M_{\odot}.$ Zanisi et al. (2020) (hereafter  Z20), found that, in the
MGs regime and at $z\sim 0.1$, the $R_{e}-R_{h}$ relation must be very tight
with a total scatter (inclusive of observational statistical uncertainties) of
only $\approx 0.1$ dex222This was necessary to reproduce the full size
distribution of early- and late-type MGs in the local Universe as measured in
the Sloan digital Sky Survey (SDSS). Stringer et al. (2014) (hereafter S14)
combined the Moster et al. (2013) SMHM relation and a constant $R_{e}-R_{h}$
relation to build a semi-empirical model which proved capable to reproduce the
size evolution of MGs in the COSMOS field (Huertas-Company et al., 2013). S14
attributed the mean size growth of the population of MGs at $z\lesssim 2$ to a
cosmological effect for which (i) the size of the host dark matter haloes of a
given mass become larger as the Universe expands and its density decreases and
(ii) MGs of the same mass form in more massive, extended dark matter haloes at
lower redshift. While these results are encouraging, the effects of assuming a
different SMHM relation in the framework outlined by S14 has remained
relatively unexplored. The only notable exception is the work by Mowla et al.
(2019), where however only the shape of the SMHM was considered and not its
scatter, $\sigma_{\rm SMHM}$. The parameter $\sigma_{\rm SMHM}$ is notoriously
degenerate with the high-mass slope of the SMHM relation in retrieving the
number density of massive galaxies (Shankar et al., 2014; Wechsler & Tinker,
2018). As we will show, instead, $\sigma_{\rm SMHM}$ has a seizable effect on
the number density of the population of compact galaxies, and therefore it is
a significant novelty that we include in the semi-empirical framework set out
in S14.

In this paper we put forward a phenomenological, transparent methodology to
probe the size evolution of massive galaxies that expands on the works by S14,
Somerville et al. 2018 and Mowla et al. (2019). Following the approach of
Grylls et al. (2020), we build mock catalogs of galaxies in dark matter haloes
using data-driven toy models where SMHM relations with different shapes and
dispersion are coupled with linear $R_{e}-R_{h}$ relations characterized by
different normalizations. The main objective of this paper is to probe the
impact of varying the input SMHM relation and its dispersion $\sigma_{\rm
SMHM}$ on: 1) the mean size evolution of MGs, 2) the full size function of MGs
across cosmic time, and 3) the number density of compact MGs. The latter point
is particularly original and powerful as the time dependence of the number
density of compact MGs is closely linked to progenitor bias: less compact
galaxies at fixed stellar mass would be observed at later epochs if they grow
in size via, e.g., mergers. We will show, in particular, that the scatter in
the SMHM relation plays a major role in setting the number density of compact
MGs, allowing to break the degeneracy between the scatter and the high-mass
slope of the SMHM (see, e.g., Figure 1). Our present work lays out an
effective strategy to unveil the evolutionary pathways of MGs by exploiting
the increased statistics of MGs that will become available from future
observations. Data for MGs are in fact at present quite sparse and uncertain
at $z\gtrsim 1$ (e.g., Kawinwanichakij et al. 2020), and effective radii have
been measured for only a handful of MGs at $z\gtrsim 2$ (e.g., Kubo et al.
2017; Patel et al. 2017; Faisst et al. 2017; Mowla et al. 2018; Stockmann et
al. 2020; Lustig et al. 2021).

The outline of the paper is as follows. In Section 2 we provide the backbone
of our framework. In Section 3.1 we explore the role of the scatter in the
SMHM relation and we present toy models inspired to these findings in Section
3.2. In Sections 3.3, 3.4 and 3.5 we show how the different SMHM implied by
the toy models results in a range of possible determinations of the size
distributions of MGs. We give an interpretation of the evolutionary pathways
of MGs implied by the different toy models in Section 4.1, and we show how our
framework can be used to constrain the shape and scatter of the SMHM relation
in Section 4.2. In Sections 4.3 we discuss the possibility that the SMHM and
$R_{e}-R_{h}$ relations are correlated, and in Section 4.4 we discuss an
extension of the model of the $R_{e}-R_{h}$ connection based on halo
concentration, and discuss its limitations and strengths. Finally, we draw our
conclusions in Section 5. Further material is available in the Appendices.

Figure 1:  The effect of different SMHM relations on the halo occupation
distribution (HOD) of MGs and the size functions implied by a linear
$R_{e}-R_{h}$ relation (eq. 4). We show that SMHM relations with different
high-mass slope $\delta$ and scatter $\sigma_{\rm SMHM}$ (shown as shaded
areas in the left panel) can produce the same number density for MGs (the mass
threshold for MGs is shown as a red horizontal line). However, the halo
occupation distribution for the two models is remarkably different (central
panel). This translates in very different size functions. In particular, the
number density of compact galaxies differs by almost a factor of two (we use
the Cassata et al. 2013 definition of compactness, i.e. 0.4 dex below the z=0
mean size, against which we calibrated the two models following Appendix B ).
The model indicated with dashed gray lines is shown to help appreciate the
effect of a lower $\sigma_{SMHM}$, at fixed $\delta$, on the halo occupation
distribution and on the size function (compare to the dot-dashed model).

Although we show results only at a given redshift as an example, the same
arguments apply at any epoch for some choices of $\delta$ and $\sigma_{\rm
SMHM}$.

Figure 2:  _Left column:_ The SMHM relation of the four models outlined at the
beginning of Section 3.2. The red line indicates the stellar mass selection
for MGs. The green shaded regions indicate the scatter of the SMHM, which
increases at higher redshift for models 2 and 4. The gray line in the center-
bottom and bottom panels indicates the SMHM for Models 1 and 2 as a reference.
_Central and Right column:_ The redshift evolution of the halo occupation
distribution $\phi(M_{h}|M_{\rm star}>10^{11.2}M_{\odot})$ and the implied
size functions $\phi(R_{e}|M_{\rm star}>10^{11.2}M_{\odot})$ of MGs for the
four models. We display results for z=0,0.5,1,1.5,2,2.5,3. Darker colours
indicate higher redshift. The gray band in the right column shows the Cassata
et al. (2013) definition for compact galaxies. It can be seen that the
increasing $\sigma_{\rm SMHM}$ of Models 2 and 4 results in broader
distributions, which have a median lower $M_{h}$ and normalised $R_{e}$
compared to Models 1 and 3, where $\sigma_{\rm SMHM}=0.15$ dex at all times.
An evolving $\sigma_{\rm SMHM}$ also results in a higher number density of MGs
at earlier times. Contrariwise, the flatter high-mass-end slope of the SMHM in
Model 3 results in overall fewer MGs and slightly larger median halo masses
compared to Model 1.

##  2 Methods

###  2.1 The galaxy-halo connection

To study the size distribution and evolution of MGs ($M_{\rm
star}>10^{11.2}M_{\odot}$), at any redshift of interest we build a catalogue
of dark matter haloes with mass $M_{h}$ and size $R_{h}$, to which we assign a
stellar mass $M_{star}$ (via an input SMHM relation) and an effective radius
$R_{e}$ (via a $R_{e}$-$R_{h}$ relation). The detailed modelling approach that
we use here closely follows the one outlined in Zanisi et al. (2020), which we
briefly summarise below for convenience333We make extensive use of the
COLOSSUS Python package (Diemer, 2017).:

  1. 1.

We use the Despali et al. (2016) halo mass function to obtain large catalogues
of dark matter haloes at all the redshifts of interest. Note that the Despali
et al. (2016) halo mass function is defined for central galaxies only, as it
does not include subhaloes (e.g. van den Bosch et al. 2014). Thus, all MGs in
this study are modelled as central galaxies.444In our model all MGs are
considered central galaxies, since satellites are negligible in this extreme
mass range (e.g., Peng et al. 2010). Using the Statistical SemiEmpirical Model
STEEL (Grylls et al., 2019b), we find that the satellite contribution to MGs
is less than 20% at z$\sim 0.1$, and it declines steeply at earlier times. The
subdominant population of satellite MGs is studied in a companion paper
(Zanisi et al. 2020).. In this paper we model dark matter haloes with a
Navarro et al. (1996) density profile with scale radius $R_{s}$,

$\rho(r)\propto\frac{1}{\frac{r}{R_{s}}\bigl{[}1+\frac{r}{R_{s}}\bigr{]}^{2}},$
(2)

and with $R_{h}=cR_{s}$ defining the _concentration_ parameter $c$. $R_{h}$ is
the dark matter halo radius,

$R_{h}=\Bigl{(}\frac{3M_{h}}{4\pi\cdot\Delta\rho_{\Delta}}\Bigr{)}^{\frac{1}{3}}$
(3)

where $\Delta$ is the virial overdensity with respect to the cosmological
critical density (Bryan & Norman, 1998). Both $\Delta$ and $\rho_{\Delta}$ are
decreasing functions of cosmic time (e.g, Mo et al. 2010). Thus, a dark matter
halo of a given mass has a smaller size at higher redshift, owing to the
higher density of the Universe.

  2. 2.

We model the link between galaxies and dark matter via the stellar mass-halo
mass relation (SMHM). The mean SMHM$\equiv M_{\rm star}(M_{\rm h})$ is a
monotonically increasing function of halo mass. We include a lognormal scatter
$\sigma_{\rm SMHM}$ at fixed halo mass that takes into account both the
intrinsic dispersion $\sigma_{\rm intr}$ in the relation and the uncertainty
in stellar mass estimates $\sigma_{*}$ (i.e., $\sigma_{\rm
SMHM}^{2}=\sigma_{\rm intr}^{2}+\sigma_{*}^{2}$, see Rodríguez-Puebla et al.
2017, Behroozi et al. 2013, Shankar et al. 2014, Tinker et al. 2017). In the
next Section we will present “toy” SMHM relations, which vary in both shape
and dispersion, to probe their impact on our galaxy mocks and on their size
distributions at different epochs. In particular, we will focus on the slope
above the knee of the SMHM relation, $\delta$ (see Figure 1), which is the
parameter in the SMHM relation controlling the number density of MGs at a
given dispersion. The precise value of $\delta$, or better of the underlying
abundances of MGs in the local and high redshift Universe, still suffer from
substantial systematic uncertainties and that are hotly debated in the
literature (see Section 1).

  3. 3.

We assign a half light radius $R_{e}$ to each galaxy according to the ansatz:

$R_{e}=A_{K}R_{h},$ (4)

which is based on the empirical findings by Kravtsov (2013), and that we call
the _K13 model_. Here $A_{K}$ is the normalization which in principle may vary
with halo mass, galaxy stellar mass and/or star formation activity (e.g. Huang
et al. 2017; Somerville et al. 2018. Z20, B). We add to eq. 4 an intrinsic
log-normal scatter $\sigma_{K}$, which, as $A_{K}$, is a free parameter. We
stress that, effectively, $R_{e}$ is a function of $M_{h}$, since there is a
direct proportionality between halo mass and halo radius (see eq. 3). While in
the remainder of the paper we will mostly comment on the K13 model, in Section
4.4 we discuss another model of galaxy sizes in which the relation between
virial radius and galaxy size is also mediated by the halo concentration
(e.g., Desmond 2017; Jiang et al. 2019; Zanisi et al. 2020). Moreover, other
definitions of galaxy sizes, such as $R_{80}$ (Miller et al., 2019) or $R_{1}$
(Trujillo et al., 2020), have been proposed to correlate to $R_{h}$ equally
well or even better than effective radius. We will discuss these models in
Appendix E.

###  2.2 Quenching

To provide a fair comparison to observations, which have so far always
distinguished between starforming and quiescent MGs (e.g., Mowla et al.,
2018), we need to include a recipe for quiescence in our galaxy mocks. To this
purpose, following the empirical calibration of Rodríguez-Puebla et al. (2015)
at $z\sim 0.1$, we assume that the probability of a galaxy being quenched in a
dark matter halo of mass $M_{h}$ is given by the fraction

$f_{Quench}(M_{h})=\frac{1}{b_{0}+[\mathcal{M}_{0}\times
10^{12}/M_{h}(M_{\odot})]}$ (5)

with $b_{0}\sim 1$ and $\mathcal{M}_{0}\sim 0.68$ at $z\sim 0.1$. $f_{Quench}$
is a monotonically increasing function of halo mass, with a characteristic
mass scale $\mathcal{M}_{0}$ above (below) which more (less) than 50% of
galaxies are quiescent (starforming).

The fraction of quenched MGs is observed to evolve with redshift (e.g.
Huertas-Company et al. 2016; Mowla et al. 2018). While it is beyond the scope
of this work to set specific constraints on the physical processes that drive
quenching (see Somerville & Davé 2015 for a review), we note that quenching is
thought to be more likely to occur in more massive haloes at higher redshift
(e.g., see the empirical models by Rodríguez-Puebla et al. (2017) and Behroozi
et al. (2019)). In our model this is achieved by replacing $\mathcal{M}_{0}$
with

$\mathcal{M}(z)=\mathcal{M}_{0}+(1+z)^{\mu},$ (6)

where $\mu>0$ is a free parameter, that regulates the increase in
characteristic quenching halo mass in the younger Universe. Figure 8 shows
examples of the evolution in $f_{Quench}$ for $\mu=1,3,5$. We note that
quiescence is defined in the literature according to different methods (e.g. 1
$\sigma$ below the main sequence, different cuts in the color-color planes, a
hard cut in specific star formation rate, see e.g., Donnari et al. 2019) which
can lead to different results (Sherman et al., 2020). Therefore, the value of
$\mu$ will depend on the method assumed. For this reason, in the following we
simply show different values of $\mu$, which we will adapt to the specific
method used once the comparison data is fixed.

Quiescent and starforming MGs of similar mass appear to grow in size at the
same rate with redshift, with quiescent galaxies being systematically smaller
at all times (Mowla et al. 2018). We assume that the two populations live on
two separate K13 relations. The normalizations $A_{K,SF}$ and $A_{K,Q}$, for
starforming and quenched MGs respectively, are calibrated at $z\sim 0.1$
following Appendix B. Following Z20 we also assume that the scatters in the
two K13 relations, $\sigma_{K,SF}$ and $\sigma_{K,Q}$, are equal to 0.1 dex.
In the remainder of this paper, we assume that this value of $\sigma_{K}$
holds at all times.

###  2.3 Target observables

Using the methodology outlined above, we will present the results of some toy
models (described in Section 3.2) for the following observables:

  * •

the evolution of the galaxy size distribution of MGs (i.e. the size function
$\phi(R_{e},z|M_{\rm star}>10^{11.2}M_{\odot})$) and its integral, the number
density of MGs

$n_{MGs}(z)=\int_{-\infty}^{\infty}\phi(R_{e},z|M_{\rm
star}>10^{11.2}M_{\odot})dlogR_{e};$ (7)

  * •

the mean size of the population of MGs as a function of redshift, $\langle
R_{e}(z)\rangle$;

  * •

the evolution in the number density of compact MGs $n_{compact}(z)$. A range
of definitions of compactness have been proposed in the literature (e.g.,
Saracco et al. 2010; Fang et al. 2013; Carollo et al. 2013; Barro et al. 2013;
van der Wel et al. 2014; Damjanov et al. 2015; van Dokkum et al. 2015; Barro
et al. 2017; Charbonnier et al. 2017; Tacchella et al. 2017; Buitrago et al.
2018; Tortora et al. 2018; Luo et al. 2020 amongst many others). Here we
define galaxies as compact systems if their size is 0.4 dex below the $z\sim
0$ $R_{e}-M_{\rm star}$ relation of quenched galaxies (Cassata et al., 2011,
2013),

$n_{compact}(z)=\int_{-\infty}^{-0.4}\phi(R_{e}/R_{e}(z=0),z)dlog(R_{e}/R_{e}(z=0)).$
(8)

In particular, we focus on compact quenched MGs (CQMGs) and compact
starforming MGs (CSFMGs).

Figure 9 shows that adopting other definitions of compactness based on the
effective radius yields qualitatively similar results to the Cassata et al.
(2011) definition. Other popular definitions of compactness based on, e.g.,
the stellar mass density in the central kiloparsec, would require information
on the light/mass profile of galaxies (e.g., the Sérsic index), which we are
not including here. This requires further modelling which we defer to future
work.

##  3 Results

###  3.1 At the core of the model

The methodology outlined in Section 2 makes use of only two ingredients: (i)
the K13 relation (eq. 4) and (ii) the SMHM relation (most notably the high-
mass slope $\delta$ and the scatter $\sigma_{\rm SMHM}$).

Figure 1 shows that two SMHM relation with different high-mass slope $\delta$
and scatter $\sigma_{\rm SMHM}$ are able to produce the same number density
for MGs. The degeneracy between $\delta$ and $\sigma_{\rm SMHM}$ in producing
the same abundances of massive galaxies was already identified in previous
studies (e.g., Behroozi et al. 2010; Shankar et al. 2014). What we emphasize
here, for the first time to the best of our knowledge, is that the
corresponding halo mass distributions (middle panel, see also Shankar et al.
2014), and thus the implied size functions computed via the linear
$R_{e}-R_{h}$ relation (right panel), remain however significantly distinct,
especially below the peaks of the distributions. The larger abundances of
compact MGs is mostly driven by a larger scatter in the input SMHM relation,
as can be inferred by comparing black dot-dashed and gray dashed lines in
Figure 1. Thus, the abundance of compact galaxies represents a valuable
observable to break the degeneracy between $\delta$ and $\sigma_{\rm SMHM}$,
allowing to set constraints on the degree of progenitor bias and ultimately to
discriminate between different models of galaxy formation.

###  3.2 Toy models

Motivated by the discussion above, we devise four toy models to show the
effect of varying $\delta$ and $\sigma_{\rm SMHM}$ on our target observables
(Section 2.3):

  * •

Model 1: $\delta=0.5$ (steep slope), $\sigma_{\rm SMHM}=0.15$ dex at all
redshifts;

  * •

Model 2: $\delta=0.5$ (steep slope), $\sigma_{\rm
SMHM}=\sqrt{(0.1z)^{2}+0.15^{2}}$;

  * •

Model 3: $\delta=0.35$ (shallow slope), $\sigma_{\rm SMHM}=0.15$ dex at all
redshifts;

  * •

Model 4: $\delta=0.35$ (shallow slope), $\sigma_{\rm
SMHM}=\sqrt{(0.1z)^{2}+0.15^{2}}$.

The slope of Model 1 (Model 2) is inspired to the Grylls et al. (2019a)
’PyMorph’ (’cmodel’) SMHM relation, which was obtained by fitting the Bernardi
et al. (2017) ’PyMorph’ (’cmodel’) stellar mass function (SMF) at $z\sim 0.1$
and the Davidzon et al. (2017) SMFs at $z\gtrsim 0.2$ where their masses have
been corrected by 0.15 dex to bring the two studies in agreement555This was
done only for the ’PyMorph’ SMF. The factor of 0.15 dex takes into account the
difference in $M/L$ used in the two studies. (see also Bernardi et al.
2016).666We also shift by -0.1 dex the knee of the SMF resulting from the
Grylls et al. (2019a) SMHM to better match the SDSS SMF.

Although some authors point to distinct SMHM relations for quiescent and
starforming galaxies (e.g., Rodríguez-Puebla et al. 2015; Moster et al. 2018;
Behroozi et al. 2019; Posti & Fall 2021), the relative content of stars in
starforming and quiescent galaxies at fixed halo mass is still highly debated
(e.g., Wechsler & Tinker, 2018). We here adopt throughout the simplest
assumption that quiescent and starforming galaxies share the same underlying
SMHM relation, and note that of our core results do not qualitatively depend
on this working assumption.

###  3.3 Halo occupation and implied size function

As a first step, in Figure 2 we show the SMHM relation and its scatter for the
four toy models, as well as the distribution of the host halos (i.e. the halo
occupation distribution) and the implied size functions. Figure 2 reveals that
different SMHM relations and their scatter $\sigma_{\rm SMHM}$ provide
significantly different size functions, that necessarily stem from distinct
host halo occupation distributions. Thus, the size functions are completely
regulated by the way the SMHM relation maps galaxies into haloes. In
particular, it is relevant to highlight the following features when comparing
different models for the input SMHM relation:

  * •

Model 1 vs Model 3. A change in the high-mass slope of the SMHM relation
generates an overall lower number density of MGs, but the mean of the halo
occupation distributions and related size functions are fairly similar in the
two cases.

  * •

Model 1 vs Model 2 and Model 3 vs Model 4. Even when the shape of the SMHM
relation is identical, if we allow for the scatter $\sigma_{\rm SMHM}$ to
evolve with redshift, and in particular to increase at earlier epochs, then
the implied halo occupation distribution drastically changes compared to the
case with constant scatter. In the former case, a higher proportion of small
MGs are hosted in less massive haloes at higher redshift, and the mean halo
occupation and galaxy size exhibit a stronger evolution, as quantitatively
described below.

###  3.4 Implied size evolution

S14 showed that, on the assumption that $R_{e}\propto R_{h}$ at all epochs,
the progressive increase in virial radii and in the number densities of
massive dark matter haloes, were sufficient conditions to produce, when
averaging over the full population, a strong size evolution in the sizes of
massive galaxies.

Figure 3 confirms and further extends the claim by S14. By using, for each of
our four toy models, a constant proportionality $R_{e}=A_{K}\times R_{h}$
calibrated at $z=0.1$ (see Appendix B), as labelled, we are always able to
reproduce the strong redshift evolution seen in the available data (Faisst et
al., 2017; Patel et al., 2017; Mowla et al., 2018), irrespective of the exact
input SMHM relation. Models with an evolving $\sigma_{\rm SMHM}$ tend to
predict up to less than 50% faster size evolutions, well within the variance
currently found in the data. We distinguish between starforming and quiescent
galaxies via the $f_{Quench}$ model with $\mu=2$. Varying the $\mu$ parameter
has little effect on our results, as it can be easily compensated by a
relative variation in $A_{K}$ and/or in the SMHM relation. Indeed, the $A_{K}$
retrieved for starforming and quenched MGs appear to be systematically
different and such difference persists even when adopting distinct SMHM
relations as, for example, in Moster et al. (2018), for which we find
$A_{K,SF}\approx 0.023$ and $A_{K,Q}\approx 0.016$.

Figure 3: The size evolution from the K13 model for starforming (left) and
quiescent (right) MGs. The black lines indicate the four toy models outlined
at the beginning of Section 3.2 and coupled with the $f_{Quench}$ model
introduced in Section 2.2 with $\mu=3$, as an example. Data are the sizes of
MGs from Mowla et al. (2018) (diamonds), Faisst et al. (2017) (circles), Patel
et al. (2017) (triangles). We also add SDSS estimates for the sizes of Massive
Late type and Early type galaxies (green and orange stars respectively) from
Zanisi et al., 2020. The normalization $A_{K}$ in each panel is chosen to
match SDSS observations. Notably, a constant normalization $A_{K}$ is able to
reproduce observations. Moreover, $A_{K}$ is lower for shallower high-mass-end
slopes of the SMHM (Models 1 and 2), while the opposite is true for steeper
SMHM relations. This indicates that $A_{K}$ and $\delta$ are degenerate in our
model.

###  3.5 Implied statistics of compact MGs

Figure 4: _Left_ : The number density of starforming and quenched MGs (cyan
and red line respectively) for models 1 and 2. Solid, dashed and dotted lines
are obtained adopting $\mu=1,2,4$ respectively. _Right:_ Prediction for the
number density of compact MGs for the two models and the different values of
$\mu$. Note that the fraction of compact MGs increases at early times. Filled
diamonds and crosses indicate the time where 20% and 70% of the population of
MGs (either starforming or quenched) are compact. The comparison data in the
left column are from the SDSS ’PyMorph’ photometry at z=0.1 (Meert et al.,
2015), Davidzon et al. (2017) and McLeod et al. (2020) (in the wo latter cases
the masses have been shifted by 0.15 dex to account for M/L differences with
our SDSS data, see Grylls et al. 2019a). Note that the data points were
retrieved from the Schechter fits provided in the two studies, extrapolated in
our mass range. With the caveat that different definitions of quiescence are
adopted in observations, we note that Model 1 is favoured by current data if
$\mu\approx 2-3$. Model 2 might provide a better fit to data if the number
density of starforming MGs is underestimated at high redshift (Franco et al.,
2018; Smail et al., 2021).

In Figure 2 we showed that the shape and scatter of the SMHM have a
significant impact on the number density of compact galaxies, a feature that
was not investigated by previous studies. We explore these trends more
quantitatively here for our toy models. The top and bottom panels of Figure 4
show the predictions of Model 1 (constant scatter) and Model 2 (evolving
scatter) for the number density of MGs (left panels) and for only compact MGs
(right panels), separately for quiescent (red) and starforming (cyan) galaxies
and for different values of the quenching parameter $\mu$, as labelled (the
predictions for Models 2 and 3 are very similar and reported in Appendix D).
All models predict a similarly sharp rise in the number density of compact
quiescent MGs (red lines) up to $z\sim 1.5-2$ and a subsequent more or less
fast drop depending on the exact value of $\mu$ adopted. All models also
predict the abundances of starforming compact MGs (cyan lines) to peak around
the same redshift $z\sim 2.5$ with a weak dependence on $\mu$ but a strong one
on scatter: a larger $\sigma_{\rm SMHM}$ at early epochs can increase by up to
a factor of ten the predicted number densities of starforming compact MGs
(bottom right panel). In Appendix C we show that adopting other definitions of
compactness (e.g., Barro et al., 2013; van der Wel et al., 2014; Gargiulo et
al., 2017) does not alter the main qualitative trends of Figure 4.

The evolution of $n_{compact}$ that we predict for compact quiescent MGs is in
qualitative agreement with observations of compact galaxies in a lower mass
range ($10.5<logM_{star}/M_{\odot}<11.5$, Cassata et al. 2011, 2013; van der
Wel et al. 2014; Barro et al. 2013). However, at present, current observations
provide rather uncertain constraints on $n_{MGs}$ at high redshift (see
Kawinwanichakij et al. 2020 for a detailed discussion of the systematics). In
addition, a secure determination of the number density of, especially compact,
MGs is hampered by the seizable but still unknown number of optically dark
starforming galaxies at high redshift (e.g., Franco et al. 2018; Wang et al.
2019; Zhou et al. 2020; Smail et al. 2021). Nevertheless, the results
presented in Figure 4 provide clear predictive trends for the evolution of
compact and large MGs that, when compared with data from the next generation
of observing facilities, will set tight constraints on the quenching
mechanisms ($\mu$ parameter) and on the level of progenitor bias in the size
evolution of MGs.

##  4 Discussion

###  4.1 Progenitor bias scenarios and continuity equation

Figure 5: The number density of CSFMGs (cyan) and CQMGs (red) for Model 1 (top
row) and Model 2 (bottom row). We adopt mu $\mu=2$ (dotted lines, left
columns) and $\mu=3$ (solid lines, right column). The number density of
compact starforming galaxies that would be obtained from continuity arguments
(see eq. 9) is shown for different values of the quenching timescale $\Delta
T_{quench}$. Increasingly larger values of $\Delta T_{quench}$ are shown with
increasing brightness. Model 2 disfavours a continuity scenario. In Model 1
continuity is broadly achieved if $\Delta T_{quench}\approx 200-400$ Myr for
$\mu=2$ and $\mu=3$ respectively. Results for Models 3 and 4 can also be found
in the online supporting material, and are qualitatively similar.

We have demonstrated that all our models are able to produce a strong
evolution in the average effective radius of the MG populations (Figure 3). On
the other hand, Figures 4 and 10 clearly show that in all models $n_{compact}$
decreases below $z\sim 1.5-2$. The peak of the abundance of compact quenched
MGs corresponds to compact fractions of $\sim 20-40$%. Thus, the ensuing
disappearance of compact galaxies as the Universe ages strongly suggests that
$\sim 20-40$% of the quenched MGs that were present at $z\sim 1.5-2$ have
grown in size _individually_ (e.g., Trujillo et al., 2011; Carollo et al.,
2013; van der Wel et al., 2014; Fagioli et al., 2016; Faisst et al., 2017;
Stockmann et al., 2020). However, it is worth pointing out that this
corresponds to only $\sim 10-15$% of the quenched MGs that are present today
(for the case of constant and evolving $\sigma_{\rm SMHM}$ respectively).

At $z\gtrsim 2$ all models instead predict a strong increase in the number
density of compact MGs, suggesting that, in line with a number of
observational studies (e.g., Barro et al., 2013; Cassata et al., 2013), a
significant fraction of MGs form in a compact phase at early epochs, most
probably due to gas dissipation following a merger (e.g., Sparre & Springel,
2016) or an in-situ burst of star formation (e.g., Lapi et al., 2011).

An interesting question that has been discussed in the literature is whether
compact quenched galaxies are the descendants of compact starforming galaxies
(e.g., van Dokkum et al., 2015; Barro et al., 2017; Gómez-Guijarro et al.,
2019). For example, based on number density conservation arguments, Barro et
al. (2013) proposed that compact starforming galaxies with
$10.5<logM_{star}/M_{\odot}<11.5$, passively evolve into quenched compact
galaxies in a timescale of $\sim 800$Myr. Following Barro et al. (2013), we
here develop basic continuity equation models without mergers in which compact
starforming MGs (CSFMGs) naturally evolve into compact quenched MGs (CQMGs) on
a given timescale $\Delta T_{quench}$ as

$n_{CSFMGs}(t)=n_{CQMGs}(t+\Delta T_{quench})-n_{CQMGs}(t)$ (9)

in which $\Delta T_{quench}$ is allowed to vary between 100 and 900 Myr, $t$
is the age of the Universe, and $n_{CQMGs}(t)$ and $n_{CSFMGs}(t)$ are the
cumulative number densities of quiescent and starforming compact MGs above
$M_{star}>10^{11.2}\,M_{\odot}$. Figure 5 shows the results of applying Eq 9
to the $n_{CSFMGs}$ extracted from Models 1 and 2 (see online supporting
material for Models 3 and 4) with quenching parameters $\mu=2,3$ as a
reference (the results derived for other values of $\mu$ are included in the
online supporting material.)

As reported in the bottom panel of Figure 5, models characterised by a scatter
$\sigma_{\rm SMHM}$ increasing at earlier epochs tend to disfavour a
continuity scenario in which all CSFMGs gradually transition into CQMGs, as
the number densities of CSFMGs (cyan lines) are always significantly larger
than those of compact quiescent galaxies (red lines). Instead, models with a
fixed $\sigma_{\rm SMHM}$ (Model 1, top row of Figure 5) are broadly
consistent with a progenitor-descendant scenario between CSFMGs and CQMGs for
some choices of $\Delta T_{quench}$. In the specific, we find that $\Delta
T_{quench}\approx 200,300,400,900$ Myr for $\mu=2,2.5,3,4$ (data shown only
for $\mu=2,3$, see online supporting material for $\mu=2.5,4$). Thomas et al.
(2005) estimated an upper limit to the main star formation episode of local
MGs around $\Delta T_{quench}\lesssim 300$Myr (see their eq. 5), which would
be consistent, at face value, with continuity in our constant $\sigma_{\rm
SMHM}$ models with $2\lesssim\mu\lesssim 3$, in line with our preferred values
of $\mu$ adopted in Figures 4 and 10. We note that continuity arguments
applied to Models 3 and 4 (see online supporting material) yield results that
are qualitatively similar to Models 1 and 2 respectively.

In a continuity scenario between CSFMGs and CQMGs (which can be produced by
Models 1 and 3), little or no size evolution occurs during quenching. This
conflicts with theoretical models where both size evolution and quenching
occur almost simultaneously as a result of AGN activity, with a predicted
expansion in size of a factor of $\gtrsim 2$ over very short timescales (i.e.,
$50-100$ Myr, Ragone-Figueroa & Granato 2011; Lapi et al. 2018a). In other
words, assuming a constant scatter $\sigma_{\rm SMHM}$ in the input SMHM
relation, would be consistent with a two-stage formation scenario in which
galaxies first quench and then grow via stochastic mergers (e.g., Hopkins et
al., 2009; Oser et al., 2012). Alternatively, an increasing $\sigma_{\rm
SMHM}$ at earlier epochs would necessarily require within our framework that
only a relatively minor fraction of the CSFMGs quench during their compact
phase, a scenario more consistent with an AGN-driven size evolution. We note
that an unbiased view of the size growth of MGs requires both optical-NIR
observations as well as FIR-submm observations (e.g., Barro et al., 2016;
Tadaki et al., 2020; Sun et al., 2021). Compact dust-enshrouded star formation
activity can in fact occur over spatial scales a factor of $\sim$3 smaller
that the $R_{e}$ measured in optical-NIR (e.g., Puglisi et al., 2019; Jiménez-
Andrade et al., 2019). AGN activity in these galaxies might cause, along with
quenching, a considerable evolution in size in a very short timescale (e.g.,
Lapi et al., 2018a).

We conclude this Section by stressing the fact that our continuity models
strictly apply to _compact_ MGs, which amount to a substantial fraction of the
total population of quiescent MGs only at $z\gtrsim 2$ (see diamonds and
crosses on Figure 4). It is evident from Figure 4, that all models predict an
increase in the number density of the overall quiescent population at $z<2$ by
up to an order of magnitude, a trend that cannot be driven by solely quenching
of the starforming MGs as the number density of the latter is always
significantly lower that those of quenched MGs at late epochs. Additional
physical processes must be at play at $z<2$ in regulating the formation and
sustenance of non-compact starforming MGs as well as the appearance of a large
population of non-compact quenched MGs.

###  4.2 The sizes of MGs as effective constraints to the galaxy-halo
connection

Providing firm constraints to the SMHM relation at different epochs can yield
invaluable information on, e.g., the merger rates of MGs (Grylls et al.,
2020), the interplay between dark matter and baryonic physics (Gu et al.,
2016; Matthee et al., 2017), the physical processes behind galaxy quenching
(Tinker, 2017). Unfortunately, the shape and scatter of the SMHM relation are
still highly debated (e.g.,Bernardi et al. 2017). In particular, there is a
well-known degeneracy between the high-mass slope, $\delta$, and the
dispersion, $\sigma_{\rm SMHM}$, of the SMHM relation (e.g. Shankar et al.
2014). Similarly to Grylls et al. (2020), in the previous Sections we made use
of toy models where only these two parameters are changed to explore their
impact on the sizes of MGs. As shown above, SMHM relations with different
values of $\delta$ and $\sigma_{\rm SMHM}$ result in distinct rates of size
increase (Figures 3) and number density evolution of compact MGs (Figures 4
and Figure 10), which are ultimately a consequence of the different implied
halo occupation distribution (Figure 1). Our results therefore suggest that
the $\delta-\sigma_{\rm SMHM}$ degeneracy may be broken by simultaneously
fitting the size growth of MGs, the redshift evolution of the number density
of compact MGs, and the number density evolution of the overall population of
MGs, in other words by an accurate measurement of the full size function of
MGs at different epochs, a goal that should be achievable with the aid of the
next-generation observational facilities such as Euclid and LSST. We note that
several previous semi-empirical studies aimed at probing the size evolution of
galaxies (e.g., Rodríguez-Puebla et al. 2017; Hearin et al. 2019; Behroozi et
al. 2021). However, they were are all limited by the use of only one SMHM
relation and dispersion, which instead, if allowed to vary, can provide
distinct structural evolutionary tracks for MGs.

It is important to highlight that the systematic uncertainties in measuring
stellar masses and number densities of galaxies substantially affect the
determination of the stellar mass function, and thus of the SMHM relation and
size distributions of MGs at different epochs, possibly explaining at least
part of the observational discrepancies in the numbers of compact galaxies
reported in the literature (e.g., Poggianti et al. 2013).

###  4.3 Covariance between the $R_{e}-R_{h}$ and the SMHM relations?

A further interesting issue that warrants more exploration is that of a
possible covariance between the $R_{e}-R_{h}$ and the SMHM relation. In our
framweork, the $R_{e}-R_{h}$ and the SMHM relations are closely intertwined.
However, we did not consider an explicit correlation between the two relations
which, instead, may be possible. For example, SDSS observations (Bernardi et
al., 2014) have shown that at fixed velocity dispersion (which is a proxy of
halo mass, e.g. Sohn et al. 2020) brighter (i.e. more massive) galaxies have
larger $R_{e}$ (Shankar & Bernardi 2009). Such a trend may be captured by
introducing a positive covariance between the scatters of the SMHM and of the
K13 relation. We ran a few simple tests and found that adding this ingredient
to our framework does not significantly affect the implied size evolution of
MGs. When including a positive covariance, less massive galaxies tend to also
naturally be the smallest galaxies. The covariance thus ultimately generates
narrower size function at fixed stellar mass, where the abundance of compact
galaxies is now only controlled by the dispersion in size at fixed halo radius
$\sigma_{K}$. Therefore, a higher $\sigma_{K}$, we find, can produce the same
amount of compact galaxies as in a model without covariance but with a
proportionally lower value of $\sigma_{K}$. The degeneracy between
$\sigma_{K}$ and a covariance between the SMHM and the K13 relations may be
broken by probing the environmental dependence of galaxy size at fixed stellar
mass, a task that is beyond the scope of the present work.

###  4.4 Including concentration in the K13 model

Some authors have argued that galaxy sizes may be regulated also by halo
concentration (Jiang et al., 2019; Desmond, 2017; Desmond et al., 2018).
Essentially, this "concentration model" is a modified version of the K13 model
where an inverse proportionality between galaxy size and halo concentration is
also considered,

$\displaystyle R_{e}$
$\displaystyle=A_{c}\bigl{(}\frac{c}{10}\bigr{)}^{\gamma}R_{h}$
$\displaystyle=f(c)R_{h}$ (10)

where we define $f(c)=A_{c}\bigl{(}\frac{c}{10}\bigr{)}^{\gamma}$ and
$\gamma<0$. Here we will leave $\gamma$ as a redshift-independent free
parameter. For the concentration, we adopt the concentration-mass relation by
Dutton & Macciò (2014),

$\log{c}=a+b\log{M_{h}[M\odot]/10^{12}/h}$ (11)

with $a(z)=0.537+(1.025-0.537)exp(-0.718z^{1.08})$ and $b(z)=-0.097+0.024z$.
Dutton & Macciò (2014) report a log-normal scatter of about $\sim 0.11$ dex
which is independent on halo mass.

The results from the concentration model with $\gamma=-0.4,-0.6,-0.8$ are
reported in Figure 6 for the four toy models explored in this paper (see
Section 3.2). The most important feature of this Figure is that all models
struggle to reproduce the size evolution of MGs, except for Model 4
characterised by a flat high-mass SMHM slope $\delta$ and an evolving scatter
$\sigma_{\rm SMHM}$. All models predict an increase in size at fixed stellar
mass, with higher (absolute) values of $\gamma$ generating a shallower
evolution. As $\gamma$ approaches zero, the trend tends to reduce to that of
the K13 model, as expected. The departure from the K13 model is explained by
the evolution of the factor $f(c)=\bigl{(}c/10\bigr{)}^{\gamma}$ (equation 10,
Figure 7), which has the effect of slowing down the evolution with respect to
the K13 model. The predicted relatively slower size evolution in the
concentration model is roughly independent of the input SMHM due the shallow
correlation between halo mass and concentration (see eq. 11). Although the
concentration model struggles to reproduce a strong size evolution, as already
noted by Jiang et al. (2019), it cannot still be ruled out as current data may
be underestimating galaxy sizes at high redshift due to surface brightness
(e.g., Ribeiro et al., 2016; Whitney et al., 2019) and/or colour gradients
effects (e.g., van der Wel et al., 2014; Mosleh et al., 2017; Suess et al.,
2019, 2020).

Figure 6: Size evolution inferred from the concentration model for
$\gamma=-0.4,-0.6,-0.8$ (dot-dashed,dotted and dashed lines respectively). Top
left is for Model 1, top right for Model 2, bottom left for Model 3 and bottom
right for Model 4. In all panels, solid lines indicate the K13 model). The
cyan shaded area broadly indicates the range of observational constraints
allowed by current data ( $R_{e}\propto(1+z)^{-\alpha}$ with
$-1.2<\alpha<-0.8$, see Faisst et al. (2017); Patel et al. (2017); Mowla et
al. (2018)). All models struggle to reproduce the observed size evolution.
Model 4, which has a shallow high-mass slope in the SMHM and for which an
evolving $\sigma_{\rm SMHM}$ is implemented, provides a better match to the
observed trend for some values of $\gamma$. Figure 7: The redshift dependence
of the factor $f(c)$ in the concentration model for $\gamma=-0.4,-0.6,-0.8$
(dot-dashed,dotted and dashed lines respectively), for Model 1 (teal lines)
and Model 3 (magenta lines). $f(c)$ increases at earlier times, and it is
weakly dependent on the SMHM relation.

##  5 Conclusions

In this work we developed accurate and transparent semi-empirical models to
study the evolution of the size (effective radius $R_{e}$) function of massive
galaxies (MGs, $M_{star}>10^{11.2}M_{\odot}$). We assumed an input stellar
mass-halo mass (SMHM) relation to populate dark matter haloes with galaxies,
and then assigned sizes to galaxies via an empirically-motivated linear and
tight relation between $R_{e}$ and the host halo virial radius $R_{h}$. We
varied the input SMHM relation to reflect the still substantial systematic
uncertainties in the stellar mass function at both low and high redshift (see
Footnote 1 in the Introduction). More specifically, we devised four toy models
with different high-mass slopes and/or dispersions at fixed halo mass,
$\sigma_{\rm SMHM}$, to probe their impact on the size function of MGs. In
particular, we focused on the mean size growth and number density evolution
since $z\sim 3$ of compact starforming and quiescent MGs, distinguished in the
mocks via a simple halo mass-dependent quenching model with only one parameter
(Rodríguez-Puebla et al., 2015). Our main results can be summarised as
follows:

  * •

The shape and evolution of the size function is completely determined by the
halo occupation distribution implied by each model. In particular, the number
density of compact galaxies, $n_{compact}$, is a strong function of the
scatter $\sigma_{\rm SMHM}$ (Figures 1 and 2).

  * •

All models are able to broadly reproduce the fast size growth of starforming
and quiescent MGs by simply assuming a redshift-independent $R_{e}-R_{h}$
relation with a different zero point for the two populations (Figure 3) and in
ways largely independent of the shape of the input SMHM relation and of its
scatter.

  * •

In all models, the number density of compact starforming MGs peaks at around
$z\sim 2.5$ and sharply declines at later times, while the peak in the number
density of compact quiescent MGs is always delayed by a characteristic
timescale which depends on the specific model (Figure 4). Our findings thus
suggest a size growth driven by newly formed MGs at $z\gtrsim 1.5-2$, e.g.
“progenitor bias”, which plays a gradually lesser (but still important) role
at $z\lesssim 1.5$.

  * •

In models in which the scatter $\sigma_{\rm SMHM}$ is strictly constant in
time, we find that our predictions are consistent with a two-phase evolution
scenario, in which compact starforming MGs first quench into compact quiescent
MGs on a timescale of a few hundred Myr (Figure 5), and then grow in size
(possibly via dry mergers). In models in which $\sigma_{\rm SMHM}$ is instead
allowed to increase at earlier epochs, a significant proportion of quiescent
MGs must increase their sizes before final quenching as in, e.g., AGN-driven
size growth.

  * •

We also implemented another variant of the models in which $R_{e}$ is
proportional to virial radius via a halo concentration-dependent factor $f(c)$
(Jiang et al., 2019). We find that, at face value, this model struggles at
reproducing the fast size growth of the population of MGs (Figure 6), although
the data may be underestimating galaxy sizes at high redshift.

All in all, our results support the view that an accurate measurement of the
full size function of MGs, which will become available with the next
generation of observing facilities such as EUCLID and the Nancy Grace Roman
Space Telescope, will be able to set constraints on: i) the high-mass slope
and scatter of the SMHM relation, ii) the rate of evolution of the number
density of compact quiescent and starforming MGs and the related degree of
progenitor bias, iii) the quenching timescales of starforming MGs, and iv) the
evolutionary processes (mergers versus AGN feedback) driving the structural
evolution of MGs.

## Acknowledgements

We thank the anonymous referee for helping us improve the presentation of the
results. L.Z. acknowledges funding from the Science and Technologies
Facilities Council, which funded his PhD through the DISCnet Center for
Doctoral Training. F.S. acknowledges partial support from a Leverhulme Trust
Research Fellowship. We thank Ignacio Trujillo, Simona Mei and Adriana
Gargiulo for useful comments and Nushkia Chamba for providing data tables for
the measurements of $R_{1}$.

## Data Availability

All data will be shared upon reasonable request to the authors. We provide our
codes at https://github.com/lorenzozanisi/SizeModels.

##  Appendix A The $f_{Quench}(z)$ relation

In Figure 8 we show the evolution of $f_{Quench}$ from eqs. 5 and 6 for
$\mu=1,3,5$. It can be seen that in models with a higher $\mu$ the halo mass
scale above which galaxies are statistically quenched evolves much faster with
redshift, and is higher at earlier cosmic times.

Figure 8:  The evolution of $f_{Quench}$ from eqs. 5 and 6 for $\mu=5$
(magenta lines) $\mu=3$ (red lines) and for $\mu=1$ (blue lines). In both
cases, dotted lines, dashed lines and solid lines are for z=3,2,1
respectively. It can be seen that in models with a higher $\mu$ the halo mass
scale above which galaxies are statistically quenched evolves much faster with
redshift, and is higher at earlier cosmic times.

##  Appendix B Calibration of $A_{K}$

We use Sloan Digital Sky Survey DR7 data (Abazajian et al., 2009) to calibrate
$A_{K,SF}$ and $A_{K,Q}$ at z$\sim$0.1. We create a mock catalog of MGs at
z=0.1 as detailed in Section 2. We then constrain the normalization $A_{K}$
for starforming and quiescent galaxies by matching the mean size of the MGs in
our mock catalog to the mean observed semi-major axis effective radius that
best fits the light profiles of Massive Late Type Galaxies (LTGs) and Early
Type Galaxies (ETGs) from the Meert et al. (2015, 2016); Domínguez Sánchez et
al. (2018) $r-$band photometric and morphological catalogues. We define LTGs
as those objects for which $TType>0$ and ETGs those that have $TType\leq 0$.
We assume that all massive LTGs are starforming and all massive ETGs are
quiescent. The light profile is truncated as in Fischer et al. (2017).

##  Appendix C Definitions of compactness

In Figure 9 we show the number density evolution of compact quenched and
compact starforming MGs, for different definitions of compactness (including
the one adopted in the main text of this paper, that is, that of Cassata et
al. 2011), and for Model 1.

Figure 9: The evolution of the number density of compact quenched (red dashed
lines) and compact starforming (blue dotted lines) MGs, for Model 1 and
$\mu=3$, as in Figure 4. Cmpactness is defined as in Cassata et al. (2011)
(left), van der Wel et al. (2014) (center left), Gargiulo et al. (2017)
(center right) and Barro et al. (2013) (right). Distinct definitions of
compactness yield qualitatively very similar results, although quantitatively
different.

##  Appendix D Compacts in models 3 and 4

Figure 10: _Left_ : The number density of starforming and quenched MGs (cyan
and red line respectively) for models 3 and 4. Solid, dashed and dotted lines
are obtained adopting $\mu=1,2,4$ respectively. _Right:_ Prediction for the
number density of compact MGs for the two models and the different values of
$\mu$. The comparison data in the left column are from the SDSS ’cmodel’
photometry at z=0.1, and Davidzon et al. (2017); McLeod et al. (2020) (not
corrected for the M/L as it was done in Figure 4, see Grylls et al. 2020) at
higher redshift.

Figure 10 shows the evolution of the number density of MGs and the
corresponding predictions for the numbers of CQMGs and CSFMGs for Models 3 and
4. The results are qualitatively (but not quantitatively) similar to Models 1
and 2 (see Figure 4. In particular, most quenching models seem to disagree
with the current data for starforming galaxies.

##  Appendix E Using other size estimators

Recent works have explored different definitions for the size of a galaxy. For
instance, Mowla et al. (2019) and Miller et al. (2019) put forward the idea
that the radius that encloses 80% of the light, $R_{80}$, might be more
fundamental than $R_{e}$. This claim is based on the observation that (i) the
size distributions of starforming and quenched galaxies are almost identical
when using $R_{80}$ as opposed to the use of $R_{e}$, and (ii) that the shape
and evolution of the $R_{80}-M_{\rm star}$ is reminiscent of the SMHM
relation. Although we make explicit mention of $R_{e}$ throughout this paper,
our model can be used to make predictions for the size evolution and the
number density evolution of compact MGs, regardless of star formation
activity, which may be interpreted in terms of $R_{80}$, rather than $R_{e}$.
In particular, the K13 model would read

$R_{80}=A_{K,80}R_{h}.$ (12)

Using the values of $R_{80}$ for SDSS that we computed in Z20, we find that
that for MGs $\langle R_{80}\rangle\approx 4.2\langle R_{e}\rangle$.
Therefore, we expect $A_{K,80}\approx 4.2A_{K}$.

Trujillo et al. (2020) and Chamba et al. (2020) used deep imaging of local
galaxies to define $R_{1}$ as the radius that encloses the region within a
physically-motivated mass surface density of $1\rm M_{\odot}\ pc^{-2}$ (see
also Sánchez Almeida 2020). Trujillo et al. (2020) found that the scatter in
the $R_{1}-M_{\rm star}$ relation is of the order of only $\approx 0.06$ dex
across five orders of magnitude, including the regime of MGs for which the
relation, which is linear at lower masses, breaks. The $R_{1}-R_{h}$ relation
would read

$R_{1}=A_{K,1}R_{h}.$ (13)

Using the publicly available catalog of $R_{1}$ measurements from Trujillo et
al. (2020) and Chamba et al. (2020) we find that $\langle R_{1}\rangle\approx
7.8\langle R_{e}\rangle$, which implies that the normalization of the K13
relation in eq. 13, $A_{K,1}\approx 7.8A_{K}$. In Section 3.4 (Figure 3) we
have shown that using a constant value of $A_{K}$ works remarkably well to
describe the evolution of $R_{e}$. Whether this will be the case also for
$R_{1}$ will be revealed by future deep high-redshift observations.

Lastly, we would like to highlight the fact that different size definitions
provide different pieces of information: while $R_{e}$ is tight to the
concentration of the light profile (see Chamba et al. 2020), $R_{80}$ probes
also the outer regions of the galaxy. Likewise, $R_{1}$ has been proposed
based on the gas mass density threshold required to initiate star formation.
We believe that this does not necessarily make a definition of size more
fundamental than another. Thus, it is possible that distinct definitions of
galaxy sizes may be related to different physical processes generated by
distinct galaxy-halo coevolution paths.

## References

  * Abazajian et al. (2009) Abazajian K. N., et al., 2009, ApJS, 182, 543
  * Andreon (2013) Andreon S., 2013, A&A, 554, A79
  * Barro et al. (2013) Barro G., et al., 2013, ApJ, 765, 104
  * Barro et al. (2016) Barro G., et al., 2016, ApJ, 827, L32
  * Barro et al. (2017) Barro G., et al., 2017, ApJ, 840, 47
  * Behroozi et al. (2010) Behroozi P. S., Conroy C., Wechsler R. H., 2010, ApJ, 717, 379
  * Behroozi et al. (2013) Behroozi P. S., Wechsler R. H., Conroy C., 2013, ApJ, 770, 57
  * Behroozi et al. (2019) Behroozi P., Wechsler R. H., Hearin A. P., Conroy C., 2019, MNRAS, 488, 3143
  * Behroozi et al. (2021) Behroozi P., Hearin A., Moster B. P., 2021, arXiv e-prints, p. arXiv:2101.05280
  * Bernardi et al. (2010) Bernardi M., Shankar F., Hyde J. B., Mei S., Marulli F., Sheth R. K., 2010, MNRAS, 404, 2087
  * Bernardi et al. (2011a) Bernardi M., Roche N., Shankar F., Sheth R. K., 2011a, MNRAS, 412, L6
  * Bernardi et al. (2011b) Bernardi M., Roche N., Shankar F., Sheth R. K., 2011b, MNRAS, 412, 684
  * Bernardi et al. (2013) Bernardi M., Meert A., Sheth R. K., Vikram V., Huertas-Company M., Mei S., Shankar F., 2013, MNRAS, 436, 697
  * Bernardi et al. (2014) Bernardi M., Meert A., Vikram V., Huertas-Company M., Mei S., Shankar F., Sheth R. K., 2014, MNRAS, 443, 874
  * Bernardi et al. (2016) Bernardi M., Meert A., Sheth R. K., Huertas-Company M., Maraston C., Shankar F., Vikram V., 2016, MNRAS, 455, 4122
  * Bernardi et al. (2017) Bernardi M., Meert A., Sheth R. K., Fischer J. L., Huertas-Company M., Maraston C., Shankar F., Vikram V., 2017, MNRAS, 467, 2217
  * Bryan & Norman (1998) Bryan G. L., Norman M. L., 1998, ApJ, 495, 80
  * Buitrago et al. (2008) Buitrago F., Trujillo I., Conselice C. J., Bouwens R. J., Dickinson M., Yan H., 2008, ApJ, 687, L61
  * Buitrago et al. (2018) Buitrago F., et al., 2018, A&A, 619, A137
  * Cappellari (2016) Cappellari M., 2016, ARA&A, 54, 597
  * Carollo et al. (2013) Carollo C. M., et al., 2013, ApJ, 773, 112
  * Cassata et al. (2011) Cassata P., et al., 2011, ApJ, 743, 96
  * Cassata et al. (2013) Cassata P., et al., 2013, ApJ, 775, 106
  * Chamba et al. (2020) Chamba N., Trujillo I., Knapen J. H., 2020, A&A, 633, L3
  * Charbonnier et al. (2017) Charbonnier A., et al., 2017, MNRAS, 469, 4523
  * Chen et al. (2020) Chen Z., et al., 2020, ApJ, 897, 102
  * Cimatti et al. (2012) Cimatti A., Nipoti C., Cassata P., 2012, MNRAS, 422, L62
  * Conroy & Wechsler (2009) Conroy C., Wechsler R. H., 2009, ApJ, 696, 620
  * Daddi et al. (2005) Daddi E., et al., 2005, ApJ, 626, 680
  * Damjanov et al. (2011) Damjanov I., et al., 2011, ApJ, 739, L44
  * Damjanov et al. (2015) Damjanov I., Geller M. J., Zahid H. J., Hwang H. S., 2015, ApJ, 806, 158
  * Davidzon et al. (2017) Davidzon I., et al., 2017, A&A, 605, A70
  * De Lucia & Blaizot (2007) De Lucia G., Blaizot J., 2007, MNRAS, 375, 2
  * Desmond (2017) Desmond H., 2017, MNRAS, 464, 4160
  * Desmond et al. (2017) Desmond H., Mao Y.-Y., Wechsler R. H., Crain R. A., Schaye J., 2017, MNRAS, 471, L11
  * Desmond et al. (2018) Desmond H., Katz H., Lelli F., McGaugh S., 2018, preprint,  (arXiv:1808.00271)
  * Despali et al. (2016) Despali G., Giocoli C., Angulo R. E., Tormen G., Sheth R. K., Baso G., Moscardini L., 2016, MNRAS, 456, 2486
  * Diemer (2017) Diemer B., 2017, preprint,  (arXiv:1712.04512)
  * Domínguez Sánchez et al. (2018) Domínguez Sánchez H., Huertas-Company M., Bernardi M., Tuccillo D., Fischer J. L., 2018, MNRAS, 476, 3661
  * Donnari et al. (2019) Donnari M., et al., 2019, MNRAS, 485, 4817
  * Duncan et al. (2019) Duncan K., et al., 2019, ApJ, 876, 110
  * Dutton & Macciò (2014) Dutton A. A., Macciò A. V., 2014, MNRAS, 441, 3359
  * Dutton et al. (2010) Dutton A. A., Conroy C., van den Bosch F. C., Prada F., More S., 2010, MNRAS, 407, 2
  * Erfanianfar et al. (2019) Erfanianfar G., et al., 2019, arXiv e-prints, p. arXiv:1908.01559
  * Fagioli et al. (2016) Fagioli M., Carollo C. M., Renzini A., Lilly S. J., Onodera M., Tacchella S., 2016, ApJ, 831, 173
  * Faisst et al. (2017) Faisst A. L., Carollo C. M., Capak P. L., Tacchella S., Renzini A., Ilbert O., McCracken H. J., Scoville N. Z., 2017, ApJ, 839, 71
  * Fall (1983) Fall S. M., 1983, in Athanassoula E., ed., IAU Symposium Vol. 100, Internal Kinematics and Dynamics of Galaxies. pp 391–398
  * Fan et al. (2008) Fan L., Lapi A., De Zotti G., Danese L., 2008, ApJ, 689, L101
  * Fan et al. (2010) Fan L., Lapi A., Bressan A., Bernardi M., De Zotti G., Danese L., 2010, ApJ, 718, 1460
  * Fang et al. (2013) Fang J. J., Faber S. M., Koo D. C., Dekel A., 2013, ApJ, 776, 63
  * Firmani & Avila-Reese (2010) Firmani C., Avila-Reese V., 2010, ApJ, 723, 755
  * Fischer et al. (2017) Fischer J. L., Bernardi M., Meert A., 2017, MNRAS, 467, 490
  * Franco et al. (2018) Franco M., et al., 2018, A&A, 620, A152
  * Gargiulo et al. (2016) Gargiulo A., Saracco P., Tamburri S., Lonoce I., Ciocca F., 2016, A&A, 592, A132
  * Gargiulo et al. (2017) Gargiulo A., et al., 2017, A&A, 606, A113
  * Gómez-Guijarro et al. (2019) Gómez-Guijarro C., et al., 2019, ApJ, 886, 88
  * Grylls et al. (2019a) Grylls P. J., Shankar F., Leja J., Menci N., Moster B., Behroozi P., Zanisi L., 2019a, MNRAS, p. 2560
  * Grylls et al. (2019b) Grylls P. J., Shankar F., Zanisi L., Bernardi M., 2019b, MNRAS, 483, 2506
  * Grylls et al. (2020) Grylls P. J., Shankar F., Conselice C., 2020, arXiv e-prints, p. arXiv:2001.06017
  * Gu et al. (2016) Gu M., Conroy C., Behroozi P., 2016, ApJ, 833, 2
  * Guarnieri et al. (2019) Guarnieri P., et al., 2019, MNRAS, 483, 3060
  * Guo et al. (2011) Guo Q., et al., 2011, MNRAS, 413, 101
  * Hearin et al. (2019) Hearin A., Behroozi P., Kravtsov A., Moster B., 2019, MNRAS, 489, 1805
  * Hopkins et al. (2009) Hopkins P. F., Hernquist L., Cox T. J., Keres D., Wuyts S., 2009, ApJ, 691, 1424
  * Hopkins et al. (2010a) Hopkins P. F., et al., 2010a, ApJ, 715, 202
  * Hopkins et al. (2010b) Hopkins P. F., et al., 2010b, ApJ, 724, 915
  * Huang et al. (2017) Huang K.-H., et al., 2017, ApJ, 838, 6
  * Huertas-Company et al. (2013) Huertas-Company M., et al., 2013, MNRAS, 428, 1715
  * Huertas-Company et al. (2016) Huertas-Company M., et al., 2016, MNRAS, 462, 4495
  * Jiang et al. (2019) Jiang F., et al., 2019, MNRAS, 488, 4801
  * Jiménez-Andrade et al. (2019) Jiménez-Andrade E. F., et al., 2019, A&A, 625, A114
  * Kawamata et al. (2015) Kawamata R., Ishigaki M., Shimasaku K., Oguri M., Ouchi M., 2015, ApJ, 804, 103
  * Kawinwanichakij et al. (2020) Kawinwanichakij L., et al., 2020, ApJ, 892, 7
  * Kocevski et al. (2017) Kocevski D. D., et al., 2017, ApJ, 846, 112
  * Kravtsov (2013) Kravtsov A. V., 2013, ApJ, 764, L31
  * Kravtsov et al. (2018) Kravtsov A. V., Vikhlinin A. A., Meshcheryakov A. V., 2018, Astronomy Letters, 44, 8
  * Kubo et al. (2017) Kubo M., Yamada T., Ichikawa T., Kajisawa M., Matsuda Y., Tanaka I., Umehata H., 2017, MNRAS, 469, 2235
  * Lapi et al. (2011) Lapi A., et al., 2011, ApJ, 742, 24
  * Lapi et al. (2018a) Lapi A., et al., 2018a, ApJ, 857, 22
  * Lapi et al. (2018b) Lapi A., Salucci P., Danese L., 2018b, ApJ, 859, 2
  * Leauthaud et al. (2012) Leauthaud A., et al., 2012, ApJ, 744, 159
  * Leja et al. (2020) Leja J., Speagle J. S., Johnson B. D., Conroy C., van Dokkum P., Franx M., 2020, ApJ, 893, 111
  * Lower et al. (2020) Lower S., Narayanan D., Leja J., Johnson B. D., Conroy C., Davé R., 2020, ApJ, 904, 33
  * Luo et al. (2020) Luo Y., et al., 2020, MNRAS, 493, 1686
  * Lustig et al. (2021) Lustig P., et al., 2021, MNRAS, 501, 2659
  * Man et al. (2016) Man A. W. S., Zirm A. W., Toft S., 2016, ApJ, 830, 89
  * Mandelbaum et al. (2016) Mandelbaum R., Wang W., Zu Y., White S., Henriques B., More S., 2016, MNRAS, 457, 3200
  * Mantha et al. (2018) Mantha K. B., et al., 2018, MNRAS, 475, 1549
  * Martig et al. (2009) Martig M., Bournaud F., Teyssier R., Dekel A., 2009, ApJ, 707, 250
  * Matthee et al. (2017) Matthee J., Schaye J., Crain R. A., Schaller M., Bower R., Theuns T., 2017, MNRAS, 465, 2381
  * McDermid et al. (2015) McDermid R. M., et al., 2015, MNRAS, 448, 3484
  * McLeod et al. (2020) McLeod D. J., McLure R. J., Dunlop J. S., Cullen F., Carnall A. C., Duncan K., 2020, arXiv e-prints, p. arXiv:2009.03176
  * Meert et al. (2015) Meert A., Vikram V., Bernardi M., 2015, MNRAS, 446, 3943
  * Meert et al. (2016) Meert A., Vikram V., Bernardi M., 2016, MNRAS, 455, 2440
  * Miller et al. (2019) Miller T. B., van Dokkum P., Mowla L., van der Wel A., 2019, ApJ, 872, L14
  * Mo et al. (1998) Mo H. J., Mao S., White S. D. M., 1998, MNRAS, 295, 319
  * Mo et al. (2010) Mo H., van den Bosch F. C., White S., 2010, Galaxy Formation and Evolution
  * Mosleh et al. (2017) Mosleh M., Tacchella S., Renzini A., Carollo C. M., Molaeinezhad A., Onodera M., Khosroshahi H. G., Lilly S., 2017, ApJ, 837, 2
  * Moster et al. (2013) Moster B. P., Naab T., White S. D. M., 2013, MNRAS, 428, 3121
  * Moster et al. (2018) Moster B. P., Naab T., White S. D. M., 2018, MNRAS, 477, 1822
  * Mowla et al. (2018) Mowla L., et al., 2018, preprint,  (arXiv:1808.04379)
  * Mowla et al. (2019) Mowla L., van der Wel A., van Dokkum P., Miller T. B., 2019, ApJ, 872, L13
  * Mundy et al. (2017) Mundy C. J., Conselice C. J., Duncan K. J., Almaini O., Häußler B., Hartley W. G., 2017, MNRAS, 470, 3507
  * Muzzin et al. (2013) Muzzin A., et al., 2013, ApJ, 777, 18
  * Naab et al. (2009) Naab T., Johansson P. H., Ostriker J. P., 2009, ApJ, 699, L178
  * Navarro et al. (1996) Navarro J. F., Frenk C. S., White S. D. M., 1996, ApJ, 462, 563
  * Newman et al. (2012) Newman A. B., Ellis R. S., Bundy K., Treu T., 2012, ApJ, 746, 162
  * Nipoti et al. (2009) Nipoti C., Treu T., Auger M. W., Bolton A. S., 2009, ApJ, 706, L86
  * Nipoti et al. (2012) Nipoti C., Treu T., Leauthaud A., Bundy K., Newman A. B., Auger M. W., 2012, MNRAS, 422, 1714
  * O’Leary et al. (2020) O’Leary J. A., Moster B. P., Naab T., Somerville R. S., 2020, arXiv e-prints, p. arXiv:2001.02687
  * O’Leary et al. (2021) O’Leary J. A., Moster B. P., Naab T., Somerville R. S., 2021, MNRAS, 501, 3215
  * Oser et al. (2010) Oser L., Ostriker J. P., Naab T., Johansson P. H., Burkert A., 2010, ApJ, 725, 2312
  * Oser et al. (2012) Oser L., Naab T., Ostriker J. P., Johansson P. H., 2012, ApJ, 744, 63
  * Patel et al. (2017) Patel S. G., Hong Y. X., Quadri R. F., Holden B. P., Williams R. J., 2017, ApJ, 839, 127
  * Peng et al. (2010) Peng Y.-j., et al., 2010, ApJ, 721, 193
  * Pillepich et al. (2018) Pillepich A., et al., 2018, MNRAS, 473, 4077
  * Poggianti et al. (2013) Poggianti B. M., Moretti A., Calvi R., D’Onofrio M., Valentinuzzi T., Fritz J., Renzini A., 2013, ApJ, 777, 125
  * Posti & Fall (2021) Posti L., Fall S. M., 2021, arXiv e-prints, p. arXiv:2102.11282
  * Puglisi et al. (2019) Puglisi A., et al., 2019, ApJ, 877, L23
  * Puglisi et al. (2021) Puglisi A., et al., 2021, arXiv e-prints, p. arXiv:2103.12035
  * Qu et al. (2017) Qu Y., et al., 2017, MNRAS, 464, 1659
  * Ragone-Figueroa & Granato (2011) Ragone-Figueroa C., Granato G. L., 2011, MNRAS, 414, 3690
  * Ribeiro et al. (2016) Ribeiro B., et al., 2016, A&A, 593, A22
  * Rodriguez-Gomez et al. (2015) Rodriguez-Gomez V., et al., 2015, MNRAS, 449, 49
  * Rodríguez-Puebla et al. (2013) Rodríguez-Puebla A., Avila-Reese V., Drory N., 2013, ApJ, 767, 92
  * Rodríguez-Puebla et al. (2015) Rodríguez-Puebla A., Avila-Reese V., Yang X., Foucaud S., Drory N., Jing Y. P., 2015, ApJ, 799, 130
  * Rodríguez-Puebla et al. (2017) Rodríguez-Puebla A., Primack J. R., Avila-Reese V., Faber S. M., 2017, MNRAS, 470, 651
  * Saglia et al. (2010) Saglia R. P., et al., 2010, A&A, 524, A6
  * Sánchez Almeida (2020) Sánchez Almeida J., 2020, MNRAS, 495, 78
  * Saracco et al. (2010) Saracco P., Longhetti M., Gargiulo A., 2010, MNRAS, 408, L21
  * Shankar & Bernardi (2009) Shankar F., Bernardi M., 2009, MNRAS, 396, L76
  * Shankar et al. (2006) Shankar F., Lapi A., Salucci P., De Zotti G., Danese L., 2006, ApJ, 643, 14
  * Shankar et al. (2010) Shankar F., Marulli F., Bernardi M., Boylan-Kolchin M., Dai X., Khochfar S., 2010, MNRAS, 405, 948
  * Shankar et al. (2013) Shankar F., Marulli F., Bernardi M., Mei S., Meert A., Vikram V., 2013, MNRAS, 428, 109
  * Shankar et al. (2014) Shankar F., et al., 2014, ApJ, 797, L27
  * Shankar et al. (2015) Shankar F., et al., 2015, ApJ, 802, 73
  * Shankar et al. (2017) Shankar F., et al., 2017, ApJ, 840, 34
  * Sherman et al. (2020) Sherman S., et al., 2020, MNRAS, 499, 4239
  * Shibuya et al. (2015) Shibuya T., Ouchi M., Harikane Y., 2015, ApJS, 219, 15
  * Smail et al. (2021) Smail I., et al., 2021, MNRAS,
  * Sohn et al. (2020) Sohn J., Geller M. J., Diaferio A., Rines K. J., 2020, ApJ, 891, 129
  * Somerville & Davé (2015) Somerville R. S., Davé R., 2015, ARA&A, 53, 51
  * Somerville et al. (2018) Somerville R. S., et al., 2018, MNRAS, 473, 2714
  * Sparre & Springel (2016) Sparre M., Springel V., 2016, MNRAS, 462, 2418
  * Stockmann et al. (2020) Stockmann M., et al., 2020, arXiv e-prints, p. arXiv:2012.05935
  * Stringer et al. (2014) Stringer M. J., Shankar F., Novak G. S., Huertas-Company M., Combes F., Moster B. P., 2014, MNRAS, 441, 1570
  * Suess et al. (2019) Suess K. A., Kriek M., Price S. H., Barro G., 2019, ApJ, 877, 103
  * Suess et al. (2020) Suess K. A., Kriek M., Price S. H., Barro G., 2020, arXiv e-prints, p. arXiv:2008.02817
  * Sun et al. (2021) Sun F., et al., 2021, ApJ, 908, 192
  * Tacchella et al. (2016) Tacchella S., Dekel A., Carollo C. M., Ceverino D., DeGraf C., Lapiner S., Mandelker N., Primack J. R., 2016, MNRAS, 458, 242
  * Tacchella et al. (2017) Tacchella S., Carollo C. M., Faber S. M., Cibinel A., Dekel A., Koo D. C., Renzini A., Woo J., 2017, ApJ, 844, L1
  * Tacchella et al. (2019) Tacchella S., et al., 2019, MNRAS, 487, 5416
  * Tadaki et al. (2020) Tadaki K.-i., et al., 2020, ApJ, 901, 74
  * Thomas et al. (2005) Thomas D., Maraston C., Bender R., Mendes de Oliveira C., 2005, ApJ, 621, 673
  * Tinker (2017) Tinker J. L., 2017, MNRAS, 467, 3533
  * Tinker et al. (2017) Tinker J. L., et al., 2017, ApJ, 839, 121
  * Tortora et al. (2018) Tortora C., et al., 2018, MNRAS, 481, 4728
  * Trujillo et al. (2007) Trujillo I., Conselice C. J., Bundy K., Cooper M. C., Eisenhardt P., Ellis R. S., 2007, MNRAS, 382, 109
  * Trujillo et al. (2011) Trujillo I., Ferreras I., de La Rosa I. G., 2011, MNRAS, 415, 3903
  * Trujillo et al. (2020) Trujillo I., Chamba N., Knapen J. H., 2020, MNRAS, 493, 87
  * Vale & Ostriker (2006) Vale A., Ostriker J. P., 2006, MNRAS, 371, 1173
  * Wang et al. (2019) Wang T., et al., 2019, Nature, 572, 211
  * Wechsler & Tinker (2018) Wechsler R. H., Tinker J. L., 2018, preprint,  (arXiv:1804.03097)
  * Whitney et al. (2019) Whitney A., Conselice C. J., Bhatawdekar R., Duncan K., 2019, ApJ, 887, 113
  * Wilman et al. (2013) Wilman D. J., Fontanot F., De Lucia G., Erwin P., Monaco P., 2013, MNRAS, 433, 2986
  * Zanisi et al. (2020) Zanisi L., et al., 2020, MNRAS, 492, 1671
  * Zhou et al. (2020) Zhou L., et al., 2020, arXiv e-prints, p. arXiv:2008.08518
  * Zolotov et al. (2015) Zolotov A., et al., 2015, MNRAS, 450, 2327
  * van Dokkum & Franx (1996) van Dokkum P. G., Franx M., 1996, MNRAS, 281, 985
  * van Dokkum et al. (2010) van Dokkum P. G., et al., 2010, ApJ, 709, 1018
  * van Dokkum et al. (2015) van Dokkum P. G., et al., 2015, ApJ, 813, 23
  * van den Bosch et al. (2014) van den Bosch F. C., Jiang F., Hearin A., Campbell D., Watson D., Padmanabhan N., 2014, MNRAS, 445, 1713
  * van der Vlugt & Costa (2019) van der Vlugt D., Costa T., 2019, MNRAS, 490, 4918
  * van der Wel et al. (2014) van der Wel A., et al., 2014, ApJ, 788, 28