# Born to be Starless: Revisiting the Missing Satellite Problem

Seyoung Jeon, Sukyoung K Yi, Emanuele Contini, Yohan Dubois, San Han, Katarina Kraljic, Sebastien Peirani, Christophe Pichon, Jinsu Rhee

► To cite this version:

Seyoung Jeon, Sukyoung K Yi, Emanuele Contini, Yohan Dubois, San Han, et al. Born to be Starless: Revisiting the Missing Satellite Problem. The Astrophysical Journal, 2025, 988 (1), pp.136. 10.3847/1538-4357/ade2e0. hal-05115939

HAL Id: hal-05115939

https://hal.science/hal-05115939v1

Submitted on 29 Jul 2025

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

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

CrossMark

# Born to be Starless: Revisiting the Missing Satellite Problem

Seyoung Jeon<fnref n="1" />, Sukyoung K. Yi<fnref n="1" />, Emanuele Contini<fnref n="1" />, Yohan Dubois<fnref n="2" />, San Han<fnref n="2" />, Katarina Kraljic<fnref n="3" />, Sebastien Peirani<fnref n="2" />,<fnref n="4" />,⁵, Christophe Pichon<fnref n="2" />,<fnref n="6" />, and Jinsu Rhee<fnref n="1" />,<fnref n="2" />,<fnref n="7" />

Received 2025 April 8; revised 2025 June 2; accepted 2025 June 6; published 2025 July 18

## Abstract

The two massive galaxies in the Local Group both host substantially fewer satellites than the subhalos expected from the cold dark matter paradigm, and recent investigations have highlighted the interplay between baryons and dark matter. We investigate the processes that make subhalos starless, using high-resolution cosmological simulations. We find that the number of satellites around Milky Way analogs closely aligns with observations, which accords with recent studies. In our simulations, the majority of subhalos are devoid of stars, i.e., “starless.” We first examined supernova feedback and the environmental effects associated with subhalos’ orbital motion as candidates for their origin. However, neither seems to be the main driver. Supernova feedback causes a reduction in cold gas in “starred” subhalos, but its impact is not significant. In the case of starless subhalos, supernova feedback is irrelevant because most of them do not have in situ star formation in the first place. The orbital motion in dense environments causes the removal of gas in all subhalos, but it is not enough to remove preexisting stars. The key is found to be the effect of reionization instead. Starless subhalos are initially born in regions that are less efficient in accreting matter. This makes them lack sufficiently dense gas to self-shield from UV background heating, preventing their gas from cooling below the star formation threshold. This indicates that starless subhalos are not made but born.

*Unified Astronomy Thesaurus concepts: Disk galaxies (391); the Milky Way (1054); Galaxy dark matter halos (1880); Hydrodynamical simulations (767)*

## 1. Introduction

In recent decades, simulations based on the Lambda cold dark matter (ΛCDM; D. N. Spergel et al. 2003; Planck Collaboration et al. 2020) cosmology have successfully reproduced large-scale structures of the Universe, enhancing our understanding of structure formation and galaxy evolution (Y. Dubois et al. 2014, 2021; M. Vogelsberger et al. 2014; J. Schaye et al. 2015; I. G. McCarthy et al. 2017; A. Pillepich et al. 2018; R. Davé et al. 2019; R. Pakmor et al. 2023). However, this model has historically faced several small-scale tensions: e.g., the missing satellite problem (A. Klypin et al. 1999; B. Moore et al. 1999; M. Boylan-Kolchin et al. 2011), the cusp-core problem (R. A. Flores & J. R. Primack 1994; B. Moore 1994), and the plane of satellites (D. Lynden-Bell 1976; M. S. Pawlowski et al. 2012). The missing satellite problem traditionally refers to the finding that the number of predicted dark matter (DM) subhalos in DM-only simulations is far larger than the number of observed satellite galaxies around the Milky Way (A. Klypin et al. 1999) and similar galaxies (M. Tanaka et al. 2018). Thanks to advanced observational instruments and techniques, the number of observed satellite galaxies has been increasing (J. Grcevich & M. E. Putman 2009; A. W. McConnachie 2012; Y.-Y. Mao et al. 2021; M. Geha et al. 2024). This progress, combined

with improved completeness corrections and forward modeling techniques (S. Y. Kim et al. 2018; D. Homma et al. 2019; E. O. Nadler et al. 2020), has significantly alleviated much of the original discrepancy for luminous dwarf galaxies; however, the discrepancy still remains, depending on the modeling. Many studies have attempted to resolve this issue by adopting different approaches (see, e.g., J. I. Read & D. Erkal 2019, and references therein).

One deals with alternative DM models by directly reducing the number of DM halos or modifying internal halo properties. Warm DM (WDM) and fuzzy DM suppress the formation of low-mass halos through free-streaming or quantum pressure, respectively, leading to a reduced number of small-scale substructures (W. Hu et al. 2000; M. Götz & J. Sommer-Larsen 2002; E. Polisensky & M. Ricotti 2011; M. R. Lovell et al. 2012; D. J. E. Marsh 2016). For instance, E. Polisensky & M. Ricotti (2011) used N-body simulations to show that WDM can reduce the number of low-mass halos and suggested a lower limit to the WDM particle mass of 2.3 keV to match satellite counts. However, their analysis excluded baryonic effects, and the suggested mass is now disfavored by Lyα forest data (A. Schneider et al. 2014). More recently, E. O. Nadler et al. (2021) suggested that the forward modeling framework with $m_{WDM} > 6.5$ keV can match the fitting results from the Dark Energy Survey (DES) and Pan-STARRS1, and V. Iršič et al. (2024) updated the constraints of $m_{WDM} > 5.7$ keV from quasar spectra. Moreover, O. Garcia-Gallego et al. (2025) showed that the extreme WDM of $m_{WDM} = 1$ keV can be viable with an upper limit to the WDM fraction of 0.16. Self-interacting DM,

<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">¹ Department of Astronomy and Yonsei University Observatory, Yonsei University, 50 Yonsei-ro, Seodaemun-gu, Seoul 03722, Republic of Korea</footnote>

<footnote n="2">² Institut d'Astrophysique de Paris, Sorbonne Universités, et CNRS, UMP 7095, 98 bis bd Arago, 75014 Paris, France</footnote>

<footnote n="3">³ Observatoire Astronomique de Strasbourg, Université de Strasbourg, CNRS, UMR 7550, Strasbourg 67000, France</footnote>

<footnote n="4">⁴ ILANCE² CNRS—University of Tokyo International Research Laboratory, Kashiwa, Chiba 277-8582, Japan</footnote>

<footnote n="3">³ Kavli IPMU (WPI), UTIAS, The University of Tokyo, Kashiwa, Chiba 277-8583, Japan</footnote>

<footnote n="6">⁶ Kyung Hee University, Dept. of Astronomy & Space Science, Yongin-shi, Gyeonggi-do 17104, Republic of Korea</footnote>

<footnote n="7">⁷ Korea Astronomy and Space Science Institute, 776 Daedeokdae-ro, Yuseong-gu, Daejeon 34055, Republic of Korea</footnote>
on the other hand, does not significantly alter halo abundance but transforms the internal structure of halos by flattening central density profiles and adjusting maximum circular velocities, which gives better comparison with observed dynamical masses and profiles (D. N. Spergel & P. J. Steinhardt 2000; L. G. van den Aarssen et al. 2012; M. Vogelsberger et al. 2016).

Alternatively, the effect of the baryonic components could result in a reduced number of satellite galaxies without modifying the nature of ΛCDM. The proper inclusion of baryons influences the DM distribution (A. Zolotov et al. 2012; A. R. Wetzel et al. 2016; S. Garrison-Kimmel et al. 2019) and star formation rates in dwarf galaxies through feedback processes, such as from supernovae (A. M. Brooks et al. 2013; S. Geen et al. 2013; R. Hazenfratz et al. 2024), active galactic nuclei (AGN; G. Dashyan et al. 2018; E. Arjo-
na-Gálvez et al. 2024), and ionizing UV radiation (J. C. Forbes et al. 2016; T. Sawala et al. 2016; A. Emerick et al. 2019). Numerous studies have explored baryonic solutions in simulations and other simplistic models and have succeeded in alleviating the classical missing satellite problem, but they often lack resolutions sufficient to resolve the problem entirely in a fully cosmological framework.

Moreover, these studies have not fully addressed the mechanisms responsible for the widespread transformation of subhalos into starless systems. Reionization, in particular, has long been recognized as a key mechanism in suppressing galaxy formation in low-mass halos. Early analytic and semianalytic works (G. Efstathiou 1992; J. S. Bullock et al. 2000; A. J. Benson et al. 2002; R. S. Somerville 2002) demonstrated that heating of the intergalactic medium after reionization raises the Jeans mass, preventing gas accretion into halos below a critical virial temperature (e.g., ~2 × 10⁴ K). A comprehensive review by R. S. Somerville & R. Davé (2015) further highlights how this mechanism naturally limits star formation in halos with $M_{vir} \lesssim 10^9 M_\odot$. However, the interplay of reionization, internal feedback, and environmental effects, especially within Local Group-like systems, still needs to be quantified using cosmological simulations with sufficiently high resolution. For example, the timing of a halo's mass assembly, which is greatly affected by cosmological environments, is known to strongly influence its ability to retain gas and form stars, not just its mass at $z=0$ (T. Okamoto et al. 2020). Understanding this process is crucial for explaining the population of subhalos, and it constitutes the central focus of this work. We hereby focus on the baryonic effects without touching the cosmological framework.

From the observational perspective, it is difficult to address this issue empirically because starless dark halos are virtually unobservable. Furthermore, analytical studies have shown that the gas content of halos is highly sensitive to the timing of reionization (Y. Noh & M. McQuinn 2014). These motivate the need for high-resolution simulations that can track gas retention and star formation across diverse environments. However, as modern hydrodynamical simulations require significantly more computational resources than DM-only simulations, achieving a balance between large cosmological volumes and the high resolution needed to resolve faint dwarf satellites remains challenging. This compromise often limits our ability to simultaneously capture both the development of galaxies in the cosmological context and the detailed structure of dwarf galaxies. This is why previous simulation-based studies have often focused on isolated dwarf galaxies or a

limited number of Milky Way analogs, restricting their ability to explore environmental diversity and statistically significant subhalo populations (F. Munshi et al. 2017; E. Applebaum et al. 2021).

To address this limitation, we employ high-resolution cosmological simulations, NEWHORIZON (Y. Dubois et al. 2021) and NEWHORIZON2 (S. K. Yi et al. 2024), whose details are given in Section 2.1. High resolution is essential to resolve dense star-forming regions, and this is required for studying the formation of faint objects. Our simulations successfully reproduce the morphology of observed spiral galaxies with thin disk structures and small star-forming clumps (S. K. Yi et al. 2024). Also, these cosmological zoom-in simulations cover a spherical volume with a diameter of 20 Mpc, which is an exceptional scale in simulations under 100 pc resolution. This ensures a relatively larger sample size evolving in the cosmological context than in idealized simulations that target individual galaxies. Thus, we can effectively capture the formation and evolution of small-scale dwarf galaxies. Furthermore, these simulations incorporate a comprehensive treatment of baryon physics, which is crucial for our analysis (see Y. Dubois et al. 2021 for details).

The remainder of the paper is structured as follows. In Section 2, we provide a detailed description of our simulations, including the methods used for halo identification and subhalo classification, incorporating information from the merger trees. Section 3 presents our main findings, which are further discussed in Section 4, where we compare our results with previous studies and explore their broader implications. Finally, in Section 5, we summarize our key conclusions and outline potential directions for future research. Throughout this work, stellar masses have been computed assuming a Chabrier (G. Chabrier 2003) initial mass function, and all masses are $h$-corrected.

## 2. Methods

In the following, we provide a more detailed description of our set of simulations, outlining their key characteristics. Given the critical role of halo and subhalo identification in our analysis, we also present a comprehensive explanation of the method used to detect them. Additionally, we describe the approach adopted to establish a connection between halos, subhalos, and their associated galaxies (if any), highlighting any matching criteria applied in our study.

### 2.1. *Numerical Simulations*

Here we describe detailed information about our simula-
tions: NEWHORIZON (NH; Y. Dubois et al. 2021) and
NEWHORIZON2 (NH2; S. K. Yi et al. 2024). NEWHORIZON
is a follow-up to the cosmological hydrodynamical simulation
Horizon-AGN (Y. Dubois et al. 2014), which was run using
the adaptive mesh refinement code, RAMSES (R. Teyssier
2002). Horizon-AGN provided numerous galaxy samples over
a large volume of (140 Mpc)³, but its kiloparsec-scale
resolution was insufficient to resolve faint dwarf galaxies. To
achieve higher resolution, NEWHORIZON targets a smaller
spherical volume of 20 Mpc diameter using the zoom-in
technique. The best spatial resolution is 34 pc, and the stellar
(DM) mass resolution is $10^4 M_\odot$ ($10^6 M_\odot$). With this high
resolution and modern astrophysical prescriptions, NEWHOR-
IZON successfully resolves small-scale galactic features such
as thin disks, star-forming clumps, and dwarf galaxies (Y. Dubois et al. 2021; G. Martin et al. 2022; J. Rhee et al. 2024; S. K. Yi et al. 2024). Although the detailed subgrid astrophysics is given by Y. Dubois et al. (2021), we want to introduce key physics very briefly. For the star formation, each leaf cell should exceed the gas density threshold of 5 H cm⁻³ (NH2) or 10 H cm⁻³ (NH). Then, stars are finally formed if the cell satisfies the gravoturbulent condition (i.e., local gravity and thermal/turbulent pressure; C. Federrath & R. S. Kles sen 2012). Our simulations follow the mechanical feedback scheme of T. Kimm et al. (2015) for Type II supernovae, which considers the energy release and momentum transfer from massive stars. We also have the radiative heating source from the uniform UV background following the F. Haardt & P. Madau (1996) model, which is turned on at z = 10.

To enhance the statistics in our main work, we also adopt a NEWHORIZON2 simulation that is a twin simulation of NEWHORIZON with the same initial condition but with poorer spatial (68 pc) and stellar mass (2 × 10⁴ M☉) resolutions. In addition to its lower resolution, NEWHORIZON2 has several minor differences. First, NEWHORIZON2 has a 50% boosted stellar feedback efficiency relative to NEWHORIZON, in the expectation that this improves the agreement of star formation in low-mass halos. Second, NEWHORIZON2 employs a more sophisticated method for selecting the zoomed-in region. Like NEWHORIZON, if we choose a simple spherical zoom-in region, the outer boundary of this region suffers from artificial boundary effects caused by low-resolution DM particles (i.e., so-called “contamination”). Since even a few low-resolution particles could artificially distort the gravitational potential in halos, we had to discard many samples due to contamination in NEWHORIZON. Instead of a simple spherical region, NEW- HORIZON2 selects the zoom-in region using the convex hull of the DM particles traced back to the initial condition. This method significantly reduces contamination, which is espe cially crucial for small-scale structures. Lastly, NEWHORI- ZON2 tracks the evolution of nine chemical elements (H, C, N, O, Mg, Si, S, Fe, and D). Although chemical evolution is not the main topic of this work, this information has been used in galactic studies, such as the decomposition of the thin and thick disks (S. K. Yi et al. 2024). Despite sacrificing resolution, many advanced techniques in NEWHORIZON2 require significantly more computational resources. We over come this issue by using the OpenMP and Message Passing Interface hybrid parallelization version of the code, RAMSES- yOMP (S. Han et al. 2025).

We also want to emphasize that we have a dense output cadence of ~1000 snapshots with an average time interval of 15 Myr for both simulations. This enables us to capture short- timescale events (e.g., the reaction of surrounding gas to stellar and AGN feedback) and detailed evolutionary tracks even in violent interactions with neighboring objects.

## 2.2. Halo (Galaxy) Identification and Sampling

The identification of galaxies and their associated halos is an important task in the context of this work. Therefore, we provide here a detailed description of our identification method. We separately identify DM halos and galaxies using the AdaptaHOP algorithm based on the particle density (D. Aubert et al. 2004; D. Tweed et al. 2009). The minimum particle number is set to 100, corresponding to a mass of ~10⁶ M☉ for galaxies and ~10⁸ M☉ for halos. The numbers of

<figure><img src="image_4.png" /><figcaption>**Figure 1.** Mass functions of identified galaxies (blue) and halos (gray) in NEWHORIZON and NEWHORIZON2 at z ≈ 0.17. The observed galaxy stellar mass functions are represented by black circles and diamonds (C. Li & S. D. M. White 2009; I. K. Baldry et al. 2012). Theoretical halo mass functions are shown as black solid and dashed lines (R. K. Sheth & G. Tormen 1999; A. Lapi et al. 2013).</figcaption></figure>

identified galaxies and halos are 4070 (1528) and 29,333 (29,176) for NEWHORIZON (NEWHORIZON2), respectively. The difference in the number of galaxies between the two simulations originates from the stellar mass resolution, which significantly affects the lowest mass of galaxies. We find that the number of DM halos is preserved due to the same DM mass resolution. The large discrepancy between the numbers of galaxies and halos arises not only from the regulated star formation but also from the mechanism of the zoom-in technique. By definition, zoom-in simulations prohibit star formation outside the zoom-in region, leaving only low-resolution DM and gas cells in those outer regions. We exclude contaminated halos from our analysis. We also present the mass function of the halos and galaxies with the comparison of observational results (C. Li & S. D. M. White 2009; I. K. Baldry et al. 2012) in Figure 1. The galaxy stellar mass functions from our simulations align well with observations, though the highest mass bins show fluctuations due to a lack of samples.

We then selected Milky Way analogs. Even though many different criteria are exercised to define Milky Way analogs (L. A. Bignone et al. 2019; J. Van Nest et al. 2023; A. Pillepich et al. 2024), we apply a simple and generous mass criterion since our simulations do not contain many massive halos. We adopt the mass estimation from the observations: the central galaxy mass of (6.08 ± 1.14) × 10¹⁰ M☉ (T. C. Licquia & J. A. Newman 2015). We double the mass range (i.e., 2σ) to secure more samples. To ensure that the host halo mass is around ~10¹² M☉ (T. M. Callingham et al. 2019), we exclude samples whose host halo mass lies outside the logarithmic range [11.65, 12.35]. In total, 26 Milky Way analogs are selected in NEWHORIZON and NEWHORIZON2 (13 in each simulation) at the last snapshot (z = 0.17). We define Milky Way analog systems (MWA systems) as the DM halos hosting the selected Milky Way analogs.
<figure><img src="image_5.png" /><figcaption>**Figure 2.** Overview of MWA systems and satellite statistics in our simulations. Panel (a) shows the large-scale distribution of DM (grayscale) and stars (reddish) in part of NEWHORIZON. The identified MWA systems are marked with yellow circles, indicating the virial radius $R_{200,c}$ and mass $M_{200,c}$. Panel (b) zooms in to a region within $R_{200,c}$ of one example system, illustrating the distribution of DM and subhalos. Yellow circles are identified subhalos in this system, with sizes corresponding to their virial radii. Solid circles represent more massive subhalos ($M_{vir} > 5 \times 10^8 M_{\odot}$), while dotted circles represent less massive ones. Panel (c) shows the stellar distribution for the same system, with blue circles marking the identified satellite galaxies. Similarly, the line style distinguishes satellites with stellar masses above (solid) or below (dotted) $5 \times 10^6 M_{\odot}$. Panel (d) displays its mock SDSS gri composite image viewed face-on.</figcaption></figure>

To investigate the subsystems (i.e., subhalos and satellite galaxies), we find all subhalos and satellite galaxies within $1.5R_{200,c}$ in the MWA systems, where $R_{200,c}$ is the radius within which the enclosed density is 200 times the critical density. From them, we want to exclude any misidentified substructures that are likely transient, such as star-forming clumps or fragmented structures. To do this, we track the progenitors for the past 1 Gyr and calculate a score based on the number of shared member particles. Using this merit score, we exclude all galaxies or halos that suddenly appear (fragmented) or disappear (transient).

Figure 2 shows an example image of MWA systems in NEWHORIZON. We emphasize that these systems evolve in cosmological environments along large-scale filaments (panel (a)), and the subkiloparsec-scale structures, such as spiral arms or dwarf satellites, are resolved for individual galaxies (panel (d)). In panels (b) and (c), we show the zoom-in images of one example system with a box size of $2R_{200,c}$ and mark identified subhalos and satellite galaxies on the DM and stellar density maps. In line with recent studies (M. Wadepuhl & V. Springel 2011; A. R. Wetzel et al. 2016; G. Lee et al. 2024), the satellite galaxies in panel (c) are far fewer than the subhalos in panel (b). In panel (d), we also present the mock Sloan Digital Sky Survey (SDSS) gri composite image of the central MWA using the Monte Carlo radiative transfer code, SKIRT (M. Baes & P. Camps 2015; P. Camps & M. Baes 2020).

## 2.3. Classification of Subhalos

We now classify subhalos into two populations: “starless” and “starred.” Starless subhalos are subhalos without a galaxy, and starred subhalos are those with one. However, since we identified the halos and the galaxies separately, we must match them. First, we identify the galaxy candidates hosted by each subhalo. Ideally, we expect that the galaxy should be located at the center of the host halo, which is generally true in most cases, mainly for massive and isolated galaxies. However, especially for low-mass galaxies in dense environments, which are the main targets of this work, the galaxy could be offset from the center of the host halo or even outside the virial radius. For instance, when a small subhalo with a galaxy orbits in a larger halo potential, the distribution of DM can get diffused, and its virialized region can shrink, while the galaxy maintains its distribution. This could lead to an offset between the centers of the subhalo and the galaxy, which makes it hard to detect their connection simply based on catalog data. Because of these offsets, we use the following criteria to ensure a robust association between halos and galaxies:

1. The center of the halo is in the galaxy: ($d < R_{gal}$) ∧ ($R_{gal} < R_{hal}$).

2. The galaxy is in the halo: $d < R_{hal, vir} - R_{gal}$.

3. The galaxy and the halo are overlapping ($d < 0.5 \times (R_{gal} + R_{hal, vir})$) ∧ ($R_{gal} < R_{hal, vir}$),
<figure><img src="image_6.png" /><figcaption>**Figure 3.** Example starless and starred subhalos in our systems at $z = 0.17$, the final snapshot. Panels (a) and (b) show the images of sample starred and starless subhalos, respectively. The DM particles are shown as a smoothed density map in the background, and the stellar particles are overlaid in orange. White circles indicate $R_{vir}$ computed from the halo finder.</figcaption></figure>

where *d* is the distance between the center of the halo and the galaxy, $R_{gal}$ and $R_{hal}$ are the radii of the galaxy and the halo, respectively, and $R_{hal,vir}$ is the virial radius of the halo. We note that our halo finder, HaloMaker, provides a quantity “r” which represents the maximum distance between the center and the farthest member particle. Another quantity, the virial radius “rvir,” is determined by satisfying the virial theorem through direct computation. $R_{gal}$ and $R_{hal}$ indicate an *r* to represent spatial occupation based on the maximum distance of member particles. $R_{hal,vir}$ is an rvir that is a boundary of virialization. Stellar, cold gas, and total gas masses in the main text are measured within rvir. If a halo does not satisfy the first criterion (i.e., no matched galaxy), we sequentially check the second and third criteria. We then classify subhalos as starless or starred based on the presence of a matched galaxy. If there are multiple galaxy candidates, we select the galaxy closest to the center of the halo.

After this procedure, we track their merger trees to validate their classified labels. We describe the detailed process of building the merger tree and further classification in

Appendices A and B, respectively. Finally, we have 2032 starless subhalos and 416 starred subhalos. Representative images of the two types of subhalos at $z = 0.17$ are shown in Figure 3. Despite their similar virial masses, only the starred subhalo hosts stars (orange).

### 3. Results

In this section, we present the main results of our analysis, demonstrating that the so-called missing satellite problem does not arise in our simulations, as the majority of subhalos do not host a galaxy. Building on this finding, we then further explore the underlying reasons why this population of starless subhalos remains devoid of galaxies, examining the physical processes that may be responsible for their evolution and fate within the host halos.

#### 3.1. Satellite Abundance and Star Formation

We present the cumulative number of satellite galaxies and subhalos as a function of mass in Figure 4(a). Statistically, one can see clearly in panel (a) that the sequences of satellite galaxies (blue) and subhalos (gray) are well separated. Especially for our satellites, the cumulative numbers (blue line) match well with the observed numbers of Milky Way and M31 satellites (black lines) of J. A. O’Leary et al. (2023). We note that recent surveys, such as DES, Pan-STARRS, and SAGA, have extended detection to fainter satellites or to other galaxy systems (A. Smercina et al. 2018; E. O. Nadler et al. 2020; Y.-Y. Mao et al. 2021; M. Nashimoto et al. 2022). We decided to use only the Local Group sample. This is, first of all, because different galaxy data show different star formation properties (A. Smercina et al. 2018; Y.-Y. Mao et al. 2021). Meanwhile, the Local Group sample is sufficient to tell us that there are far more starless subhalos than starred ones in the real Universe, which is the main motivation for our study. On a different account, we did not use fainter satellite observations than those used in this study, because our simulations do not reach so faint anyway. This reproduction of satellite abundance aligns with results reported previously from other models (T. Sawala et al. 2016; A. R. Wetzel et al. 2016; C. Engler et al. 2021; M. Jung et al. 2024; G. Lee et al. 2024), further supporting the alleviation of the classical missing satellite problem for this stellar mass range.

The star formation properties of our satellite samples are also aligned with the observational data. In Figure 4(b), we measure the star formation rate of our samples and observational results and classify satellite galaxies as either star-forming or quenched using the cut in specific star formation rate (S. Tacchella et al. 2019). We find that the star-forming satellites (blue circles) are well located on the star-forming main sequence (gray lines; D. Elbaz et al. 2007; T. H. Jarrett et al. 2017; A. Saintonge & B. Catinella 2022). The trend that a significant fraction of the satellites are quenched (red circles) is also well matched with the observed quenched fraction in the Local Group (A. R. Wetzel et al. 2015), although it is much higher than the result of the SAGA survey (Y.-Y. Mao et al. 2021).

Since our samples seem to agree largely with the general understanding of not only abundance but also star formation properties, we now want to focus on the individual subhalos to understand their stellar contents.
## 3.2. Distinct Subhalo Populations: Starless and Starred

The effect of baryons seems to play a crucial role in resolving the missing satellite problem. Their impact extends beyond just reducing the stellar mass of satellite galaxies and even influences their existence. From the large discrepancy in mass functions of subhalos and satellite galaxies, we see that only a few subhalos host galaxies. We call the subhalos without significant stellar components “starless,” and those with prominent satellite galaxies “starred,” as mentioned in Section 2.3.<fnref n="8" />

In Figure 5, the statistical comparison of the two subhalo types reveals a clear difference in their overall mass distribution (panel (a)). In general, starred subhalos have larger masses than starless subhalos. More massive halos tend to attract more gas, leading to increased star formation and the development of starred subhalos. For the least massive subhalos, some of them barely exceed the minimum threshold for identification, making them unreliable. Nevertheless, there is a range of common masses between the two populations. We select the subhalos of $M_{vir} = 10^{8.4-8.9} M_{\odot}$ where the $1\sigma$ ranges of the two mass functions overlap. We now examine the subhalos in this common mass range to find what makes a subhalo starred or starless, setting aside the effect of halo mass.

We measure the fractional mass distribution of subhalo components, such as stars, total gas, and cold gas (below $10^4$ K, regardless of hydrogen species). Panel (b) shows the dramatic difference in stellar content between the two subhalo types. A small fraction of starless halos have stars (red peak near $f > 0$). In most cases, these stars are interlopers from larger neighboring galaxies that temporarily pollute these subhalos. Despite this large difference, the two populations contain similar amounts of gas (panel (c)), which is very amusing. The two populations are well separated in their cold gas content (panel (d)). Starless subhalos completely lack cold gas, indicating that star formation is totally absent. Then, the key question is whether starless subhalos have ever hosted galaxies or whether they have lost them during their evolution.

We examine two of the most widely discussed mechanisms that may explain the lack of cold gas during the evolution of starless subhalos: supernova feedback and environmental effects associated with the orbital motion of subhalos. These two are representative processes for low-mass galaxy evolution and quenching, namely “nature or nurture” (Y. J. Peng et al. 2010; M. L. Balogh et al. 2016; B. Darvish et al. 2016; L. Kawinwanichakij et al. 2017; E. Contini et al. 2020; J. Rhee et al. 2020). The final mass of a subhalo does not necessarily reflect its entire growth history. Therefore, to minimize the influence of different assembly histories, we sample subhalos with similar evolution using peak masses.<fnref n="9" /> Using these subsamples, we first investigate the significance of the two candidate processes mentioned above for removing cold gas.

<figure><img src="image_7.png" /><figcaption>**Figure 4.** Number and star formation rates (SFRs) of satellite galaxy samples. (a) The cumulative number of satellite galaxies (blue) and subhalos (gray) in our MWA systems in NEWHORIZON and NEWHORIZON2, plotted with respect to the stellar mass of satellites and the virial mass of subhalos, respectively. Arrow markers at $5 \times 10^6 M_{\odot}$ (galaxy) and $5 \times 10^8 M_{\odot}$ (halo) correspond to ~500 particles. Thicker lines are the median, and shaded regions are the 16th–84th percentiles. Black lines represent observed values for the Milky Way and M31 (stellar mass extracted from Table A1 of J. A. O’Leary et al. 2023). (b) The star formation rate and the stellar mass of the satellite galaxies. The blue circles exhibit the star-forming satellites, the red circles are the quenched satellites, the red triangles are those with zero SFR, and the black dots are the Milky Way analogs in our simulations. The gray lines are the observational star-forming main sequence (D. Elbaz et al. 2007; T. H. Jarrett et al. 2017; A. Saintonge & B. Catinella 2022), and the gray squares are the observed values of the Milky Way and M31 (S. Rahmani et al. 2016; G. Natale et al. 2022).</figcaption></figure>

### 3.2.1. Supernova Feedback

First, supernova feedback is expected to be effective for removing cold gas through energy and momentum injection. The surrounding gas can be heated by the release of supernova energy and, if loosely bound, expelled from the subhalo, resulting in the suppression of star formation (R. B. Larson 1974; A. Dekel & J. Silk 1986; Y. Dubois & R. Teyssier 2008; T. Kimm et al. 2018). To investigate whether supernova feedback had a more significant impact on starless subhalos, we measure $P_{SN}$, the expected value of supernova energy for any time step when supernovae occur, in Figure 6(a). The

<footnote n="8">⁸ A main criterion to separate them is sustainability to host stars. Even if there are some stellar particles within the subhalo in question, we classify it as a starred subhalo only when it continues to contain the stars for enough time, excluding a temporary invasion from neighboring galaxies. The detailed procedure is explained in Appendix B.</footnote>

<footnote n="9">⁹ We use the same method previously described for final mass binning (Figure 5(a)) to determine the common mass range between starred and starless subhalos based on their peak masses. We apply a mass range between the 8th and 92nd percentiles to ensure a sufficient number of samples.</footnote>
<figure><img src="image_8.png" /><figcaption>**Figure 5.** Structural and compositional differences between starless and starred subhalos in our systems at $z = 0.17$, the final snapshot. (a) Their mass functions. Arrows in the distribution indicate the $1\sigma$ ranges. The overlapping range (gray band) of $10^{8.4} \lesssim M_{\text{vir}}/M_{\odot} \lesssim 10^{8.9}$ is used to select subhalos for analysis in panels (b), (c), and (d). The black dotted vertical line indicates the minimum threshold for halo identification ($1.3 \times 10^8 M_{\odot}$). Panels (b), (c), and (d) display the cumulative fractional distributions of stars, total gas, and cold gas, respectively. The fraction of subhalos containing each component is also shown. Note that the starless subhalos have comparable amounts of gas (panel (c)) but lack stars (panel (b)), probably because their gas was not cold and dense enough (panel (d)).</figcaption></figure>

detailed calculation is explained in Appendix C, but for instance, a high value of $P_{\text{SN}}$ indicates that a burst of supernovae occurred over a short period of time.

As is visible in Figure 6(a), all the starred subhalos have large values of $P_{\text{SN}}$. In order to check its impact on cold gas content, we show the cold gas fraction before and after the supernova events (vertical line) in panel (b). The cold gas mass indeed slightly decreases after the supernova explosion in the starred subhalos (blue), as expected. Nevertheless, these subhalos succeed in retaining a significant amount of cold gas afterward. Supernova feedback probably does not affect existing stars of starred subhalos, but it may affect their growth by reducing future star formation rates.

Another striking feature in panel (a) is that most (93.8%) of the starless subhalos do not experience any supernova events. In fact, they do not have any in situ star formation. Therefore, supernova feedback is irrelevant as a process for removing cold gas in most starless subhalos.

Some starless subhalos have comparable values of $P_{\text{SN}}$ to starred subhalos in panel (a). Before examining them, we would like to remind the reader that some subhalos can be influenced by supernovae from neighboring halos. To account for this, we classify starless subhalos with supernova events into internal and external cases based on the true ownership of young stars. When inspecting the star formation histories, we find that starless subhalos with external supernovae alone
<figure><img src="image_9.png" /><figcaption>**Figure 6.** The first main possible scenario responsible for the removal of the cold gas in certain halos: supernova feedback. (a) The distribution of supernova energy rate from the starred (blue) and starless (red) subhalos. We also present the percentage of starless subhalos that do not have any supernova events (93.8%). (b) The time evolution of the cold gas mass fraction, with the x-axis representing the time offset from the supernova explosion. The black vertical line indicates the epoch of the supernova explosion. Thin solid lines are the evolution of individual subhalos, while thick lines and shaded regions represent the median values and the 16th–84th percentiles, respectively. Colors represent different populations of supernova-hosting subhalos: blue for starred subhalos, green for starless subhalos with internal supernovae, and red for starless subhalos with external supernovae from neighboring halos.</figcaption></figure>

(N = 14) have indeed not formed their own stars, like most starless subhalos. Panel (b) shows that they (red line at y = 0) had no cold gas, regardless of supernova events. Thus, the case of external supernovae is also irrelevant as a test for the effect of supernova feedback.

A small number (N = 8) of starless subhalos have internal supernova events. As shown in panel (b), they lose cold gas after their supernova events, and their cold gas decays faster than in their starred counterparts. These eight out of 372 starless subhalos (2%) could be said to have become starless due to supernova feedback, but are too few to represent the starless subhalos.

In summary, supernova feedback is not the main driver of the removal of cold gas from starless subhalos.

### 3.2.2. Environmental Effects Associated with Orbital Motion

We now investigate the orbital motion-associated environmental effects on subhalos. Specifically, Local Group environments could remove cold gas from satellite subhalos through ram pressure stripping or even directly strip stars through tidal stripping (J. E. Gunn et al. 1972; D. O. Richstone 1976; Y. M. Bahé & I. G. McCarthy 2015; E. Contini et al. 2024; J. Rhee et al. 2024), potentially driving their transition to starless subhalos. To test this, we examine the evolution of their components using the time offset from the first infall normalized by the crossing time. In this normalized unit, key orbital stages, such as first infall (at 0), pericenter passage (at 1), and apocenter passage (at 2), are conspicuous. Figure 7(a) shows the evolution of DM mass, which is similar across the two populations, as expected from similar assembly histories of the sampled subhalos. Before infalling into the MWA system, both populations exhibit a similar flat trend. Following the first infall, both experience a decrease in DM due to tidal stripping. We present the stellar and cold gas fractions relative to the total mass in panels (b) and (c), respectively.

It is inappropriate to discuss the environmental effects on starless subhalos because they hardly had any cold gas or stars

in the first place. The peak at the first pericenter passage (x ≈ 1) in Figure 7(b) shows a transient feature of having interloper stars from neighboring halos, primarily from their central Milky Way analogs. In the case of starred subhalos, their stellar fractions are not significantly affected by the orbital motion. Aside from the transient bump near the first pericenter passage, the stellar fraction appears similar or even elevated. This can be understood by the fact that DM stripping precedes stellar stripping due to differences in spatial distribution and concentration (E. Contini et al. 2014, 2023; R. Smith et al. 2016; J. K. Jang et al. 2024). So, the environmental effects on stars are negligible. On the other hand, the cold gas of starred subhalos is dramatically reduced through the infall due to ram pressure stripping (panel (c)). This explains why starred subhalos are largely deficient in cold gas in the final stage (Figure 5(d)). In conclusion, for the case of starred subhalos, environmental effects are effective on cold gas, removing it quickly and quenching star formation, but they do not remove preexisting stars. As a result, they cannot transform starred halos into starless ones.

We conclude that environmental effects associated with a subhalo's orbital motions are not responsible for the absence of stellar components in starless subhalos either. What, then, caused starless halos to be starless?

### 3.2.3. Born to be Starless

The prescription for star formation in our simulations follows the gravoturbulent scheme (C. Federrath & R. S. Klessen 2012) that considers local gravity and thermal/turbulent pressure (Y. Dubois et al. 2021). These criteria are applied to dense gas cells with a density threshold of 10 H cm⁻³ (5 H cm⁻³ for NEWHORIZON2). However, in most starless subhalos, these dense cells rarely exist for the following reasons.

UV background is the dominant mechanism preventing cooling and subsequent star formation in starless subhalos. Since most starless subhalos are very low in (virial) mass, they are sensitive to UV heating (J. F. Navarro & M. Steinmetz 1997;
<figure><img src="image_10.png" /><figcaption>**Figure 7.** The second main possible scenario responsible for the removal of the cold gas from certain halos: the environmental effects associated with the orbital motion. (a) The evolution of DM mass around the first infall into the proto-MWA systems. The x-axis is the time offset from the first infall, normalized by the crossing time. The black vertical line marks the first infall epoch. (b) The stellar mass fraction. (c) The cold gas fraction.</figcaption></figure>

N. Y. Gnedin 2000; M. Pereira-Wilson et al. 2023). The uniform UV background in our simulations after z=10 heats and evaporates gas, therefore suppressing star formation (see the visualized temperature map in Figure D1 in Appendix D). To overcome this heating, halos must, in the first place, have a sufficiently deep potential and host gas with a density of at least 0.01 H cm⁻³, the self-shielding threshold (J. Rosdahl & J. Blaizot 2012; Y. Dubois et al. 2021). Subhalos that sustain dense gas can enrich their gas metallicity through star formation, enabling additional metal-cooling. While most starred subhalos contain dense and cold gas before reionization, starless subhalos fail to reach this density threshold. Thus, their gas is rapidly heated from early times.

The failure of starless subhalos to have dense gas may be related to the conditions of their birthplace before reionization. To verify this, we investigate the differences in the long-range environment. In Figure 8, we measure various initial evolutions within a box of 100 comoving kpc centered on the birthplace of each subhalo. The masses inside the box are nearly identical at birth (panel (a)), so both populations appear to form in similar environments with comparable matter content. However, the mass growth rates of starless subhalos are slightly lower than those of starred subhalos, leading to a divergence in total mass over time. The birth sites of starred subhalos exhibit higher matter accretion rates (panels (b) and (c)), indicating that these regions gather matter more efficiently.

Higher accretion rates in starred subhalos also drive rapid growth in gas density (Figure 9(a)). Similar to bulk gas growth shown in Figure 8(a), the average gas densities (dotted lines) of the two populations in Figure 9(a) start at similar levels and slowly diverge over time due to different growth rates. The maximum gas densities (solid lines) show more dramatic divergence early on. After the initial density decreases, both the maximum density and the average density in the virialized region (dashed lines) increase in starred subhalos. Thanks to more efficient accretion, gas in starred subhalos can exceed the self-shielding density (i.e., 0.01 H cm⁻³) before the reionization becomes prevalent at z ~ 7. This dense gas can resist UV heating, and further cool and collapse, eventually reaching the

star-forming density (i.e., 5–10 H cm⁻³). By contrast, starless subhalos show a stagnant level of gas density far below the self-shielding threshold. They do not have dense enough gas to self-shield even through the completion of reionization.

Reionization is completed at z ~ 7 in our simulations, which is aligned with the sudden decrease in the cold gas mass in Figure 8(a). During the reionization transition period (i.e., 7 < z < 10), the directional difference between the two populations becomes noticeable. The temperature (Figure 9(b)) of the densest gas in starred subhalos peaks at ~10,000 K in this period, and it then cools down marginally but steadily. In contrast, starless subhalos show a gradual increase in temperature, and even for the densest cells, the trend of temperature evolution mirrors the average temperature in the box. In Figure 9(c), we can see the impact of reionization more explicitly by measuring the change in gas temperature caused only by radiative cooling and heating, excluding adiabatic effects. The key difference begins near the onset of reionization (z = 10). Before the completion of reionization, gas in starred subhalos becomes dense enough to self-shield, suppressing radiative heating and enabling continuous cooling. In contrast, gas in starless subhalos remains too diffuse (panel (a)) to self-shield, leading to inefficient radiative cooling. After the UV background becomes fully prevalent, starless subhalos undergo a dominantly heating phase, while starred subhalos exhibit a strong net cooling phase (panels (b) and (c)).

The fate of subhalos seems to be determined at birth. The difference in birth sites may reflect variations in the subhalo's connection to the cosmic web (e.g., M. Musso et al. 2018). While further analysis is required to clarify this point, the initial difference appears to be amplified and sealed during reionization. Starless subhalos live their starless fate thereafter.

Building upon the analysis presented in this section, the next section will be dedicated to a comprehensive discussion of our results, placing them in the broader context of previous key studies on the same topic. We will highlight the main similarities and differences between our findings and those reported in the literature, addressing potential reasons for any discrepancies.
<figure><img src="image_11.png" /><figcaption>**Figure 8.** Comparison of the initial evolution between starred (blue) and starless (red) subhalos, measured within a box of 100 comoving kpc centered on their birthplace. (a) The median mass evolution of DM (dashed), gas (solid), and cold gas (dotted). Shaded regions are the 16th–84th percentiles for each ingredient. The gray vertical line indicates the onset of reionization in our models. To visualize the end of reionization, we color the background of panels in grayscale (e.g., a white background indicates a fully ionized state). We extract the hydrogen ionization fraction for a representative gas cell ($10^4$ K and $0.001$ H cm$^{-3}$) to illustrate the progression of reionization. (b) The median evolution of DM accretion rates in the *box*, which indicates the amount of matter gathering on a large scale. (c) Similar to panel (b) but showing baryon accretion rates in *virialized regions*.</figcaption></figure>

## 4. Discussion

The results obtained in this analysis provide a detailed picture of the relationship between subhalo mass and the ability to host visible galaxies, contributing to our understanding of the properties of dwarf galaxies in the Local Group and the role of baryonic feedback in suppressing star formation in low-mass subhalos. Our data suggest that the abundance and distribution of satellite galaxies are strongly influenced by feedback processes and tidal stripping, in agreement with the findings of T. Sawala et al. (2016). In that study, the authors demonstrated, through hydrodynamic simulations, that baryonic heating and reionization significantly reduce the final mass of subhalos, making them more vulnerable to tidal destruction. This insight builds on earlier works, such as those by G. Efstathiou (1992), A. A. Thoul & D. H. Weinberg (1996), and A. J. Benson et al. (2002), which proposed that reionization prevents gas from cooling and condensing in halos below a mass threshold determined by the virial temperature of the heated gas. Our results confirm that including baryonic effects can explain the relatively low number of observed satellites in the Milky Way compared to the predictions of purely gravitational simulations. However, as noted by

<figure><img src="image_12.png" /><figcaption>**Figure 9.** Similar evolutionary information to Figure 8, but focused on gas. (a) The evolution of the gas density for each population. Solid lines represent the maximum density, while dashed (dotted) lines are the average density within the virialized region (box). We also indicate the density thresholds for star formation ($10 \text{ H cm}^{-3}$ for NEWHORIZON and $5 \text{ H cm}^{-3}$ for NEWHORIZON2; yellow shading) and self-shielding ($0.01 \text{ H cm}^{-3}$; black horizontal line). (b) Similar to panel (a) but depicting the evolution of temperature. Solid lines are the temperature of the densest gas cells, while dashed (dotted) lines are the average temperature in the virial region (box). (c) The net change in temperature of the top 1% densest gas cells, derived only from the *radiative cooling and heating* rates.</figcaption></figure>

A. Benítez-Llambay et al. (2019) and A. Benitez-Llambay & C. Frenk (2020), this effect is compounded by gas cooling and dynamical instabilities, which further influence subhalo survival.

One of the most significant aspects emerging from our study concerns the mass distribution of subhalos that host luminous galaxies. While previous work, such as that by C. M. Simpson et al. (2018)—based on hydrodynamic simulations—suggested that the visibility of subhalos is determined by reionization and star formation feedback, we find that the picture is more nuanced. Thanks to the higher mass and spatial resolution of our simulations, we observe that subhalos with and without stars occupy broad and overlapping mass ranges below the suggested critical mass threshold of $M_{\text{halo}} \sim 10^9 M_{\odot}$ (Figure 5). This suggests that satellite survival is not governed by a sharp mass threshold. Instead, the formation of luminous galaxies in low-mass subhalos is likely influenced by their early assembly history, not merely their halo mass at $z=0$. This behavior has also been reported in zoom-in simulations by A. Fitts et al. (2017), where late-forming and low-concentration halos exhibit a total absence of stars despite reaching masses comparable to luminous counterparts.
In this context, it is important to note the seminal
contributions from earlier studies on the effect of modeled
reionization on star formation in halos. For example, A. Bení-
tez-Llambay et al. (2015) and A. Benitez-Llambay & C. Frenk
(2020) showed that halos just above the cooling threshold can
retain gas in equilibrium with the dark matter potential and UV
background. Those halos exhibit significantly delayed star
formation histories, depending on their trajectory of mass
growth. Recent models have also reported that the halo
occupation fraction depends sensitively on the timing of
reionization (A. Kravtsov & V. Manwadkar 2022) and on the
availability of molecular cooling in early-forming halos
(E. O. Nadler 2025). Our study aligns with those findings,
suggesting that the role of reionization is crucial in determin-
ing whether subhalos can host galaxies. Beyond this, we also
emphasize the importance of the subhalo's birthplace and early
mass assembly history well before reionization.

These conclusions connect with prior research, and we
acknowledge that the key claims of this paper—(1) that
reionization plays a crucial role in the “missing satellites”
problem and (2) that a subhalo’s birth environment and mass
accretion history determine whether it remains dark—are not
entirely novel, even if approached from a different angle. The
depth of our analysis, while significant, could benefit from
statistically robust results across a larger sample of MWAs,
allowing for a more detailed exploration of how self-shielding
and gas properties influence the outcome from diverse
environments.

Despite these connections with previous literature, this
study presents some methodological differences and new
perspectives. For example, our simulations include more
detailed modeling of stellar feedback than previous works,
which could explain some discrepancies in the relative
abundances of satellite galaxies. However, the limited
resolution of our simulations—particularly the absence of
radiative transfer modeling for ionizing photons—prevents a
more accurate treatment of self-shielding and high-density
regions. This issue, along with the neglect of local stellar
radiation, suggests that further exploration of the self-shielding
mechanisms and their impact on star formation would
strengthen the findings.

Furthermore, while previous studies have primarily focused
on isolated dwarf galaxies or a limited number of Milky Way-
like simulations (F. Munshi et al. 2017; E. Applebaum et al.
2021), our approach has been applied to a larger sample of
halos with similar masses, allowing us to evaluate object-to-
object variations and better understand the role of the local
environment alongside reionization in satellite survival within
the cosmological context. A more detailed exploration of how
these environmental factors, combined with the modeling
choices (e.g., properties of the interstellar medium and star
formation), influence the results would be a beneficial next
step for future work.

Let us discuss the possible shortcomings of the present
study. One of the main constraints is the resolution of the
simulations, which, although sufficient to resolve the internal
structure of the most massive subhalos, may underestimate the
effect of feedback in smaller galaxies. Additionally, the
adopted star formation model might not fully capture the
stochastic variations in star formation within marginal-mass
subhalos, suggesting that the results should be interpreted with
caution in this regime. However, these limitations do not

invalidate the general conclusions but rather highlight areas
where future research could improve our understanding of the
dynamical evolution that affects the occupation fraction of
subhalos and the faint end of luminous satellites.

Given these results, several directions for future studies
emerge. First, extending the analysis to higher-resolution
simulations would be useful in order to achieve a better
characterization of the population of lower-mass satellites.
Furthermore, a direct comparison with observational data, such
as those from the Dark Energy Survey 5 yr supernova analysis
(M. Vincenzi et al. 2024) or the Large Synoptic Survey
Telescope (Ž. Ivezić et al. 2019) surveys, could provide even
tighter constraints on feedback models and satellite survival.
Finally, future studies could explore the role of nonthermal
processes, such as cosmic-ray heating or the influence of the
galactic magnetic field, which might have a nonnegligible
impact on the dynamics of subhalos and their ability to host
visible galaxies.

We also acknowledge that our results depend on the
specifics of the subgrid physics used in the simulation. For
example, our simulations lack full radiative transfer, which is
crucial for tracing the local ionization and, in turn, gas
temperature. Recent studies based on radiative hydrodynamic
simulations have indeed suggested that the reionization
suppresses the growth of gas mass and subsequent star
formation in low-mass halos (J. Oñorbe et al. 2015; P. Ocvirk
et al. 2020; M. P. Rey et al. 2020). Specifically, direct
computation of ionizing photons could modify our results.
H. Katz et al. (2020) suggested that the hydrogen neutral
fraction is affected by the large-scale structure (i.e., filaments),
and A. Bhagwat et al. (2024) showed that a different scheme of
supernova feedback, incorporating the radiative transfer, can
change the distribution and timing of reionization, and also the
final stellar mass. Using zoom-in simulation, O. Zier et al.
(2025) compared the uniform UV background and the
radiation field extracted from the parent full-box model and
found that the difference in ionization is much more severe in
low-mass halos. Moreover, given that low-mass halos
dominantly contribute to the ionizing photon budget
(J. S. W. Lewis et al. 2020), the ionization state of halos
likely depends not only on the large-scale structure but also on
the local neighborhood, adding complexity to the picture. It
would be exciting if such simulations could reach close to the
local Universe and cover a large cosmological volume in the
near future. Considering that NEWHORIZON consumed nearly
100 million CPU hours of computing time, it is not feasible to
run a similarly high-resolution simulation with radiative
transfer. However, if what we found in our study is largely
true, that is, if the fate of starred and starless subhalos is
determined so early, future studies could extend this work by
exploring the epoch of reionization more thoroughly and
monitoring galaxies during the first few gigayears.

5. Conclusions

Recent hydrodynamical simulations have shown that the use
of updated baryonic physics suppresses the number of
luminous satellite galaxies, tending in the direction of solving
the classical missing satellite problem. In this work, we further
investigate the origin of the observed dichotomy between
luminous and dark subhalos using the NEWHORIZON and
NEWHORIZON2 simulations. The high-resolution simulations
indeed roughly reproduce the number of satellite galaxies in
the Local Group. We further explored the detailed origin of “starless” subhalos. Following the dark halo mass function, less massive halos are more abundant. In our simulations, starless subhalos are substantially less massive than starred ones. We selected a sample of subhalos in a narrow mass range that contains both starred and starless subhalos and attempted to pin down what makes them starred or starless.

The key difference lies in their birthplace and UV background heating, probably associated with reionization. When small subhalos are born in less gas-attracting regions, the gas within them struggles to condense. Their fate is sealed by reionization, as only starred subhalos can maintain dense gas, self-shield, and continue to form stars. The two most widely suspected candidate processes, i.e., supernova feedback and the environmental effects associated with orbital motion, were found to be irrelevant, primarily because starless subhalos never experienced in situ star formation in the first place. Moreover, these processes remove gas but not stars, and thus cannot transform star-forming subhalos into starless ones. In conclusion, starless subhalos are born to be starless, not made.

This study contributes to the existing literature by providing new insights into the relationship between subhalo mass and the survival of satellite galaxies, with direct implications for the missing satellites problem and the modeling of low-mass satellite galaxy formation. The results obtained confirm the trends observed in the works of T. Sawala et al. (2016), C. M. Simpson et al. (2018), A. Benítez-Llambay et al. (2019), and A. Benitez-Llambay & C. Frenk (2020), while also suggesting new perspectives that warrant further investigation. Hence, our results should trigger future research aimed at improving our understanding of the interplay between baryonic feedback, subhalo dynamical evolution, and star formation in Milky Way satellites.

## Acknowledgments

We are particularly grateful to the referee for pointing us to numerous previous studies that were relevant to our investigation. This work was granted access to the HPC resources of CINES under the allocations c2016047637, A0020407637, and A0070402192 by Genci, KSC-2017-G2-0003, KSC-2020-CRE-0055, and KSC-2020-CRE-0280 by KISTI, and as a “Grand Challenge” project granted by GENCI on the AMD Rome extension of the Joliot Curie supercomputer at TGCC. The large data transfer was supported by KREONET, which is managed and operated by KISTI. S.K.Y. acknowledges support from the Korean National Research Foundation (RS-2025-00514475 and RS-2022-NR070872). E.C. acknowledges support from the Korean National Research Foundation (RS-2023-00241934). J.R. was supported by the KASI-Yonsei Postdoctoral Fellowship and by the Korea Astronomy and Space Science Institute under the R&D program (Project No. 2025-1-831-02), supervised by the Korea AeroSpace Administration. This work was partially supported by the Institut de physique des Deux Infinis of Sorbonne Université and by the ANR grant ANR-19-CE31-0017 of the French Agence Nationale de la Recherche.

We thank Maarten Baes for providing useful information for the radiative transfer code, SKIRT, which was used for generating mock images of our simulated galaxies.

## Appendix A
Merger Tree

High resolutions alone do not guarantee accurate detection and classification of halos. Building a merger tree for small halos, especially inside another halo, is not trivial because even very sophisticated merger tree algorithms cannot connect halos if they are not detected in the first place. Besides, when a halo temporarily disappears or gets distorted due to a close encounter or merger, it may be linked to an incorrect progenitor.

Our main results, discussed in Section 3, track the very early stage of the halo evolution, so we need to be particularly careful with the merger tree. First, we utilize the in-house merger tree, PhantomTree (S. Han 2025, private communication), based on the particle membership. We generate a directed graph structure by treating halos as nodes and connections as edges. For each halo (node), we find all descendant candidate halos where the member particles are at the next four snapshots, and then we calculate the score based on the fraction of the received and sent particles:

$$S_{i,j} = \frac{n(i \cap j)}{n(j)} \times \frac{n(i \cap j)}{n(i)}, \qquad (A1)$$

where $S_{i,j}$ is the score between the target halo i and the descendant candidate halo j, $n(i \cap j)$ is the number of member particles shared between them, $n(j)$ is the number of member particles of the descendant candidate halo, and $n(i)$ is the number of member particles of the target halo. We directly connect two halos if they are mutually linked. Otherwise, we connect them to the highest-scoring candidate node, which would be the merger or fragmentation cases.

Second, we rebuild the merger tree using “stable” member particles that are defined when their occurrence exceeds the threshold in the main branch of the tree. This occurrence threshold is initially set to the length of the branch but is iteratively reduced until stable member particles constitute more than 10% of the number of members at the last snapshot. To ensure more reliable connections in the tree, we re-identify the main progenitor, which has the most stable particles in each halo. Based on these stable members, we also recalculate the give and take scores:

$$S_{\text{give},i} = \frac{N_i}{N}, \quad S_{\text{take},i} = \frac{N_i}{n(i)}, \qquad (A2)$$

where $S_{\text{give},i}$ is the “give” score of the target halo i, $S_{\text{take},i}$ is the “take” score of the target halo i, $N_i$ is the number of stable member particles of the target halo i, $N$ is the total number of stable member particles in the main branch, and $n(i)$ is the number of particles in the target halo. We discard the halo if $S_{\text{take}}$ is less than 5% or $S_{\text{give}} \times S_{\text{take}}$ is less than 10%.

We further exclude nodes with masses that deviate significantly from the local and global trends, as described below. We calculate the mean and standard deviation of the mass using the previous 100 nodes and the next 100 nodes in the main branch. If the mass of the node is more than 4σ from the mean or 1000 times the last mass, we tag it as a locally outlying node. Also, if the mass of the node is more than 3σ from the median of the main branch, we tag it as a globally outlying node. We exclude both locally and globally outlying nodes from the main branch.
We then check the position of each node in the main branch. By taking advantage of the dense cadence of the stored snapshots, we can estimate the expected position at the next snapshot. If the node is located beyond the expected position by more than the sum of the virial radii of the two adjacent nodes, we further investigate. For those nodes, we compare the raw positions of the nodes and the polynomial fitted positions. If the offset between the two positions is larger than the standard deviation of the residuals and exceeds 10 comoving kpc, we exclude the node from the main branch. This process results in the reliable merger tree used in this work.

We also extend the merger tree to the first snapshot where no halo is detected by tracking the position of stable particles, which is used in Figures 8 and 9. We acknowledge that this method has inherent uncertainties and cannot entirely avoid incorrect detections or broken tree issues, but repeated checks and strict criteria help minimize unwanted artifacts.

## Appendix B
### Validity of Subhalo Classification

We introduce the classification of subhalos into “starred” and “starless” in Section 2.3 based on matching galaxies and halos. However, this kind of automated method could lead to some misidentification. Therefore, we further checked using their merger tree and star formation history to identify and correct false classifications. We found false starred subhalos that contain galaxies inside but show no star formation during their evolution. These subhalos contain interloper stars from larger neighbors, and thus, they are numerically misidentified as starred subhalos. We also identify false starless subhalos that experienced significant star formation during their history. This case is when subhalos do not have enough stars to be

identified as galaxies. We reclassify false starred subhalos as starless subhalos and false starless subhalos as starred subhalos.

To demonstrate the reliability of our classification, we show the mass evolution of the four populations, (true, false) × (starred, starless), in Figure B1. We show the median and 16th–84th percentiles of mass evolution for starred (blue) and starless (red) subhalos, distinguishing between true (solid lines and shaded areas) and false (dashed lines and hatched areas) populations. It is clearly visible that the true populations have trends in mass evolution distinct from the false populations. For instance, in the stellar mass evolution (panel (c)) of the false starless subhalos, there is significant star formation at the early epoch. However, their stellar mass barely reaches the minimum required to be detected as a galaxy, leading to their classification as starless subhalos. Also, the false starred subhalos have a lot of stars at the last snapshot, but we can infer that those are the interloper stars in the group systems. The overall trends of the false starless closely resemble those of the true starred, and vice versa, which is the reason why we modified the original classification.

One might be concerned about the possibility that subhalos capture field stars or stars from neighboring systems. J. Peñarrubia et al. (2024), using analytic and N-body methods, demonstrated that subhalos in a galactic environment do capture field stars and could keep them permanently under certain conditions. However, our subhalo samples have various epochs of infall to the main system, and before that, we find that capture events rarely happen. Truly, most cases of false starless and true starred subhalos have in situ star formation, which contributes most of their final stellar components.

<figure><img src="image_13.png" /><figcaption>**Figure B1.** Comparisons of the mass evolution of four populations with similar final DM masses. Starless subhalos are shown in red, and starred subhalos in blue. Solid lines and shading represent the median and 16th–84th percentiles of the true populations. Dashed lines and hatched areas are those of the false populations. Panels (a)–(d) show the mass evolution of DM, gas, stars, and cold gas, respectively.</figcaption></figure>
## Appendix C

### Supernova Energy and Rates

In the main text, we analyze the impact of supernova feedback with supernova power, $P_{SN}$. In this appendix, we want to define it. As prescribed in the subgrid model of supernova feedback (see the detailed process in Section 2.4 of Y. Dubois et al. 2021), a young stellar particle releases energy and mass during a supernova explosion, following the Chabrier initial mass function (G. Chabrier 2003). At each time step, we calculate the total supernova energy, $E_{SN}$, by summing the energy released from young stars within the boundary of each subhalo. The supernova energy rate, $\epsilon_{SN}$, at each time step ($t_i$) is calculated as

$$ \epsilon_{SN}(t_i) = \frac{E_{SN}(t_i)}{\Delta t_i}, \qquad (C1) $$

where $\Delta t_i$ is the time interval between the ith and (i - 1)th snapshots. To quantify both the amplitude and duration of supernova events, we define the “supernova power,” $P_{SN}$, by averaging these energy rates over time steps where $E_{SN}(t_i) > 0$. This physically indicates an expected value of supernova energy for a given time interval if supernovae occur. For instance, a high value of $P_{SN}$ indicates that a burst of supernovae occurred over a short period of time.

## Appendix D

### Radiative Cooling and Heating

As we have demonstrated the suppression of radiative cooling due to UV background heating in starless subhalos in the main text, we hereby provide a detailed description of the radiative cooling and heating prescription. Since radiative transfer is computationally expensive and not applied in NEWHORIZON and NEWHORIZON2, our code computes the density evolution of six species (e⁻, H I, H II, He I, He II, and He III) and then constructs the cooling and heating tables, taking into account the bremsstrahlung, ionization, recombination, dielectric recombination, line cooling, radiative heating, and Compton cooling and heating (model described by S. Courty & J. M. Alimi 2004). We reconstruct the cooling and heating rates based on the density and temperature of each

cell using the stored tables:

$$ \Lambda_{\text{net}} = \Lambda_{\text{cool}} + (Z/Z_{\odot})\Lambda_{\text{metal}} - \Gamma_{\text{heat}} + (\Lambda_{\text{Comp}} - \Gamma_{\text{Comp}})/n_{\text{H}}, \qquad (D1) $$

$$ \Lambda'_{\text{net}}(\frac{d\Lambda_{\text{net}}}{dT}) = \Lambda'_{\text{cool}} + (Z/Z_{\odot})\Lambda'_{\text{metal}} - \Gamma'_{\text{heat}} + (\Lambda'_{\text{Comp}} - \Gamma'_{\text{Comp}})/n_{\text{H}}, \qquad (D2) $$

where $\Lambda_{net}$ and $\Lambda'_{net}$ are the net cooling rate and its derivative. $\Lambda_{cool}$ and $\Lambda_{metal}$ are the cooling rates of the gas and the metal, $\Gamma_{heat}$ is the heating rate, $\Lambda_{Comp}$ and $\Gamma_{Comp}$ are the Compton cooling and heating rates, and $n_{H}$ is the number density of hydrogen.

The uniform UV background heating after $z_{reion} = 10$ is implicitly incorporated into the radiative cooling and heating tables (J. Rosdahl & J. Blaizot 2012), rather than through direct radiative transfer. The self-shielding effect is included as a boosting factor as follows:

$$ n_{\text{H,corr}} = n_{\text{H}}/e^{-n_{\text{H}}/(0.01 \text{ H cm}^{-3})}, \qquad (D3) $$

where $n_{H,corr}$ is the corrected number density of hydrogen that is practically used to compute the cooling and heating rates, and the self-shielding threshold of $0.01 \text{ H cm}^{-3}$ is motivated by the result of J. Rosdahl & J. Blaizot (2012). In Figure D1, we display the sequential temperature map for high-redshift snapshots. Although the uniform UV background turns on at $z = 10$, it is prominent after $z \sim 7$.

Figure D2 shows the reconstructed change in temperature based on the above procedures, excluding adiabatic processes. We compute the net change in temperature for different gas densities and temperatures. Bluer colors represent net cooling, and redder colors represent net heating. The three panels show three different epochs around reionization. The heating zones are most sensitive to UV background heating. A sharp transition from cooling to heating occurs around a density of $0.01 \text{ H cm}^{-3}$, corresponding to the self-shielding threshold. This effect becomes more pronounced around redshift 7 when the epoch of reionization is almost completed. The overall distribution of dense cells in starless (gray and magenta) and starred (black and cyan) subhalos is also shown with contours and markers. Around $z = 7$, the divergence is noticeable, and the starless subhalos still remain in the heating zones, while starred subhalos have already entered the cooling zones.
<figure><img src="image_14.png" /><figcaption>**Figure D1.** Sequential temperature map of the NEWHORIZON simulation, as it visually demonstrates the stage of reionization. From left to right, each slice indicates a different redshift from 11 to 5. Redder colors represent higher temperatures.</figcaption></figure>

<figure><img src="image_15.png" /><figcaption>**Figure D2.** The temperature changes within a single computational time step, following the radiative cooling and heating prescriptions. The metallicity is assumed to be 0.001 times the solar metallicity. The three panels are for the three different epochs marked at the top. We color-coded the net temperature change across the parameter space (cell density and temperature). Contours represent the 0.5σ, 1σ, 2σ, and 3σ distributions of dense cells of starless (gray) and starred (black) subhalos. The dense cells are defined as the top 25 densest within the virial radius. Markers represent peak values in gas phase space.</figcaption></figure>

## ORCID iDs

Seyoung Jeon https://orcid.org/0000-0002-1270-4465

Sukyoung K. Yi https://orcid.org/0000-0002-4556-2619

Emanuele Contini https://orcid.org/0000-0002-2873-8598

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

San Han https://orcid.org/0000-0001-9939-713X

Katarina Kraljic https://orcid.org/0000-0001-6180-0245

Sebastien Peirani https://orcid.org/0000-0001-6902-2898

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

Jinsu Rhee https://orcid.org/0000-0002-0184-9589

## References

* Applebaum, E., Brooks, A. M., Christensen, C. R., et al. 2021, ApJ, 906, 96

* Arjona-Gálvez, E., Di Cintio, A., & Grand, R. J. J. 2024, A&A, 690, A286

* Aubert, D., Pichon, C., & Colombi, S. 2004, MNRAS, 352, 376

* Baes, M., & Camps, P. 2015, A&C, 12, 33

* Bahé, Y. M., & McCarthy, I. G. 2015, MNRAS, 447, 969

* Baldry, I. K., Driver, S. P., Loveday, J., et al. 2012, MNRAS, 421, 621

* Balogh, M. L., McGee, S. L., Mok, A., et al. 2016, MNRAS, 456, 4364

* Benítez-Llambay, A., & Frenk, C. 2020, MNRAS, 498, 4887

* Benítez-Llambay, A., Frenk, C. S., Ludlow, A. D., & Navarro, J. F. 2019, MNRAS, 488, 2387

* Benítez-Llambay, A., Navarro, J. F., Abadi, M. G., et al. 2015, MNRAS, 450, 4207

* Benson, A. J., Frenk, C. S., Lacey, C. G., Baugh, C. M., & Cole, S. 2002, MNRAS, 333, 177

* Bhagwat, A., Costa, T., Ciardi, B., Pakmor, R., & Garaldi, E. 2024, MNRAS, 531, 3406

* Bignone, L. A., Helmi, A., & Tissera, P. B. 2019, ApJL, 883, L5

* Boylan-Kolchin, M., Bullock, J. S., & Kaplinghat, M. 2011, MNRAS, 415, L40

* Brooks, A. M., Kuhlen, M., Zolotov, A., & Hooper, D. 2013, ApJ, 765, 22

* Bullock, J. S., Kravtsov, A. V., & Weinberg, D. H. 2000, ApJ, 539, 517

* Callingham, T. M., Cautun, M., Deason, A. J., et al. 2019, MNRAS, 484, 5453

* Camps, P., & Baes, M. 2020, A&C, 31, 100381

* Chabrier, G. 2003, PASP, 115, 763

* Contini, E., De Lucia, G., Villalobos, Á., & Borgani, S. 2014, MNRAS, 437, 3787

* Contini, E., Gu, Q., Ge, X., et al. 2020, ApJ, 889, 156
Contini, E., Jeon, S., Rhee, J., Han, S., & Yi, S. K. 2023, ApJ, 958, 72

Contini, E., Rhee, J., Han, S., Jeon, S., & Yi, S. K. 2024, AJ, 167, 7

Courty, S., & Alimi, J. M. 2004, A&A, 416, 875

Darvish, B., Mobasher, B., Sobral, D., et al. 2016, ApJ, 825, 113

Dashyan, G., Silk, J., Mamon, G. A., Dubois, Y., & Hartwig, T. 2018, MNRAS, 473, 5698

Davé, R., Anglés-Alcázar, D., Narayanan, D., et al. 2019, MNRAS, 486, 2827

Dekel, A., & Silk, J. 1986, ApJ, 303, 39

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

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

Dubois, Y., & Teyssier, R. 2008, A&A, 477, 79

Efstathiou, G. 1992, MNRAS, 256, 43P

Elbaz, D., Daddi, E., Le Borgne, D., et al. 2007, A&A, 468, 33

Emerick, A., Bryan, G. L., & Mac Low, M.-M. 2019, MNRAS, 482, 1304

Engler, C., Pillepich, A., Pasquali, A., et al. 2021, MNRAS, 507, 4211

Federrath, C., & Klessen, R. S. 2012, ApJ, 761, 156

Fitts, A., Boylan-Kolchin, M., Elbert, O. D., et al. 2017, MNRAS, 471, 3547

Flores, R. A., & Primack, J. R. 1994, ApJL, 427, L1

Forbes, J. C., Krumholz, M. R., Goldbaum, N. J., & Dekel, A. 2016, Natur, 535, 523

García-Gallego, O., Iršič, V., Haehnelt, M. G., Viel, M., & Bolton, J. S. 2025, arXiv:2504.06367

Garrison-Kimmel, S., Hopkins, P. F., Wetzel, A., et al. 2019, MNRAS, 487, 1380

Geen, S., Slyz, A., & Devriendt, J. 2013, MNRAS, 429, 633

Geha, M., Mao, Y.-Y., Wechsler, R. H., et al. 2024, ApJ, 976, 118

Gnedin, N. Y. 2000, ApJ, 542, 535

Götz, M., & Sommer-Larsen, J. 2002, in The Evolution of Galaxies (Dordrecht: Springer), 415

Grcevich, J., & Putman, M. E. 2009, ApJ, 696, 385

Gunn, J. E., Gott, I., & J. Richard, I. 1972, ApJ, 176, 1

Haardt, F., & Madau, P. 1996, ApJ, 461, 20

Han, S., Dubois, Y., Lee, J., et al. 2025, ApJ, 978, 96

Hazenfratz, R., Barai, P., Lanfranchi, G. A., & Caproni, A. 2024, ApJ, 969, 65

Homma, D., Chiba, M., Komiyama, Y., et al. 2019, PASJ, 71, 94

Hu, W., Barkana, R., & Gruzinov, A. 2000, PhRvL, 85, 1158

Iršič, V., Viel, M., Haehnelt, M. G., et al. 2024, PhRvD, 109, 043511

Ivezić, Ž., Kahn, S. M., Tyson, J. A., et al. 2019, ApJ, 873, 111

Jang, J. K., Yi, S. K., Rey, S.-C., et al. 2024, ApJ, 969, 59

Jarrett, T. H., Cluver, M. E., Magoulas, C., et al. 2017, ApJ, 836, 182

Jung, M., Roca-Fàbrega, S., Kim, J.-H., et al. 2024, ApJ, 964, 123

Katz, H., Ramsoy, M., Rosdahl, J., et al. 2020, MNRAS, 494, 2200

Kawinwanichakij, L., Papovich, C., Quadri, R. F., et al. 2017, ApJ, 847, 134

Kim, S. Y., Peter, A. H. G., & Hargis, J. R. 2018, PhRvL, 121, 211302

Kimm, T., Cen, R., Devriendt, J., Dubois, Y., & Slyz, A. 2015, MNRAS, 451, 2900

Kimm, T., Haehnelt, M., Blaizot, J., et al. 2018, MNRAS, 475, 4617

Klypin, A., Kravtsov, A. V., Valenzuela, O., & Prada, F. 1999, ApJ, 522, 82

Kravtsov, A., & Manwadkar, V. 2022, MNRAS, 514, 2667

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

Larson, R. B. 1974, MNRAS, 169, 229

Lee, G., Hwang, H. S., Lee, J., Shin, J., & Song, H. 2024, ApJ, 962, 129

Lewis, J. S. W., Ocvirk, P., Aubert, D., et al. 2020, MNRAS, 496, 4342

Li, C., & White, S. D. M. 2009, MNRAS, 398, 2177

Licquia, T. C., & Newman, J. A. 2015, ApJ, 806, 96

Lovell, M. R., Eke, V., Frenk, C. S., et al. 2012, MNRAS, 420, 2318

Lynden-Bell, D. 1976, MNRAS, 174, 695

Mao, Y.-Y., Geha, M., Wechsler, R. H., et al. 2021, ApJ, 907, 85

Marsh, D. J. E. 2016, PhR, 643, 1

Martin, G., Bazkiaei, A. E., Spavone, M., et al. 2022, MNRAS, 513, 1459

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

McConnachie, A. W. 2012, ApJ, 144, 4

Moore, B. 1994, Natur, 370, 629

Moore, B., Ghigna, S., Governato, F., et al. 1999, ApJL, 524, L19

Munshi, F., Brooks, A. M., Applebaum, E., et al. 2017, arXiv:1705.06286

Musso, M., Cadiou, C., Pichon, C., et al. 2018, MNRAS, 476, 4877

Nadler, E. O. 2025, ApJL, 983, L23

Nadler, E. O., Drlica-Wagner, A., Bechtol, K., et al. 2021, PhRvL, 126, 091101

Nadler, E. O., Wechsler, R. H., Bechtol, K., et al. 2020, ApJ, 893, 48

Nashimoto, M., Tanaka, M., Chiba, M., et al. 2022, ApJ, 936, 38

Natale, G., Popescu, C. C., Rushton, M., et al. 2022, MNRAS, 509, 2339

Navarro, J. F., & Steinmetz, M. 1997, ApJ, 478, 13

Noh, Y., & McQuinn, M. 2014, MNRAS, 444, 503

Ocvirk, P., Aubert, D., Sorce, J. G., et al. 2020, MNRAS, 496, 4087

Okamoto, T., Gao, L., & Theuns, T. 2020, MNRAS, 390, 920

O'Leary, J. A., Steinwandel, U. P., Moster, B. P., Martin, N., & Naab, T. 2023, MNRAS, 520, 897

Oñorbe, J., Boylan-Kolchin, M., Bullock, J. S., et al. 2015, MNRAS, 454, 2092

Pakmor, R., Springel, V., Coles, J. P., et al. 2023, MNRAS, 524, 2539

Pawlowski, M. S., Kroupa, P., Angus, G., et al. 2012, MNRAS, 424, 80

Peñarrubia, J., Errani, R., Walker, M. G., Gieles, M., & Boekholt, T. C. N. 2024, MNRAS, 533, 3263

Peng, Y. J., Lilly, S. J., Kovac, K., et al. 2010, ApJ, 721, 193

Pereira-Wilson, M., Navarro, J. F., Benítez-Llambay, A., & Santos-Santos, I. 2023, MNRAS, 519, 1425

Pillepich, A., Sotillo-Ramos, D., Ramesh, R., et al. 2024, MNRAS, 535, 1721

Pillepich, A., Springel, V., Nelson, D., et al. 2018, MNRAS, 473, 4077

Planck Collaboration, Aghanim, N., Akrami, Y., et al. 2020, A&A, 641, A6

Polisensky, E., & Ricotti, M. 2011, PhRvD, 83, 043506

Rahmani, S., Lianou, S., & Barmby, P. 2016, MNRAS, 456, 4128

Read, J. I., & Erkal, D. 2019, MNRAS, 487, 5799

Rey, M. P., Pontzen, A., Agertz, O., et al. 2020, MNRAS, 497, 1508

Rhee, J., Smith, R., Choi, H., et al. 2020, ApJS, 247, 45

Rhee, J., Yi, S. K., Ko, J., et al. 2024, ApJ, 971, 111

Richstone, D. O. 1976, ApJ, 204, 642

Rosdahl, J., & Blaizot, J. 2012, MNRAS, 423, 344

Saintonge, A., & Catinella, B. 2022, ARA&A, 60, 319

Sawala, T., Frenk, C. S., Fattahi, A., et al. 2016, MNRAS, 457, 1931

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

Schneider, A., Anderhalden, D., Maccio, A. V., & Diemand, J. 2014, MNRAS, 441, L6

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

Simpson, C. M., Grand, R. J. J., Gómez, F. A., et al. 2018, MNRAS, 478, 548

Smercinina, A., Bell, E. F., Price, P. A., et al. 2018, ApJ, 863, 152

Smith, R., Choi, H., Lee, J., et al. 2016, ApJ, 833, 109

Somerville, R. S. 2002, ApJL, 572, L23

Somerville, R. S., & Davé, R. 2015, ARA&A, 53, 51

Spergel, D. N., & Steinhardt, P. J. 2000, PhRvL, 84, 3760

Spergel, D. N., Verde, L., Peiris, H. V., et al. 2003, ApJS, 148, 175

Tacchella, S., Diemer, B., Hernquist, L., et al. 2019, MNRAS, 487, 5416

Tanaka, M., Chiba, M., Hayashi, K., et al. 2018, ApJ, 865, 125

Teyssier, R. 2002, A&A, 385, 337

Thoul, A. A., & Weinberg, D. H. 1996, ApJ, 465, 608

Tweed, D., Devriendt, J., Blaizot, J., Colombi, S., & Slyz, A. 2009, A&A, 506, 647

van den Aarssen, L. G., Bringmann, T., & Pfrommer, C. 2012, PhRvL, 109, 231301

Van Nest, J., Munshi, F., Christensen, C., et al. 2023, ApJ, 956, 96

Vincenzi, M., Brout, D., Armstrong, P., et al. 2024, ApJ, 975, 86

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

Vogelsberger, M., Zavala, J., Cyr-Racine, F.-Y., et al. 2016, MNRAS, 460, 1399

Wadepuhl, M., & Springel, V. 2011, MNRAS, 410, 1975

Wetzel, A. R., Hopkins, P. F., Kim, J.-h., et al. 2016, ApJL, 827, L23

Wetzel, A. R., Tollerud, E. J., & Weisz, D. R. 2015, ApJL, 808, L27

Yi, S. K., Jang, J. K., Devriendt, J., et al. 2024, ApJS, 271, 1

Zier, O., Kannan, R., Smith, A., et al. 2025, arXiv:2503.02927

Zolotov, A., Brooks, A. M., Willman, B., et al. 2012, ApJ, 761, 71