# Low-surface-brightness Galaxies are Missing in the Observed Stellar Mass Function

Juhan Kim, Jaehyun Lee, Clotilde Laigle, Yohan Dubois, Yonghwi Kim, Changbom Park, Christophe Pichon, Brad K. Gibson, C. Gareth Few, Jihye Shin, et al.

► To cite this version:

Juhan Kim, Jaehyun Lee, Clotilde Laigle, Yohan Dubois, Yonghwi Kim, et al. Low-surface-brightness Galaxies are Missing in the Observed Stellar Mass Function. The Astrophysical Journal, 2023, 951, 10.3847/1538-4357/acd251. insu-04234588

HAL Id: insu-04234588

https://insu.hal.science/insu-04234588

Submitted on 13 Oct 2023

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
<figure><img src="image_2.png" /></figure>

# Low-surface-brightness Galaxies are Missing in the Observed Stellar Mass Function

Juhan Kim<fnref n="1" />, Jaehyun Lee<fnref n="2" />, Clotilde Laigle<fnref n="3" />, Yohan Dubois<fnref n="3" />, Yonghwi Kim<fnref n="2" />,<fnref n="4" />, Changbom Park<fnref n="2" />, Christophe Pichon<fnref n="2" />,<fnref n="3" />,<fnref n="5" />, Brad K. Gibson<fnref n="6" />, C. Gareth Few<fnref n="6" />, Jihye Shin<fnref n="7" />, and Owain Snaith<fnref n="8" />

Received 2022 December 29; revised 2023 April 27; accepted 2023 May 2; published 2023 July 11

## Abstract

We investigate the impact of the surface-brightness (SB) limit on the galaxy stellar mass functions (GSMFs) using galaxy catalogs 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 observational 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 the low-mass end. This approach can ease the discrepancy between the GSMFs obtained from simulations and observations at $0.6 \lesssim z \lesssim 2$. We also find that a redshift survey with an SB selection limit of $\langle \mu_r \rangle^e = 25 \ \text{mag arcsec}^{-2}$ will miss 20% of galaxies with $M_*^g = 10^9 M_\odot$ at $z = 0.625$. The missing fraction of low-surface-brightness galaxies increases to 35%, 55%, and 80% at $z = 0.9, 1.1$, and 1.9, respectively, at the same SB limit.

Unified Astronomy Thesaurus concepts: Hydrodynamical simulations (767)

## 1. Introduction

The observed galaxy stellar mass function (GSMF; for a noncomplete 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 (Conroy et al. 2007; Puchwein & Springel 2013; Madau & Dickinson 2014; Grazian et al. 2015; Katsianis et al. 2015; Adams et al. 2021; Boco et al. 2021). 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; Panter et al. 2007; Bolzonella et al. 2010; Conroy 2013) where the galaxy history of star formation activities (Madau & Dickinson 2014; Lower et al. 2020) is 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 and halo relation, see Wechsler & Tinker 2018). The GSMFs are accordingly expected to follow the halo-mass functions (HMFs), but their shapes are quite different to each other. The HMF rises more steeply in the low-mass tail (see Kim et al. 2015 for a noncomplete list of various fitting functions) while the GSMF has a much shallower (power law) tail below the knee of the Schechter function (Song et al. 2016; Weigel et al. 2016; Davidzon et al. 2017; Adams et al. 2021; McLeod et al.

2021). On the low-mass scale, this discrepancy is mainly attributed to stellar feedback (to name a few, Benson et al. 2003; Baldry et al. 2008; Silk & Mamon 2012; Puchwein & Springel 2013). Supernova (SN) explosions and stellar winds heat and eject the interstellar medium from galaxies in low-mass halos, accordingly suppressing the growth of galaxies (e.g., Silk & Mamon 2012).

Motivated by this, some efforts have been made to reproduce the observed GSMFs by implementing stellar feedback into gravitohydrodynamics simulations (Vogelsberger et al. 2013; Dubois et al. 2014a, 2016; Katsianis et al. 2015; Schaye et al. 2015; Pillepich et al. 2018). The IllustrisTNG 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 Marinacci et al. 2018; Naiman et al. 2018; Nelson et al. 2018; Pillepich et al. 2018; Springel et al. 2018). The EAGLE simulation focused on stochastic thermal feedback to prevent the overcooling problem, without turning off radiative cooling (for details, see Crain et al. 2015; Schaye 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. 2014a) 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; Valls-Gabaud & MESSIER Collaboration 2017; Martin et al. 2019; Tang et al. 2021; Greene et al. 2022). A well-known example of this issue in low-redshift surveys (typically with a surface brightness of $\langle \mu_R \rangle^e \ge 21.82 \ \text{mag arcsec}^{-2}$) is given in Figure 10 of Geller et al. (2012). It shows that the LSBGs dominate the

<figure><img src="image_3.png" /></figure>

Original content from this work may be used under the terms of the Creative Commons Attribution 4.0 licence. Any further distribution of this work must maintain attribution to the author(s) and the title of the work, journal citation and DOI.

<footnote n="1">¹Center for Advanced Computation, Korea Institute for Advanced Study, 85 Hoegiro, Dongdaemun-gu, Seoul 02455, Republic of Korea</footnote>

<footnote n="2">²Korea Institute for Advanced Study, 85 Hoegiro, Dongdaemun-gu, Seoul 02455, Republic of Korea; syncphy@gmail.com</footnote>

<footnote n="3">³CNRS and Sorbonne Université, UMR 7095, Institut d'Astrophysique de Paris, 98 bis, Boulevard Arago, F-75014 Paris, France</footnote>

<footnote n="4">⁴Department of Astronomy, Yonsei University, 50 Yonsei-ro, Seodaemun-gu, Seoul 03722, Republic of Korea</footnote>

<footnote n="5">⁵IPhT, DRF-INP, UMR 3680, CEA, L'Orme des Merisiers, Bât 774, F-91191 Gif-sur-Yvette, France</footnote>

<footnote n="6">⁶E.A. Milne Centre for Astrophysics, University of Hull, Hull, HU6 7RX, UK</footnote>

<footnote n="7">⁷Korea Astronomy and Space Science Institute, 776 Daedeokdae-ro, Yuseong-gu, Daejeon 34055, Republic of Korea</footnote>

<footnote n="8">⁸University of Exeter, School of Physics and Astronomy, Stocker Road, Exeter, EX4 4QL, UK</footnote>
faint-end regime of the galaxy luminosity function (LFs). Therefore, a somewhat shallow surface brightness (SB) cut probably misses the majority of LSBGs at the low-brightness end of the galaxy LFs.

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 Hubble Space Telescope (HST) Cosmic Assembly Near-infrared Deep Extragalactic Legacy Survey (CANDELS) archive image data with the SB limits of 28.3 mag arcsec⁻² ≲ ⟨μ<sub>F160W</sub>⟩ ≲ 28.6 mag arcsec⁻² corresponding to the 3σ limiting magnitude in the Cosmic Evolution Survey (COSMOS), Great Observatories Origins Deep Survey-North (GOODS-N), UltraDeep Survey (UDS), and Extended Groth Strip fields for an aperture size of 4" (main) and 1" (small galaxies). Grazian et al. (2015) adopted 27 mag arcsec⁻² ≲ ⟨μ<sub>F160W</sub>⟩ ≲ 28.5 mag arcsec⁻² at the 1σ limit in GOODS-S for galaxies at 3.5 ≲ z ≲ 7.5 with an aperture size 2 × the FWHM. Adams et al. (2021) derived the GSMFs from the Hyper Suprime-Cam (HSC) data of the COSMOS and XMM-LSS fields with an SB limit of 25.9 mag arcsec⁻² ≲ ⟨μ<sub>r</sub>⟩⁰ ≲ 26.7 mag arcsec⁻² at 5σ with an aperture of 2" 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σ depths of 24.6 < ⟨μ<sub>K<sub>i</sub></sub>⟩ < 25.2 in a circular diameter of 0".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 Sloan Digital Sky Survey (SDSS) has the SB selection criterion of ⟨μ<sub>r</sub>⟩<sup>e</sup> ≲ 23 mag arcsec<sup>-2</sup> (Blanton et al. 2005). The SB selection criteria are ⟨μ<sub>r</sub>⟩ ≲ 23 mag arcsec<sup>-2</sup> for the LCRS (Cross et al. 2001), ⟨μ<sub>bJ</sub>⟩<sup>e</sup> ≲ 24.7 mag arcsec<sup>-2</sup> for 2dFGRS and APM (Cross et al. 2001), ⟨μ<sub>r</sub>⟩<sup>petro</sup> ≲ 21.3 mag arcsec<sup>-2</sup> for HectoMAP (Sohn et al. 2018), and ⟨μ<sub>R</sub>⟩<sup>e</sup> ≲ 21.82 mag arcsec<sup>-2</sup> for SHELS<sub>0.1</sub> (Geller et al. 2012). Wide-field surveys enable us to obtain accurate GSMFs in the local universe, but they are still incomplete for detecting LSBGs due to their shallow detection limits. Given that the GSMF is the key statistics in calibrating physical ingredients of galaxy formation models, it necessitates 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, Lee et al. 2021; Park et al. 2022) 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 the low-redshift GSMFs are fully addressed in Sections 4 and 5, respectively. We conclude this study with discussions in Section 6. Additionally we describe the galaxy finding method used for HR5 in Appendix A. We also provide auxiliary calculations on the magnitude transform in Appendix B, the source detection criteria in Appendix C, the dependence of GSMFs on the star formation efficiency in Appendix D, and the definition of the effective radius in Appendix E.

## 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_{\text{box}} = 1049$ cMpc. A zoomed region is set to have a cuboid geometry of $L_{(x,y,z)}^{\text{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 is given in Appendix A.

### 2.2. Star Formation and Feedback

Here, we summarize the subgrid physics associated with star formation and stellar feedback adopted in HR5. The star formation rates (SFRs) are computed based on the Schmidt law (Schmidt 1959):

$$ \frac{d\rho_{*}}{dt} = \epsilon_{*} \frac{\rho_{\text{gas}}}{t_{\text{ff}}}, \qquad (1) $$

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

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

$$ \chi \equiv \left( \frac{\dot{M}_{\text{BH}}}{\dot{M}_{\text{Edd}}} \right), \qquad (2) $$

where $\dot{M}_{\text{BH}}$ is the growth rate of black hole mass and $\dot{M}_{\text{Edd}}$ is the Eddington limit. If $\chi \le \chi_{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}_{\text{BH,h}} = \epsilon_{r} \epsilon_{f,h} \dot{M}_{\text{BHL}} c^2$, where $\epsilon_{r}$ is the spin-dependent radiative efficiency, $\epsilon_{f,h}$ is the coupling efficiency of the thermal feedback, $\dot{M}_{\text{BHL}}$ is the Bondi–Hoyle–Lyttleton accretion rate, and c is the speed of light. A jet mode releases an amount of energy of $\dot{E}_{\text{BH,j}} = \epsilon_{f,j} \dot{M}_{\text{BHL}} c^2$, where $\epsilon_{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_{c,jet}$ and $\epsilon_{f,h}$. We adopt $\epsilon_{f,h} = 0.15$ to reproduce the observed $M_{\text{BH}}-M_{*}$ relation at z = 0 (Booth & Schaye 2009; Dubois et al. 2012) and $\chi_{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. (2014b, 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 and 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 the 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 (Ricciardelli et al. 2012; Vazdekis et al. 2012, 2016). 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 brightest cluster galaxy (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 an overestimated stellar mass due to its far-extended stellar envelope tangled with intracluster 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 galaxies of total stellar masses of $M \sim 10^{12} M_{\odot}$ have half their masses 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 3D 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 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 an increase in the mass-to-light ratios of galaxies (van de Ven et al. 2003; van der Wel et al. 2005).

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

$$ \log_{10} \left( \frac{M_{\star}^g}{M_{\odot}} \right) = 1.097(g - r) - 0.406 - 0.4(M_r - 4.67) + f_{\text{Bell}}(z), \quad (3) $$

where $f_{\text{Bell}}(z) = -0.19z$ reflects the redshift evolution in Two Micron All Sky Survey and 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_{\text{Bell}}$, we propose a new fitting model of

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

which is formulated to capture the coevolution of the redshift and galaxy color at $0.625 \le z \le 5$.

<figure><img src="image_4.png" /><figcaption>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 Galaxy And Mass Assembly survey (hatched area; Mahajan et al. 2018) with our fitting model simply extrapolated to $z = 0$ (colored region). (top): Fitting residuals. The filled circle is the median in each stellar mass bin, and the colored box is the 1σ distribution. The error bar marks the $2\sigma$ scattering.</figcaption></figure>

Figure 1 shows the evolution of the mass-brightness relation at a given galaxy stellar mass as a function of the redshift. In this plot, the r-band magnitude $M_r$ is given at fixed $M_{\star}^g$ because the 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 the aging of the stellar populations. The aging effect is taken into consideration by the correction term ($f_{\text{Bell}}$ or $f_{\text{HR5}}$). In the top panel of the figure, we show the residuals of the fit. Here the residual absolute magnitude ($R$) is the difference in the absolute magnitude between the simulated galaxy and the modeled one and is written as

$$ R(M_{\star}^g) \equiv M_r(\text{sim}) - M_r(\text{model}), \quad (5) $$

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

## 3.3. Environmental Dependence of Galaxy Brightness

Now, we present the environmental effects on the galaxy absolute magnitude in the SDSS r band. From the simulation catalog of HR5, we divide the entire galaxy sample ($G_{\text{all}}$) into the satellite ($G_{\text{sat}}$), central ($G_{\text{cen}}$), and isolated ($G_{\text{iso}}$) samples. The most massive galaxy in each halo is classified as the central, while the satellite sample contains the rest, less massive galaxies. Isolated galaxies are those with no companion in their friend-of-friend (FoF) halos.
<figure><img src="image_5.png" /><figcaption>**Figure 2.** Relation between the environment and the galaxy brightness as a function of the galaxy stellar mass. The differences in the average absolute magnitudes between the $G_{cen}$ and $G_{all}$, $G_{sat}$ and $G_{cen}$, $G_{cen}$ and $G_{iso}$, and $G_{sat}$ and $G_{iso}$ samples are shown clockwise starting from the bottom right panel.</figcaption></figure>

In the bottom right panel of Figure 2, $G_{cen}$ tends to be brighter than the average obtained from the $G_{all}$ sample, while this is less significant for more massive galaxies. The redshift evolution of $\Delta M_r$ is apparent at the low-mass end due to the decline in the 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 central galaxies of the same mass, and this tendency is more pronounced at the low-mass end.

The two top panels demonstrate that isolated galaxies ($G_{iso}$) are fainter than both 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. 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 galaxies (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).

### 3.4. Galaxy Luminosity Functions

Ilbert et al. (2005) and Dahlen et al. (2005) derived galaxy LFs at $z \le 2$ from the VISIBLE MultiObject Spectrograph-Very Large Telescope (VLT) Deep Survey, HST Advanced Camera for Survey (ACS), and GOODS-S in the rest-frame R band of the Cousins filter system. While the R band has similar wavelength coverage to the SDSS r band, we have applied the

$\chi^2$ fit to obtain the r-band luminosity functions (for the fitting form; see Appendix B). We also add the LFs obtained from the FOcal Reducer and low dispersion Spectrograph 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_{Johnson}$ band, we convert the SDSS r-band magnitudes of the HR5 galaxies to the $B_{Johnson}$-band magnitudes by adopting the band transforms of $B_{Johnson} = g + 0.115 + 0.370 \times 0.82$ and calculate $M_r^*$ 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 by thick black solid lines. As mentioned above, we apply the 30 pkpc aperture cut that is substantially expected to reduce the total brightness of the BCGs. In spite of the aperture cut, we still have substantial excess of the BCG-scale galaxies in HR5 at $z \le 1.13$. Another notable feature seen in Figure 3 is that the excess of the simulated LFs to the observations becomes larger with increasing $M_r$ below the knee of the Schechter functions at all 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 and luminosity scales will be quantitatively addressed in the following sections.

### 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 are presented in Figure 5. At $z=3$ and $5$, the relations obtained from HR5 well follow the observations, while at $z \sim 1.3$ and $0.6$, HR5 overproduces the stellar mass on the BCG scale ($M \ge 3 \times 10^{11} M_\odot$; see Bellstedt et al. 2016 for reference), compared to Behroozi et al. 2013 (gray 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, which are based on the local SDSS BCG catalogs ($0 \le z \le 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. 2014a), and TNG100 (Pillepich et al. 2018). At all the redshifts investigated in this study, EAGLE and TNG100 show SHM distributions that agree well with the empirical data, particularly in $M_r^* \lesssim 10^{11} M_\odot$. In H-AGN, galaxies are indeed far too massive at a fixed halo mass compared to other simulations and
<figure><img src="image_6.png" /><figcaption>**Figure 3.** Galaxy luminosity functions at $z = 0.625$, 1.13, 1.89, and 3 (from bottom to top). The thick solid lines are the mock galaxy luminosity functions from HR5. The hatched or shaded areas are the observed Schechter functions (Dahlen et al. 2005; Ilbert et al. 2005; Gabasch et al. 2006; Faber et al. 2007).</figcaption></figure>

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 the high-halo-mass end at which AGN feedback dominantly regulates star formation; the slopes of the SHM relation are 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 the high-mass end is to some degree relieved. However, the slope still does not seem to be consistent with the empirical results.

### 3.6. Effective Sizes of Galaxies

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 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

<figure><img src="image_7.png" /><figcaption>**Figure 4.** Relations between $M_h$ and $M_*^g$. The 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; gray 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. (2021; open circles, $z = 6$).</figcaption></figure>

means that the 1 pkpc resolution is probably insufficient to properly simulate early-type galaxies with masses below $M_*^g \sim 10^{11} 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 the galaxy size $R_{1/2}$ as a function of the 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
<figure><img src="image_8.png" /><figcaption>Figure 5. Same as Figure 4 but with the 30 pkpc aperture correction only to HR5.</figcaption></figure>

correlation between the galaxy size and stellar mass becomes stronger (see also Dubois et al. 2016). The scatter in the galaxy-size distribution increases with increasing mass and decreasing redshift.

### 3.7. 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

<figure><img src="image_9.png" /><figcaption>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 innermost contours, we show distributions of 10%, 68% (1σ), 95% (2σ), and 99.7% (3σ) 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).</figcaption></figure>

magnitudes in the frame as (Blanton & Roweis 2007)

$$K \equiv -2.5 \log \left[ \frac{1}{(1+z)} \frac{\int d\lambda_0 \lambda_0 L_{\lambda}(\lambda_0/1+z) R(\lambda_0)}{\int d\lambda_0 \lambda_0 L_{\lambda}(\lambda_0) R(\lambda_0)} \right], \quad (6)$$

where R is the filter response function and L_λ is the rest-frame SED of the galaxy. 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 in the r band is

$$K_r^{\text{fit}} = \{13.3(g-r)_{\text{rest}} - 0.5 - 3.5(z-1.47)^2\} \log(1+z), \quad (7)$$

for a redshift range of 0.625 ≤ z ≤ 2. Figure 8 shows the fitting results at several redshifts. The k-correction difference is defined as ΔK ≡ K_r^{sim} - K_r^{fit}, where K_r^{sim} is directly derived from the simulated SEDs. They seem to have a small nonlinearity 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)_{rest} color space.

It is worthwhile to note that Equation (7) introduces the logarithmic form of the redshift rather than polynomials of z (see O'Mill et al. 2011), to avoid a linear increase in K(z) at high redshifts. Moreover, the k-correction of Equation (7) 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 in galaxy colors (−0.2 ≲ (g - r) ≲ 2), we have −0.6 ≲ K^{fit}/log(1 + z) ≲ 6, which is similar to the observational result derived by Stabenau et al. (2008), who measured the change in the surface brightness of galaxies with various
<figure><img src="image_10.png" /><figcaption>**Figure 7.** Galaxy mass–size relation as a function of the redshift. The distribution of galactic size for each stellar mass bin is shown with a filled circle (median), a colored box (1σ), and an error bar (2σ). 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-N 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).</figcaption></figure>

<figure><img src="image_11.png" /><figcaption>**Figure 9.** Fitting model of the k-correction formulated in Equation (7). We colorized the distribution of the modeled k-correction for each ($g - r_{rest}$) as shown in the right scale bar. Note that the valid redshift range for the fit is 0.625 ≲ $z$ ≲ 2, but we extrapolate the fit down to $z = 0$, in the redshift range marked by gray.</figcaption></figure>

<figure><img src="image_12.png" /><figcaption>**Figure 8.** Fitting uncertainties of the r-band k-correction. The symbols locate the median values; the boxes show the 1σ dispersion while the error bars indicate the 2σ scatters.</figcaption></figure>

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.8. 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 LSBGs are missed due to the SB limit of observations (Driver 1999; Geller et al. 2012; Tang et al. 2021).

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

$$ \theta''_{1/2} \approx 206,265 \frac{r_{1/2}}{d_c(z)}, \qquad (8) $$

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

$$ d_c(z) = \frac{c}{H_0} \int_0^z \frac{dz}{E(z)}, \qquad (9) $$

where $E(z) \equiv \{\Omega_m^0(1+z)^3 + \Omega_k(1+z)^2 + \Omega_\Lambda\}^{1/2}$ in the ΛCDM model. 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 the effective radius. The surface
brightness of a galaxy in units of mag arcsec⁻² is

$$ \langle \mu_{1/2} \rangle = \mathcal{M}_{\odot} - \frac{5}{2} \log \left[ \frac{L_{1/2}/d_c^2(1+z)^2}{L_{\odot}/(10\text{pc})^2} \right] + \frac{5}{2} \log(\pi\theta''_{1/2}), \quad (10) $$

where $\mathcal{M}_{\odot}$ is the absolute magnitude of the Sun, $L_{\odot}$ is the solar luminosity, and $(1+z)^2$ is inserted to reflect the cosmological expansion. The first two terms in the right-hand side of Equation (10) correspond to 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 finite bandwidth, Equation (10) needs to be modified into

$$ \langle \mu_{1/2}^X \rangle \approx \mathcal{M}_{\odot}^X - \frac{5}{2} \log \left( \frac{\mathcal{L}_{1/2}}{\mathcal{R}_{1/2}^2} \right) + \frac{5}{2} \log(1+z)^3 + \mathcal{K}(z; X) + 37.82, \quad (11) $$

where $\mathcal{K}$ is the k-correction term for the redshifted 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\text{kpc}$, where $r_{1/2}$ is the half-light radius in the proper scale. Equation (11) can be simplified as

$$ \langle \mu_{1/2}^X \rangle \approx \mathcal{M}^X + 5 \log \mathcal{R}_{1/2} + \mathcal{C}(z; X) + 38.57, \quad (12) $$

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) + \frac{5}{2} \log(1 + z)^3, \quad (13) $$

which depends only on the redshift, irrespective of the cosmological model parameters (Stabenau et al. 2008). When the cosmological factor is not included, we derive the rest-frame surface brightness as

$$ \langle \mu_{1/2}^{X} \rangle \approx \mathcal{M}^{X} + 5 \log \mathcal{R}_{1/2} + 38.57. \quad (14) $$

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

## 3.9. Distribution of the Surface Brightness of Galaxies

We calculate the mean SB of the simulated galaxies of HR5 using Equation (11) 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 in the brightness or stellar mass outpaces the increase in 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 the cosmological expansion term $\mathcal{C}(z; X)$ in Equation (13). Accordingly, for instance, a galaxy redshift survey with the target SB limit of $\langle \mu_r \rangle = 26 \text{ mag arcsec}^{-2}$ will be able to build a complete catalog of galaxies more massive than about $10^{11} M_{\odot}$ in the redshift range of $0.7 \lesssim z \lesssim 2$. However, if the SB limit is lower than, e.g., $\langle \mu_r \rangle = 24 \text{ mag arcsec}^{-2}$, the galaxy survey would only cover $M_{\star} \gtrsim 10^{12} M_{\odot}$ in the redshift range.

Figure 11 demonstrates the impact of the cosmological expansion term $\mathcal{C}(z; X)$ in Equation (12) on the SB of galaxies. In this case, the SB is only governed by the evolution of

<figure><img src="image_13.png" /><figcaption>Figure 10. Distributions of mean SB of the HR5 galaxies that are assumed to be observed at $z = 0$. The filled circle with error bars marks the median value and the $2\sigma$ (95%) scatter while the colored box ranges show the $1\sigma$ (68%) distribution of the mean surface brightness.</figcaption></figure>

<figure><img src="image_14.png" /><figcaption>Figure 11. Rest frame of the HR5 galaxies obtained by using Equation (12) without the terms depending on the cosmic expansion. The mean SB vs. stellar mass relation of the New Horizon simulation (Jackson et al. 2021) at $z = 0.25$ is included for comparison.</figcaption></figure>

astrophysical parameters: the effective radius and the absolute magnitude. Because of their compact size and younger stellar populations, higher-$z$ galaxies have notably higher SB.

## 3.10. Completeness of GSMF by the SB Limit

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
<figure><img src="image_15.png" /><figcaption>**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. Note that in this plot we do not apply the 30 ckpc aperture correction to the stellar mass.</figcaption></figure>

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 the SB limit on the GSMF at $z = 0.625$. As expected, the SB limit has more impact on low-mass galaxies; more dwarf galaxies are missed with a lower SB limit. The vertical axis ($\mathcal{C}_{\mu_r}$) of the top panel is the survey completeness defined as

$$ \mathcal{C}_{\mu_r} \equiv \frac{N_{\text{obs}}}{N_{\text{all}}}, \qquad (15) $$

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

In Figure 13, we show the completeness of the GSMF in mock observations at $z = 0.875, 1.13$, and $1.89$ (from top to bottom) as a function of the SB limit. For example, we may fail to detect 35% of galaxies with $M_*^g = 10^9 M_\odot$ if the SB limit of the survey is $\langle \mu_r \rangle = 25 \text{ 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 2$ would be complete down to $M = 10^{11} M_\odot$, if a survey adopts the SB condition of $\langle \mu_r \rangle \le 24 \text{ mag arcsec}^{-2}$.

Figure 15 presents an example of the distribution of LSBGs identified around the 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 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 and 11 o’clock) of the image. The

<figure><img src="image_16.png" /><figcaption>**Figure 13.** Completeness of mock observations made with several surface-brightness cuts (values in the legend). From top to bottom, the survey completeness (curves) is shown in a mock survey of galaxies at $z = 0.875, 1.13$, and $1.89$.</figcaption></figure>

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.

## 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, which is controlled by the accretion of gas via mergers or secular evolution.

### 4.1. HR5 versus JWST at High Redshifts

Our universe still has much to reveal and remains open to new discoveries. Recently, JWST has observed galaxies in their infant era beyond $z = 10$, which will give us a hint on the early structure formation and enable us to test the standard model of cosmology and galaxy formation.

Haslbauer et al. (2022) argued that JWST observations may falsify the standard $\Lambda$CDM model due to the relatively massive galaxies found at high redshifts compared to cosmological
<figure><img src="image_17.png" /><figcaption>Figure 14. Redshift evolution of the fraction of observationally selected galaxies for the mean surface-brightness cut of $\langle \mu_r \rangle = 25$ (top), 24 (middle), and 23 mag arcsec$^{-2}$ (bottom panel).</figcaption></figure>

hydrodynamical simulations such as EAGLE and IllustrisTNG. However, the farthest galaxy candidate used in their analysis is now suspected to be a dusty galaxy at lower $z$ (Naidu et al. 2022). Therefore, we do not use the galaxy sample of Haslbauer et al. (2022) in this comparison. In Figure 16, the black downside arrows mark the upper limit of the galaxy number density at the scale of the galaxy stellar mass derived from the JWST observations (Haslbauer et al. 2022). The mass scale marked by the black error bars is lowered to the mass scale marked by the gray error bars if an environment-dependent IMF is adopted. This figure shows that HR5 is consistent with the JWST observations at $z \simeq 10$, depending on the IMF. It is noteworthy that the results are based on the photometric redshift, and thus they should be treated with caution. However, the redshifts of the galaxy sample in the JWST Advanced Deep Extragalactic Survey are obtained from the spectroscopic observations (magenta downside arrows; Curtis-Lake et al. 2023). If the simulation volume is bigger by an order of magnitude, galaxies could properly be simulated on the observed mass scale (see also Keller et al. 2023, for other possible factors). On the other hand, Donnan et al. (2023) confirmed that the galaxy UV luminosity function from JWST does not significantly exceed the $\Lambda$CDM prediction. We may conclude with HR5 that the high-$z$ observations are consistent with the standard hierarchical structure formation model, so far, as illustrated in the figure.

## 4.2. GSMF at $5 \le z \le 10$

Galaxy observational data compiled over the last two decades in many legacy surveys have enabled us to not only trace the global star formation and reionization history in the early universe but also study proto-galaxy formation with 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. 2021) that have been obtained from the GOODS and CANDELS fields observed by the HST and Spitzer IRAC to name a few.

In Figure 17, the simulated GSMFs (bluish 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 or the small star formation efficiency (SFE; see Appendix D) adopted in HR5. Even with the small volume size of the test simulations, the GSMFs at $z \ge 5$ are underestimated by 1.5–3 factors with respect to the simulation with SFE = 4%, which may partially explain the difference shown in the figure.

It may, however, also be attributed to numerical artifacts of adaptive mesh refinement (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 regions where the mass inflow becomes significant. In HR5, this global refinement is carried at the expansion factors of $a_{amr} = 0.0125, 0.25, 0.05, 0.1, 0.2, 0.4$, and $0.8$, in order to maintain a spatial resolution of $\sim 1$ pkpc. As the simulation evolves, a dense gas grid inside a massive galaxy must wait for the next $a_{amr}$ to be refined into eight one-level-finer grids; hence star formation is suppressed until next-level refinement is allowed.

At $z = 5$, the simulated GSMF at the low-mass end is slightly higher than the observations (González et al. 2011; Song et al. 2016) by a few ten percent. This mismatch could possibly be due to the limit to the observational completeness of galaxy sample. Faint-end galaxies usually lie at the observational boundaries, and we assume that faint galaxy samples could suffer from completeness (see Finkelstein et al. 2015, for the simulated completeness of observations). The same argument may apply to the observed GSMFs at $z = 2.5$ and $3.5$ in the next subsection.

## 4.3. GSMF at $2 \le z \le 4$

In order to carry out a robust comparison of GSMFs between the observational results and simulations at intermediate redshifts, we take five different observed GSMFs (Marchesini et al. 2009; Song et al. 2016; Davidzon et al. 2017; Wright et al. 2018; Thorne et al. 2021 among others; for a full list of references, see the legend in Figure 18). Figure 18 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 \le 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.
<figure><img src="image_18.png" /><figcaption>**Figure 15.** Pseudo color map of the star+gas density field of Cluster 1 of HR5 at z = 0.625. The open circles are the LSBGs located within 10 cMpc from the cluster center that have stellar mass of $10^9 M_{\odot} \le M_* \le 10^{10} M_{\odot}$. The white, green, and red circles with $r_{\theta} = 4''$ mark the LSBGs in the SB ranges of [23, 24), [24, 25), and [25, $\infty$) mag arcsec$^{-2}$, respectively. Note that most of the unmarked small faint blobs scattered around the cluster are LSBGs located beyond $10 h^{-1}$ cMpc from the cluster center. The dotted circle has a radius of R = 1 $h^{-1}$ cMpc.</figcaption></figure>

**5. Galaxy Stellar Mass Function at 0.625 $\le z \le 2$**

**5.1. SB Limits of Low-z Observations**

In most galaxy catalogs, the galaxy faint-end magnitude ($m_{\text{ap}}^{\text{lim}}$) is typically provided, which is determined by the background level and aperture size applied to measure the galaxy total flux. However, full information on the pixel-based source detection is unavailable, which is, on the other hand, crucial to determine the SB limit ($\mu_{p}^{\text{lim}}$) of the given galaxy catalog. Due to the limited availability of direct information on SB limits, we must rely on other related information provided in the literature.

Here, we apply some assumptions to derive the SB limit hidden in the galaxy catalog. In Appendix C.1, we show the relation between $m_{\text{ap}}^{\text{lim}}$ and $\mu_{p}^{\text{lim}}$ as a function of the source detection level ($f_{\text{p}}$) above the background noise ($\sigma_{\text{p}}$) on the pixel scale. According to Figure 27, the limited SB of galaxies is about 1–2 mag brighter than $m_{\text{ap}}^{\text{lim}}$ for the same pixel scale ($l_{\text{p}} = 0.1''15$) as of COSMOS (McCracken et al. 2012) between $f_{\text{p}} = 3$ and 5.

**5.2. Seeing Effects on Observational SB Limits**

There is a serious impact of the seeing conditions on the effective SB measurement when the effective radius of a galaxy is comparable to or smaller than the FWHM of the point-spread function (PSF; see Appendix C.2 for angular-size comparisons and Figure 6 of Trujillo et al. 2001, for the seeing effect). The seeing effect depends on the galaxy profile, ellipticity, and specific PSF models. However, the overall observed SB becomes fainter than the true SB of a galaxy, and this difference becomes more severe by $(\Delta\mu)^{\text{seeing}} \gtrsim 0.5$ mag arcsec$^{-2}$ when the FWHM of the PSF is comparable to or larger than the angular scale of the galaxy of interest.

Although we may, as a result, be able to roughly estimate the SB limit of each survey catalog, we prefer to use the derived value as a reference to fit the observed mass function and provide the best-fitting SB limit in the subsequent subsections.

**5.3. How to Fit to Observations**

To fit the observed GSMFs, we generate the simulation GSMFs for a range of 23 mag arcsec$^{-2} \le \mu^{\text{lim}} \le 25$ mag
<figure><img src="image_19.png" /><figcaption>**Figure 16.** Comparison with JWST observations. From top to bottom, we show the HR5 mass functions (thick line) at $z=13.6, 11.7$, and 10. (observations) Haslbauer et al. (2022): Ⓐ is ID 1514 ($z=9.9$), Ⓑ is GL-z11 ($z=11$), and Ⓒ is GL-z13 ($z=13$). Curtis-Lake et al. (2023): Ⓐ is GS-z10-0 ($z=10.38$), Ⓒ is GS-z11-0 ($z=11.5$), Ⓑ is GS-z12-0 ($z=12.63$), and Ⓒ is GS-z13-0 ($z=13.2$). The gray downside arrow marks the galaxy stellar mass calculated when an environment-dependent IMF is applied (for details, see Haslbauer et al. 2022).</figcaption></figure>

$\text{arcsec}^{-2}$ with a step size of $\Delta\mu = 0.1$ mag $\text{arcsec}^{-2}$ and apply $\chi^2$ fitting for data points at $5 \times 10^9 M_{\odot} \le M_*^g \le 10^{11} M_{\odot}$.

For the survey catalogs used in this subsection, Table 1 lists the fitting results (fourth column) with the catalog information on the magnitude limit (second column) and the seeing of the image survey in the r band (fifth column). For $\mu_{\text{lim}}$ we simply subtract a factor of 2 from $m_{\text{ap}}^{\text{lim}}$ assuming $f_{\text{p}} = 5$ and $f_{\text{ap}} = 5$. The most probable cause of the difference between $\mu_{\text{fit}}^{\text{lim}}$ (third column) and $\mu_{\text{p}}^{\text{lim}}$ (second column) may come from the seeing effect of $(\Delta\mu)^{\text{seeing}} = 0.7-1.1$ mag $\text{arcsec}^{-2}$.

5.4. Results

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 \le z \le 2$ grows quickly with a nearly fixed shape, the slope of the GSMF on the 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 as 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 the stellar mass from the brightness of galaxies are essentially different to simulations in which the galaxy stellar mass is rather an intrinsic variable given by the stellar particle mass. On the BCG scale, we measure the galaxy stellar mass using the 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 (Williams et al. 2016; Martin et al. 2019).

Using the galaxy samples in the COSMOS and XMM-LSS fields, together with Spitzer IRAC observations, Adams et al. (2021) fitted the observed stellar mass functions with 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 the redshift evolution in the double Schechter function into account. They found the best-fit models by investigating a wider range of parameters. The resulting fits of the GSMFs turned 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.

The left panel of Figure 19 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 overpopulations 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 cutoff 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 observational data have a flatter slope below the knees of the GSMFs than those of HR5. Second, the simulation predicts an overpopulation of low-mass galaxies ($M_*^g \le 10^{10} M_{\odot}$) compared to observations. On the right panel of the figure, we show the effect of an SB limit of $\langle \mu_r \rangle \le 23.8$ mag $\text{arcsec}^{-2}$. Note that the prominent knee of the observed mass function is not reproduced even with the SB limit. This
<figure><img src="image_20.png" /><figcaption>**Figure 17.** Comparison of GSMFs at high redshifts. The 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. (2021).*</figcaption></figure>

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 the
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.

In Figure 20 we compare the simulated GSMFs with more
observations: Leja et al. (2020; *upper left*), McLeod et al.
(2021; *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).

We measure the amplitude differences of the GSMFs
between observations and simulations. Figure 21 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; McLeod et al. 2021; Thorne et al. 2021). We note
that the aperture correction and SB limit are applied only to
HR5. Among the three simulations, EAGLE is the least deviated
from the three observations over the redshift range of
z ~ 0.5 - 2. The GSMFs of HR5 shows considerable scatters

around the observations while the TNG100 has a systematic
overestimation on all mass scales. To reduce this systematic
bias of TNG100, Tang et al. (2021) showed that the systematic
bias of TNG100 can be significantly reduced by applying the
SB limit and an aperture cut to the mock surveys in a similar
way to ours.

**6. Discussions and Conclusions**

One of the key purposes of numerical simulations is to
reproduce observations by tuning the model parameters,
to understand underlying physics encoded in the observed
properties. This is valid only when comprehensive under-
standing precedes 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 an SB limit of ⟨μ_r⟩ ≤ 25
mag arcsec⁻², our mock surveys reach about 65% complete-
ness at the stellar mass of M_s^g = 10⁹ M☉ at z = 0.875, while the
completeness of the same mass drops sharply down to ~20% at
z = 1.89. Even with ⟨μ_r⟩ ≤ 25 mag arcsec⁻², 80% of galaxies
of M_s^g = 10⁹ M☉ 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 a 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.
<figure><img src="image_21.png" /><figcaption>**Figure 18.** Comparison of GSMFs at intermediate redshifts ($z = 2.5, 2.5, \text{ and } 4$ from left to right). 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), and Grazian et al. (2015).</figcaption></figure>

<table><caption>**Table 1**

Best-fitting SB Limit ($\mu_p^{\text{lim}}$) for each Galaxy Catalog in the r Band for $f_p = 5$ and $f_{ap} = 5$</caption><thead><tr><th>Name</th><th><i>m</i><sub>ap</sub><sup>lim</sup></th><th><i>μ</i><sub>p</sub><sup>lim</sup></th><th><i>μ</i><sub>fit</sub><sup>lim</sup></th><th>Seeing<sup>a</sup></th></tr></thead><tbody><tr><td>Adams et al. (2021)</td><td>26.5</td><td>24.5</td><td>23.8</td><td>0.′′6–0.′′8<sup>b</sup></td></tr><tr><td>McLeod et al. (2021)</td><td>26.4</td><td>24.4</td><td>23.8</td><td>0.′′9–1.′′<sup>c</sup></td></tr><tr><td>Leja et al. (2020)</td><td>27.0<sup>f</sup></td><td>24.8</td><td>24.1</td><td>≤1.′′3</td></tr><tr><td>Wright et al. (2018)</td><td>(23.5)<sup>g</sup></td><td>...</td><td>23.7</td><td>1.′′05<sup>d</sup></td></tr><tr><td>Thorne et al. (2021)</td><td>26.9</td><td>24.9</td><td>23.8</td><td>0.′′74–0.′′76<sup>e</sup></td></tr></tbody></table>

**Notes.**

a FWHM.

b UltraVISTA; McCracken et al. (2012).

c UDS, UltraVISTA, and CFHTLS.

d Andrews et al. (2017).

e HSC-SSP (Aihara et al. 2018).

f $f_{ap} = 3$ (Laigle et al. 2016).

g Magnitude limit in the i band (Andrews et al. 2017).

Tantalizing evidence has been reported by a number of literature works that are consistent with our claim of a substantial amount of missing LSBGs in galaxy surveys. Greene et al. (2022) performed image simulations to obtain 80% detection completeness at an SB limit of $\langle\mu_r\rangle^e \sim 28$ mag arcsec$^{-2}$. van der Burg et al. (2017) also reported a ~80% completeness of LSBGs from the ESO Kilo-Degree Survey with an SB limit of $\langle\mu_r\rangle^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 an SB limit of about $\langle\mu_r\rangle^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 as the missing LSBGs insignificantly contribute to the global SFR in the low-mass range (Davies et al. 2016).

As of today, the ultrafaint field at $\langle\mu_r\rangle \gtrsim 28$ mag arcsec$^{-2}$ is an almost uncharted territory. The Dark Energy Camera (DECam) reaches $\langle\mu_r\rangle = 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 $\langle\mu_g\rangle \sim 31$ mag arcsec$^{-2}$ (Danieli et al. 2020), and many extragalactic survey projects (HST: Borlaff et al. 2019; LSST: Brough et al. 2020; Euclid: Euclid Collaboration et al. 2022) are ongoing or scheduled for ultrafaint-field observations of LSBGs. Recently, a space-born observatory, the JWST reached faint surface-brightness limits of $\langle\mu_{\text{NIRCam}}\rangle = 31.1–31.3$ mag arcsec$^{-2}$ at $3\sigma$ fluctuations in $10'' \times 10''$ size with about 24 hr exposure (Montes & Trujillo 2022).

Even with the correction for the low-surface-brightness effects, the GSMFs of HR5 still has a knee around the characteristic mass less distinct than that of observed GSMFs. It may come from an insufficient star formation efficiency ($\epsilon_* = 2\%$) 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 mass scales.

Theoretical studies have suggested, based on controlled hydrodynamical simulations, that SFE ranges from 2% to 9.5% (e.g., Matzner 2002; Semenov et al. 2016; Kim et al. 2021). We adopt $\epsilon_* = 2\%$ in HR5, which is at the lower end of the SFE range. From a pair of simulations with two different SFEs ($\epsilon_* = 2\%$ and 4%), we found that the simulation for $\epsilon_* = 4\%$ produces 1.5–3 times more galaxies than the simulation for $\epsilon_* = 2\%$ does at $z \simeq 10$. The ratio however drops to ~1.3 for $M \gtrsim 10^9 M_\odot$ at $z = 1.8$ (see Appendix D). This indicates that higher SFEs would not seriously change our fitting to $\mu_{\text{fit}}^{\text{lim}}$ at $z \lesssim 5$, even though the SFE is still a substantial factor regulating the growth of galaxies. As seen in Figure 16, HR5 is also consistent with the current observations made by the JWST for the number density of massive protogalaxies at high $z$. Further observations and simulations are, however, required to statistically confirm the populations of massive protogalaxies in the current $\Lambda$CDM paradigm.

Also, as shown in Appendix E, the definition of the effective radius has a significant impact on the GSMF measurement in simulations. As the 3D effective radius is usually bigger than that measured in two dimensions, the galaxy SB is lower and, consequently, the GSMFs are more seriously affected by the SB limit. Likewise, the SB limit has more significant impact on
<figure><img src="image_22.png" /><figcaption>**Figure 19.** GSMFs from the HR5 (*symbols*) and observations (Adams et al. 2021; *lines*) at several low redshifts. The left panel shows the uncorrected HR5 GSMF while the middle panel shows the simulation results after applying the 30 kpc aperture correction only. The right panel shows the HR5 GSMF after applying an additional constraint on the surface brightness of the galaxies. In this plot we assume that the galaxy data of Adams et al. (2021) have a fitting limit SB of about 23.8 $mag\ arcsec^{-2}$ (see text for details).</figcaption></figure>

<figure><img src="image_23.png" /><figcaption>**Figure 20.** Comparisons of GSMFs between HR5 (*symbols*) and observations. Counterclockwise from the upper left panel, the lines are the observed GSMFs given by Leja et al. (2020), McLeod et al. (2021), Wright et al. (2018), and Thorne et al. (2021) at several low redshifts. Each color indicates the same redshift as used in Figure 19. 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.</figcaption></figure>

the GSMFs when the effective radius is measured in the face-on plane, compared with that measured in the plane projected along the *x*-axis of the simulation box.

The simulation resolution, the seeing effect, and the source detection criteria are perhaps other important factors that influence the completeness of surveys. HR5 is not able to resolve galaxies in scales smaller than 1 pkpc. The resolution limit may severely affect early-type galaxies that have radii typically smaller than those of late types at a given stellar mass (e.g., Nair & Abraham 2010). However, because of the relatively small fraction of early-type galaxies at $z \simeq 1$ (~15%; Hwang & Park 2009), the resolution limit of HR5 would not

<figure><img src="image_24.png" /><figcaption>**Figure 21.** 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 top to bottom. The base GSMFs of HR5 are corrected with a 3D aperture of 30 pkpc and an SB limit of 23.8 $mag\ arcsec^{-2}$.</figcaption></figure>

significantly affect our results. The seeing effect can also influence the galaxy SB in observations (see Appendix C.2). It reduces the observed SB of galaxies by more than 0.5 $mag$
$arcsec^{-2}$ when the effective radius of a galaxy is not sufficiently larger than the FWHM of the PSF. This may lead to the substantial uncertainty in deriving the SB limit in each survey catalog. Also, it is worthwhile to note that the pixel-level source detection limit may be another key factor in measuring the completeness of a galaxy catalog based on photometry (see Appendix C.1). If proper information on the source detection is provided, it would be more helpful in fine-tuning the GSMFs. These three factors are closely related to each other, but disentangling them is beyond the scope of this paper.

## Acknowledgments

The authors thank the Korea Institute for Advanced Study for providing computing resources (KIAS Center for Advanced Computation Linux Cluster System) for this work. A special thank should be given to Profs. Ho Seong Hwang and Jihoon Kim who commented on the photometry of extragalactic images. J.K. was supported by a KIAS Individual Grant (KG039603) via the Center for Advanced Computation at Korea Institute for Advanced Study. J.L. is supported by the National Research Foundation of Korea (NRF-2021R1C1C2011626). This work benefited from the out-standing support provided by the KISTI National Super-computing Center and its Nurion Supercomputer through the Grand Challenge Program (KSC-2018-CHA-0003). A large data transfer was supported by KREONET, which is managed and operated by KISTI. B.K.G. 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 program (ChETEC-INFRA—Project No. 101008324). This research was also partially supported by the ANR-19-CE31-0017; http://www.secular-evolution.org. Y.K. is supported by NRF-2022M3K3A1097100. J.H.S. acknowledges support from the Korea Astronomy and Space Science Institute under the R&D program (Project No. 2023-1-830-00) and the National Research Foundation of Korea grant (2021R1C1C1003785) funded by the Ministry of Science, ICT & Future Planning. High-performance computing resources for this research were partially supported by the Research Solution Center and the National Supercomputing Center for Astrophysics and Space Sciences in the Institute for Basic Science. This research was supported in part by the National Science Foundation under Grant No. NSF PHY-1748958.

# Appendix A

## Halo and Galaxy Findings

### A.1. AMR and Unified Data Types

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 data type for the four different matter

species. The unified data type 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 data type 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

To find virialized structures in cosmological N-body simulations, a percolation method like the FoF algorithm is frequently employed. One fifth of 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 ~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_{\text{link}} = \frac{1}{5} \left( \frac{m_p}{\Omega_{m0} \rho_c} \right)^{1/3}, \qquad (\text{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_{\text{comb}} = \frac{1}{2}(\ell_1 + \ell_2), \qquad (\text{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 pairwise 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 22 shows the images of three representative FoF halos identified at $z=0.625$. In this figure, a clear difference in the density distributions is observed between the dark matter (left) and stellar components (middle column). The stellar distribution is more compact than the dark matter and gas distributions. 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.
<figure><img src="image_25.png" /><figcaption>**Figure 22.** Three projected density maps of massive halos identified with the variable length in the FoF at $z = 0.625$. From left to right, we show the dark matter, stellar, and gas density fields. In the middle panels, stellar particles with higher metallicity have a more pinkish color. The white bar at the bottom in each panel marks the scale of $l_{bar} = 0.5 h^{-1} \text{cMpc}$. From top to bottom, the total masses of the FoF halo are $1.09 \times 10^{13} h^{-1} M_{\odot}$, $1.42 \times 10^{13} h^{-1} M_{\odot}$, and $2.63 \times 10^{13} h^{-1} M_{\odot}$, respectively.</figcaption></figure>

### A.3. Halo Mass Function of HR5

The HMF is one of the fundamental statistics that closely relates to the cosmological models (Press & Schechter 1974; Sheth & Tormen 1999; Jenkins et al. 2001; Ganeshaiah et al. 2018). It has been used as the fiducial statistics method 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 percent deviations.

Throughout this paper, the mass function is defined as

$$ \Phi(M) = \frac{\Delta N}{\Delta \log_{10} M}, \qquad (A3) $$

where $\Delta N$ is the number of galaxies in a mass bin of size $\Delta \log_{10} M$. In Figure 23, we show the HMFs at $z = 0.625, 1, 2, 3$, and $4$. In the mass range of $10^{11} M_{\odot} \le M_{\text{FoF}} \le 5 \times 10^{14} M_{\odot}$, the HR5 HMFs are well described by the reference model of the Sheth & Tormen functions (Sheth & Tormen 1999).
<figure><img src="image_26.png" /><figcaption>Figure 23. FoF HMFs at different redshifts. The symbols with Poisson error bars are simulated HMFs while the lines are the reference mass functions (Sheth & Tormen 1999).</figcaption></figure>

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

We developed the PSB-based Galaxy Finder (pGalF) based on the 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 subhalo, 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 2008; Renaud et al. 2016; van den Bosch et al. 2018). 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( \frac{m/M}{\alpha + \beta} \right)^{1/3}, \qquad (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, α, to account for the mass distribution of the host at the satellite position as

$$\alpha(R) \equiv 2 - \left( \frac{d \ln M}{d \ln R} \right) |_{R}. \qquad (A5)$$

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

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

where Ω(≡ V_c/R) is the angular velocity (King 1962) for a circular velocity V_c. From this equation, one may easily derive that β = −2K_t/W, where K_t is the rotational kinetic energy of the satellite and W is its potential energy. For a circular orbit, β = 1 but β = 0 for a radial motion. If the satellite is not bound to the host, then we neglect this effect and set β = 0. From this parameterization we recover nearly all the modeled radii of the satellite. For example, we recover the Jacobi radius with α + β = 3 for a satellite in a circular orbit around a noncontacting host, or we recover the Roche limit (when α + β = 2) for a radial orbit around the host.

Now, we discuss the α parameter in more detail for the cosmological extended object. Assuming a Navarro–Frenk–White 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) - \frac{cs}{1 + cs} \right], \qquad (A7)$$

where c is the concentration index, s is the scaled radius to the virial radius as s ≡ R/R_v, and

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

Then, we finally get

$$\alpha(s) \equiv 2 - \left( \frac{cs}{1 + cs} \right)^2. \qquad (A9)$$
The concentration index is empirically given as (Comerford & Natarajan 2007)

$$c = 14.5 \left( \frac{M_v}{1.3 \times 10^{13} h^{-1} M_{\odot}} \right)^{-0.15} (1+z)^{-1}, \quad (A10)$$

where z is the redshift and h is the Hubble expansion rate divided by 100 km s⁻¹ Mpc⁻¹.

### A.4.3. Hierarchical Membership Determination

Figure 24 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 innermost 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 watershed technique. It is a percolation method that grows a volume of interest by lowering the density threshold around the target density peak. Therefore, 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 α 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 α 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 β 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 reapply the membership determination for member particles to confirm whether they could still hold their membership under the updated circumstances.

<figure><img src="image_27.png" /><figcaption>Figure 24. 2D example of the density field and isodensity contours, which split and arrange areas in a hierarchical manner. The density field is made using a virialized FoF halos using its particles.</figcaption></figure>

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

Like the PSB method, 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 a 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 Figure 25) due to the dynamical frictions and tidal stripping in a group or cluster region.
<figure><img src="image_28.png" /><figcaption>**Figure 25.** Galaxy finding results applied to the FoF halos as given in Figure 22. Contours are shown after fitting the density distribution of each component of galaxies.</figcaption></figure>

## Appendix B

### Difference between the Absolute Magnitudes in the SDSS r and Cousin R bands

We have examined the difference between the absolute magnitudes of the HR5 galaxies of $M_{\star} > 10^{9} M_{\odot}$ in the SDSS r

and Cousin R- ands at $z=0.625, 1.13$, and $1.89$ (see Figure 26). The R-band absolute magnitude $M_{R}$ is slightly lower than $M_{r}$, and the difference becomes smaller with decreasing redshifts. We present the $\chi^{2}$-minimization fitting result with a magenta solid curve in each panel.

<figure><img src="image_29.png" /><figcaption>**Figure 26.** Difference Between the Absolute Magnitudes in the SDSS r and Cousin R bands. The magenta solid lines present the fitting functions derived from the $M_{r} - M_{R}$ distribution using the $\chi^{2}$ minimization.</figcaption></figure>
<figure><img src="image_30.png" /><figcaption>**Figure 27.** Derived SB limit as a function of the pixel scale for several $f_p$. The solid, dotted, short-dashed, and long-dashed vertical lines mark the pixel scales of Hyper Suprime-Cam (Aihara et al. 2018), CFHT MegaCam, UltraVISTA/COSMOS (McCracken et al. 2012), and ACS of HST, respectively.</figcaption></figure>

# Appendix C
## Observational Constraints on the Source Detection

### C.1. Pixel-scale Source Detection Criteria

The source detection in the photometry of extragalactic images usually begins with searching image pixels that have values larger than $f_p\sigma_p$, where $f_p$ is the multiplication factor above the background noises, ($\sigma_p$). The pixels are grouped by a given connection criterion, being regarded as a galaxy candidate. Then, an aperture with the solid angle $A_{ap}$ is located at the center of the candidate to calculate the total flux by integrating the pixel values enclosed in the aperture. Usually, there is a criterion to limit the faint-end magnitude ($m_{ap}^{lim}$) of galaxy candidates. The average pixel value should be $f_{ap}$ times larger than the background fluctuations on the aperture scale ($\sigma_{ap}$). For the Gaussian fluctuations of the background, we may assume the scaling relation between these two noises as, $\sigma_p = \sigma_{ap}(A_{ap}/A_p)^{1/2}$, where $A_p$ is the solid angle of the pixel. Finally we may formulate the SB limit of the source detection as (also see Román et al. 2020),

$$ \mu_{p}^{\text{lim}} = -2.5 \log \left[ \frac{f_{p}}{f_{\text{ap}}} \left( \frac{A_{\text{ap}}}{A_{p}} \right)^{1/2} \right] + \mu_{\text{ap}}^{\text{lim}}, \quad (\text{C1}) $$

where the average SB of the aperture is related to the galaxy magnitude limit as $\mu_{ap}^{lim} = m_{ap}^{lim} + 2.5 \log(A_{ap}/\text{arcsec}^2)$.

Figure 27 shows the offset between $\mu_{p}^{lim}$ and $m_{ap}^{lim}$ as a function of the pixel scale with the aperture size ($r_{ap} = 1''$) and the source detection level ($f_{ap} = 5$), which are recommended by SExtractor (Bertin & Arnouts 1996). For the source detection criterion of $3\sigma_p \le \mu_{p}^{lim} \le 5\sigma_p$, the SB limit would be about 1–2 mag brighter than the aperture magnitude limit for the faint galaxy photometry of the COSMOS field.

<figure><img src="image_31.png" /><figcaption>**Figure 28.** Angular scales of simulated galaxies (candle bars), pixel (horizontal bar), and average seeing at low $z$. The gray region delineates the distribution of FWHM of PSF in the UltraVISTA observations (McCracken et al. 2012).</figcaption></figure>

### C.2. Seeing Effects

In this subsection, we briefly discuss the seeing effects on the galaxy SB. First, in Figure 28, we compare the angular scales of simulated galaxies (candle bars), COSMOS pixel (0"15; a solid horizontal bar), and the range of the observed FWHM of a PSF (0"4–0"8; a gray region). Even though the pixel scale of, for example, the COSMOS field is sufficiently smaller than the galaxy angular size, the galaxy shape and SB are severely contaminated by the seeing effect (McCracken et al. 2012).

By simulating the seeing effect on the galaxy models with Sérsic profiles, Trujillo et al. (2001) showed that the observed SB becomes fainter than the input value by $\Delta\mu^{\text{seeing}} \gtrsim 0.5$ when the galaxy angular size is comparable to or smaller than the FWHM of the PSF.

# Appendix D
## Effect of Star Formation Efficiency on the GSMF

In this section, we examine the effect of the star formation efficiency $\epsilon_*$ on the GSMF using two auxiliary simulations. The first simulation (M2) shares the same simulation parameters with HR5 except for the smaller box size, $L_{\text{box}}^{\text{aux}} = 128$ cMpc. The second simulation (M4) is intended to study the role of SFE having the same parameters as of M2 except for $\epsilon_* = 4\%$ without any significant systematic bias. Both simulations are run down to $z = 1.7$.

The global SFR is directly regulated by $\epsilon_*$, producing, as a result, different GSMFs. The top panel of Figure 29 shows the ratio of the global SFRs between M2 and M4 as a function of the redshift. The ratio decreases over time, but it is still higher than unity at $z \sim 2$. We also compare the GSMFs of the two simulations in the bottom panel. The difference gradually decreases with decreasing redshift at least until $z = 1.8$. This is because higher SFEs essentially accompany more energetic stellar feedback, which blows gas from galaxies while baryon reservoirs are the same between the two. Note that the volume
<figure><img src="image_32.png" /><figcaption>**Figure 29.** Ratios of the global SFR (top) and GSMFs (bottom panel) between M2 and M4. The spurious spike at $z=4$ is caused by the discrete global mesh refinement scheme of RAMSES in the global SFR ratio.</figcaption></figure>

of the two simulations is not large enough to make the statistic of GSMFs at $z > 10$.

# Appendix E

## Dependence of GSMFs on the Effective Radius

The effective radius of a galaxy can be defined in various ways, and thus we examine in this section how the different definition may influence the resulting GSMFs. Figure 30 shows

<figure><img src="image_33.png" /><figcaption>**Figure 30.** Simulated GSMFs for different definitions of the galaxy radius. On the bottom panel we show the GSMFs for the cases of effective radii measured in the projected distributions of mass (dotted) and *r*-band brightness (long dotted–dashed). On the other hand, the GSMFs with the radii measured in the face-on view of a galaxy are shown with long-dashed (mass) and short dotted–dashed (*r*-band) lines. The short-dashed line is for the case of radii measured in three dimensions. The true GSMF is shown with the thick solid line (all). In the top panel, we show the differences in the GSMFs with respect to the true GSMF. In this plot, we use an *SB* cut of $\mu_r^{\text{lim}} = 25 \ \text{mag arcsec}^{-2}$.</figcaption></figure>

the GSMFs derived based on five different definitions of the effective radius at $z = 0.625$ (left column) and $z = 1.89$ (right column). Because the effective radius is generally larger when it is measured in three dimensions than in two dimensions, the GSMFs are most significantly lowered with $R_{1/2}$ in 3D (dashed purple). Likewise, the GSMFs are diminished more due to the *SB* limit in the face-on projection than in the case of projection along the *x*-axis, which is effectively a random projection. On the other hand, there is no substantial difference between the radii measured from the stellar mass and the *r*-band luminosity distribution.
ORCID iDs

Juhan Kim https://orcid.org/0000-0002-4391-2275

Jaehyun Lee https://orcid.org/0000-0002-6810-1778

Clotilde Laigle https://orcid.org/0000-0003-4446-3130

Yohan Dubois https://orcid.org/0000-0003-0225-6387

Yonghwi Kim https://orcid.org/0000-0003-4164-5414

Changbom Park https://orcid.org/0000-0001-9521-6397

Christophe Pichon https://orcid.org/0000-0003-0695-6735

Brad K. Gibson https://orcid.org/0000-0003-4446-3130

Jihye Shin https://orcid.org/0000-0001-5135-1693

References

Adams, N. J., Bowler, R. A. A., Jarvis, M. J., Häußler, B., & Lagos, C. D. P. 2021, MNRAS, 506, 4933

Aihara, H., Arimoto, N., Armstrong, R., et al. 2018, PASJ, 70, S4

Andrews, S. K., Driver, S. P., Davies, L. J. M., et al. 2017, MNRAS, 464, 1569

Baldry, I. K., Glazebrook, K., & Driver, S. P. 2008, MNRAS, 388, 945

Behroozi, P. S., Wechsler, R. H., & Wu, H.-Y. 2013, ApJ, 762, 109

Bell, E. F., & de Jong, R. S. 2001, ApJ, 550, 212

Bell, E. F., McIntosh, D. H., Katz, N., & Weinberg, M. D. 2003, ApJS, 149, 289

Bellstedt, S., Lidman, C., Muzzin, A., et al. 2016, MNRAS, 460, 2862

Benson, A. J., Bower, R. G., Frenk, C. S., et al. 2003, ApJ, 599, 38

Bertin, E., & Arnouts, S. 1996, A&AS, 117, 393

Binney, J., & Tremaine, S. 2008, Galactic Dynamics (2nd ed.; Princeton, NJ: Princeton Univ. Press)

Blanton, M. R., Lupton, R. H., Schlegel, D. J., et al. 2005, ApJ, 631, 208

Blanton, M. R., & Roweis, S. 2007, AJ, 133, 734

Boco, L., Lapi, A., Chruslinska, M., et al. 2021, ApJ, 907, 110

Bolzonella, M., Kovač, K., Pozzetti, L., et al. 2010, A&A, 524, A76

Booth, C. M., & Schaye, J. 2009, MNRAS, 398, 53

Borlaff, A., Trujillo, I., Román, J., et al. 2019, A&A, 621, A133

Boylan-Kolchin, M., Ma, C.-P., & Quataert, E. 2008, MNRAS, 383, 93

Brough, S., Collins, C., Demarco, R., et al. 2020, arXiv:2001.11067

Chabrier, G. 2003, PASP, 115, 763

Comerford, J. M., & Natarajan, P. 2007, MNRAS, 379, 190

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

Conroy, C., Wechsler, R. H., & Kravtsov, A. V. 2007, ApJ, 668, 826

Crain, R. A., Schaye, J., Bower, R. G., et al. 2015, MNRAS, 450, 1937

Cross, N., Driver, S. P., Couch, W., et al. 2001, MNRAS, 324, 825

Curtis-Lake, E., Carniani, S., Cameron, A., et al. 2023, NatAs, 7, 622

Dabringhausen, J., & Fellhauer, M. 2016, MNRAS, 460, 4492

Dahlen, T., Mobasher, B., Somerville, R. S., et al. 2005, ApJ, 631, 126

Danieli, S., Lokhorst, D., Zhang, J., et al. 2020, ApJ, 894, 119

Davidzon, I., Ilbert, O., Faisst, A. L., Sparre, M., & Capak, P. L. 2018, ApJ, 852, 107

Davidzon, I., Ilbert, O., Laigle, C., et al. 2017, A&A, 605, A70

Davies, J. I., Davies, L. J. M., & Keenan, O. C. 2016, MNRAS, 456, 1607

Donnan, C. T., McLeod, D. J., Dunlop, J. S., et al. 2023, MNRAS, 518, 6011

Donnari, M., Pillepich, A., Nelson, D., et al. 2019, MNRAS, 485, 4817

dos Reis, S. N., Buitrago, F., Papaderos, P., et al. 2020, A&A, 634, A11

Driver, S. P. 1999, ApJL, 526, L69

Dubois, Y., Beckmann, R., Bournaud, F., et al. 2021, A&A, 651, A109

Dubois, Y., Devriendt, J., Slyz, A., & Teyssier, R. 2012, MNRAS, 420, 2662

Dubois, Y., Peirani, S., Pichon, C., et al. 2016, MNRAS, 463, 3948

Dubois, Y., Pichon, C., Welker, C., et al. 2014a, MNRAS, 444, 1453

Dubois, Y., Volonteri, M., & Silk, J. 2014b, MNRAS, 440, 1590

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

Durkalec, A., Le Fèvre, O., de la Torre, S., et al. 2015, A&A, 576, L7

Euclid Collaboration, Borlaff, A. S., Gómez-Alvarez, P., et al. 2022, A&A, 657, A92

Faber, S. M., Willmer, C. N. A., Wolf, C., et al. 2007, ApJ, 665, 265

Fernández, Lorenzo, M., Sulentic, J., Verdes-Montenegro, L., & Argudo-Fernández, M. 2013, MNRAS, 434, 325

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

Foucaud, S., Conselice, C. J., Hartley, W. G., et al. 2010, MNRAS, 406, 147

Furlong, M., Bower, R. G., Theuns, T., et al. 2015, MNRAS, 450, 4486

Gabasch, A., Hopp, U., Feulner, G., et al. 2006, A&A, 448, 101

Ganeshaiha, V. P., Cautun, M., van de Weygaert, R., et al. 2018, MNRAS, 481, 414

Geller, M. J., Diaferio, A., Kurtz, M. J., Dell'Antonio, I. P., & Fabricant, D. G. 2012, AJ, 143, 102

Golden-Marx, J. B., Miller, C. J., Zhang, Y., et al. 2022, ApJ, 928, 28

Gonzalez, A. H., Zaritsky, D., & Zabludoff, A. I. 2007, ApJ, 666, 147

González, V., Labbé, I., Bouwens, R. J., et al. 2011, ApJL, 735, L34

Graham, A. W., Driver, S. P., Petrosian, V., et al. 2005, AJ, 130, 1535

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

Greene, J. E., Greco, J. P., Goulding, A. D., et al. 2022, ApJ, 933, 150

Guennou, L., Adami, C., Da Rocha, C., et al. 2012, A&A, 537, A64

Haslbauer, M., Kroupa, P., Zonoozi, A. H., & Haghi, H. 2022, ApJL, 939, L31

Hopkins, A. M. 2004, ApJ, 615, 209

Hwang, H. S., & Park, C. 2009, ApJ, 700, 791

Ichikawa, T., Kajisawa, M., & Akhlaghi, M. 2012, MNRAS, 422, 1014

Ilbert, O., Tresse, L., Zucca, E., et al. 2005, A&A, 439, 863

Jackson, R. A., Martin, G., Kaviraj, S., et al. 2021, MNRAS, 502, 4262

Jenkins, A., Frenk, C. S., White, S. D. M., et al. 2001, MNRAS, 321, 372

Katsianis, A., Tescari, E., & Wyithe, J. S. B. 2015, MNRAS, 448, 3001

Kaviraj, S., Laigle, C., Kimm, T., et al. 2017, MNRAS, 467, 4739

Keller, B. W., Munshi, F., Trebitsch, M., & Tremmel, M. 2023, ApJL, 943, L28

Kim, J., & Park, C. 2006, ApJ, 639, 600

Kim, J., Park, C., Gott, J., Richard, I., & Dubinski, J. 2009, ApJ, 701, 1547

Kim, J., Park, C., L'Huillier, B., & Hong, S. E. 2015, JKAS, 48, 213

Kim, J.-G., Ostriker, E. C., & Filippova, N. 2021, ApJ, 911, 128

King, I. 1962, AJ, 67, 471

Kravtsov, A. V., Vikhlinin, A. A., & Meshcheryakov, A. V. 2018, AstL, 44, 8

Laigle, C., McCracken, H. J., Ilbert, O., et al. 2016, ApJS, 224, 24

Lee, J., Kim, S., Jeong, H., et al. 2018, ApJ, 864, 69

Lee, J., Shin, J., Snaith, O. N., et al. 2021, ApJ, 908, 11

Lee, K.-S., Ferguson, H. C., Wiklind, T., et al. 2012, ApJ, 752, 66

Legrand, L., McCracken, H. J., Davidzon, I., et al. 2019, MNRAS, 486, 5468

Leja, J., Speagle, J. S., Johnson, B. D., et al. 2020, ApJS, 893, 111

Łokas, E. L., & Mamon, G. A. 2001, MNRAS, 321, 155

Lower, S., Narayanan, D., Leja, J., et al. 2020, ApJ, 904, 33

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

Mahajan, S., Drinkwater, M. J., Driver, S., et al. 2018, MNRAS, 475, 788

Marchesini, D., van Dokkum, P. G., Schreiber, N. M. F., et al. 2009, ApJ, 701, 1765

Marinacci, F., Vogelsberger, M., Pakmor, R., et al. 2018, MNRAS, 480, 5113

Martin, G., Kaviraj, S., Laigle, C., et al. 2019, MNRAS, 485, 796

Matzner, C. D. 2002, ApJ, 566, 302

McCarthy, I. G., Schaye, J., Bird, S., & Le Brun, A. M. C. 2017, MNRAS, 465, 2936

McCracken, H. J., Milvang-Jensen, B., Dunlop, J., et al. 2012, A&A, 544, A156

McKinney, J. C., Tchekhovskoy, A., & Blandford, R. D. 2012, MNRAS, 423, 3083

McLeod, D. J., McLure, R. J., Dunlop, J. S., et al. 2021, MNRAS, 503, 4413

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

Mihos, J. C. 2019, arXiv:1909.09456

Montes, M. 2022, NatAs, 6, 308

Montes, M., & Trujillo, I. 2022, ApJL, 940, L51

Naidu, R. P., Oesch, P. A., Setton, D. J., et al. 2022, arXiv:2208.02794

Naiman, J. P., Pillepich, A., Springel, V., et al. 2018, MNRAS, 477, 1206

Nair, P. B., & Abraham, R. G. 2010, ApJS, 186, 427

Navarro, J. F., Frenk, C. S., & White, S. D. M. 1997, ApJ, 490, 493

Nedkova, K. V., Häußler, B., Marchesini, D., et al. 2021, MNRAS, 506, 928

Nelson, D., Pillepich, A., Springel, V., et al. 2018, MNRAS, 475, 624

O'Mill, A. L., Duplancic, F., García Lambas, D., & Sodré, L. 2011, MNRAS, 413, 1395

Panter, B., Jimenez, R., Heavens, A. F., & Charlot, S. 2007, MNRAS, 378, 1550

Park, C., Lee, J., Kim, J., et al. 2022, ApJ, 937, 15

Pillepich, A., Nelson, D., Hernquist, L., et al. 2018, MNRAS, 475, 648

Press, W. H., & Schechter, P. 1974, ApJ, 187, 425

Puchwein, E., & Springel, V. 2013, MNRAS, 428, 2966

Renaud, F., Agertz, O., & Gieles, M. 2016, MNRAS, 465, 3622

Ricciardelli, E., Vazdekis, A., Cenarro, A. J., & Falcón-Barroso, J. 2012, MNRAS, 424, 172

Román, J., Castilla, A., & Pascual-Granado, J. 2021, A&A, 656, A44

Román, J., Trujillo, I., & Montes, M. 2020, A&A, 644, A42

Roy, N., Napolitano, N. R., La Barbera, F., et al. 2018, MNRAS, 480, 1057

Samuel, J., Wetzel, A., Santistevan, I., et al. 2022, MNRAS, 514, 5276

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

Schmidt, M. 1959, ApJ, 129, 243

Semenov, V. A., Kravtsov, A. V., & Gnedin, N. Y. 2016, ApJ, 826, 200

Sheth, R. K., & Tormen, G. 1999, MNRAS, 308, 119

Shuntov, M., McCracken, H. J., Gavazzi, R., et al. 2022, A&A, 664, A61
Silk, J., & Mamon, G. A. 2012, RAA, 12, 917

Sohn, J., Geller, M. J., Rines, K. J., et al. 2018, ApJ, 856, 172

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

Springel, V., Pakmor, R., Pillepich, A., et al. 2018, MNRAS, 475, 676

Stabenau, H. F., Connolly, A., & Jain, B. 2008, MNRAS, 387, 1215

Stefanon, M., Bouwens, R. J., Labbé, I., et al. 2021, ApJ, 922, 29

Suto, D., Kitayama, T., Osato, K., Sasaki, S., & Suto, Y. 2016, PASJ, 68, 14

Tang, L., Lin, W., Wang, Y., & Napolitano, N. R. 2021, MNRAS, 508, 3321

Thorne, J. E., Robotham, A. S. G., Davies, L. J. M., et al. 2021, MNRAS, 505, 540

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

Trujillo, I., Aguerri, J. A. L., Cepa, J., & Gutiérrez, C. M. 2001, MNRAS, 328, 977

Trujillo, I., Chamba, N., & Knapen, J. H. 2020, MNRAS, 493, 87

Valls-Gabaud, D. & MESSIER Collaboration 2017, in IAU Symp. 321, Formation and Evolution of Galaxy Outskirts, ed. A. Gil de Paz, J. H. Knapen, & J. C. Lee (Cambridge: Cambridge Univ. Press), 199

van de Ven, G., van Dokkum, P. G., & Franx, M. 2003, MNRAS, 344, 924

van den Bosch, F. C., Ogiya, G., Hahn, O., & Burkert, A. 2018, MNRAS, 474, 3043

van der Burg, R. F. J., Hoekstra, H., Muzzin, A., et al. 2017, A&A, 607, A79

van der Wel, A., Franx, M., van Dokkum, P. G., et al. 2005, ApJ, 631, 145

Vazdekis, A., Koleva, M., Ricciardelli, E., Röck, B., & Falcón-Barroso, J. 2016, MNRAS, 463, 3409

Vazdekis, A., Ricciardelli, E., Cenarro, A. J., et al. 2012, MNRAS, 424, 157

Vogelsberger, M., Genel, S., Sijacki, D., et al. 2013, MNRAS, 436, 3031

Vogelsberger, M., Genel, S., Springel, V., et al. 2014, MNRAS, 444, 1518

Weaver, J. R., Davidzon, I., Toft, S., et al. 2022, arXiv:2212.02512

Wechsler, R. H., & Tinker, J. L. 2018, ARA&A, 56, 435

Weigel, A. K., Schawinski, K., & Bruderer, C. 2016, MNRAS, 459, 2150

Williams, R. P., Baldry, I. K., Kelvin, L. S., et al. 2016, MNRAS, 463, 2746

Wright, A. H., Driver, S. P., & Robotham, A. S. G. 2018, MNRAS, 480, 3491

Yoo, J., Ko, J., Kim, J.-W., & Kim, H. 2021, MNRAS, 508, 2634

Zwicky, F. 1951, PASP, 63, 61