# Low-Surface-Brightness Galaxies are missing in the observed Stellar Mass
Function

Juhan Kim Center for Advanced Computation, Korea Institute for Advanced Study,
85 Hoegiro, Dongdaemun-gu, Seoul 02455, Republic of Korea  Jaehyun Lee Korea
Institute for Advanced Study, 85 Hoegi-ro, Dongdaemun-gu, Seoul 02455,
Republic of Korea  Clotilde Laigle  Yohan Dubois CNRS and Sorbonne Université,
UMR 7095, Institut d’Astrophysique de Paris, 98 bis, Boulevard Arago, F-75014
Paris, France  Yonghwi Kim  Korea Institute for Advanced Study, 85 Hoegi-ro,
Dongdaemun-gu, Seoul 02455, Republic of Korea  Department of Astronomy, Yonsei
University, 50 Yonsei-ro, Seodaemun-gu, Seoul 03722, Republic of Korea
Changbom Park Korea Institute for Advanced Study, 85 Hoegi-ro, Dongdaemun-gu,
Seoul 02455, Republic of Korea  Christophe Pichon CNRS and Sorbonne
Université, UMR 7095, Institut d’Astrophysique de Paris, 98 bis, Boulevard
Arago, F-75014 Paris, France  IPhT, DRF-INP, UMR 3680, CEA, L’Orme des
Merisiers, Bât 774, 91191 Gif-sur-Yvette, France  Korea Institute for Advanced
Study, 85 Hoegi-ro, Dongdaemun-gu, Seoul 02455, Republic of Korea  Brad K.
Gibson C. Gareth Few  E.A. Milne Centre for Astrophysics, University of Hull,
Hull, HU6 7RX, United Kingdom  Jihye Shin Korea Astronomy and Space Science
Institute, 776 Daedeokdae-ro, Yuseong-gu, Daejeon 34055, Republic of Korea
Owain Snaith  GEPI, Observatoire de Paris, PSL Research University, CNRS, 5
Place Jules Janssen, 92190, Meudon, France Jaehyun Lee<EMAIL_ADDRESS>
###### Abstract

We investigate the impact of the surface brightness (SB) limit on the galaxy
stellar mass functions (GSMFs) using mock surveys generated from the Horizon
Run 5 (HR5) simulation. We compare the stellar-to-halo-mass relation, GSMF,
and size-stellar mass relation of the HR5 galaxies with empirical data and
other cosmological simulations. The mean SB of simulated galaxies are computed
using their effective radii, luminosities, and colors. To examine the cosmic
SB dimming effect, we compute $k$–corrections from the spectral energy
distributions of individual simulated galaxy at each redshift, apply the
$k$–corrections to the galaxies, and conduct mock surveys based on the various
SB limits. We find that the GSMFs are significantly affected by the SB limits
at a low-mass end. This approach can ease the discrepancy between the GSMFs
obtained from simulations and observations at $0.625\leq z\leq 2$. We also
find that a redshift survey with a SB selection limit of
$\left<\mu_{r}\right>^{e}=$ 28 mag arcsec-2 will miss 20 % of galaxies with
$M_{\star}^{g}=10^{9}~{}{\rm M_{\odot}}$ at $z=0.625$. The missing fraction of
low-surface-brightness galaxies increases to 50 %, 70 %, and 98 % at $z=0.9$,
1.1, and 1.9, respectively, at the SB limit.

Hydrodynamical simulations (767), Galaxy abundances (574), Galaxy
evolution(594), Redshift surveys(1378), Low surface brightness galaxies (940)

##  1 Introduction

The observed galaxy stellar mass function (GSMF; for a non-complete list of
references, see Bell et al., 2003; Baldry et al., 2008; Weigel et al., 2016;
Davidzon et al., 2017; Weaver et al., 2022) has usually been regarded as an
essential tool for understanding galaxy evolution and star formation history
(Madau & Dickinson, 2014; Katsianis et al., 2015; Boco et al., 2021; Adams et
al., 2021; Puchwein & Springel, 2013; Conroy et al., 2007; Grazian et al.,
2015). The galaxy stellar mass is derived from the color-dependent mass-to-
light ratio (Bell & de Jong, 2001; Madau & Dickinson, 2014) or the galaxy
spectral energy distribution (SED; Conroy, 2013; Panter et al., 2007;
Bolzonella et al., 2010) where the galaxy history of star formation activities
(Lower et al., 2020; Madau & Dickinson, 2014) are fully encoded. We are able
to study the evolution of the baryonic content in galaxies using this
encrypted global star formation history or the GSMF.

In the cosmological context, galaxies are believed to form inside dark matter
halos (for a recent review on the galaxy & halo relation, see Wechsler &
Tinker, 2018). The GSMFs are accordingly expected to follow the halo mass
functions (HMFs), but their shapes are largely different to each other. The
HMF rises more steeply in the low-mass tail (see Kim et al. 2015 for a non-
complete list of various fitting functions) while the GSMF has a much
shallower (power-law) tail below the knee of the Schechter function (Davidzon
et al., 2017; Adams et al., 2021; Song et al., 2016; McLeod et al., 2021;
Weigel et al., 2016). On the low-mass scale, this discrepancy is mainly
attributed to stellar feedback (Benson et al., 2003; Baldry et al., 2008; Silk
& Mamon, 2012; Puchwein & Springel, 2013, to name a few). Supernova (SN)
explosions and stellar winds heat and eject the interstellar medium from
galaxies in low mass halos, and accordingly suppressing the growth of galaxy
(e.g. Silk & Mamon, 2012).

Motivated by this, some efforts have been make to reproduce the observed GSMFs
by implementing stellar feedback into gravito-hydrodynamics simulation
(Vogelsberger et al., 2013; Dubois et al., 2014b; Katsianis et al., 2015;
Schaye et al., 2015; Dubois et al., 2016; Pillepich et al., 2018). The
Illustris TNG simulations are updated from the Illustris simulations
(Vogelsberger et al., 2014) to better fit the low-$z$ GSMFs. They reduced the
star formation in low-mass galaxies by adopting a redshift-dependent wind
velocity floor in the mass loading factor (for details, see Springel et al.,
2018; Pillepich et al., 2018; Naiman et al., 2018; Marinacci et al., 2018;
Nelson et al., 2018). The EAGLE simulation focused on stochastic thermal
feedback to prevent the overcooling problem, without turning-off the radiative
cooling (for details, see Schaye et al., 2015; Crain et al., 2015). They
showed that EAGLE could well reproduce the GSMF for small-mass galaxies. The
Horizon-AGN (hereafter H-AGN, Dubois et al., 2014b) simulation adopts a dual
mode of SN feedback (kinetic and thermal).

The low-mass end slope of a GSMF has been debated due to incomplete
observations for low-surface-brightness galaxies (LSBGs; Baldry et al., 2008;
Tang et al., 2021; Martin et al., 2019; Valls-Gabaud & MESSIER Collaboration,
2017; Greene et al., 2022). A well known example of this issue in low redshift
surveys (typically with the surface brightness of
$\left<\mu_{R}\right>^{e}\geq 21.82$ mag arcsec-2) is given in the figure 10
of Geller et al. (2012). It shows that the LSBGs dominate the faint-end regime
of the galaxy luminosity functions (LFs). Therefore, a somewhat shallow
surface brightness (SB) cut probably misses the majority of LSBGs at the low
brightness end of the galaxy LFs.

The pencil-beam surveys have achieved deep SB limits above the background
fluctuations, at the cost of the area of fields (see Lee et al., 2012). For
example, dos Reis et al. (2020) identified high-$z$ galaxies from the HST
CANDELS archive image data with the SB limits of 28.3 mag arcsec-2
$\leq\left<\mu_{\rm F160W}\right>\leq 28.6$ mag arcsec-2 corresponding to the
$3\sigma$ limiting magnitude in the COSMOS, GOODS-N, UDS, and EGS fields.
Grazian et al. (2015) adopt 27 mag arcsec${}^{-2}\leq\left<\mu_{\rm
F160W}\right>\leq 28.5$ mag arcsec-2 at 1$\sigma$ limit in GOODS-S. Adams et
al. (2021) derived the GSMFs from the HSC (Hyper Suprime-Cam) data of the
COSMOS & XMM-LSS fields with the SB limit of 25.9 mag arcsec-2
$\leq\left<\mu_{r}\right>^{0}\leq$ 26.7 mag arcsec-2 at 5$\sigma$ with the
aperture of $2^{\arcsec}$ diameter. Tomczak et al. (2014) measured the GSMF
from NIR imaging of the CDFS, COSMOS, and UDS fields with Magellan Baade
telescope, which reaches the $5\sigma$ depths of
$24.6<\left<\mu_{K_{s}}\right><25.2$ in a circular diameter of $0\arcsec.6$.

On the other hand, wide-field surveys impose relatively shallow SB limits to
the source detection criteria, to achieve high statistical significance in
galaxy sampling. To list a few, the SDSS has the SB selection criterion of
$\left<\mu_{r}\right>^{e}\leq 23$ mag arcsec-2 (Blanton et al., 2005). The SB
selection criteria are $\left<\mu_{r}\right>\leq 23$ mag arcsec-2 for the LCRS
(Cross et al., 2001), $\left<\mu_{bJ}\right>^{e}\leq 24.7$ mag arcsec-2 for
2dFGRS & APM (Cross et al., 2001), $\left<\mu_{r}\right>^{\rm petro}\leq 21.3$
mag arcsec-2 for HectoMAP (Sohn et al., 2018), and
$\left<\mu_{R}\right>^{e}\leq 21.82$ mag arcsec-2 for SHELS0.1 (Geller et al.,
2012). The wide field surveys enable us to obtain accurate GSMFs in the local
universe, but they are still incomplete to detect LSBGs with high completeness
due to their shallow detection limits. Given that the GSMF is the key
statistics in calibrating physical ingredients of galaxy formation models, it
is necessitated to evaluate the missing LSBG populations in the observed
GSMFs. This study aims at quantifying the fraction of missing LSBGs as a
function of the SB detection limit using the Horizon Run 5 (hereafter HR5)
simulation.

This paper is organized as follows. In Section 2, we describe the HR5
simulation and its output galaxy catalog. The galaxy properties of HR5 are
statistically analyzed in Section 3, and the high-redshift and low-redshift
GSMFs are fully addressed in Sections 4 and 5, respectively. We close this
study with discussions in Section 6. Additionally we describe the galaxy
finding method used for HR5 in A.

##  2 Simulation

###  2.1 Horizon Run 5

HR5 is a cosmological hydrodynamic simulation aiming at studying galaxy
formation and evolution from high- to intermediate-redshift (down to
$z=0.625$) in a cubic volume with a side length of
$L_{\mathrm{box}}=1049~{}\mathrm{cMpc}$. A zoomed region is set to have a
cuboid geometry of $L_{\mathrm{(x,y,z)}}^{\mathrm{zoom}}=(1049,119,127)$ cMpc
crossing the central region of the simulation box. Detailed information of the
adopted cosmology and simulation setup can be found in Lee et al. (2021). We
identify halos and galaxies from the entire snapshots of HR5 using the
Physically Self-Bound (PSB)-based galaxy finder (pGalF). A detailed
description of pGalF are given in A.

###  2.2 Star Formation and Feedback

We summarize the sub-grid physics associated with star formation and stellar
feedback adopted in HR5 here. Star formation rates are computed based on a
Schmidt law (Schmidt, 1959):

$\frac{d\rho_{\star}}{dt}=\epsilon_{\star}\frac{\rho_{\rm gas}}{t_{\rm ff}},$
(1)

where $\rho_{\rm gas}$ is the gas density in a cell, $\epsilon_{\star}$ is the
star formation efficiency per freefall time, and $t_{\rm
ff}=\sqrt{3\pi/32G\rho_{\rm gas}}$ is the freefall time of a gas cell, where
$G$ is the gravitational constant. We adopt a constant star formation
efficiency of $\epsilon_{\star}=0.02$ to reproduce the global SFR evolution.

The version of RAMSES adopted for HR5 has two modes of AGN feedback switched
by the Eddington ratio (Dubois et al., 2012),

$\chi\equiv\left({\dot{M}_{\rm BH}\over\dot{M}_{\mathrm{Edd}}}\right),$ (2)

where $\dot{M}_{\rm BH}$ is the growth rate of black hole mass and
$\dot{M}_{\mathrm{Edd}}$ is the Eddington limit. If $\chi\leq\chi_{\rm
c,jet}=0.01$, the radio (dual jet–heating) mode is turned on, and otherwise,
the quasar (thermal) mode operates. The total energy released in the thermal
form is $\dot{E}_{\mathrm{BH,h}}=\epsilon_{\rm r}\epsilon_{\rm
f,h}\dot{M}_{\rm BHL}c^{2}$ where $\epsilon_{\rm r}$ is the spin-dependent
radiative efficiency, $\epsilon_{\rm f,h}$ is the coupling efficiency of the
thermal feedback, $\dot{M}_{\rm BHL}$ is the Bondi-Hoyle-Lyttleton accretion
rate, and $c$ is the speed of light. A jet mode releases the amount of energy
of $\dot{E}_{\rm BH,j}=\epsilon_{\rm f,j}\dot{M}_{\rm BHL}c^{2}$, where
$\epsilon_{\rm f,j}$ is the spin-dependent coupling efficiency of the jet mode
(McKinney et al., 2012). Therefore, AGN feedback is controlled by the two free
parameters, $\chi_{\rm c,jet}$ and $\epsilon_{\rm f,h}$. We adopt
$\epsilon_{\rm f,h}=0.15$ to reproduce the observed $M_{\rm BH}$–$M_{\star}$
relation at $z=0$ (Booth & Schaye, 2009; Dubois et al., 2012) and $\chi_{\rm
c,jet}=0.01$ (Merloni & Heinz, 2008). The gas accretion rate is capped by the
Eddington limit in HR5. For more details of AGN feedback in HR5, see Dubois et
al. (2014a) and Dubois et al. (2021).

In HR5, we assume stellar feedback operated by the winds from asymptotic giant
branch (AGB) stars and Type Ia and Type II supernovae (SNe). The amount of
energy released from an SN is set to $2\times 10^{51}$ erg. When stellar
particles are sufficiently young, we deposit 30% of the SNII energy as kinetic
energy and the rest as thermal. Old stellar populations release their energy
in the SNIa and AGB wind modes, deposit the energy in its thermal form, and
eject their mass to the nearest grid cells. The initial mass function (IMF) of
stellar populations is one of primary parameters governing stellar feedback.
We adopt a Chabrier IMF (Chabrier, 2003) for all stellar populations in HR5.

###  2.3 Size & Brightness of Simulated Galaxies

We define the galaxy size as the half mass radius of stellar components
projected on an $XY$ plane of the simulation coordinate. We assume the density
peak of stellar mass distribution as the center of a galaxy. The rest-frame
SEDs of galaxies are modeled using the masses, ages and metallicities of
stellar particles based on the E-MILES single stellar population SEDs
(Vazdekis et al., 2012, 2016; Ricciardelli et al., 2012). We assume a Chabrier
IMF in this calculation, for consistency with the stellar population model
adopted in HR5. We also compute the luminosities of galaxies using the
photometric predictions of E-MILES in the Johnson and SDSS filter systems.

##  3 Distributions of Galaxy Properties

###  3.1 Aperture Correction

Numerical simulations have suffered from the excess of BCG-scale galaxies,
compared to the empirical GSMFs, and thus the idea of contained mass has been
proposed to reconcile the discrepancy(Schaye et al. 2015; McCarthy et al.
2017; Adams et al. 2021; Tang et al. 2021, and for an extensive discussion on
the choice of the aperture size, see Donnari et al. 2019). The BCG of a
cluster can have stellar mass overestimated due to its far-extended stellar
envelope tangled with the intra-cluster light (Zwicky, 1951; Gonzalez et al.,
2007; Guennou et al., 2012; Mihos, 2019; Yoo et al., 2021; Montes, 2022). For
instance, Pillepich et al. (2018) demonstrated that the galaxies of total
stellar mass of $M\sim 10^{12}\,\mathrm{M}_{\odot}$ have half their mass
outside 30 pkpc. Observationally, however, the galaxy stellar mass is
generally measured in the Petrosian aperture or by assuming a Sérsic-model
profile (Graham et al., 2005), both of which may neglect the extended stellar
components of the BCGs. Following the idea of the aperture correction, we also
adopt the 3-D aperture of a 30 pkpc radius, to alleviate the mass
overestimation of BCG-scale galaxies.

###  3.2 Relation between Stellar Brightness and Mass

The global star formation rate (SFR) of the universe peaks at $z\sim$ 2–3 and
declines after the cosmic noon (Hopkins, 2004; Behroozi et al., 2013; Madau &
Dickinson, 2014), resulting in the increase of the mass-to-light ratios of
galaxies (van de Ven et al., 2003; van der Wel et al., 2005).

Figure 1: (bottom): Redshift evolution of the galaxy stellar-mass function
with respect to the absolute magnitude in the rest-frame $r$ band. For
comparison, we overlap the observed relation for blue spheroids from the GAMA
survey (the hatched area, Mahajan et al., 2018) with our fitting model simply
extrapolated to $z=0$ (the colored region). (top): Fitting residuals. The
filled circle is the median in each stellar-mass bin and the colored box is
the 1–$\sigma$ distribution. The error bar marks the 2–$\sigma$ scattering.

Bell et al. (2003) presented the relation between the stellar mass and the
$r$-band absolute magnitude ($\mathcal{M}_{r}$) as a function of redshift and
rest-frame galaxy color ($g-r$) as

$\displaystyle\log_{10}\\!\left({M_{\star}^{g}\over M_{\rm\odot}}\right)$
$\displaystyle=$ $\displaystyle
1.097(g-r)\\!-\\!0.406\\!-\\!0.4(\mathcal{M}_{r}\\!-\\!4.67)$ (3)
$\displaystyle+f_{\rm Bell}(z),$

where $f_{\rm Bell}(z)=-0.19z$ reflects the redshift evolution in the
2MASS/SDSS galaxies. However, this relation is based on the sample of the
local universe and only applicable to low redshift galaxies. Hence, instead of
$f_{\rm Bell}$, we propose a new fitting model of

$f_{\rm HR5}(z)=-\left\\{0.02+4.96(g-r)\right\\}\log_{10}(1+z),$ (4)

which is formulated to capture the co-evolution of redshift and galaxy color
at $0.625\leq z\leq 5$.

Figure 1 shows the evolution of the mass-to-brightness relation at given
galaxy stellar mass as a function of redshift. In this plot, the r-band
magnitude $\mathcal{M}_{r}$ is given at fixed $M^{*}_{g}$ because galaxy
stellar mass is intrinsic in simulations while the brightness is derived based
on assumed stellar models. Throughout this paper, we maintain this convention.

At $z=0.625$, galaxies tend to be fainter than their higher-redshift
counterparts at a fixed mass due to aging of the stellar populations. The
aging effect is taken into consideration by the correction term ($f_{\rm
Bell}$ or $f_{\rm HR5}$). In the top panel of the figure, we show the
residuals of the fit. Here the residual absolute magnitude ($\mathcal{R}$) is
the difference of the absolute magnitude between the simulated galaxy and the
modeled one and is written as

$\mathcal{R}(M^{g}_{\star})\equiv\mathcal{M}_{r}({\rm
sim})-\mathcal{M}_{r}({\rm model}),$ (5)

for a given galaxy stellar mass. The majority of the 1–$\sigma$ scatterings
are within the 0.2 magnitude error.

###  3.3 Environmental Dependence of Galaxy Brightness

Now, we present the environmental effects on galaxy absolute magnitude in the
SDSS $r$-band. From the simulation catalog of HR5, we divide entire galaxy
sample ($G_{\texttt{all}}$) into the satellite ($G_{\texttt{sat}}$), central
($G_{\texttt{cen}}$), and isolated ($G_{\texttt{iso}}$) samples. The most
massive galaxy in each halo is classified as centrals, while the satellite
sample contains the rest, less-massive galaxies. The isolated galaxies are
those with no companion in their FoF halos.

In the bottom right panel of Figure 2, central galaxies ($G_{\texttt{cen}}$)
tend to be brighter than the average obtained from the $G_{\texttt{all}}$
sample, while this is less significant for more massive galaxies. Redshift-
evolution of $\Delta\mathcal{M}_{\rm r}$ is apparent at the low-mass end due
to the decline of SFRs in small satellite galaxies and the late formation of
small central galaxies. The bottom left panel shows this more clearly.
Satellite galaxies are generally fainter than the central galaxies of the same
mass, and this tendency is more pronounced at a low-mass end.

The two top panels demonstrate that the isolated galaxies ($G_{\texttt{iso}}$)
are fainter than both the central and satellite galaxies. Given that a younger
stellar population is brighter at a fixed mass, these results indicate that
star formation activities are lowered in the isolated galaxies earlier than
any other galaxy groups in HR5. The massive satellite galaxies are more likely
to have been the centrals of groups or clusters until recently because their
merging timescales are shorter than those of smaller ones (e.g., Binney &
Tremaine, 2008; Boylan-Kolchin et al., 2008). This is consistent with the
small magnitude difference between the centrals and satellites at the massive
end. On the other hand, small satellites are not only vulnerable to
environmental effects (Samuel et al., 2022) but also likely to have
experienced environmental effects in their hosts for a long period of time due
to long merging timescales (Lee et al., 2018).

Figure 2: Relation between environments and galaxy brightness as a function of
galaxy stellar mass. The differences of average absolute magnitudes are shown
between $G_{\texttt{cen}}$ & $G_{\texttt{all}}$, $G_{\texttt{sat}}$ &
$G_{\texttt{cen}}$, $G_{\texttt{cen}}$ & $G_{\texttt{iso}}$, and
$G_{\texttt{sat}}$ & $G_{\texttt{sio}}$ samples clockwise from the bottom-
right panel, respectively.

###  3.4 Galaxy Luminosity Functions

Ilbert et al. (2005) and Dahlen et al. (2005) derived galaxy LFs at $z\leq 2$
from the VIMOS-VLT Deep Survey, the HST ACS and GOODS South in the rest frame
of R-band of the Cousins filter system. Given the similar wavelength coverage
to that of the SDSS $r$-band, we simply use their LFs for a direct comparison
with the LFs of HR5. We also add the LFs obtained from the FORS Deep Field
data observed with the ESO VLT by Gabasch et al. (2006), who derived the
characteristic parameters of the Schechter function in the rest-frame $r$-band
from the observations. Since the characteristic absolute magnitude of the
Schechter function in Faber et al. (2007) are given in the $B_{\rm Johnson}$
band, we convert the SDSS $r$-band magnitudes of the HR5 galaxies to the the
$B_{\rm Johnson}$ band magnitudes by adopting the band transforms of $B_{\rm
Johnson}=g+0.115+0.370\times 0.82$ and calculate $\mathcal{M}_{r}^{\star}$ by
assuming the average color of $(g-r)=0.82$ (for details, see Faber et al.
2007).

We summarize the comparison between HR5 and the observations in Figure 3, in
which the LFs of HR5 in the rest-frame $r$-band are illustrated in the thick
black solid lines. As mentioned above, we apply the 30 pkpc aperture cut that
is expected to reduce the total brightness of the BCGs substantially. In spite
of the aperture cut, we still have substantial excess of the BCG-scale
galaxies in HR5 at $z\leq 1.13$. Another notably feature seen in Figure 3 is
that the excess of the simulated LFs to the observations becomes larger with
decreasing $\mathcal{M}_{\rm r}$ below the knee of the Schechter functions at
all the redshifts. This is also seen in the galaxy LFs or GSMFs of H-AGN (see
Kaviraj et al., 2017). Kaviraj et al. (2017) speculated that the excess of
low-mass galaxies is attributed to insufficient stellar feedback in H-AGN.
This systematic excess on both mass or luminosity scales will be addressed
quantitatively in the following sections.

Figure 3: Galaxy luminosity functions at $z=0.625$, 1.13, 1.89, and 3 (from
the bottom panel). Thick solid lines are the mock galaxy luminosity functions
from HR5. The hatched or shaded areas are observed Schechter functions (Ilbert
et al. 2005; Dahlen et al. 2005; Gabasch et al. 2006; Faber et al. 2007).

###  3.5 Stellar Mass versus Halo Mass Relation

The stellar-to-halo-mass (SHM) relation is one of key properties that
cosmological hydrodynamical simulations should be able to reproduce. The SHM
relation connects galaxy observations to $N$-body simulations, which enables
one to pin down the cosmological models and to study the coupled evolution of
halos and galaxies.

Figure 4 shows the SHM relation of four different simulations (filled circles
with error bars) and several observations at $z\sim 0$, 0.6, 1.3, 3, and 5. In
this plot, no aperture cut is applied to HR5, for a fair comparison with other
simulations. The results with an aperture cut is presented in Figure 5. At
$z=3$ and 5, the relations obtained from HR5 well recover the observations,
while at $z\sim 1.3$ and 0.6, HR5 overproduces the stellar mass on the BCG
scale, compared to Behroozi et al. (2013) (the grey shaded regions). HR5 seems
to be comparable with the SHM relations of Kravtsov et al. (2018) (open blue
stars) and Golden-Marx et al. (2022) (pink hatched) at the massive end that
are based on the local SDSS BCG catalogs ($0\leq z\leq 0.15$). The redshift
evolution below $z\sim 1$ is known to be insignificant in the SHM relation
(Behroozi et al., 2013; Legrand et al., 2019; Shuntov et al., 2022), and thus
the comparison between HR5 at $z=0.625$ and the local observations is viable
to some degree.

We also overplot the SHMs of other simulations, EAGLE (Schaye et al., 2015),
H-AGN (Dubois et al., 2014b), and TNG100 (Pillepich et al., 2018). At all the
redshifts investigated in this study, EAGLE and TNG100 show the SHM
distributions well agreeing with the empirical data, particularly in
$M_{\star}^{g}\lesssim 10^{11}\,\mathrm{M}_{\odot}$.In H-AGN, galaxies are
indeed by far too massive at fixed halo mass compared to other simulations and
observations, especially at lower $z$. As mentioned above, this may be due to
insufficient SN feedback (and consequently higher global star formation
history) at high redshift (Kaviraj et al., 2017). On the other hand, all the
simulations produce more stellar mass than the empirical relation at high
halo-mass end at which AGN feedback dominantly regulates star formation; the
slopes of SHM relation is not reproduced at the massive end at low $z$. In
Figure 5 we show the results of the 30 pkpc aperture correction applied only
to HR5, where the discrepancy between simulations and observations at high-
mass end is to some degree relieved. However, the slope still does not seem to
be consistent with the empirical results.

Figure 4: Relations between $M_{h}$ and $M_{*}^{g}$. Symbols with error bars
are the medians and 1–$\sigma$ scatters in each galaxy mass bin of the HR5(red
filled circles), EAGLE (blue filled squared), TNG100 (blue open squares), and
H-AGN(yellow filled hexagons) at $z\sim$ 0.6, 1.3, 3, and 5 (from bottom to
top). In this plot, we do not apply any aperture corrections to HR5, for a
fair comparison with other simulations. Other simulation results are those
provided by Shuntov et al. (2022). Observations: Golden-Marx et al. (2022)
(pink hatched, $z\sim 0$), Behroozi et al. (2013) (grey shades), Foucaud et
al. (2010) (black triangles, $z~{}1.3$), Durkalec et al. (2015) (open navy
squares, $z\sim 3$ ), Kravtsov et al. (2018)(open stars, $z\sim 0$), Shuntov
et al. (2022)(green hatched), Legrand et al. (2019) (vertically hatched), and
Stefanon et al. (2021a)(open circles, $z=6$). Figure 5: Same as Fig. 4 but
with the 30 pkpc aperture correction only to HR5.

###  3.6 Effective Sizes of Galaxies

Figure 6: Comparison of galaxy size distributions between observations and HR5
at $z=0.625$. The background contour regions are the galaxy distributions of
HR5 in this plane. From the inner most contours, we show distribution of 10 %,
68% (1–$\sigma$), 95% (2–$\sigma$), and 99.7% (3–$\sigma$) around the median
effective radius. Observations: Nedkova et al. (2021, lines), Trujillo et al.
(2020), Fernández Lorenzo et al. (2013), and Dabringhausen & Fellhauer (2016).

We look into the size distribution of HR5 galaxies compared with observations
in Figure 6, where observations are plotted with various symbols. The
background colored contour regions are obtained from the HR5 galaxies at
$z=0.625$. The lower outline of the size distributions of the HR5 galaxies are
above 1 pkpc, which is the limit imposed by the simulation resolution. This
means that the 1 pkpc resolution is probably insufficient to properly simulate
the early-type galaxies with mass below $M_{\star}^{g}\sim 10^{11}~{}{\rm
M_{\odot}}$ (small red and gray circles in the figure) at $z=0.625$,
suggesting the necessity of higher-resolution simulations

Figure 7 shows the time evolution of galaxy size $R_{1/2}$ as a function of
stellar mass. At high redshift, galaxies are compact and the stellar-mass
dependency is low, compared to those at lower redshift. However, as time goes
by, the correlation between galaxy size and stellar mass becomes stronger (see
also Dubois et al., 2016). The scatter of the galaxy-size distribution
increases with increasing mass and decreasing redshift.

Figure 7: Galaxy mass-size relation as a function of redshift. The
distribution of galactic size for each stellar mass bin is shown with a filled
circle (median), a colored box (1–$\sigma$), and an error bar (2–$\sigma$).
The gray region is from the observations of Roy et al. (2018) at $0.4<z<0.6$.
and the hatch areas are obtained from the star-forming galaxy samples in the
GOODS-north region by Ichikawa et al. (2012). For simulation comparisons, we
add the Illustris (thick-dashed) and TNG100 (thick-solid line) results
obtained from Pillepich et al. (2018).

###  3.7 Cosmological Surface Brightness Dimming Effect

The stellar mass of observed galaxies can be derived from their luminosities
and colors based on the stellar population models and IMFs, while it is
intrinsically given from the sum of stellar particle masses in simulations.
Thus, it has been suggested that the discrepancy between the GSMFs of
simulations and observations in a low-mass end can be eased with the
assumption that a substantial fraction of small low-surface-brightness (LSB)
galaxies are missed due to the SB limit of observations (Geller et al., 2012;
Driver, 1999; Tang et al., 2021).

In a flat universe, the half-mass angular radius of a galaxy is given by

$\theta^{\arcsec}_{1/2}\simeq 206,265{r_{1/2}\over d_{c}(z)},$ (6)

where $r_{1/2}$ is the half-mass radius in a proper scale and $d_{c}$ is the
comoving distance to the galaxy, which is a function of the cosmological
parameters via

$d_{c}(z)={c\over H_{0}}\int_{0}^{z}{dz\over E(z)},$ (7)

where
$E(z)\equiv\\{{\Omega_{m}^{0}(1+z)^{3}+\Omega_{k}(1+z)^{2}+\Omega_{\Lambda}}\\}^{1/2}$.
The luminosity distance is related to the comoving distance as $d_{\rm
L}=d_{c}(1+z)$. Now we set the average SB of a galaxy to be the brightness
averaged within its half-mass radius ($r_{1/2}$). Note that $r_{1/2}$ is often
called as the effective radius. The surface brightness of a galaxy in the unit
of mag arcsec-2 is

$\hskip-5.69046pt\left<\mu_{1/2}\right>\\!=\\!\mathcal{M}_{\odot}\\!-\\!{5\over
2}\log\left[{L_{1/2}/d_{\rm L}^{2}\over L_{\odot}/(10{\rm
pc})^{2}}\right]\\!+\\!{5\over 2}\log\left(\pi\theta_{1/2}^{\arcsec
2}\right),$ (8)

where $\mathcal{M}_{\odot}$ is the absolute magnitude of the Sun and
$L_{\odot}$ is the solar luminosity. The first two terms in the right-hand
side of equation (8) is the observed magnitude of the galaxy, while the last
term is added to convert the observed magnitude into average surface
brightness per square arcsecond.

For an observation adopting a filter ($X$) of a finite bandwidth, Equation (8)
needs to be modified into

$\displaystyle\left<\mu_{1/2}^{X}\right>$ $\displaystyle\simeq$
$\displaystyle\mathcal{M}_{\odot}^{X}-{5\over
2}\log\left({\mathcal{L}_{1/2}\over\mathcal{R}_{1/2}^{2}}\right)+5\log(1+z)$
(9) $\displaystyle+\mathcal{K}(z;X)+40.33,$

where $\mathcal{K}$ is the $k$-correction term for the redshifted spectral
energy distribution (SED) of the galaxy, $\mathcal{L}_{1/2}\equiv
L_{1/2}^{X}/L_{\odot}^{X}$, and $\mathcal{R}_{1/2}\equiv r_{1/2}/1{\rm kpc}$,
where $r_{1/2}$ is the half-light radius in a proper scale. Equation 9 can be
simplified as

$\left<\mu_{1/2}^{X}\right>\simeq\mathcal{M}^{X}+5\log\mathcal{R}_{1/2}+\mathcal{C}(z;X)+41.08,$
(10)

where $\mathcal{M}^{X}$ is the absolute magnitude of a galaxy and
$\mathcal{C}(z;X)$ is the cosmological expansion term:

$\mathcal{C}(z;X)\equiv\mathcal{K}(z,X)+5\log(1+z),$ (11)

which depends only on redshift, irrespective of the cosmological model
parameters (Stabenau et al., 2008). The last term comes from the difference
between the luminosity distance and angular-diameter distance. When the
cosmological factor is not included, we call it the rest-frame surface
brightness as

$\left<\mu_{1/2}^{\prime
X}\right>\simeq\mathcal{M}^{X}+5\log\mathcal{R}_{1/2}+41.08.$ (12)

If we neglect the expansion effect, the surface brightness of a galaxy is
invariant irrespective of the distance from an observer.

###  3.8 A Model for k-Corrections

Observers can obtain the rest-frame magnitude of a target object by applying
the $k$-correction to the observed magnitude (Blanton et al., 2005). In the
analysis of simulated galaxies, conversely, it is required to apply the
$k$-correction in the opposite direction, to convert the rest-frame magnitude
to the observed magnitude at a given redshift. In this study, we take the
second approach to simulate the magnitude in observed frames.

For the simulated galaxies, we may directly shift their rest frame SEDs to a
target redshift and measure the absolute magnitudes in the frame. We have
found that the $k$-correction can be formulated as a function of the rest-
frame color and redshift. The best-fit model for the $k$-correction is

$\hskip-8.5359pt\mathcal{K}^{\rm fit}_{r}\\!=\\!\left\\{13.3(g\\!-\\!r)_{\rm
rest}\\!-\\!0.5\\!-\\!3.5(z\\!-\\!1.47)^{2}\right\\}\log(1\\!+\\!z),$ (13)

for the redshift range of $0.625\leq z\leq 2$. Figure 8 shows the fitting
results at several redshifts. The $k$-correction difference is defined as
$\Delta\mathcal{K}\equiv\mathcal{K}_{r}^{\rm sim}-\mathcal{K}_{r}^{\rm fit}$,
where $\mathcal{K}^{\rm sim}$ is directly derived from the simulated SEDs.
They seem to have a small non-linearity in the dependence on the color. Also
in Figure 9, we show the redshift evolution of $k$-correction as a function of
the galaxy rest-frame color. This figure shows that a smaller $k$-correction
is required for bluer galaxies in the $(g-r)_{\rm rest}$ color space.

Figure 8: Fitting uncertainties of the $r$–band $k$-correction. Symbols locate
the median values and boxes show the 1–$\sigma$ dispersion while the error
bars indicate the 2–$\sigma$ scatters. Figure 9: The fitting model of the
$k$-correction formulated in Eq. (13). We colorized the distribution of the
modeled $k$-correction for each $(g-r)_{\rm rest}$ as shown in the right scale
bar. Note that the valid redshift range for the fit is $0.625\leq z\leq 2$,
but we extrapolate the fit down to $z=0$, in the redshift range marked by
gray.

It is worthwhile to note that equation (13) introduces the logarithmic form of
redshift rather than polynomials of $z$ (c.f., O’Mill et al., 2011), to avoid
a linear increase of $\mathcal{K}(z)$ at high redshifts. Moreover, the
$k$-correction of equation (13) may be a linear function of $\log(1+z)$ with a
slope depending on galaxy colors or spectral ages (Stabenau et al., 2008). For
the typical range of galaxy colors ($-0.2\lesssim{(g-r)}\lesssim 2$), we have
$-0.6\lesssim\mathcal{K}^{\rm fit}/\log(1+z)\lesssim 6$, which is similar to
the observational result derived by Stabenau et al. (2008), who measured the
change of the surface brightness of galaxies with various types of samples up
to $z\sim 5$. We use this modeled $k$-corrections to quantify the contribution
of the $k$-correction to the GSMFs of galaxies of mock surveys.

###  3.9 Distribution of the Surface Brightness of Galaxies

We calculate the mean SB of the simulated galaxies of HR5 using equation (9)
at four redshifts of $z=0.625$, 1.13, 3, and 5. In each stellar mass bin, we
measure the median, 1–$\sigma$, and 2–$\sigma$ distributions as shown in
Figure 10. More massive galaxies tend to have higher SB because the increase
of brightness or stellar mass outpaces the increase of the projected surface
area. Although galaxies are more compact and brighter at a fixed mass at
higher $z$, the SB is higher at lower $z$ due to thecosmological expansion
term $\mathcal{C}(z;X)$ in Equation 11. Accordingly, for instance, a galaxy
redshift survey with the target SB limit of $\left<\mu_{r}\right>$=26 mag
arcsec-2 will be able to build a complete catalog of galaxies more massive
than about $1\times 10^{11}{\rm M_{\odot}}$ in the redshift range of
$0.7\lesssim z\lesssim 2$. However, if the SB limit is lower than e.g.,
$\left<\mu_{r}\right>$=24 mag arcsec-2, the galaxy survey would only cover
$M_{\star}\gtrsim 10^{12}{\rm M_{\odot}}$ in the redshift range.

Figure 10: Distributions of mean SB of simulated galaxies in HR5. The filled
circle with error bars marks the median value and the 2–$\sigma$ (95%) scatter
while the colored box ranges the 1–$\sigma$ (68%) distribution of the mean
surface brightness.

Figure 11 demonstrates the impact of the cosmological expansion term
$\mathcal{C}(z;X)$ in Equation (10) on the SB of galaxies. In this case, the
SB is only governed by the evolution of astrophysical parameters: effective
radius and the absolute magnitude. Because of their compact size and younger
stellar populations, higher-$z$ galaxies have notably higher SB.

Figure 11: Same as Fig. 10 but excluding the cosmological factor in Equation
(10). We add the mean SB-stellar mass relation found in New Horzion simulation
(Jackson et al., 2021) at $z=0.25$.

Now, we examine the effect of the SB detection limit on mock galaxy surveys.
Mock galaxy surveys are made with various SB limits without implementing any
apparent magnitude limit, to isolate the SB limit effect. Therefore, this is
an extreme case of the extragalactic surveys, where the SB of a galaxy is the
only constraint on source detection. Figure 12 shows the effect of SB limit on
the GSMF at $z=0.625$. As expected, the SB limit has more impact in low mass
galaxies; dwarf galaxies are missed more with a lower SB limit. The vertical
axis ($\mathscr{C}_{\rm sb}$) of the top panel is the survey completeness
defined as

$\mathscr{C}_{\left<\mu_{r}\right>}\equiv{N_{\rm obs}\over N_{\rm all}},$ (14)

where $N_{\rm obs}$ is the number of galaxies which satisfy the SB condition
of the mock survey and $N_{\rm all}$ is the total number of target galaxies.

Figure 12: Mock observations of the GSMFs with various surface-brightness
limits at $z=0.625$. (bottom) The black solid line is the stellar mass
functions with all galaxies while the other lines are obtained with surface
brightness cuts. Error bars mark the 1–$\sigma$ Poisson errors. (top) The
dependence of the completeness of mock samples on the input surface-brightness
cut.  Figure 13: Distributions of the observation completeness for several
surface brightness cuts (values in the legend). From the top panel, the survey
completeness (curves) is shown in a mock survey of galaxies at $z=0.875$,
1.13, and 1.89. The horizontal solid line marks the 100% completeness.

In Figure 13, we show the completeness of the GSMF in mock observations at
$z=0.875$, 1.13 and 1.89 (from the top panel) as a function of the SB limit.
For example, we may fail to detect 10 % of galaxies with
$M^{g}_{\star}=10^{9}~{}{\rm M_{\odot}}$ if the SB limit of the survey is
$\left<\mu_{r}\right>=29$ mag arcsec-2 at $z=0.875$. The fraction of missing
LSBGs inevitably increases with increasing redshifts.

It is valuable to note that, in spite of strong star formation activities at
high redshifts, the cosmological SB dimming effect becomes dominant,
eventually lowering the SB of galaxies, as shown in Figure 14. This figure
suggests that the GSMF at $z\lesssim 5$ would be complete down to
$M=10^{11}{\rm M_{\odot}}$, if a survey adopts the SB condition of
$\left<\mu_{r}\right>\leq 27$ mag arcsec-2.

A caveat of this analysis is that the size of a galaxy is limited by the
resolution of a simulation. HR5 has a spatial resolution of $\sim 1\,$pkpc,
while galaxies can have half mass radii close to the resolution limit in the
mass range of $M_{\star}^{g}\sim 10^{9}-10^{10}\mathrm{M}_{\odot}$ (see the
red solid line and the thin blue shade in Figure 6). Thus, the size low mass
galaxies can be overestimated in HR5. Figure 6 implies that low mass quiescent
galaxies in HR5 may have sizes twice larger than observations, on average,
which leads to 1.5 mag arcsec-1 fainter SB. Therefore, depending on the
morphology distribution of the simulated galaxies in HR5, the missing fraction
of LSBGs given by the SB limit of $\left<\mu_{r}\right>=26.5$ mag arcsec-2 may
be close to that given by the SB limit of $\left<\mu_{r}\right>=28$ mag
arcsec-2.

Figure 14: Redshift evolution of galaxy fractions with the mean surface
brighter than $\left<\mu_{r}\right>$=29 (top), 28 (middle), and 27 mag
arcsec-2 (bottom panel).

Figure 15 presents an example of the distribution of LSBGs identified around
Cluster 1 region in HR5 at $z=0.625$ (Lee et al., 2021). The LSBGs located
within a radius of 10 $h^{-1}$ cMpc from the cluster’s center are marked by
the open circles with colors for various SB ranges. The LSBGs are faint but
distinct with compact sizes, without prominent tidal features. There are two
groups of galaxies falling into the cluster in the upper right and left parts
(around 1 & 11 o’clock) of the image. The group also contains several LSBG
satellites. In this LSBG distribution, we can see that the cluster LSBGs
scatter around the central part of the cluster, and some of them stretch along
the local filament connecting the cluster to neighboring galaxy groups.

Figure 15: A pseudo color map of the star+gas density field of Cluster 1 of
HR5 at $z=0.625$. We mark the LSBGs using open circles if they are located
within 10 cMpc distance from the cluster center and have stellar masses of
$10^{9}{\rm M_{\odot}}\leq M_{\star}\leq 10^{10}{\rm M_{\odot}}$. The white,
blue, and red circles indicates the LSBGs in the SB ranges of [26, 27), [27,
28), and $[28,\infty)$ mag arcsec-2, respectively. Note that the majority of
unmarked small faint blobs scattered around the cluster are LSBGs located
beyond 10 $h^{-1}$cMpc from the cluster center. The dotted circle has the
radius of $R=1~{}h^{-1}$ cMpc. Figure 16: Comparison of GSMFs at high
redshifts. Symbols with error bars are the observations (see the main text for
references) at $z=$ 5, 6, 7, 8, 9, and 10 (from the top-left to the bottom-
right panel in the sinistrodextral order) while solid lines are the
corresponding simulated stellar mass functions of HR5. Observations: Song et
al. (2016), González et al. (2011), Duncan et al. (2014), Grazian et al.
(2015), and Behroozi et al. (2013), and Stefanon et al. (2021b).

##  4 Galaxy Stellar Mass Functions at $z\gtrsim 2$

The GSMF is one of the fundamental statistics that galaxy formation models aim
to reproduce. Since the cosmic dawn, the average star formation density of the
universe has increased through $z=2$–3. Then, star formation is suppressed
particularly in massive galaxies as AGNs start to emit sufficient energy that
blows away or heats up the interstellar medium. The global SFR declines until
today, with many implications imprinted on the GSMF, galaxy morphology, and
color distributions. The shape of the GSMF reflects a sophisticated history of
cosmic matter evolution. For example, the slope above the knee of a Schechter
function is known to be suppressed by energetic AGN feedback in massive
galaxies that is controlled by the accretion of gas via mergers or secular
evolution.

###  4.1 GSMF at $4\leq z\leq 10$

Galaxy samples compiled over the last two decades in many legacy surveys have
enabled us not only to trace the global star formation and reionization
history in the early universe, but also to study proto-galaxy formation with a
high precision. In this study we utilize various observational data (González
et al., 2011; Duncan et al., 2014; Grazian et al., 2015; Song et al., 2016;
Stefanon et al., 2021a) which have been obtained from the GOODS and CANDELS
fields observed by HST and Spitzer IRAC to name a few.

In Figure 16, the simulated GSMFs (blueish solid lines) are consistent with
observations, except for the high-mass end at $z<7$. This deviation partly
comes from the cosmic variance effect caused by the finite simulation box size
of HR5. It may, however, also be attributed to numerical artifacts of AMR,
which causes suppression of star formation between global mesh refinement
redshifts. AMR is a classical implementation to enhance the spatial and mass
resolution in Eulerian hydrodynamical simulations, in the regions where the
mass inflow becomes significant. In HR5, this global refinement is carried at
the expansion factors of $a_{\texttt{amr}}$ = 0.0125, 0.25, 0.05, 0.1, 0.2,
0.4, and 0.8, in order to maintain spatial resolution of $\sim 1$ pkpc. As the
simulation evolves, a dense gas grid inside a massive galaxy must wait for the
next $a_{\mathrm{amr}}$ to be refined into eight one-level-finer grids; hence
star formation is suppressed until the next-level refinement is allowed.

###  4.2 GSMF at $2\lesssim z\lesssim 4$

Figure 17: Comparison of GSMFs at intermediate redshifts ($z=$ 2.5, 2.5, and 4
from the left panel). The GSMFs of HR5 are denoted by the thick solid lines
while the observed mass functions are marked with dashed lines or symbols with
error bars. Observations: Behroozi et al. (2013), Davidzon et al. (2017),
Marchesini et al. (2009), Wright et al. (2018), Thorne et al. (2021), González
et al. (2011), and Duncan et al. (2014); Grazian et al. (2015).

In order to carry out a robust comparisons of GSMFs between the observational
results and simulations at intermediate redshifts, we take five different
observed GSMFs (Marchesini et al., 2009; Davidzon et al., 2017; Wright et al.,
2018; Thorne et al., 2021; Song et al., 2016) among others (for a full list of
references, see the legend in Fig. 17). Figure 17 shows the GSMFs of HR5 and
observations at the three redshifts, $z=4$, 3.5, and 2.5. The HR5 predictions
become more consistent with observations at the massive end after $z=4$ at
which global refinement takes place. From this trend, it is expected that the
simulated BCG populations will overshoot the observations at lower redshifts.
Meanwhile, we begin to observe the overproduction of low-mass galaxies at
$z\leq 3.5$. This may also indicate less efficient stellar feedback in HR5,
like H-AGN (Kaviraj et al., 2017), or that a different approach is needed to
explain this low-mass discrepancy. We address this in the next section.

##  5 Galaxy Stellar Mass Function at $0.625\leq z\leq 2$

Before carrying out a direct comparison with low-$z$ observations, let us
briefly summarize the previous findings from observations. There is a
consensus that the low-$z$ GSMF evolves slowly under $z\sim 1$ (Behroozi et
al., 2013; Legrand et al., 2019; Shuntov et al., 2022). While the GSMF at
$1\leq z\leq 2$ grows quickly with a nearly fixed shape, the slope of the GSMF
on low-mass scale rises steeply with time (Marchesini et al., 2009; Katsianis
et al., 2015; Leja et al., 2020; Adams et al., 2021).

In contrast, the GSMF of HR5 is shaped on the dwarf mass scale before $z=2$
and it barely evolves in $z\sim 1-2$. On the BCG mass scale, HR5 shows a
strong number-density evolution with decreasing redshift while observations
report relatively weak evolution. The TNG100 and EAGLE simulations demonstrate
the GSMFs that agree well with the observed GSMFs at low redshift since they
employ updated versions of the stellar winds and SN feedbacks to fit the
simulation GSMF for the low-$z$ observations (see the detailed calibration
schemes of EAGLE in Schaye et al. 2015 and TNG100 in Pillepich et al. 2018).
On the other hand, HR5 is mainly calibrated to fit the global SFR (Lee et al.,
2021).

In this section, we take two approaches to alleviate these discrepancies based
on the fact that the observational techniques for deriving stellar mass from
the brightness of galaxies are essentially different to simulations in which
galaxy stellar mass is rather intrinsic variables given by stellar particle
masses. On the BCG scale, we measure the galaxy stellar mass using finite
aperture size to mimic observations. Meanwhile, mainly for the dwarf galaxy
scale, the SB of galaxies is taken into account, to understand the shallow
slope of GSMFs in a low mass scale in observations (Martin et al., 2019;
Williams et al., 2016).

Using the galaxy samples in the COSMOS and XMM-LSS fields, together with
Spitzer IRAC, Adams et al. (2021) fitted the observed stellar mass functions
with the double Schechter functions, in the redshift range of $z=0.2-2$. They
corrected the cluster stellar mass with the 30 pkpc aperture size and found
that the number density of high-mass end populations drops dramatically
compared to the uncorrected ones.

On the other hand, Leja et al. (2020) employed a single continuity model for
the fitting functions, to take into account the redshift evolution in the
double Schechter function. They find the best fit models by investigating a
wider range of parameters. The resulting fits of GSMFs turn out to relax the
factor-of-two tension (Davidzon et al., 2018) between the specific star
formation rate (sSFR) and the stellar density by lowering 0.2 dex in the sSFR
and raising 0.2 dex in the GSMF.

Figure 18: GSMFs from the HR5 (symbols) and observations (Adams et al. 2021;
lines) at several lower redshifts. The left panel shows the uncorrected HR5
GSMF while the middle panel shows the simulation results after applying the
30kpc aperture correction only. The right panel shows the HR5 GSMF after
applying an additional constraint on the surface brightness of galaxies. Note
that Adams et al. (2021)’s data have the SB limit of 25.9 mag arcsec-2
$\leq\left<\mu_{r}\right>^{0}\leq$ 26.7 mag arcsec-2.

The left panel of Figure 18 shows the empirical fits (lines; Adams et al.
2021) and the HR5 results (symbols) before applying the 30 pkpc aperture
correction. The GSMFs of HR5 notably deviate from the observations; the over-
populations of simulated galaxies are clear on both the high-mass and low-mass
ends. However, after applying the 30 pkpc-aperture correction to the simulated
galaxies (bottom panel), we obtain a better correspondence between HR5 and
observations on the BCG scale, recovering the exponential cut-off observed at
the high-mass end (also see Schaye et al., 2015). However, there are still
several differences that should be briefly mentioned here. First, the
observation data has a flatter slope below the knees of the GSMFs than those
of HR5. Second, the simulation predicts an overpopulation of low-mass galaxies
($M_{\star}^{g}\leq 10^{10}{\rm M_{\odot}}$) compared to observations. On the
right panel of the figure, we show the effect of a SB limit of
$\left<\mu_{r}\right>\leq 27$ mag arcsec-2, which best fits the observations.
Note that the prominent knee of the observed mass function is not reproduced
even with the SB limit. This implies that we need to apply more fine tuning to
the simulation parameters and the SB limit parameters. However, it is
important to acknowledge that the redshift-evolution of dwarf population
amplitude is consistent with the observed evolution, which is hardly seen in
the original simulation outputs presented in the left panel of the figure.

Figure 19: Comparisons of GSMFs between the HR5 (symbols) and observations.
Counterclockwise from the upper-left panel, the lines are the observed GSMFs
given by Leja et al. (2020), Kawinwanichakij et al. (2020), Wright et al.
(2018), and Thorne et al. (2021) at several lower redshifts. Each color
indicates the same redshift as used in Fig. 18. In this plot we have used the
modified versions of GSMFs of HR5, after applying the SB observation limit
together with the 30 pkpc aperture correction.

In Figure 19 we compare the simulated GSMFs with more observations: Leja et
al. (2020) (upper-left), Kawinwanichakij et al. (2020) (lower-left), Wright et
al. (2018) (lower-right), and Thorne et al. (2021) (upper-right panel). As
stated above, the GSMFs derived by Leja et al. (2020) have the highest
amplitudes compared to the others. The GSMFs of HR5 show that the size of the
BCG population is comparable to the observations, except for that of Leja et
al. (2020), whose BCG population size is substantially larger and, moreover,
hardly shows any time evolution, which is in tension with the other empirical
fits (for a detailed discussion, see Leja et al., 2020).

Figure 20: Amplitude difference of the GSMFs between the EAGLE (Furlong et
al., 2015), TNG100 (Pillepich et al., 2018), and HR5 (Lee et al., 2021)
simulations and observations of Adams et al. (2021), Thorne et al. (2021), and
McLeod et al. (2021) from the top to bottom, respectively. The base GSMFs of
HR5 are corrected with the 3D-aperture of 30 pkpc and the SB limit of 27 mag
arcsec-2.

We measure the amplitude differences of GSMFs between observations and
simulations. Figure 20 shows how much the GSMFs of simulations (marked by
regions for EAGLE, lines for TNG100, and symbols with error bars for HR5)
deviate from the three different observations (Adams et al., 2021; Thorne et
al., 2021; McLeod et al., 2021). We note that the aperture correction and SB
limit is applied only to HR5. Among the three simulations, EAGLE is least
deviated from the three observations over the redshift range of $z\sim 0.5-2$.
The GSMFs of HR5 shows considerable scatters around observations while the
TNG100 has a systematic overestimation on all mass scales. To reduce this
systematic bias of TNG100, Tang et al. (2021) show that the systematic bias of
TNG100 can be significantly reduced by applying the SB limit and an apperture
cut to the mock survyes in a similar way to ours.

##  6 Discussions & Conclusions

One of the key purposes of numerical simulations is to reproduce observations
by tuning model parameters to understand underlying physics encoded in
observed properties. This is valid only when comprehensive understanding
precedes for empirical data. In this study, we have found that the surface-
brightness limit of galaxy surveys is one of the important factors fixing the
survey completeness at the faint-end of GSMFs. For instance, with the SB limit
of $\left<\mu_{r}\right>\leq 28$ mag arcsec-2, our mock surveys reach 80 %
completeness at the stellar mass of $M_{\star}^{g}=10^{9}{\rm M_{\odot}}$,
while the completeness drops sharply down to $\sim 2\%$ at $=1.89$. Even with
$\left<\mu_{r}\right>\leq 29$ mag arcsec-2, 70 % of galaxies of
$M_{\star}^{g}=10^{9}{\rm M_{\odot}}$ are missed at $z=1.89$. One certainly
needs stronger stellar and SN feedback models if observed GSMFs are simply
accepted with no evaluation of the missing LSBGs at given SB limit. However,
the fraction of LSBGs can be overestimated to some extent at the low-mass end
due to the resolution limit of HR5.

Tantalizing evidences have been reported by a number of literature that are
consistent with our claim of a substantial amount of missing LSBGs in galaxy
surveys. Greene et al. (2022) performed the image simulation to obtain 80 %
detection completeness at the SB limit of $\left<\mu_{r}\right>^{e}\sim$ 28
mag arcsec-2. van der Burg et al. (2017) also showed that they reaches
$\sim$80 % completeness of LSBGs from ESO Kilo-Degree Survey with the SB limit
of $\left<\mu_{r}\right>^{e}\sim$ 25.5 mag arcsec-2. For a lengthy list of
complementary observational references, see also Greene et al. (2022). If the
cosmological surface brightness dimming effect is added to the findings of van
der Burg et al. (2017) for LSBGs at $z=0.625$ – 2, neglecting the evolution
correction, then it leads to the SB limit of about
$\left<\mu_{r}\right>^{e}\sim$ 26 – 28 mag arcsec-2. This is roughly
consistent with the arguments presented in the main part of this paper.
However, the global SFR may not change significantly after the SB limit
corrections since the missing LSBGs insignificantly contribute to the global
SFR in a low mass range (Davies et al., 2016).

As of today, the ultra faint field at $\left<\mu_{r}\right>\gtrsim$ 28 mag
arcsec-2 is an almost uncharted territory. The Dark Energy Camera (DECam)
reaches $\left<\mu_{r}\right>=28$ mag arcsec-2 in the Dark Energy Camera
Legacy Survey searching for LSBGs (Román et al., 2021). The Dragonfly Wide
Field Survey also reaches 1$\sigma$ depths of $\left<\mu_{g}\right>\sim$ 31
mag arcsec-2 (Danieli et al., 2020), and many extra-galactic survey projects
(HST: Borlaff et al. 2019; LSST: Brough et al. 2020; Euclid: Euclid
Collaboration et al. 2022) are on-going or scheduled for ultra faint-field
observations of LSBGs.

Even with the correction to the LSB effects, the GSMFs of HR5 still show a
knee around the characteristic mass less distinct than that of observed GSMFs.
It may come from an insufficient star formation efficiency
($\epsilon_{\star}=0.02$) or inefficient SN and AGN feedback. This calls for
future studies based on the simulations with more complete physics and higher
resolution to better understand the star formation and galaxy evolution on a
wide-range of a mass scale.

The authors thank the Korea Institute for Advanced Study for providing
computing resources (KIAS Center for Advanced Computation Linux Cluster
System) for this work. JK was supported by a KIAS Individual Grant (KG039603)
via the Center for Advanced Computation at Korea Institute for Advanced Study.
JL is supported by the National Research Foundation of Korea
(NRF-2021R1C1C2011626). This work benefited from the outstanding support
provided by the KISTI National Supercomputing Center and its Nurion
Supercomputer through the Grand Challenge Program (KSC-2018-CHA-0003). Large
data transfer was supported by KREONET, which is managed and operated by
KISTI. BKG acknowledges the support of STFC through the University of Hull
Consolidated Grant ST/R000840/1, access to viper, the University of Hull High
Performance Computing Facility, and the European Union’s Horizon 2020 research
and innovation programme (ChETEC-INFRA – Project no. 101008324). This research
was also partially supported by the ANR-19-CE31-0017 http://www.secular-
evolution.org.

##  Appendix A Halo and Galaxy Findings

###  A.1 AMR & Unified Datatypes

Dark matter, super massive black holes (SMBHs), and stars are represented with
point-mass particles in RAMSES, while gas is a hydrodynamical component of a
static structured mesh. Hence, in order to post-process the RAMSES
simulations, we integrate all the elements into a single unified data
structure. We treat gas cells as particles that inherit all the mesh cell
information. The center of a cell is set as the position of the cell
‘particle’. Although gas cells are hierarchically structured, we only consider
leaf cells or terminal cells.

To reduce complexity in coding and enhance performance, we introduce a unified
datatype for the four different matter species. The unified datatype has a
fixed length of the common block containing the species, mass, position, and
velocity, while extra space is arranged as a private block. The purpose of the
private blocks is to save a union datatype for the raw simulation particle
information, which has a variable length depending on the particle type
(stellar, AGN, dark matter, and gas). The galaxy finding is processed with the
data contained in the common block, while we save the finding results by
dumping the data in the private blocks.

###  A.2 Adaptive Friend-of-Friend Method

Figure 21: Three projected density maps of massive halos identified with the
variable length in the FoF at $z=0.625$. From the left column shown are the
dark matter, stellar, and gas density fields. In the middle panels, stellar
particles with higher metallicity have more pinkish color. The white bar at
the bottom in each panel marks the scale of $l_{\rm bar}=0.5~{}h^{-1}{\rm
cMpc}$. From the top-row panel, the total mass of the FoF halo is $1.09\times
10^{13}~{}h^{-1}{\rm M_{\odot}}$, $1.42\times 10^{13}~{}h^{-1}{\rm
M_{\odot}}$, and $2.63\times 10^{13}~{}h^{-1}{\rm M_{\odot}}$, respectively.

To find virialized structures in cosmological $N$-body simulations, a
percolation method like the Friend-of-Friend (FoF) algorithm is frequently
employed. One fifth the mean particle separation is adopted as the standard
linking length to identify structures from the distribution of uniform-mass
particles. The linking length is known to find structures with a mean density
of $\sim 178$ times the critical density satisfying the cosmological
virialization condition, according to the spherical top-hat collapse model
(Suto et al., 2016). The standard linking length is formulated as

$\ell_{\mathrm{link}}\equiv{1\over
5}\left({m_{p}\over\Omega_{m0}\rho_{c}}\right)^{1/3},$ (A1)

where $m_{p}$ is the particle mass and $\rho_{c}$ is the critical density at
the current epoch.

However, in $N$-body or hydrodynamical simulations, matter components may have
varying mass not only between different types but also in a type. In this
case, the standard linking length is needed to be modified to a general form.
To identify FoF halos with varying mass particles, we use the averaged linking
length between a pair of two particles of different types given as

$\ell_{\mathrm{comb}}={1\over 2}\left({\ell_{1}+\ell_{2}}\right),$ (A2)

, where $\ell_{1}$ and $\ell_{2}$ are the linking lengths of the pair of each
particle calculated by using Equation (A1). This combined form is commutative
or pair-wise mutual.

We apply this adaptive linking length to the HR5 data of multiple particle
types to identify FoF halos. Therefore, a FoF halo may consist of multiple
types of matter components: dark matter, stellar, gas, and SMBH particles.
Figure 21 show the images of three representative FoF halos identified at
$z=0.625$. In this figure, a clear difference in density distributions are
observed between the dark matter (left) and stellar components (middle
column). The stellar distribution is more compact than the dark matter or gas.
The tidal tails of gas (right column) in the merger remnants are most
prominent among the tree components, and stellar streams between galaxies are
also seen. The stellar density map is overlaid with the stellar metallicity
colored in pink. From this, we can clearly observe that more massive galaxies
are more metal enriched in the central regions. There is a good one-to-one
correspondence among the density maps of the three different components,
confirming that the adaptive linking length produces a consistent result from
the distribution of different matter components.

###  A.3 Halo Mass Function of HR5

The halo mass function (HMF) is one of the fundamental statistics which
closely relates to the cosmological models (Press & Schechter, 1974;
Ganeshaiah Veena et al., 2018; Jenkins et al., 2001; Sheth & Tormen, 1999). It
has been used as a fiducial statistic of $N$-body simulations, because it
converges well between the $N$-body simulations of different resolutions at
most redshifts. For example, Figure 7 of Kim et al. (2015) shows consistent
fitting functions with less than a few percents of deviations.

Throughout this paper, the mass function is defined as

$\Phi(M)\equiv{\Delta N\over\Delta\log_{10}M}\,,$ (A3)

where $\Delta N$ is the number of galaxies in a mass bin of size
$\Delta\log_{10}M$. In Figure 22, we show the HMFs at $z=0.625$, 1, 2, 3, and
4. In the mass range of $10^{11}~{}{\rm M_{\odot}}\leq M_{\rm FoF}\leq 5\times
10^{14}~{}{\rm M_{\odot}}$, the HR5 HMFs are well described by the reference
model of the Sheth & Tormen functions (Sheth & Tormen, 1999).

Figure 22: Friend-of-Friend HMFs at several redshifts. Symbols with Poisson
error bars are simulated HMFs while lines are the reference mass functions
(Sheth & Tormen, 1999).

###  A.4 PGalF: A New Galaxy Finder

We developed the PSB-based Galaxy Finder (pGalF) based on the Physically Self-
Bound (PSB) method, to identify galaxies composed of multiple matter
components in HR5. We revisit PSB and describe the details of pGalF in this
subsection.

####  A.4.1 Revisiting the PSB Method

Kim & Park (2006) first introduced PSB to identify self-bound structures (or
subhalos) in their high-resolution $N$-body simulations (Kim et al., 2015,
2009). PSB operates structure identification based on the coordinate-free
density map and a web of network made of neighbor links of particles. From the
input number of neighbors, we build a network of neighboring. Densities are
measured at all particle positions, with an adaptive smoothing kernel. This
coordinate-free density field is advantageous over the regular density grid in
several aspects. The result is independent of the size of the grids and it
does not rely on a preferential shape of the particle distribution.

On the coordinate-free density map, density peaks at particle positions are
identified on the condition that all the neighboring particles should have
lower density values. We define a core-density region around a peak by
applying a water shedding method, to find the lowest threshold density for
which no other density peaks are surrounded by the threshold contour. Then,
each density peak may be extended to its core-density region, and we assume
all particles in the region become core members of the subhalo candidate.
After extracting core particles, we apply the hierarchical density contours to
separate the remaining particles into multiple hierarchical sets of particles.

The aforementioned subhalo membership is mainly controlled by two conditions:
the tidal boundary and total energy. If a particle is bound to a subhalo
candidate and it is within the tidal boundary of the subahlo, it becomes the
member of it. Unless it satisfies these two conditions, we then check the
membership with respect to other subhalos. A particle may pass the membership
check for two subhalos. Then, we set the particle to be a member of the less-
massive subhalo.

####  A.4.2 Tidal Radius

There have been various definitions for the boundary in a binary system.
(Binney & Tremaine, 1987; van den Bosch et al., 2018; Renaud et al., 2016).
They depend on the assumptions made for the simplified model, such as the
circularity of the orbital motion and point-mass or extended body of the host.

We parameterize the various definitions of the tidal radius to account for
various situations. The tidal radius of a satellite can be generalized as (Kim
& Park, 2006; van den Bosch et al., 2018)

$r_{t}=R\left({{m/M\over\alpha+\beta}}\right)^{1/3},$ (A4)

where $m$ is the satellite mass, $R$ is that distance from the host to the
satellite, and $M$ is the host mass contained within $R$. We introduce a
dimensionless parameter, $\alpha$, to account for the mass distribution of the
host at the satellite position as

$\alpha(R)\equiv 2-\left({d\ln M\over d\ln R}\right)\bigg{\rvert}_{R}.$ (A5)

Therefore, if the host is a point mass, we get $\alpha=2$. In Equation (A4),
$\beta$ is a function of orbital motion reflecting the effect of centrifugal
force on the tidal radius, and can be formulated as (van den Bosch et al.,
2018),

$\beta\equiv{\Omega^{2}R^{3}\over GM},$ (A6)

where $\Omega(\equiv V_{c}/R)$ is the angular velocity (King, 1962) for the
circular velocity, $V_{c}$. From this equation, one may easily derive that
$\beta=-2\mathcal{K}_{t}/\mathcal{W}$ where $\mathcal{K}_{t}$ is the
rotational kinetic energy of the satellite and $\mathcal{W}$ is its potential
energy. For a circular orbit, $\beta=1$ but $\beta=0$ for a radial motion. If
the satellite is not bound to the host, then we neglect this effect and set
$\beta=0$. From this parameterization we recover nearly all the modeled radii
of the satellite. For example, we recover the Jacobi radius with
$\alpha+\beta=3$ for a satellite in a circular orbit around a non-contacting
host, or we recover the Roche limit (when $\alpha+\beta=2$) for a radial orbit
around the host.

Now, we discuss the $\alpha$ parameter in more detail for the cosmological
extended object. Assuming the NFW density profile (Navarro et al., 1997) for a
host with virial mass, $M_{v}$, the enclosed mass, $M(s)$, is modeled as
(Łokas & Mamon, 2001),

$M(s)=M_{v}g(c)\left[\ln(1+cs)-{cs\over 1+cs}\right],$ (A7)

where $c$ is the concentration index, $s$ is the scaled radius to the virial
radius as $s\equiv R/R_{v}$, and

$g(c)\equiv\left[\ln(1+c)-c/(1+c)\right]^{-1}.$ (A8)

Then, we finally get

$\alpha(s)\equiv 2-\left({cs\over 1+cs}\right)^{2}.$ (A9)

The concentration index is empirically given as (Comerford & Natarajan, 2007)

$c=14.5\left({M_{v}\over 1.3\times 10^{13}~{}h^{-1}{\rm
M_{\odot}}}\right)^{-0.15}(1+z)^{-1}\,,$ (A10)

where $z$ is the redshift and $h$ is the Hubble expansion rate divided by 100
km/second/Mpc.

####  A.4.3 Hierarchical Membership Determination

Figure 23 exemplifies how the subhalo finding of PSB is working for a given
density field. In this density map, there are three density peaks and each
peak has the corresponding core region marked with A, B, and C. First, the A
and B regions are found separated by the inner most contour while the C area
is delineated by another contour that separates the C region from the other
two core regions. This separation is done via a water-shedding technique. It
is a percolation method which grows a volume of interest by lowering the
density threshold around the target density peak. Therefore, the regions A and
B are specified by the density threshold given by the saddle point between
those nearby density peaks. Particles in the core regions are named core
particles and they are assumed to be members of the subhalo candidate. From
now on, we mix the use of A, B and C, to denote core regions or corresponding
subhalo candidates depending on the meaning in the text.

Particles in the shell region of $\boldsymbol{\alpha}$ are supposed to be
members of one of two halo candidates, A and B. We calculate the tidal radius
of the less massive subhalo and check whether a particle in the
$\boldsymbol{\alpha}$ shell region lies within the tidal radius and its total
energy with respect to the subhalo is negative. If not, we move on to the next
more massive subhalo and judge its membership in the same manner. At last,
particles in the $\boldsymbol{\beta}$ region are picked up to check membership
with respect to the subhalos, A, B, and C.

Whenever we move our focus to the next shell particles, we update the tidal
radius of subhalos and re-apply the membership determination for member
particles to confirm whether they could still hold their membership under the
updated circumstances.

Figure 23: Two-dimensional example of the density field and iso-density
contours which split and arrange areas in the hierarchical manner. The density
field is made using a virialized FoF halos using its particles.

####  A.4.4 Stellar Density Field in pGalF

Like the PSB method, the pGalF is based on the coordinate-free density map and
a network of particle neighbors. In pGalF, however, the galaxy finding is
based on the stellar mass density, instead of the dark-matter mass density
used in the analysis for pure $N$-body simulations. The density kernel, W4, is
also used for a fixed number of nearest stellar particles. However, a
neighboring link can be made between different matter species in pGalF. Except
for this, the overall scheme is similar to the PSB method.

The stellar density-based approach has several advantages. First, stellar
distributions are more compact than those of dark matter, which helps us
clearly detect the boundary of a galaxy system. Also, stellar components of
galaxies tend to be at the near bottom of potential well, and this helps to
easily add member particles starting from the core region. Also the
substantial parts of the stellar component are not bound to galaxies (see the
panels in the middle column of Fig. 24) due to the dynamical frictions and
tidal stripping in a group or cluster region.

Figure 24: Galaxy finding results applied to the FoF halos as given in Fig.
21. Contours are shown after fitting the density distribution of each
component of galaxies.

## References

  * Adams et al. (2021) Adams, N. J., Bowler, R. A. A., Jarvis, M. J., Häußler, B., & Lagos, C. D. P. 2021, MNRAS, 506, 4933
  * Baldry et al. (2008) Baldry, I. K., Glazebrook, K., & Driver, S. P. 2008, MNRAS, 388, 945
  * Behroozi et al. (2013) Behroozi, P. S., Wechsler, R. H., & Wu, H.-Y. 2013, ApJ, 762, 109
  * Bell & de Jong (2001) Bell, E. F., & de Jong, R. S. 2001, ApJ, 550, 212
  * Bell et al. (2003) Bell, E. F., McIntosh, D. H., Katz, N., & Weinberg, M. D. 2003, ApJS, 149, 289
  * Benson et al. (2003) Benson, A. J., Bower, R. G., Frenk, C. S., et al. 2003, ApJ, 599, 38
  * Binney & Tremaine (1987) Binney, J., & Tremaine, S. 1987, Galactic dynamics
  * Binney & Tremaine (2008) —. 2008, Galactic Dynamics: Second Edition
  * Blanton et al. (2005) Blanton, M. R., Lupton, R. H., Schlegel, D. J., et al. 2005, The Astrophysical Journal, 631, 208
  * Boco et al. (2021) Boco, L., Lapi, A., Chruslinska, M., et al. 2021, The Astrophysical Journal, 907, 110
  * Bolzonella et al. (2010) Bolzonella, M., Kovač, K., Pozzetti, L., et al. 2010, A&A, 524, A76
  * Booth & Schaye (2009) Booth, C. M., & Schaye, J. 2009, MNRAS, 398, 53
  * Borlaff et al. (2019) Borlaff, A., Trujillo, I., Román, J., et al. 2019, A&A, 621, A133
  * Boylan-Kolchin et al. (2008) Boylan-Kolchin, M., Ma, C.-P., & Quataert, E. 2008, MNRAS, 383, 93
  * Brough et al. (2020) Brough, S., Collins, C., Demarco, R., et al. 2020, arXiv e-prints, arXiv:2001.11067
  * Chabrier (2003) Chabrier, G. 2003, PASP, 115, 763
  * Comerford & Natarajan (2007) Comerford, J. M., & Natarajan, P. 2007, MNRAS, 379, 190
  * Conroy (2013) Conroy, C. 2013, ARA&A, 51, 393
  * Conroy et al. (2007) Conroy, C., Wechsler, R. H., & Kravtsov, A. V. 2007, ApJ, 668, 826
  * Crain et al. (2015) Crain, R. A., Schaye, J., Bower, R. G., et al. 2015, MNRAS, 450, 1937
  * Cross et al. (2001) Cross, N., Driver, S. P., Couch, W., et al. 2001, MNRAS, 324, 825
  * Dabringhausen & Fellhauer (2016) Dabringhausen, J., & Fellhauer, M. 2016, MNRAS, 460, 4492
  * Dahlen et al. (2005) Dahlen, T., Mobasher, B., Somerville, R. S., et al. 2005, The Astrophysical Journal, 631, 126
  * Danieli et al. (2020) Danieli, S., Lokhorst, D., Zhang, J., et al. 2020, The Astrophysical Journal, 894, 119
  * Davidzon et al. (2018) Davidzon, I., Ilbert, O., Faisst, A. L., Sparre, M., & Capak, P. L. 2018, ApJ, 852, 107
  * Davidzon et al. (2017) Davidzon, I., Ilbert, O., Laigle, C., et al. 2017, A&A, 605, A70
  * Davies et al. (2016) Davies, J. I., Davies, L. J. M., & Keenan, O. C. 2016, MNRAS, 456, 1607
  * Donnari et al. (2019) Donnari, M., Pillepich, A., Nelson, D., et al. 2019, MNRAS, 485, 4817
  * dos Reis et al. (2020) dos Reis, S. N., Buitrago, F., Papaderos, P., et al. 2020, A&A, 634, A11
  * Driver (1999) Driver, S. P. 1999, ApJ, 526, L69
  * Dubois et al. (2012) Dubois, Y., Devriendt, J., Slyz, A., & Teyssier, R. 2012, MNRAS, 420, 2662
  * Dubois et al. (2016) Dubois, Y., Peirani, S., Pichon, C., et al. 2016, MNRAS, 463, 3948
  * Dubois et al. (2014a) Dubois, Y., Volonteri, M., & Silk, J. 2014a, MNRAS, 440, 1590
  * Dubois et al. (2014b) Dubois, Y., Pichon, C., Welker, C., et al. 2014b, MNRAS, 444, 1453
  * Dubois et al. (2021) Dubois, Y., Beckmann, R., Bournaud, F., et al. 2021, A&A, 651, A109
  * Duncan et al. (2014) Duncan, K., Conselice, C. J., Mortlock, A., et al. 2014, MNRAS, 444, 2960
  * Durkalec et al. (2015) Durkalec, A., Le Fèvre, O., de la Torre, S., et al. 2015, A&A, 576, L7
  * Euclid Collaboration et al. (2022) Euclid Collaboration, Borlaff, A. S., Gómez-Alvarez, P., et al. 2022, A&A, 657, A92
  * Faber et al. (2007) Faber, S. M., Willmer, C. N. A., Wolf, C., et al. 2007, ApJ, 665, 265
  * Fernández Lorenzo et al. (2013) Fernández Lorenzo, M., Sulentic, J., Verdes-Montenegro, L., & Argudo-Fernández, M. 2013, MNRAS, 434, 325
  * Foucaud et al. (2010) Foucaud, S., Conselice, C. J., Hartley, W. G., et al. 2010, MNRAS, 406, 147
  * Furlong et al. (2015) Furlong, M., Bower, R. G., Theuns, T., et al. 2015, MNRAS, 450, 4486
  * Gabasch et al. (2006) Gabasch, A., Hopp, U., Feulner, G., et al. 2006, A&A, 448, 101
  * Ganeshaiah Veena et al. (2018) Ganeshaiah Veena, P., Cautun, M., van de Weygaert, R., et al. 2018, Monthly Notices of the Royal Astronomical Society, 481, 414
  * Geller et al. (2012) Geller, M. J., Diaferio, A., Kurtz, M. J., Dell’Antonio, I. P., & Fabricant, D. G. 2012, AJ, 143, 102
  * Golden-Marx et al. (2022) Golden-Marx, J. B., Miller, C. J., Zhang, Y., et al. 2022, The Astrophysical Journal, 928, 28
  * Gonzalez et al. (2007) Gonzalez, A. H., Zaritsky, D., & Zabludoff, A. I. 2007, ApJ, 666, 147
  * González et al. (2011) González, V., Labbé, I., Bouwens, R. J., et al. 2011, ApJ, 735, L34
  * Graham et al. (2005) Graham, A. W., Driver, S. P., Petrosian, V., et al. 2005, The Astronomical Journal, 130, 1535
  * Grazian et al. (2015) Grazian, A., Fontana, A., Santini, P., et al. 2015, A&A, 575, A96
  * Greene et al. (2022) Greene, J. E., Greco, J. P., Goulding, A. D., et al. 2022, arXiv e-prints, arXiv:2204.11883
  * Guennou et al. (2012) Guennou, L., Adami, C., Da Rocha, C., et al. 2012, A&A, 537, A64
  * Hopkins (2004) Hopkins, A. M. 2004, ApJ, 615, 209
  * Ichikawa et al. (2012) Ichikawa, T., Kajisawa, M., & Akhlaghi, M. 2012, MNRAS, 422, 1014
  * Ilbert et al. (2005) Ilbert, O., Tresse, L., Zucca, E., et al. 2005, A&A, 439, 863
  * Jackson et al. (2021) Jackson, R. A., Martin, G., Kaviraj, S., et al. 2021, MNRAS, 502, 4262
  * Jenkins et al. (2001) Jenkins, A., Frenk, C. S., White, S. D. M., et al. 2001, Monthly Notices of the Royal Astronomical Society, 321, 372
  * Katsianis et al. (2015) Katsianis, A., Tescari, E., & Wyithe, J. S. B. 2015, MNRAS, 448, 3001
  * Kaviraj et al. (2017) Kaviraj, S., Laigle, C., Kimm, T., et al. 2017, Monthly Notices of the Royal Astronomical Society, 467, 4739
  * Kaviraj et al. (2017) Kaviraj, S., Laigle, C., Kimm, T., et al. 2017, MNRAS, 467, 4739
  * Kawinwanichakij et al. (2020) Kawinwanichakij, L., Papovich, C., Ciardullo, R., et al. 2020, ApJ, 892, 7
  * Kim & Park (2006) Kim, J., & Park, C. 2006, ApJ, 639, 600
  * Kim & Park (2006) Kim, J., & Park, C. 2006, The Astrophysical Journal, 639, 600
  * Kim et al. (2009) Kim, J., Park, C., Gott, J. Richard, I., & Dubinski, J. 2009, ApJ, 701, 1547
  * Kim et al. (2015) Kim, J., Park, C., L’Huillier, B., & Hong, S. E. 2015, Journal of Korean Astronomical Society, 48, 213
  * King (1962) King, I. 1962, AJ, 67, 471
  * Kravtsov et al. (2018) Kravtsov, A. V., Vikhlinin, A. A., & Meshcheryakov, A. V. 2018, Astronomy Letters, 44, 8
  * Lee et al. (2018) Lee, J., Kim, S., Jeong, H., et al. 2018, ApJ, 864, 69
  * Lee et al. (2021) Lee, J., Shin, J., Snaith, O. N., et al. 2021, ApJ, 908, 11
  * Lee et al. (2012) Lee, K.-S., Ferguson, H. C., Wiklind, T., et al. 2012, The Astrophysical Journal, 752, 66
  * Legrand et al. (2019) Legrand, L., McCracken, H. J., Davidzon, I., et al. 2019, MNRAS, 486, 5468
  * Leja et al. (2020) Leja, J., Speagle, J. S., Johnson, B. D., et al. 2020, ApJ, 893, 111
  * Łokas & Mamon (2001) Łokas, E. L., & Mamon, G. A. 2001, MNRAS, 321, 155
  * Lower et al. (2020) Lower, S., Narayanan, D., Leja, J., et al. 2020, ApJ, 904, 33
  * Madau & Dickinson (2014) Madau, P., & Dickinson, M. 2014, ARA&A, 52, 415
  * Madau & Dickinson (2014) Madau, P., & Dickinson, M. 2014, Annual Review of Astronomy and Astrophysics, 52, 415
  * Mahajan et al. (2018) Mahajan, S., Drinkwater, M. J., Driver, S., et al. 2018, MNRAS, 475, 788
  * Marchesini et al. (2009) Marchesini, D., van Dokkum, P. G., Schreiber, N. M. F., et al. 2009, ApJ, 701, 1765
  * Marinacci et al. (2018) Marinacci, F., Vogelsberger, M., Pakmor, R., et al. 2018, MNRAS, 480, 5113
  * Martin et al. (2019) Martin, G., Kaviraj, S., Laigle, C., et al. 2019, MNRAS, 485, 796
  * McCarthy et al. (2017) McCarthy, I. G., Schaye, J., Bird, S., & Le Brun, A. M. C. 2017, MNRAS, 465, 2936
  * McKinney et al. (2012) McKinney, J. C., Tchekhovskoy, A., & Blandford, R. D. 2012, MNRAS, 423, 3083
  * McLeod et al. (2021) McLeod, D. J., McLure, R. J., Dunlop, J. S., et al. 2021, MNRAS, 503, 4413
  * McLeod et al. (2021) McLeod, D. J., McLure, R. J., Dunlop, J. S., et al. 2021, MNRAS, 503, 4413
  * Merloni & Heinz (2008) Merloni, A., & Heinz, S. 2008, MNRAS, 388, 1011
  * Mihos (2019) Mihos, J. C. 2019, arXiv e-prints, arXiv:1909.09456
  * Montes (2022) Montes, M. 2022, Nature Astronomy, 6, 308
  * Naiman et al. (2018) Naiman, J. P., Pillepich, A., Springel, V., et al. 2018, MNRAS, 477, 1206
  * Navarro et al. (1997) Navarro, J. F., Frenk, C. S., & White, S. D. M. 1997, ApJ, 490, 493
  * Nedkova et al. (2021) Nedkova, K. V., Häußler, B., Marchesini, D., et al. 2021, MNRAS, 506, 928
  * Nelson et al. (2018) Nelson, D., Pillepich, A., Springel, V., et al. 2018, MNRAS, 475, 624
  * O’Mill et al. (2011) O’Mill, A. L., Duplancic, F., García Lambas, D., & Sodré, L. 2011, Monthly Notices of the Royal Astronomical Society, 413, 1395
  * Panter et al. (2007) Panter, B., Jimenez, R., Heavens, A. F., & Charlot, S. 2007, MNRAS, 378, 1550
  * Pillepich et al. (2018) Pillepich, A., Nelson, D., Hernquist, L., et al. 2018, MNRAS, 475, 648
  * Press & Schechter (1974) Press, W. H., & Schechter, P. 1974, ApJ, 187, 425
  * Puchwein & Springel (2013) Puchwein, E., & Springel, V. 2013, MNRAS, 428, 2966
  * Renaud et al. (2016) Renaud, F., Agertz, O., & Gieles, M. 2016, Monthly Notices of the Royal Astronomical Society, 465, 3622
  * Ricciardelli et al. (2012) Ricciardelli, E., Vazdekis, A., Cenarro, A. J., & Falcón-Barroso, J. 2012, MNRAS, 424, 172
  * Román et al. (2021) Román, J., Castilla, A., & Pascual-Granado, J. 2021, A&A, 656, A44
  * Roy et al. (2018) Roy, N., Napolitano, N. R., La Barbera, F., et al. 2018, MNRAS, 480, 1057
  * Samuel et al. (2022) Samuel, J., Wetzel, A., Santistevan, I., et al. 2022, MNRAS, 514, 5276
  * Schaye et al. (2015) Schaye, J., Crain, R. A., Bower, R. G., et al. 2015, MNRAS, 446, 521
  * Schmidt (1959) Schmidt, M. 1959, ApJ, 129, 243
  * Sheth & Tormen (1999) Sheth, R. K., & Tormen, G. 1999, MNRAS, 308, 119
  * Shuntov et al. (2022) Shuntov, M., McCracken, H. J., Gavazzi, R., et al. 2022, arXiv e-prints, arXiv:2203.10895
  * Silk & Mamon (2012) Silk, J., & Mamon, G. A. 2012, Research in Astronomy and Astrophysics, 12, 917
  * Sohn et al. (2018) Sohn, J., Geller, M. J., Rines, K. J., et al. 2018, ApJ, 856, 172
  * Song et al. (2016) Song, M., Finkelstein, S. L., Ashby, M. L. N., et al. 2016, The Astrophysical Journal, 825, 5
  * Springel et al. (2018) Springel, V., Pakmor, R., Pillepich, A., et al. 2018, MNRAS, 475, 676
  * Stabenau et al. (2008) Stabenau, H. F., Connolly, A., & Jain, B. 2008, MNRAS, 387, 1215
  * Stefanon et al. (2021a) Stefanon, M., Bouwens, R. J., Labbé, I., et al. 2021a, ApJ, 922, 29
  * Stefanon et al. (2021b) —. 2021b, ApJ, 922, 29
  * Suto et al. (2016) Suto, D., Kitayama, T., Osato, K., Sasaki, S., & Suto, Y. 2016, PASJ, 68, 14
  * Tang et al. (2021) Tang, L., Lin, W., Wang, Y., & Napolitano, N. R. 2021, MNRAS, 508, 3321
  * Thorne et al. (2021) Thorne, J. E., Robotham, A. S. G., Davies, L. J. M., et al. 2021, MNRAS, 505, 540
  * Tomczak et al. (2014) Tomczak, A. R., Quadri, R. F., Tran, K.-V. H., et al. 2014, ApJ, 783, 85
  * Trujillo et al. (2020) Trujillo, I., Chamba, N., & Knapen, J. H. 2020, Monthly Notices of the Royal Astronomical Society, 493, 87
  * Valls-Gabaud & MESSIER Collaboration (2017) Valls-Gabaud, D., & MESSIER Collaboration. 2017, in Formation and Evolution of Galaxy Outskirts, ed. A. Gil de Paz, J. H. Knapen, & J. C. Lee, Vol. 321, 199–201
  * van de Ven et al. (2003) van de Ven, G., van Dokkum, P. G., & Franx, M. 2003, MNRAS, 344, 924
  * van den Bosch et al. (2018) van den Bosch, F. C., Ogiya, G., Hahn, O., & Burkert, A. 2018, MNRAS, 474, 3043
  * van der Burg et al. (2017) van der Burg, R. F. J., Hoekstra, H., Muzzin, A., et al. 2017, A&A, 607, A79
  * van der Wel et al. (2005) van der Wel, A., Franx, M., van Dokkum, P. G., et al. 2005, ApJ, 631, 145
  * Vazdekis et al. (2016) Vazdekis, A., Koleva, M., Ricciardelli, E., Röck, B., & Falcón-Barroso, J. 2016, MNRAS, 463, 3409
  * Vazdekis et al. (2012) Vazdekis, A., Ricciardelli, E., Cenarro, A. J., et al. 2012, MNRAS, 424, 157
  * Vogelsberger et al. (2013) Vogelsberger, M., Genel, S., Sijacki, D., et al. 2013, MNRAS, 436, 3031
  * Vogelsberger et al. (2014) Vogelsberger, M., Genel, S., Springel, V., et al. 2014, Monthly Notices of the Royal Astronomical Society, 444, 1518
  * Weaver et al. (2022) Weaver, J. R., Davidzon, I., Toft, S., et al. 2022, arXiv e-prints, arXiv:2212.02512
  * Wechsler & Tinker (2018) Wechsler, R. H., & Tinker, J. L. 2018, ARA&A, 56, 435
  * Weigel et al. (2016) Weigel, A. K., Schawinski, K., & Bruderer, C. 2016, MNRAS, 459, 2150
  * Williams et al. (2016) Williams, R. P., Baldry, I. K., Kelvin, L. S., et al. 2016, Monthly Notices of the Royal Astronomical Society, 463, 2746
  * Wright et al. (2018) Wright, A. H., Driver, S. P., & Robotham, A. S. G. 2018, MNRAS, 480, 3491
  * Yoo et al. (2021) Yoo, J., Ko, J., Kim, J.-W., & Kim, H. 2021, MNRAS, 508, 2634
  * Zwicky (1951) Zwicky, F. 1951, PASP, 63, 61