# On the Origin of Intracluster Light based on the High-resolution Simulation, NewCluster

Seyoung Jeon, Emanuele Contini, San Han, Jinsu Rhee, Garreth Martin, Juhan Kim, Jaehyun Lee, Taysun Kimm, Christophe Pichon, Gyeong-Hwan Byun, et al.

► To cite this version:

Seyoung Jeon, Emanuele Contini, San Han, Jinsu Rhee, Garreth Martin, et al. On the Origin of Intracluster Light based on the High-resolution Simulation, NewCluster. The Astrophysical Journal, 2026, 998 (1), pp.30. 10.3847/1538-4357/ae2eaa. hal-05500412

HAL Id: hal-05500412

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

Submitted on 9 Feb 2026

**HAL** is a multi-disciplinary open access archive for the deposit and dissemination of scientific research documents, whether they are published 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 CC BY 4.0 - Attribution - International License
On the Origin of Intracluster Light based on the High-resolution Simulation, NewCluster

SEYOUNG JEON,¹ EMANUELE CONTINI,¹ SAN HAN,² JINSU RHEE,² GARRETH MARTIN,³ JUHAN KIM,⁴
JAEHYUN LEE,⁵ TAYSUN KIMM,¹ CHRISTOPHE PICHON,²,⁶ GYEONG-HWAN BYUN,¹ YOHAN DUBOIS,²
CORENTIN CADIOU,² J. K. JANG,¹ AND SUKYOUNG K. YI,¹

¹Department of Astronomy and Yonsei University Observatory, Yonsei University, 50 Yonsei-ro, Seodaemun-gu, Seoul 03722, Republic of Korea

²Institut d'Astrophysique de Paris, CNRS and Sorbonne Université, UMR 7095, 98 bis Boulevard Arago, F-75014 Paris, France

³School of Physics and Astronomy, University of Nottingham, University Park, Nottingham NG7 2RD, UK

⁴Center for Advanced Computation, Korea Institute for Advanced Study, 85 Hoegi-ro, Dongdaemun-gu, Seoul 02455, Republic of Korea

⁵Korea Astronomy and Space Science Institute, 776, Daedeokdae-ro, Yuseong-gu, Daejeon 34055, Republic of Korea

⁶Kyung Hee University, Dept. of Astronomy & Space Science, Yongin-shi, Gyeonggi-do 17104, Republic of Korea

ABSTRACT

Intracluster light (ICL) is a key component of galaxy clusters, with the potential to trace their dynamical assembly histories and the underlying dark matter distribution. Despite these prospects, its faint nature makes a consensus on its origin or population properties difficult to achieve, both in observations and simulations. In the hope of finding a breakthrough, we utilize the ongoing high-resolution cluster simulation, NEWCLUSTER. By classifying billions of particles in and around the cluster with a rigorous tracking procedure, we find that the majority of the ICL originates from satellites, including surviving and disrupted galaxies. Another notable finding is that the *preprocessed* component follows the density profile of dark matter better than the other components and has distinctive properties: old age, low metallicity, and enhanced α-element abundance. We further investigate the orbital dynamics, and our results demonstrate that the stripped fraction of satellites is primarily determined by the *time since infall* and the pericenter distance. By linking the demographic, chemical, and orbital properties of ICL stars to their origins, this work proposes a quantitative approach for tracing the assembly history of galaxy clusters from the ICL.

Keywords: galaxies: clusters: general

1. INTRODUCTION

Intracluster light (ICL) is a diffuse stellar component residing in the gravitational potential of a galaxy cluster, without being bound to the potential of individual satellite galaxies (Mihos 2016; Montes & Trujillo 2019; Contini 2021; Montes 2022). Since its first detection by Zwicky (1937, 1951) in the Coma cluster, it has been recognized as an important cluster component and visible evidence of the hierarchical structure formation via numerous mergers and tidal stripping (Willman et al. 2004; Mihos et al. 2005; Zibetti et al. 2005; Monaco et al. 2006; Murante et al. 2007; Contini et al. 2014; Montes 2022; Brough et al. 2024; Jiménez-Teja et al. 2025). Therefore, studying the ICL is anticipated to reveal the past assembly history of galaxy clusters. Furthermore, the ICL is expected to follow the cluster's total dynamical mass distribution (Montes & Trujillo 2019; Contreras-Santos et al. 2024; Yoo et al. 2025), which is dominated

by dark matter. As a consequence, its observation can
 unlock an independent way to trace dark matter, com-
plementary to X-ray emission (Borgani & Guzzo 2001)
or gravitational lensing analysis (Kochanek 2006; Clowe
et al. 2006; Cha & Jee 2023).

Despite these promising prospects, the depth of analy-
sis is relatively limited due to its faint and diffuse nature.
The typical surface brightness of the ICL is fainter than
~ 26 mag arcsec⁻² in the optical range, and can reach
~ 30 mag arcsec⁻² or beyond, which requires substan-
tial sensitivity and resolution (Mihos et al. 2005; Zibetti
et al. 2005; Rudick et al. 2011). In recent decades, mod-
ern deep imaging surveys have reached this limit (Tru-
jillo & Fliri 2016; Montes et al. 2021; Gonzalez et al.
2021; Montes & Trujillo 2022; Martínez-Lombilla et al.
2023; Joo & Jee 2023; Kluge et al. 2025), enabling sta-
tistical studies. However, the uncertainty in detection
is still huge, and many properties obtained by differ-
ent studies have not yet converged. For example, the
primary property of the ICL is the mass fraction rel-
ative to the total stellar mass or the brightest cluster
galaxy (BCG) mass in a galaxy cluster. Current studies
suggest the fraction measured is 5–50 percent (Contini
et al. 2024c, see references therein), but it is still un-
clear whether this large range of fractions comes from
cluster-to-cluster variation or methodological differences
(Rudick et al. 2011; Kimmig et al. 2025). Moreover, how
the fraction changes with redshifts or halo properties
(e.g., mass, concentration, or dynamical states) is more
complex (Golden-Marx et al. 2025), or even whether it
shows little variation, remains a debated topic (Ko &
Jee 2018; Joo & Jee 2023; Contini et al. 2024b).

Another important issue is the origin of the ICL. Tidal stripping of satellite galaxies, by-products of mergers with the BCG, or preprocessing are considered important channels (Mihos 2003; Chun et al. 2024; Contini et al. 2024b); however, which process has a dominant role or how these different origins affect the stellar populations are still under debate. Perhaps the most difficult issue is the absence of consensus on definitions. Specifically, as recently introduced in Brough et al. (2024), there are various ways of determining the ICL itself, such as using a surface brightness threshold (Montes & Trujillo 2014; Martínez-Lombilla et al. 2023), masking galactic light (Annunziatella et al. 2016; Ahad et al. 2023; Euclid Collaboration et al. 2025), or a wavelet technique (Adami et al. 2005; Ellien et al. 2021). This lack of agreement fundamentally stems from the non-trivial ICL detection.

Numerical simulations, in parallel, have provided a
powerful tool to study the formation of the ICL. Early
on, Murante et al. (2004) identified the diffuse light in
the simulated galaxy cluster and found that the dif-
fuse light is more centrally concentrated and older than
the galaxy light. Using the TNG300 simulation (Nelson
et al. 2019), Montenegro-Taborda et al. (2023) showed
that 90% of the ICL mass is contributed by “accreted”
stars. Recently, Kimmig et al. (2025) compared several
independent cosmological simulations, including MAG-
NETICUM PATHFINDER (Dolag et al. 2025), TNG100
(Nelson et al. 2019), HORIZON-AGN (Dubois et al.
2014b), and HYDRANGEA (Bahé et al. 2017); they found
that the formation redshift of galaxy clusters is tightly
correlated with the mass fraction of the sum of BCG and
ICL (hereafter, BCG+ICL) relative to the cluster’s total
stellar mass. Complementing these hydrodynamical ef-
forts, semi-analytic models (SAMs) have also provided
key insights. For example, since Monaco et al. (2006)
predicted a total amount of diffuse stellar components
with the MORGANA model, Contini et al. (2024d) in-

troduced the FEGA model that provides a large statisti-
cal sample of clusters, a task that remains computation-
ally too expensive for fully-hydrodynamical simulations.
Despite the merits of using simulations, the faint nature
of the ICL still remains a major challenge for theoretical
studies as well. A galaxy cluster is a very massive ob-
ject, and it requires much more computational resources
than field environments to maintain particle mass and
spatial resolutions. Given the limited resolution of con-
ventional cosmological simulations, their resolution has
been insufficient to fully resolve the detailed structure
and stellar populations of faint tidal features and merger
relics (Martin et al. 2024; Brown et al. 2024; Kimmig
et al. 2025). Although the semi-analytic approach gives
us significant benefits, particularly in statistical trends
(Guo et al. 2011; Contini et al. 2023, 2024a), it is also
inherently limited by its nature of solving analytic equa-
tions for galaxy properties rather than employing indi-
vidual stellar particles. For a detailed analysis of the
ICL, it is thus important to resolve the dynamics of stars
finely to track the dynamical stripping of stars.

In response to this demand, we use NEWCLUSTER, a high-resolution cosmological hydrodynamical simulation targeting a galaxy cluster (Han et al. 2025b). As we will describe the details of NEWCLUSTER in Section 2.1, we hereby highlight the strengths of NEWCLUSTER in studying the ICL. NEWCLUSTER achieves substantially improved resolutions compared to other cluster-scale simulations (cf. Tremmel et al. 2019, Table 1), specifically a stellar mass resolution of 2×10⁴ M☉ and a spatial resolution of 68 pc, allowing studies down to a limiting r-band surface brightness of 33.5 mag arcsec⁻² at z = 0.8. We also have an order of magnitude higher output cidence of ~15 Myr (i.e., ~1000 snapshots for redshifts 50 to 0) than the conventional setting of ~100 snapshots. This enables better tracking of galaxies across snapshots and capturing the orbital motion in dense environments.

Throughout this work, utilizing NEWCLUSTER, we
study the origin of the ICL from the detailed orbital
history of individual stellar particles, which has never
been attempted at this resolution. In Section 2, we de-
scribe information on NEWCLUSTER, the halo/galaxy
identification, the merger tree construction, and the def-
inition of satellite galaxies. In Section 3, we describe
comprehensively the procedure that defines BCG+ICL
stars and further divides them into subsamples by their
origins. The overall results and their demographics are
presented in Section 4. Section 5 presents the analysis
of the properties of their stellar populations. In Sec-
tion 6, we further investigate the role of the orbital mo-
tion in the formation of the ICL. We then discuss the
current limitations and the future outlook in Section 7, and finally, we summarize our findings and conclude in Section 8.

## 2. METHODS

In this section, we describe the details of the simulation. Building merger trees and classifying satellites in the galaxy cluster are critical factors of this work, so we give particular attention to them.

### 2.1. *Numerical simulations*

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

**Figure 1.** Overview of the NEWCLUSTER simulation centered on the primary cluster halo at $z = 0.79$. Panel (a) presents a dark matter density map (grayscale), with the virial radii of the two cluster halos (magenta and cyan circles). The dotted line shows the virial radius ($R_{vir}$) from ADAPTAHOP, and the dashed line is the $R_{200}$, which are described in Section 2.2. Panel (b) shows the corresponding $r$-band surface brightness map derived from the stellar density. Note that faint but complex tidal features are resolved in the surface brightness map.

NEWCLUSTER is a cosmological zoom-in simulation targeting a galaxy cluster with a virial mass comparable to that of the Virgo cluster at $z = 0$ (Han et al. 2025b).

The simulation is conducted using the RAMSES-yOMP code (Han et al. 2025a), the hybrid parallelized version of the adaptive mesh refinement (AMR) code, RAMSES (Teyssier 2002), assuming the Chabrier (2003) initial mass function and the WMAP7 cosmology (Kobayashi et al. 2011). The initial volume of the targeted zoom-in region is ~23,300 Mpc$^3$ at the center of the total (100 $h^{-1}$ Mpc)$^3$ box. In this zoom-in region, each cell can form stars following the gravo-thermo-turbulent model (Kimm et al. 2017) when its density exceeds 5 Hcm$^{-3}$. The subsequent stellar feedback process is composed of Type II and Ia supernovae with prescriptions from Kimm & Cen (2014). NEWCLUSTER also tracks the evolution of chemical elements (H, D, He, C, N, O, Mg, Fe, Si, and S) ejected from stellar wind (Schaller et al. 1992; Maeder & Meynet 2000), Type Ia supernovae (Kobayashi et al. 2006), and Type II supernovae (Iwamoto et al. 1999). The dual-mode feedback (i.e., radio and quasar modes by the accretion rates) from active galactic nuclei (AGN) is also embedded in NEWCLUSTER (Dubois et al. 2014a,b). For other detailed prescriptions, we refer readers to the introductory paper of NEWCLUSTER (Han et al. 2025b).

In Figure 1, we present an overview of NEWCLUSTER at $z = 0.79$. Panel (a) shows the dark matter density map and the primary and secondary halos. As the two halos are about to merge, the primary halo is expected to be highly perturbed in the near future. To avoid such unrelaxed dynamical states, we limit our study to the snapshot at $z = 0.79$. This snapshot is also the latest one that has been fully post-processed for analysis.<fnref n="7" />

NEWCLUSTER has many updated astrophysical prescriptions compared to its predecessors—HORIZON-AGN (Dubois et al. 2014a) and NEWHORIZON (Dubois et al. 2021)—, such as chemical elements, on-the-fly dust evolution, Lagrangian gas tracer particles (Genel et al. 2013; Cadiou et al. 2019), or dynamical friction of a supermassive black hole (Dubois et al. 2014a; Pfister et al. 2019; Han et al. 2025b). However, the most important strength related to this work is the resolution. The spatial resolution is 68 pc in the most refined cell, and the stellar mass resolution is ~$2 \times 10^4$ M$_\odot$. This allows us to resolve the primary cluster region into ~200 million stellar particles, which is an unprecedented scale. Low surface brightness features, including the ICL, require fine resolutions to better reproduce their diffuse nature and complicated dynamics. Using a suite of TNG simulations, Lovell et al. (2025) assessed the numerical effects and demonstrated that the resolution highly affects the

<footnote n="7">⁷ NEWCLUSTER is still ongoing and has currently passed $z \sim 0.56$ as of Dec. 2025.</footnote>
stripped mass of satellites. Although they argued that the stripping times of the dark matter are converged in their flagship runs—TNG50-1 ($m_{DM} = 4 \times 10^5$ M$_{\odot}$ and $\epsilon_{DM} = 0.29$ kpc) to TNG300-1 ($m_{DM} = 5.9 \times 10^7$ M$_{\odot}$ and $\epsilon_{DM} = 1.48$ kpc), where $m_{DM}$ and $\epsilon_{DM}$ are the dark matter particle mass resolution and softening length, respectively—, the stripped mass and stripping times of stellar particles are not converged in their resolutions. This magnifies the strength of our spatial resolution of 68 pc in minimizing the numerical artifacts on resolving tidal features. For example, in Figure 1(b), we present the r-band surface brightness map of NEWCLUSTER. The faint but highly disturbed features are prominent in the surface brightness map.

## 2.2. Halo (galaxy) identification

We utilize the ADAPTAHOP halo-finding algorithm to identify galaxies (GALAXYMAKER) and dark matter halos (HALOMAKER) from stellar and dark matter particles, respectively. As a detailed description of the algorithm can be found in Aubert et al. (2004); Tweed et al. (2009), we briefly outline the main procedure here. This method finds local density maxima and saddle points that lie above a minimum density threshold—conventionally $80 \times \rho_{crit}$ for halos and $178 \times \rho_{crit}$ for galaxies, where $\rho_{crit}$ is the critical density of the Universe—, which are regarded as initial guesses of the center and boundary of structures. After identifying the topological structures, ADAPTAHOP iteratively eliminates structures that do not physically satisfy the input conditions, such as the density threshold or the gradient of the substructure density. It also builds the structure-substructure hierarchy based on spatial occupation. For each ground-level structure, the structures hierarchically connected to the most massive substructure are integrated into a single structure, which is designated as the main structure, while the remaining structures are tagged as substructures (i.e., subhalos or satellite galaxies). To robustly identify the center, ADAPTAHOP iteratively shrinks concentric shells toward the densest region rather than simply selecting the densest particle or the center of mass. In this work, we use the term $R_{max}$, which is the maximum distance from the galaxy/halo center to any of its member particles, indicating the spatial occupation of the structure.

With these structures, ADAPTAHOP calculates the physical properties, such as the virial radius, the virial mass, or the spin parameters. For the virial radius of a dark matter halo, ADAPTAHOP directly calculates the virial condition using the member dark matter particles:

$$ \eta_{vir} = \left| \frac{2T + W}{T + W} \right|, \qquad (1) $$

where $\eta_{vir}$ is the virial ratio between the virial residual and the total energy, T is the kinetic energy, and W is the potential energy. By iteratively shrinking the concentric ellipsoidal radius from where $T + W < 0$ and computing $\eta_{vir}$ at each step, the virial radius ($R_{vir}$) is defined at the radius where $\eta_{vir}$ first drops below $0.2$.<fnref n="8" /> If a halo is too unrelaxed for this criterion to be met—i.e., $\eta_{vir}$ never falls below 0.2 over the inner sweep—, ADAPTAHOP regards the radius within which the density is greater than $\Delta_{vir}(z) \times \rho_{mean}$ as the virial radius, where $\Delta_{vir}(z)$ is the virial overdensity relative to the mean matter density computed from a spherical top hat collapse calculation, and $\rho_{mean}$ is the mean density of the Universe. The virial mass ($M_{vir}$) is defined as the dark matter mass within an ellipsoid corresponding to $R_{vir}$. However, because this calculation considers only the member dark matter particles, it can be offset from the total dynamical mass. Therefore, we further calculate $R_{200}$ as the radius within which the mean enclosed density is 200 times the critical density, taking baryon components into account, and the corresponding mass is $M_{200}$. In this work, we adopt $R_{200}$ and $M_{200}$ as the halo virial radius and mass. For reference, the $M_{vir}$ and $M_{200}$ of the primary halo are $1.09 \times 10^{14}$ M$_{\odot}$ and $7.76 \times 10^{13}$ M$_{\odot}$, and those of the secondary halo are $1.52 \times 10^{13}$ M$_{\odot}$ and $4.38 \times 10^{13}$ M$_{\odot}$ at $z = 0.79$.

## 2.3. Building merger trees

Since halo finders provide halo/galaxy catalogs only for individual snapshots, constructing a merger tree is essential for studying the time-series evolution of galaxies. Importantly, the merger tree must properly track low-mass galaxies in very dense environments, which is particularly crucial in galaxy cluster studies. The simplest and conventional method is to calculate the shared member particles between two adjacent snapshots and connect progenitors and descendants with the highest number of shared particles. However, this method can lead to a “broken tree” issue (see e.g., Srisawat et al. 2013; more recently Poole et al. 2017; Park et al. 2022; Chandro-Gómez et al. 2025). Commonly, when small galaxies closely pass through a larger galaxy, most particles could be falsely allocated to the larger one, or worse, a halo-finder temporarily fails to identify small structures embedded in a high-density region. Studies with high-resolution simulations have more technical problems due to, ironically, their good output resolutions. The standard parameter set of halo-finders often

<footnote n="8">⁸ This threshold is a default setting in AdaptaHOP, which is a small and non-zero threshold to account for minor deviations from a perfectly relaxed state.</footnote>
<figure><img src="image_3.png" /><figcaption>**Figure 2.** Schematic diagram of the merger tree. The spiral icons indicate all leaf galaxies connected to the target galaxy (z) across different snapshots, marked with alphabetical IDs. Each line represents a connected branch based on the score (described in the text). The main branch (black) has the same IDs for FIRST and HEAD, indicating a non-fragmented branch. Likewise, the same IDs for TAIL and FINAL mean that this branch is not merged into any other branch.</figcaption></figure>

overpredicts highly fragmented substructures in high-resolution simulations, requiring manual removal or adjustment.

To overcome this issue, we hereby introduce our in-house merger tree algorithm.<fnref n="9" /> To robustly measure the connectivity among galaxies across different snapshots, we measure the combined scores of 1) the bidirectional match score using shared particles, 2) the velocity offset of shared particles, and 3) the stellar mass differences. We first assign a proximity-based weight in phase space to each member particle. Consider two galaxies, Aᵢ and Bⱼ, at snapshots i and j, with Nₐᵢ and Nₐⱼ member particles, respectively. The weight (w̃ₖ|Aᵢ) of the k-th member particle of the Aᵢ is defined as:

$$\tilde{w}_{k|A_i} \equiv \frac{m_k}{\sqrt{\left(\frac{\|r_k - r_{A_i}\|}{\sigma_{r,A_i}}\right)^2 + \left(\frac{\|v_k - v_{A_i}\|}{\sigma_{v,A_i}}\right)^2}}, \quad (2)$$

where mₖ is the k-th particle mass, rₖ and vₖ are the position and velocity vectors of the k-th particle, and rₐᵢ and vₐᵢ are the position and velocity of the host Aᵢ. σₖ,Aᵢ and σv,Aᵢ are the standard deviations of distances and relative velocities of all members in Aᵢ, respectively. We then divide the weights by the sum of them, wₖ|Aᵢ = w̃ₖ|Aᵢ / Σₚ=1 w̃ₚ|Aᵢ, which ensures that Σₖ=1 wₖ|Aᵢ = 1.

Using the weights, we define the match rate as

$$S_{A_i}(B_j) \equiv \sum_{k \in A_i \cap B_j} w_{k|A_i}, \quad (3)$$

where Aᵢ ∩ Bⱼ denotes the shared particles between Aᵢ and Bⱼ. We define the bidirectional match score,

$$q_{A_i, B_j}^{(1)} = S_{A_i}(B_j) + S_{B_j}(A_i). \quad (4)$$

If Aᵢ and Bⱼ have identical members, q(Aᵢ∩Bⱼ) = 2.

We also incorporate kinematic information. For a single galaxy across different snapshots, the velocity of its shared particles relative to the bulk motion should not vary significantly. We therefore define a relative velocity of the shared particles to the host as

$$\Delta \mathbf{v}_{A_i}(B_j) \equiv \mathbf{v}_{A_i \cap B_j} - \mathbf{v}_{A_i}, \quad (5)$$

where vₐᵢ∩Bⱼ is the mean velocity of the shared particles, and vₐᵢ is the bulk velocity of Aᵢ. vₐᵢ∩Bⱼ and vₐᵢ are weighted averages using the normalized weights, wₖ|Aᵢ. We then define the second score of velocity offset,

$$q_{A_i, B_j}^{(2)} = 1 - \frac{\|\Delta \mathbf{v}_{B_j}(A_i) - \Delta \mathbf{v}_{A_i}(B_j)\|}{\|\Delta \mathbf{v}_{B_j}(A_i)\| + \|\Delta \mathbf{v}_{A_i}(B_j)\|}. \quad (6)$$

The score q(Aᵢ, Bⱼ) is unity for identical kinematics and zero for opposite.

As a final criterion, we also consider the stellar mass difference,

$$q_{A_i, B_j}^{(3)} = \left( \frac{\min(m_{A_i}, m_{B_j})}{\max(m_{A_i}, m_{B_j})} \right)^C, \quad (7)$$

<footnote n="9">⁹ The source code is available on GitHub (YoungTree: https://github.com/svj3514/YoungTree.git) under a 3-Clause BSD License, and version 1.0.0 is archived in Zenodo (Jeon 2025).</footnote>
where $m_{A_i}$ and $m_{B_j}$ are the stellar masses of $A_i$ and $B_j$, respectively, and $C$ is a tunable exponent. If $A_i$ and $B_j$ have similar masses, $q_{A_i,B_j}^{(3)}$ is close to unity. Finally, we calculate the combined score as

$$q_{A_i,B_j} = q_{A_i,B_j}^{(1)} + q_{A_i,B_j}^{(2)} + q_{A_i,B_j}^{(3)}, \quad (8)$$

for all galaxy pairs if they have at least one shared member particle.

Furthermore, since small galaxies are often missed in detection, we evaluate the score mentioned above not only between adjacent snapshots, but also across up to four snapshots (i.e., roughly ±60 Myr) before and after. After scoring, all galaxies have their progenitor and descendant candidates in previous and subsequent snapshots. If two galaxies mutually point at each other, we directly connect them into a single branch. Otherwise, we connect two branches, but instead record them as merger or fragmented events. As a result, the merger tree contains branches classified as main, merged, or fragmented. We formalize this classification by tagging the FIRST, HEAD, TAIL, and FINAL parameters:

* FIRST: first leaf ID of the tree (same FIRST for the same origin),

* HEAD: first leaf ID of the independent branch (same HEAD for the same branch),

* TAIL: last leaf ID of the independent branch (same TAIL for the same branch),

* FINAL: final leaf ID of the tree (same FINAL for the same galaxy at the final snapshot).

From these parameters, we define a “branch” as a sequence of leaves sharing the same HEAD and TAIL, while a “tree” is defined as the leaves connected with the same FIRST or FINAL.

We present the schematic diagram of an example tree of the target galaxy (z) in Figure 2. The main branch (black) represents a simple case with no fragmentation (FIRST = HEAD = (b)) or merger (TAIL = FINAL = (z)). The purple line is a merged branch. Initially, it has a different FIRST ((a)) from the main branch (black), but it merged later ((t) to (s)), resulting in a different TAIL ((t)) from FINAL ((z)). In addition, the branch is not broken even when (k) is not detected by ADAPTAHOP, because, as we mentioned, our tree algorithm checks the progenitor/descendant candidates in multiple snapshots. The red and blue lines are fragmented branches from the main branch. In most cases, they are star-forming clumps that are misidentified as individual galaxies by ADAPTAHOP. Our tree algorithm can capture these structures by checking FIRST and HEAD.

Our method builds on the standard technique of connecting galaxies via shared particles (Srisawat et al. 2013) by incorporating kinematic and mass information, conceptually similar to other modern algorithms that use velocity information to improve accuracy (Behroozi et al. 2013; Elahi et al. 2019). In a galaxy cluster, there are numerous galaxy-galaxy interactions or mergers with the BCG. The key strength of our tree is its stability in connecting small galaxies in crowded environments. We now robustly identify satellite galaxies that are still surviving or already disrupted based on our tree, as will be shown in Section 2.4.

## 2.4. Satellite galaxies

We identify the BCG from the galaxy catalogs as the most massive galaxy located at the center of the primary cluster halo. For satellite galaxies, we first identify all satellite galaxies within $1.5 \times R_{200}$, adopting a generous range. We exclude satellites that have the same origin (FIRST) as the BCG merger tree, which are likely fragmented structures of the BCG. We further check the merger tree of the remaining satellites and exclude those originating from other satellites, because they are transient fragments of real satellites. This selection process yields a final sample of 403 galaxies, which we define as “surviving” satellites.<fnref n="10" /> In Figure 3, we show their orbital tracks in blue lines, but only for galaxies more massive than $10^9 M_{\odot}$ in panel (b).

Since BCGs experience numerous merger events, the contribution of merged galaxies to the ICL is likely non-negligible (Murante et al. 2007; Brown et al. 2024; Mayes et al. 2025; Joo et al. 2025), although the relative contributions of merger and stripping depend on their definition (cf. Contini et al. 2018). To investigate this contribution, we also identify galaxies that have already merged and been disrupted. We select all galaxies that have the same destination (FINAL) as the BCG tree but different origins (FIRST). We then apply two cleaning steps. First, we exclude structures fragmented from the BCG if the maximum distance from the BCG is less than $1.5 \times R_{200}$. Second, we exclude fragments of surviving satellites. This complex case occurs when any of the progenitors of a disrupted candidate is temporarily linked to the branches of surviving satellites. If the origin (FIRST) of this candidate is the same as that of any surviving satellites, we regard it as a fragment and exclude it. After this, 67 galaxies are identified as “dis-

<footnote n="10">¹⁰ Note that we do not apply a simple mass cut for sampling but select galaxies that have a persistent merger tree across time (e.g., 441 progenitors for each satellite on average), which ensures more robust identification by excluding transient structures.</footnote>
<figure><img src="image_4.png" /></figure>

**Figure 3.** Orbital tracks of galaxies derived from the merger tree in the NEWCLUSTER cluster. Panel (a) shows the comoving-scale tracks in the cosmological box, while panel (b) shows the physical-scale tracks centered on the main cluster at $z = 0.79$ overlaid on a stellar density map (grayscale). The black dashed circles mark the $R_{200}$ of the main cluster halo. The red line in panel (a) traces the BCG, and the red marker in panel (b) indicates its final position. The green lines show satellites that have completely merged into the BCG. The blue circles indicate $R_{90}$, the radius enclosing 90% of its stellar mass, and the blue lines show orbits of surviving satellites. Unlike panel (a), we only present the tracks of surviving satellites more massive than $10^9 M_{\odot}$ in panel (b), for clarity.

rupted (merged)” satellites, and their tracks are shown in Figure 3 (green).

### 3. CLASSIFICATION OF BCG+ICL STARS

The ICL is conventionally defined as the diffuse stellar component bound to the cluster potential but unbound

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

**Figure 4.** Illustration of the definition of the ICL within $1.5 \times R_{200}$. Panel (a) shows all stars bound to satellite galaxies at $z = 0.79$, while panel (b) presents the density map of the final BCG+ICL component, after excluding the satellite stars shown in panel (a). As described in the main text, the tight (red) and loose (blue) member stars are identified and follow the motion of their host satellite galaxy.

from individual satellite galaxies. However, it is not

<table><caption>**Table 1.** Classification of BCG+ICL subsamples.</caption><thead><tr><td>Component</td><td>Definition</td></tr></thead><tbody><tr><td>stripped</td><td>stripped from surviving satellites</td></tr><tr><td>disrupted</td><td>stripped and disintegrated from disrupted satellites</td></tr><tr><td>in-situ</td><td>born in the BCG</td></tr><tr><td>preprocessed</td><td>from satellites but unbound before infall by group environments or neighboring galaxies</td></tr></tbody></table>
trivial to determine the true stellar membership of an individual galaxy. For example, the stars in the galactic outskirts are barely bound to the host galaxy and are often not identified as members and instead classified as ICL. However, when they are still kinematically aligned to the host and follow the motion of the host, this method potentially leads to an overestimation of the ICL. Alternatively, direct binding-energy calculations combined with kinetic information are more robust (Contreras-Santos et al. 2024; Joo et al. 2025), but uncertainties in center finding and, more practically, the computational cost are prohibitive for high-resolution simulations (Rudick et al. 2009; Puchwein et al. 2010).

Furthermore, decomposing the BCG and ICL components is one of the main challenges (Jiménez-Teja & Dupke 2016). It is also physically ambiguous whether the BCG and ICL are truly separated populations or the ICL is just an extended feature of BCGs (Donzelli et al. 2011; Mihos 2016; Pillepich et al. 2018; Montes & Trujillo 2019), apart from the technical difficulty in the separation (Contini et al. 2022). In observations, two- (or three-) component fitting and fixed-aperture masking are widely used to distinguish those (Gonzalez et al. 2005; Annunziatella et al. 2016; Kravtsov et al. 2018; Spavone et al. 2022; Ahad et al. 2023; Euclid Collaboration et al. 2025), acknowledging the overlap of populations. On the contrary, most theoretical studies adopt their own criteria; conventionally, particle membership from halo finder results (Dolag et al. 2010; Pillepich et al. 2018; Alonso Asensio et al. 2020; Yoo et al. 2022; Montenegro-Taborda et al. 2023; Brown et al. 2024; Kimmig et al. 2025). This method has the drawback of missing stars in the outskirts (i.e., stellar halos), which are supposed to contribute significantly to the ICL because of their weak boundness (Mayes et al. 2025). Also, similar to observational limits, a halo finder relying on spatial distribution cannot distinguish true galaxy members from ICL stars temporarily overlapping with galaxies. To make matters worse, since most halo finders rely on single-snapshot data, stellar membership—and even galaxy detection itself—can vary strongly from snapshot to snapshot (Behroozi et al. 2013; Srisawat et al. 2013).

Considering these, to isolate the BCG+ICL component from all stellar populations in the cluster, we employ the time-series membership approach, utilizing the merger tree. The membership from ADAPTAHOP, at the very least, indicates that the stars are temporarily in the same spatial occupation of the galaxy. We evaluate their kinematic alignment to the host (i.e., continue to follow the host's motion) by checking whether they consistently appear in the membership. In particular, we classify the member stars of the surviving

satellites within a 500 Myr time window before the final snapshot. During this time, stars identified as members of the satellite galaxy in every post-birth snapshot are defined as “tight,” while those present in more than 10% of their lifetime as “loose.” Although this threshold is somewhat arbitrary, it approximates binding without explicit potential calculations and is more robust than single-snapshot assignments.

Figure 4(a) presents the resulting tight (red) and loose (blue) stellar members of surviving satellites. In this work, to ensure a conservative definition of the BCG+ICL components, we include the loose members in the satellite membership, which can be interpreted as the upper limit of the bound stars.<fnref n="11" /> By subtracting these satellite members from the total stellar population, we define our final sample, BCG+ICL, shown in Figure 4(b). The mass fraction of BCG+ICL ($f_{BCG+ICL}$) is ~ 41%, which is broadly consistent with other studies (Pillepich et al. 2018; Brough et al. 2024; Kluge et al. 2025; Ellien et al. 2025; Mayes et al. 2025). As previous studies generally predict an increasing ICL fraction toward $z=0$ but no significant change for $z > 0.5$ (Burke et al. 2015; Joo & Jee 2023; Contini et al. 2024b), this provides another benchmark for studies of the ICL fraction’s evolution at intermediate and high redshifts.

The final subsamples of BCG+ICL are classified into four components based on their origin: stripped, disrupted, in-situ, and preprocessed. The definitions of each component are summarized in Table 1, and the detailed procedures are described in the following subsections.

### 3.1. BCG+ICL from surviving satellites

For each surviving satellite, we trace its full progenitor history using the merger tree (sharing the same FIRST or FINAL) and then extract all member particles of any progenitors within a window from 500 Myr before infall to the first infall (entering $R_{200}$ of the cluster, $t_{inf}$). As in Section 3, we classify stars as tight (> 90%), loose (10–90%), or non-members (< 10%) by their occurrence across this window; for stars born within it, by the fraction of their lifetime. We define the “stripped” components as stars that are initially tight or loose members of surviving satellites but are later stripped into the BCG+ICL. Inspecting orbits, we find that a small fraction of loose members become unbound even before cluster infall, which arises from uncertainties in the binding criterion. We therefore reclassify them as preprocessed

<footnote n="11">¹¹ The mass fraction of the loose members in the total satellite membership is ~ 24%, a proxy for the uncertainty in our methodology.</footnote>
<figure><img src="image_6.png" /></figure>

**Figure 5.** Top (a)-(g): example of the infall member stars carried by a target surviving satellite galaxy. Panel (a) shows the overall projected smoothed trajectory (green) of the target satellite. Each position is normalized by the $R_{200}$ of the main cluster progenitors. The purple line highlights the time interval (500 Myr before infall) used to identify member stars. Panels (b)-(g) present zoom-in views of the target satellite at each stage, showing the tight (red) and loose (blue) member stars along with other stars (grayscale). We also indicate the direction of the BCG with arrows and markers. Bottom (h)-(n): example of infall stars associated with a target disrupted galaxy. Similar to the top panels, panel (h) shows the overall projected trajectory, and the other panels present zoom-in views. Since the target does not survive to the final snapshot, we include $t_{disrupted}$ (red), the epoch when the satellite is last detected.

instead of stripped components. The preprocessed components are discussed in Section 3.4.

In the top panels of Figure 5, we present the track and stellar distribution of a sample satellite. Panel (a) shows the satellite's main-branch track (green line),
normalized by $R_{200}$ of the cluster progenitor at each time, and other panels present zoom-in views exhibiting the stripped stellar distribution (red and blue) over a grayscale background of other stars. The membership remains well sustained before infall—panels (b)–(c), whereas after infall the stripped components are removed by the cluster tidal field, with the effect especially severe for the loose members—panels (d)–(g). For completeness, we also include stars formed within 500 Myr after $t_{inf}$ in the satellite, although they are not shown in Figure 5 for clarity.

## 3.2. BCG+ICL from disrupted satellites

We identify disrupted (or merged) galaxies as those sharing the same final destination (FINAL) as the BCG but with a different origin (FIRST), and determine their stellar membership using the same criteria as for the surviving satellites. The subset of these members later found in the BCG+ICL is defined as the “disrupted” component, after reclassifying any pre-infall unbound stars as preprocessed.

The bottom panels of Figure 5 display the track and distribution of the disrupted components for a sample galaxy, in the same format as the top panels. Since the target is disrupted, we highlight in red the snapshot when this satellite was last detected ($t_{disrupted}$). As with surviving satellites, the disrupted components remain well bound prior to infall. However, after infall, they are strongly stripped, even more severely than for surviving satellites. At the final snapshot (Figure 5(n)), both the tight (red) and loose (blue) members are highly scattered along their orbital path, although the loose members show a more extended distribution.

## 3.3. BCG+ICL born in BCG

Since the BCG has resided at the center of the overdensity, it likely formed a substantial number of in-situ stars through continuous gas accretion, in addition to those contributed by direct mergers (Bonaventura et al. 2017; Montenegro-Taborda et al. 2023; Mayes et al. 2025). To identify these in-situ stars, we extract the full merger tree of the BCG (same FIRST), including its fragmented structures—mostly internal star-forming clumps—. We classify stars associated with the extracted BCG tree based on the fraction of their lifetime spent as members, defining them as tight if the fraction is > 90% and loose if it is ≥ 10%. These stars are collectively referred to as the “in-situ” components.

## 3.4. BCG+ICL from preprocessing

Meanwhile, the environmental effects below the cluster scale or before the cluster infall (i.e., preprocessing)

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

**Figure 6.** Example of stellar trajectories in a target satellite galaxy. The x-axis shows the time relative to infall ($t_{inf}$; black dashed vertical line), and the y-axis is the moving distance of stars from the satellite center. The black lines denote “stable” stars that remain tightly bound to the satellite throughout the simulation. The blue lines indicate “stripped” stars that are initially bound, but stripped after infall. The red lines represent “preprocessed” stars that are already unstable or unbound well before infall. The satellite radius is also shown at each time as a guideline for boundness (green).

are now actively studied, mainly on the theoretical side. For instance, Contini et al. (2024a) revealed the non-negligible contribution of preprocessing to the diffuse light fractions of low-mass ($M_{halo} < 10^{13} M_{\odot}$) halos. Conventionally in the ICL field, “preprocess” refers to the ICL stars that initially reside in their host satellites but become unbound before cluster infall (Fujita 2004; Han et al. 2018; Jung et al. 2018; Contini 2021; Contini et al. 2024b; Jiménez-Teja et al. 2025). However, in this work, we regard all processes that unbind previously bound stars before cluster infall (predominantly tidal stripping) as preprocessing, because tidal heating/stripping are largely self-similar across mass scales in the hierarchical structure formation. To finalize the classification of the preprocessed candidates mentioned in Sections 3.1 and 3.2, we evaluate their boundness by examining their orbital motions. The kinematic signature of a preprocessed component is generally either continued recession from its host or a sudden increase in apocenter distance well before infall. In practice, stars showing clear evidence of prior unbinding are classified as “preprocessed” components. A visual justification is provided in Figure 6.
To justify and illustrate the classification of the preprocessed components, we present the distances of member stars of a surviving satellite in Figure 6. We display three populations: stable stars that consistently reside in the host to the end (black; satellite stars), stars tidally stripped after infall (blue; stripped), and the preprocessed stars that become unbound before infall (red). The preprocessed stars are evidently unbound and exhibit more extended orbits. The main cause of preprocessing for this satellite is a major merger at 2 Gyr before $t_{inf}$, which produces the corresponding bump in its radius (green).

In summary, having classified the BCG+ICL stars into four subsamples by their origin, we now proceed to investigate their properties. We note that while a very recent work (Mayes et al. 2025) applied a similar classification scheme, our detailed treatment of the preprocessed components represents a key novelty.

## 4. DEMOGRAPHICS OF THE ICL IN NEWCLUSTER

In this section, we present a demographic analysis of the BCG+ICL stars by origin. We note that the results represent a case study rather than a statistical sample, as our analysis is based on a single galaxy cluster. Nevertheless, we emphasize that our results remain informative due to the unprecedented spatial and mass resolutions of stellar particles at the galaxy cluster scale.

First, we present the classification results. Figure 7 illustrates the stellar density maps of the classified populations. Panel (a) shows the overall BCG+ICL sample, and panels (b) to (e) display the subsamples of the BCG+ICL: (b) those from surviving satellites, (c) from disrupted satellites, (d) born in-situ, and (e) from preprocessing.

The first noticeable feature is the tidal tail features stripped from surviving satellites (stripped, blue) in Figure 7(b). These stars are stripped inside the cluster and account for the largest fraction (33.4%) of the BCG+ICL components. In terms of morphological features, they exhibit prominent, elongated tail-like structures, as their host satellites are still alive at the final snapshot. We remind readers that these stripped components were relatively tightly bound to the host satellite until infall, as we separate the preprocessed populations that are already unbound before infall. The stripped mass of $4.8 \times 10^{11} M_{\odot}$ remains small compared to the stellar mass of the surviving satellites ($2.1 \times 10^{12} M_{\odot}$; Figure 4).

Figure 7(c) shows the stripped stars from satellites that have been disrupted (disrupted, green). Most disrupted satellites are stripped in the tidal field of the cluster or merge into the BCG, rather than into other satellites in the cluster, due to their higher relative velocity (Ghigna et al. 1998; Oh et al. 2018; Pearson et al.

2024). Moreover, these satellites typically enter the cluster earlier than the surviving ones, which gives them more time to reach dynamical relaxation, in addition to the smaller size of the cluster when they arrive. As a result, this population exhibits a smoother and less extended distribution, compared to the stripped ones. Likewise, we can expect that the distribution of stripped components would be relaxed toward lower redshift, as many currently surviving satellites complete mergers.

Figure 7(d) displays the in-situ components (red) that are predominantly born in the BCG, but also include a small minority of stars born in the ICL region. While BCGs are generally thought to acquire the majority of their stellar mass from mergers (Aragon-Salamanca et al. 1998; De Lucia & Blaizot 2007), we find that a comparable amount of stellar mass forms in-situ, which might be due to the current redshift ($z = 0.79$), with an increased ex-situ fraction expected toward $z = 0$. We also note that there can be merger-induced star formation when satellites pass through the BCG (Guo et al. 2011; Martin et al. 2017, 2021; Contreras-Santos et al. 2022; Aldás et al. 2025). This temporarily enhanced star formation can occur in both the BCG and its satellites; however, our method cannot separate them. Nevertheless, despite this potential contamination, we observe that the distribution of the majority of the in-situ components is highly concentrated, suggesting that they have kinematics similar to that of the BCG. Furthermore, we find that the majority of these stars formed deep within the BCG's potential well and do not migrate to the outskirts throughout their evolution. As these in-situ stars occupy the central region, this population can provide a clue to the physical demarcation between the BCG and the ICL, a long-standing challenge in the field, as mentioned in Section 1.

The final population, preprocessed components (yellow), is shown in Figure 7(e). We find that the spatial distribution closely resembles the overall distribution of the BCG+ICL stars. As the preprocessed components have experienced dynamical heating at earlier times, they exhibit a smoother and more extended distribution than the stripped components that have been tightly bound to host satellites until recently. In other words, the kinematics of preprocessed components are more strongly governed by the cluster potential than by the potential of their original host satellites. Given that the ICL is well regarded as a prospective and luminous tracer of dark matter (Yoo et al. 2024), we anticipate that preprocessed better follows the dark matter distribution than the other subsamples—Figure 7 shows stellar maps only, and we return to this point below.
<figure><img src="image_8.png" /><figcaption>**Figure 7.** Stellar density maps of classified BCG+ICL stars. Panel (a) presents all BCG+ICL stars within $1.5 R_{200}$. Panels (b)-(e) show the subsamples of the BCG+ICL stars. Each subsample is defined as follows in Section 3 and displayed in a different color. Panel (f) shows the mass fractions of the subsamples.</figcaption></figure>

We also present the mass fractions of the subsamples in Figure 7(f), with the color scheme of the subcomponents in panels (b)-(e). The most noticeable feature is that all channels contribute comparable fractions to the total BCG+ICL. The preprocessed component, for instance, accounts for ~12%, which is a substantial fraction, given that preprocessing is generally regarded as a secondary channel for ICL formation. The dominant origin of BCG+ICL is stars from satellites (surviving or disrupted), which account for more than half (~55%) of the total BCG+ICL. Montenegro-Taborda et al. (2023, 2025) have similarly investigated the demographics of the BCG+ICL using the TNG300 simulation, broadly consistent with our results. However, the detailed numbers slightly differ from ours, such as a low (~19%) in-situ fraction or a larger merger fraction than the stripped fraction from surviving satellites. This gap would be narrowed at lower redshifts as continued accretion increases the ex-situ fraction, and some of our currently ($z=0.79$) surviving satellites will eventually merge with the BCG. In terms of spatial extent, the outer parts of the BCG+ICL are mainly composed of accreted stars from satellites. This indicates that the

ICL, apart from the BCG, is predominantly formed by stellar stripping. Stars formed in situ comprise a third of the total BCG+ICL, but they rarely migrate to the outskirts. Contreras-Santos et al. (2024) also exhibit similar results with us, using GADGET-X and GIZMO-SIMBA codes based on THE THREE HUNDRED catalogs. The GADGET-X simulation indicated a low fraction of in-situ stars (12.5%) and a dominant fraction (75.5%) of “post-stripped” components. This may appear different from our results. However, the apparent difference is mostly coming from the difference in definitions. They define ICL as isolated from the BCG, which corresponds to the outer part of our BCG+ICL.

To quantify the similarity of spatial distributions, we present the density profile of dark matter and stellar components in Figure 8(a). The x-axis is the radial distance from the cluster center, normalized by $R_{200}$, and the y-axis is the three-dimensional density of each shell. The dark matter profile of the main cluster is shown in the black dashed line, while the profiles of stellar components are represented in solid lines, colored according to the legend. We also exhibit the stellar density profile divided by the dark matter density in Figure 8(b).
<figure><img src="image_9.png" /><figcaption>**Figure 8.** Density profiles of the BCG+ICL stars and dark matter. Panel (a) shows the density profile of dark matter (black dashed line) and all BCG+ICL stars (black solid line). The subsamples of BCG+ICL stars are shown in different colors as indicated in the legend. The gray dashed line indicates a slope of −2 for reference. Panel (b) presents the residuals with respect to the dark matter density profile.</figcaption></figure>

Again, the horizontal dashed line at unity in panel (b) indicates the dark matter density for reference.

The dark matter profile has a much shallower slope than the stellar components, and it shows a slight break near $r_s \simeq 0.15 R_{200}$, where $r_s$ is the scale radius of the NFW (Navarro, Frenk, & White 1997) profile formula, exhibiting an indication of a weak core in the central region. The overall profiles of dark matter and the BCG+ICL are consistent with those in other studies: dark matter profile being shallower than that of BCG+ICL (Sampaio-Santos et al. 2021; Alonso Asensio et al. 2020; Diego et al. 2023; Butler et al. 2025), and a roughly constant gradient in the ratio of the two density profiles (Alonso Asensio et al. 2020; Reina-Campos et al. 2022; Diego et al. 2023; Butler et al. 2025; Manuwal et al. 2025). As expected, the preprocessed population (yellow) has the shallowest slope in the density profile, making it the most similar to the dark matter profile of all the stellar components. It is not necessarily expected that the stellar distribution should follow that of dark matter. Most accreted stars inherit the biased phase-space of their parent satellites, whose orbits are reshaped by dynamical friction (Chandrasekhar 1943; Binney & Tremaine 2008). In particular, recently formed stars originate in dissipative, shock-triggered star formation and are born on a biased subset of orbits, with the gas itself following trajectories distinct from the collisionless

dark matter. Hence, the closer similarity of the preprocessed component to the dark matter can be understood as a mixing-time effect. Stars unbound well before cluster infall have completed multiple orbital periods and experienced substantial phase-angle mixing, erasing coherent tidal features (Helmi & White 1999).

## 5. STELLAR POPULATION STUDY

<figure><img src="image_10.png" /><figcaption>**Figure 9.** Physical properties of the BCG+ICL stars. The left panel shows the radial profile, and the right panel shows the probability density function. From top to bottom, the properties shown are stellar age, metallicity, and specific angular momentum. All properties are mass-weighted. In the radial profiles, the 16th–84th percentile ranges of the in-situ (red) and preprocessed components (yellow) are shown as shaded regions.</figcaption></figure>

Thanks to advanced observational techniques, our understanding of the ICL is expanding beyond simple abundance measurements. Many recent works have investigated detailed photometric properties. Montes et al. (2021) found changes in the slopes of the surface brightness and color profiles in Abell 85 using Subaru Hyper Suprime-Cam, and more recently, Ellien et al. (2025) and Englert et al. (2025) also investigated the
color profiles of the ICL of Abell 2390 and Abell 3667, respectively, finding that the profile can be varied by modeling methods or dynamical states of clusters. However, the general stellar population of the ICL has not yet been systematically examined (Montes & Trujillo 2018; Contini 2021). Having analyzed the origins and spatial distributions of the BCG+ICL subsamples, we now turn to their stellar population properties. Figure 9 displays the analysis of the stellar age (top), metallicity (middle), and specific angular momentum (bottom). The left (right) panel shows the mass-weighted radial profile (the probability density function). The y-axis indicates the values of the corresponding properties. The x-axis of the profile panels is the clustocentric distance normalized by $R_{200}$, and we also show the 16th–84th percentiles of the in-situ (red) and preprocessed components (yellow) as shaded regions to highlight the most distinctive populations.

First, we examine the stellar age (first row). At $z = 0.79$, the age distribution peaks at $\sim 4$ Gyr, which corresponds to $z \sim 2.3$ or the cosmic noon (Madau & Dickinson 2014). As seen in the radial profile, the central region ($< 0.3 R_{200}$) exhibits younger populations. The most noticeable difference among populations is that the disrupted (green) and preprocessed (yellow) components are dominated by old populations, while the in-situ (red) and stripped (blue) components are young. Specifically, although all populations have similar age peaks near $\sim 4$ Gyr, the in-situ and stripped components show pronounced tails toward younger ages in the PDF. This reflects recent or ongoing star formation inside the BCG and infalling satellites. On the contrary, the disrupted and preprocessed components rarely include young stars. The preprocessed components become unbound before the infall of satellites, and therefore have an inherent lower age limit set by the infall epoch of their host satellites. Likewise, the disrupted components are naturally older than the stripped components because their host galaxies stopped forming stars earlier than the host satellites of the stripped components. This echoes the earlier study of Murante et al. (2004), using a hydrodynamical cosmological simulation, which found that the diffuse components are older than other cluster stars on average.

Next, we examine the metallicity (second row). From the perspective of stellar evolution, younger stars are likely more metal-enriched due to the enrichment of recycled gas. Consistent with this expectation, in-situ components exhibit a distinct distribution compared to other populations, showing a more skewed inclination towards higher metallicity. In contrast, preprocessed components display the opposite behavior, forming a

<figure><img src="image_11.png" /><figcaption>**Figure 10.** Dominant populations in the phase-space diagram. For each phase-space bin, we color the most dominant populations following the text shown. The opacity represents how dominant a population is in each bin; colors range from white (32%) to vivid (100%). The black contours indicate 0.5 to 3 $\sigma$ of the density distribution for all BCG+ICL in phase-space. The white dashed line shows the escape velocity curve of an NFW ($c = 6.5$) halo, while the white solid line (shaded region) represents the median (16th–84th) escape velocity directly derived from the gravitational potential.</figcaption></figure>

predominantly metal-poor stellar population. Although the age or metallicity of ICL by different formation channels has hardly been studied, the metallicity of the BCG+ICL component in our sample is broadly consistent with previous studies (Contini et al. 2019; Mayes et al. 2025), exhibiting the drop of $\log(Z/Z_\odot)$ toward the outskirts. For the overall metallicity profile, the radial variation in metallicity is much more prominent than what can be inferred from the age difference.

We suggest that this difference in metallicity may provide a potential criterion for distinguishing the BCG and the ICL. By our definition, in-situ components in our study represent the in-situ stellar population, which is mainly born inside the BCG. In contrast, our preprocessed components represent the accreted stars that are already free from the potential of individual galaxies, reflecting the mass growth driven by the surrounding environment. This intrinsic difference in stellar populations can be traced by their metallicity, producing a potential demarcation of $Z \sim 0.5 Z_\odot$, where the two curves (red and yellow) intersect in the PDF panel, which roughly corresponds to a radius of $0.1 - 0.2 R_{200}$ in the profile panel.
Lastly, we examine the specific angular momentum ($j_*$) of each population to investigate the kinematic differences. The overall distributions, as shown in the PDF, are clearly distinct. For example, the preprocessed components exhibit higher $j_*$ values, peaking at $\sim 10^5$ kpc km s$^{-1}$, whereas the in-situ components peak around $10^3$ kpc km s$^{-1}$. However, this difference in $j_*$ is primarily a consequence of the different spatial distributions among populations, rather than an independent kinematic signature as shown in the radial profiles of $j_*$, which are nearly identical across all populations. This suggests that the populations do not differ intrinsically in angular momentum at a given radius, but rather in their radial distributions. Thus, $j_*$ does not provide an independent criterion for separating the BCG+ICL components.

Nevertheless, we can use the kinematic information of the BCG+ICL stars to qualitatively separate the populations. In Figure 10, we present the phase-space diagram, with each bin colored by its predominant component. For example, the blue-colored regions indicate where the stripped components (i.e., stripped stars from surviving satellites) dominate. We note that the pre-dominance does not necessarily mean that one component outnumbers the others at each bin. We therefore adjust the opacity following the fraction of the dominant population. In other words, the most white bins indicate the minimum fraction of $\sim 32\%$, while the most vivid bins represent 100%. The x-axis shows the clustocentric distance normalized by $R_{200}$, and the y-axis shows the three-dimensional velocity normalized by the maximum circular velocity ($V_{\text{cir}} = 825$ km s$^{-1}$) of the cluster. We also mark the $0.5-3\sigma$ contours (black) of the overall BCG+ICL components distribution and display the escape velocities derived from the NFW profile (white dashed line) and the gravitational potential (white solid line and shades). In the central region within $\sim 0.1 R_{200}$, most stars are in-situ (red) regardless of their velocity. The disrupted components (green) occupy the central ($d < 0.4 R_{200}$) region, but notably have lower velocities, typically less than $V_{\text{cir}}$. By contrast, the stripped components (blue) show a more extended distribution toward the outskirts and greater inhomogeneity in velocity space. The outermost region with high velocities and large distances is predominantly occupied by the preprocessed components (yellow). It is also noteworthy that most unbound BCG+ICL stars (i.e., above the escape velocity) are likely to be stripped and preprocessed components, indicating a first-infall phase. To summarize, a radius of $\sim 0.1 R_{200}$ can serve as an empirical demarcation separating the in-situ BCG stars from those

originating in satellites, which can be further divided using velocity information.

NEWCLUSTER computes the chemical evolution on-the-fly, which provides additional observable information. We analyze the chemical abundance ratios of each subsample, specifically [α/Fe] and [Fe/H]. [α/Fe] is calculated by averaging [O/Fe], [Mg/Fe], and [Si/Fe], all of which are directly available in NEWCLUSTER. We use the solar chemical abundances from Asplund et al. (2009). Figure 11 presents the chemical distributions. Panels (a) and (b) show the PDFs of [Fe/H] and [α/Fe], respectively, and panel (d) exhibits radial profiles. Panel (c) displays the [α/Fe]-[Fe/H] plane, and different colors of contours indicate the BCG+ICL subsamples following the legend. We also show the distribution of all components as a grayscale background. Again, the preprocessed components occupy a distinct region compared to other components. They show lower [Fe/H] and higher [α/Fe], indicating older stellar ages (Tinsley 1979). Combining with the result that the pre-processed components are old and metal-poor, as seen in Figure 9, this suggests that the preprocessed components formed in relatively sparse regions where gas recycling and chemical enrichment are inefficient. The differences are still prominent in the radial profile. In Figure 11(d), unlike the results shown in Figure 9, the [α/Fe] and [Fe/H] profiles of preprocessed components show offsets from the mean profile (black) at all radii. We also note that our result—super-solar at the center and [Fe/H] < −0.5 at the outskirts—is consistent with the observed profile (Montes & Trujillo 2018). In summary, our results highlight that a detailed analysis using various chemical elements can offer a more independent way to separate the stellar populations of the ICL.

## 6. RELATION BETWEEN ORBITAL MOTIONS AND STRIPPING

As shown in previous sections, most of the BCG+ICL stars originate from satellite galaxies, both surviving and disrupted, whose spatial and velocity distributions differ. The distinct properties among different stellar populations suggest that their formation is linked to the dynamical history of their host satellites. In this section, we evaluate the relationship between the orbital parameters of the host satellites in the cluster and the efficiency of stellar stripping.

The measured orbital parameters are: 1) “time since infall” (TSI), the time since the first infall into the cluster; 2) min($r_p$), the minimum pericenter distance; 3) $N_{\text{orbit}}$, the number of orbits (i.e., half of the number of
<figure><img src="image_12.png" /><figcaption>**Figure 11.** Distributions of [α/Fe] and [Fe/H] for the BCG+ICL stars. Panel (c) shows the [α/Fe]-[Fe/H] diagram for each subsample. The contours indicate 1, 1.5, and 2 σ regions, while the shaded background shows the distribution of all BCG+ICL stars. Panels (a) and (b) display the probability density functions of [Fe/H] (left) and [α/Fe] (right). Panel (d) presents the corresponding radial profiles. The preprocessed component notably shows a distinctive distribution and profile, reflecting its different properties in the stellar population.</figcaption></figure>

peri- and apocenters passages); 4) $\beta \equiv 1 - v_t^2 / (2v_r^2)^{1/2}$, the characterized orbital anisotropy parameter (An & Evans 2006; Binney & Tremaine 2008), where $v_t$ and $v_r$ denote the tangential and radial velocities, respectively; and 5) $M_*/M_h$, the stellar-to-halo mass ratio measured at $t_{inf}$, where $M_*$ is the stellar mass of the satellite galaxy at its infall and $M_h$ is the halo mass of the main cluster progenitor at the corresponding epoch. We measure these parameters from $t_{inf}$ to the final ($z = 0.79$) snapshot. To quantitatively assess the effect of the parameters, we define the stripped fraction ($f_{strip}$) as the fraction of the satellite (both surviving and disrupted ones) member stars at infall that are found in the BCG+ICL at the final snapshot.

We present these relationships in the corner plot shown in Figure 12. The diagonal panels show the one-

dimensional relations between each orbital parameter and $f_{strip}$, with the error bars indicating the 16th–84th percentile range. The off-diagonal panels show $f_{strip}$ as a function of two orbital parameters. $f_{strip}$ is color-coded from blue to red with increasing values. The three parameters TSI, $\min(r_p)$, and $N_{orbit}$ have clear relations with $f_{strip}$, consistent with previous studies of tidal stripping in cluster simulations (Smith et al. 2016; Rhee et al. 2017; Han et al. 2018). $f_{strip}$ increases as satellites spend more time in the cluster, approach closer to the cluster center, or complete more orbits, although the scatters are large. By contrast, $\beta$ and $M_*/M_h$ show no clear dependencies. The most significant relationship is between $\min(r_p)$ and TSI, indicating that satellites are more likely to lose their member stars if they fall into the cluster earlier and pass closer to the cluster center (e.g. Taylor & Babul 2001). These two parameters consistently exhibit clear dependencies, as seen in the left two columns. We also note that, for a fixed TSI, $M_*/M_h$ at infall shows a marginal negative trend with $f_{strip}$, and satellites with higher $\beta$ (i.e., more ra-

<footnote n="12">¹² Although this parameter is originally used to analyze the ensemble anisotropy of orbital systems, we extend this to the satellite galaxy orbiting the cluster for tracing the orbital shape.</footnote>
<figure><img src="image_13.png" /><figcaption>**Figure 12.** Relationships between orbital parameters and the stripped fraction ($f_{strip}$). The parameters considered are the time since infall (TSI), the minimum pericenter distance (min($r_p$)), the number of orbits ($N_{orbit}$), the orbital anisotropy parameter ($\beta$), and the stellar-to-halo mass ratio at infall (M*/M$_h$), as described in the text. The diagonal panels show the relation between each parameter and $f_{strip}$, including the 16th–84th percentile ranges. The off-diagonal panels present comparisons between parameter pairs, with colors indicating the value of $f_{strip}$.</figcaption></figure>

dial orbits) are likely to be more stripped. These results largely agree with the previous study of Martin et al. (2024), which studied the effects of orbital motions and found that satellites with more radial orbits and more orbital cycles exhibit larger stripped fractions in idealized environments.

Taken together, this analysis suggests that the efficiency of stellar stripping from satellites, specifically into the BCG+ICL components, is primarily driven by the time since infall and the depth of penetration into the cluster potential, rather than by their orbital shapes or the stellar-to-halo mass ratio at infall. Satellites that

have resided in the cluster longer, reached closer pericenter distances, and completed more orbital passages are more likely to lose their members into the BCG+ICL. This finding, which highlights the central role of orbital dynamics in shaping the origin of the ICL, is consistent with previous theoretical studies of satellite evolution (Puchwein et al. 2010; Contini et al. 2018; Martin et al. 2024).

## 7. CAVEATS

In this section, we discuss the current caveats of our study, including a limited sample, the difficulty of direct comparison with observations, and methodological
arbitrariness. These caveats are not trivial but have the value of emphasizing the complicated nature of the ICL and guiding us in future research directions.

## 7.1. Limiting sample and dynamical state

The main limitation of our study is that it focuses on only one galaxy cluster and does not capture the cluster-to-cluster variation. Since galaxy clusters have a wide range of mass and dynamical states (Kravtsov & Borgani 2012), it would be ideal to investigate a larger sample of clusters to ensure meaningful statistics. In particular, the dynamical state of the cluster should impact the demographics and origins of the ICL because the formation and evolution of the ICL are naturally affected by the assembly history of the cluster itself (Rudick et al. 2011; Mostoghiu et al. 2019; Contini et al. 2023). For example, Contini et al. (2023) used a semi-analytic model to evaluate halo concentration, which is directly connected to the dynamical state. They found that clusters with higher halo concentrations or dynamically older have higher ICL fractions, especially for low-mass halos. Also, Yoo et al. (2022) applied the weighted overlap coefficient method that quantifies the similarity of two-dimensional spatial distributions to the Galaxy Replacement Technique (Chun et al. 2022) simulation, and found a stronger similarity in the relaxed cluster samples than dynamically younger ones. Montes & Trujillo (2019) also suggested that the baryonic distribution traced from X-rays does not match the total mass distribution in disturbed clusters, although they did not find a similar offset in the ICL distribution. Observationally, the results regarding color and metallicity have not yet converged, which may also be related to dynamical states (Jiménez-Teja & Dupke 2016; Contini 2021, see references therein). Our primary cluster has experienced recent mergers and is about to merge, indicating that it is a relatively unrelaxed system, which could lead to biased results.

## 7.2. Current limitation of direct comparison with observations

NEWCLUSTER has a high stellar mass resolution of $2 \times 10^4$ M$_\odot$ and a spatial resolution of 68 pc compared to other cluster-scale simulations—such as HYDRANGEA (Bahé et al. 2017), ROMULUSC (Tremmel et al. 2019), or TNG-CLUSTER (Nelson et al. 2024)—, which allows faint structures and outskirts to be resolved. Han et al. (2025a), the introductory work of NEWCLUSTER, showcased the mock radiative transfer images of galaxies in NEWCLUSTER, which resolve galactic structures in detail. At the same time, they also introduced the mock observations of low-surface brightness features at

$z = 0.8$, and demonstrated that modern instruments can resolve faint structures. However, NEWCLUSTER is just passing $z \sim 0.56$, i.e., relatively high in redshift for practical observations, and a direct comparison with observations at this redshift is challenging. For example, Ko & Jee (2018) used a deep Hubble Space Telescope Wide Field Camera 3 (HST/WFC3) observation, but the surface brightness limit was still $\sim 29$ mag arcsec$^{-2}$. Moreover, Dacunha et al. (2025) pointed out that the data reduction process, such as background subtraction, should be treated carefully to obtain reliable ICL measurements.

In addition, we utilize the three-dimensional position and time series data, which are fully available only in simulations. Although we separate the BCG+ICL stars into four different channels and show that they exhibit different stellar properties, this approach is technically infeasible to replicate with current observational techniques. For example, our preprocessed component has the most distinct properties, such as old age, low metallicity, and far distance from the cluster center in the phase-space. However, in observations, the quantities to estimate those require high-resolution spectroscopy, which makes the direct comparison challenging.

A fair and direct comparison requires constructing mock images under the same conditions as observations. Many recent instruments, such as the James Webb Space Telescope (JWST), the Vera Rubin Observatory's Legacy Survey of Space and Time (Vera Rubin/LSST), and the Euclid mission, are capable of observing the low-surface brightness features, and surveys are planned to study the ICL (Montes & Trujillo 2022; Martis et al. 2024; Martin et al. 2022; Brough et al. 2024; Englert et al. 2025; Kluge et al. 2025; Euclid Collaboration et al. 2025). To meet these observational capabilities, we can perform mock observations by combining radiative transfer with on-the-fly dust calculation—which is provided by NEWCLUSTER—, enabling direct and fair comparisons between simulations and future surveys. For instance, Byun et al. (2025) demonstrated that mock observations using on-the-fly dust calculations of NEWCLUSTER can reproduce the galaxy's morphology of the early Universe observed by JWST. Furthermore, as resolved spectroscopic studies of low-surface brightness are now emerging (Adami et al. 2016; Gu et al. 2020; Spavone et al. 2022), we can even build mock integral field unit (IFU) data using the simulation. While observing techniques gradually reach farther with time, simulations (NEWCLUSTER in this case) approach nearer to the local Universe. The quantities derived in our study may not be directly verified through
observations yet, but it is promising that both observations and simulations are nearing their front lines.

## 8. CONCLUSIONS

The main goal of this paper is to present the first results of an ICL population study with high resolution and to highlight the potential of next-generation cluster simulations for ICL studies. We used the high-resolution zoom-in cluster simulation, NEWCLUSTER, at $z = 0.79$ to examine the origins and stellar population properties of the BCG+ICL. To ensure a complete census of stellar particles, we developed a dedicated merger tree algorithm and tracked individual stars across closely spaced simulation outputs.

Our main results can be summarized as follows.

1. We found that tracking the full orbital history of satellite galaxies is essential for robustly determining the boundness of stellar particles. The member stellar particles of satellites are tightly bound to their host before infall, but they are then significantly stripped as the satellites orbit within the cluster. However, the survival—or the cluster infall time—of host satellites strongly influences the final distribution of stripped stars.

2. We derived the demographics of the BCG+ICL based on their origin: stripped (BCG+ICL stripped from surviving satellites), disrupted (BCG+ICL stripped and merged from disrupted satellites), in-situ (formed in situ within the BCG), and preprocessed (unbound from their host satellites prior to the cluster infall). The majority of the BCG+ICL components originate from satellite galaxies (~60%), including surviving and disrupted ones.

3. Our classification of the preprocessed component reveals its distinctive properties and shows that it comprises a significant fraction (~12%) of the total BCG+ICL. They exhibit shallower density profiles than the other BCG+ICL components and are more closely aligned with the dark matter profile. The stellar population of the preprocessed components shows more dramatic differences in age, metallicity, and alpha abundance (e.g., older, metal-poor, and enhanced [α/Fe]), suggesting this population as a key fossil record of the cluster's dynamical history. We also introduced a clear division of the subsamples in the phase-space, which provides a useful diagnostic for the ICL composition.

4. We investigated the effect of orbital parameters on the fraction of stripped stars ($f_{strip}$). We found

that the time since infall and the minimum period center passage distance are the primary drivers in determining $f_{strip}$. The satellite galaxies that join the cluster earlier or reach closer to the BCG in their orbital motion are stripped more. Also, a larger number of orbital cycles ($N_{orbit}$) increases the fraction of stripped stars. The shapes of orbits and the stellar-to-halo mass ratio at infall have only a marginal effect on $f_{strip}$.

While we focus on a single cluster and acknowledge the limitations in the investigation, these significant caveats highlight the intrinsic complexity of ICL studies rather than diminishing the significance of our results. Nevertheless, the methodological novelty of our work lies in our detailed classification of stellar populations, tracking billions of particles across dense temporal outputs to separate the ICL into four distinct origins. Our work suggests new directions for future ICL studies, particularly in understanding the role of preprocessing and in using detailed stellar properties to disentangle formation channels.

In closing, we emphasize the prospects of the upcoming low-redshift data from NEWCLUSTER. The primary and secondary halos in NEWCLUSTER are currently in the process of merging ($z \sim 0.56$), and the merged final cluster is expected to undergo violent mixing and relaxation. This stage provides an opportunity to investigate the migration of stars among the BCG, satellites, and the ICL through the cluster merger. Crucially, because we follow the same cluster across time, we can quantify the ICL's evolution within one halo rather than combining heterogeneous samples. As the evolutionary trend of the ICL fraction in a single galaxy cluster has not yet converged, we can disentangle the methodological differences and cosmological dimming from the true evolution, which provides the prediction of the later growth of the ICL as a benchmark. The forthcoming low-redshift data from NEWCLUSTER will provide a unique laboratory for bridging theoretical and observational studies of the ICL.

## ACKNOWLEDGEMENTS

This work was granted access to the HPC resources of KISTI under the allocations KSC-2021-CRE-0486, KSC-2022-CRE-0088, KSC-2022-CRE-0344, KSC-2022-CRE-0409, KSC-2023-CRE-0343, KSC-2024-CHA-0009, and KSC-2025-CRE-0031 and of GENCI under the allocations A0150414625 and A0180416216. 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). J.R. acknowledges support from the Insti-
tut de Physique des deux infinis de Sorbonne Université
and by the ANR grant ANR-19-CE31-0017 of the French
Agence Nationale de la Recherche. J.L. is supported
by the National Research Foundation (NRF) of Korea
grant funded by the Korea government (MSIT, RS-2021-
NR061998). J.L. also acknowledges the support of the
NRF of Korea grant funded by the Korea government

(MSIT, RS2022-NR068800). G.M. acknowledges sup-
port from the UK STFC under grant ST/X000982/1.
J.K. is supported by KIAS individual Grant (KG039603)
at the Korea Institute for Advanced Study.". E.C. and
S.J. acknowledge support from the Korean National Re-
search Foundation (RS-2023-00241934) T.K. was sup-
ported by the Yonsei Fellowship, funded by Lee Youn
Jae.

REFERENCES

Adami, C., Slezak, E., Durret, F., et al. 2005, A&A, 429, 39, doi: 10.1051/0004-6361:20041322

Adami, C., Pompei, E., Sadibekova, T., et al. 2016, A&A, 592, A7, doi: 10.1051/0004-6361/201526831

Ahad, S. L., Bahé, Y. M., & Hoekstra, H. 2023, MNRAS, 518, 3685, doi: 10.1093/mnras/stac3357

Aldás, F., Gómez, F. A., Vega-Martínez, C., Zenteno, A., & Carrasco, E. R. 2025, A&A, 699, A313, doi: 10.1051/0004-6361/202451801

Alonso Asensio, I., Dalla Vecchia, C., Bahé, Y. M., Barnes, D. J., & Kay, S. T. 2020, MNRAS, 494, 1859, doi: 10.1093/mnras/staa861

An, J. H., & Evans, N. W. 2006, AJ, 131, 782, doi: 10.1086/499305

Annunziatella, M., Mercurio, A., Biviano, A., et al. 2016, A&A, 585, A160, doi: 10.1051/0004-6361/201527399

Aragon-Salamanca, A., Baugh, C. M., & Kauffmann, G. 1998, MNRAS, 297, 427, doi: 10.1046/j.1365-8711.1998.01495.x

Asplund, M., Grevesse, N., Sauval, A. J., & Scott, P. 2009, ARA&A, 47, 481, doi: 10.1146/annurev.astro.46.060407.145222

Aubert, D., Pichon, C., & Colombi, S. 2004, MNRAS, 352, 376, doi: 10.1111/j.1365-2966.2004.07883.x

Bahé, Y. M., Barnes, D. J., Dalla Vecchia, C., et al. 2017, MNRAS, 470, 4186, doi: 10.1093/mnras/stx1403

Behroozi, P. S., Wechsler, R. H., Wu, H.-Y., et al. 2013, ApJ, 763, 18, doi: 10.1088/0004-637X/763/1/18

Binney, J., & Tremaine, S. 2008, Galactic Dynamics: Second Edition

Bonaventura, N. R., Webb, T. M. A., Muzzin, A., et al. 2017, MNRAS, 469, 1259, doi: 10.1093/mnras/stx722

Borgani, S., & Guzzo, L. 2001, Nature, 409, 39, doi: 10.1038/409039A010.1038/35051000

Brough, S., Ahad, S. L., Bahé, Y. M., et al. 2024, MNRAS, 528, 771, doi: 10.1093/mnras/stad3810

Brown, H. J., Martin, G., Pearce, F. R., et al. 2024, MNRAS, 534, 431, doi: 10.1093/mnras/stae2084

Burke, C., Hilton, M., & Collins, C. 2015, MNRAS, 449, 2353, doi: 10.1093/mnras/stv450

Butler, J., Martin, G., Hatch, N. A., et al. 2025, MNRAS, 539, 2279, doi: 10.1093/mnras/staf615

Byun, G.-H., Jang, J. K., Scofield, Z. P., et al. 2025, arXiv e-prints, arXiv:2508.18374. https://arxiv.org/abs/2508.18374

Cadiou, C., Dubois, Y., & Pichon, C. 2019, A&A, 621, A96, doi: 10.1051/0004-6361/201834496

Cha, S., & Jee, M. J. 2023, ApJ, 951, 140, doi: 10.3847/1538-4357/acd111

Chabrier, G. 2003, PASP, 115, 763, doi: 10.1086/376392

Chandrasekhar, S. 1943, ApJ, 97, 255, doi: 10.1086/144517

Chandro-Gómez, Á., Lagos, C. d. P., Power, C., et al. 2025, MNRAS, 539, 776, doi: 10.1093/mnras/staf519

Chun, K., Shin, J., Ko, J., Smith, R., & Yoo, J. 2024, ApJ, 969, 142, doi: 10.3847/1538-4357/ad4a52

Chun, K., Shin, J., Smith, R., Ko, J., & Yoo, J. 2022, ApJ, 925, 103, doi: 10.3847/1538-4357/ac2cbe

Clowe, D., Bradač, M., Gonzalez, A. H., et al. 2006, ApJL, 648, L109, doi: 10.1086/508162

Contini, E. 2021, Galaxies, 9, 60, doi: 10.3390/galaxies9030060

Contini, E., Chen, H. Z., & Gu, Q. 2022, ApJ, 928, 99, doi: 10.3847/1538-4357/ac57c4

Contini, E., De Lucia, G., Villalobos, Á., & Borgani, S. 2014, MNRAS, 437, 3787, doi: 10.1093/mnras/stt2174

Contini, E., Han, S., Jeon, S., Rhee, J., & Yi, S. K. 2024a, ApJL, 962, L10, doi: 10.3847/2041-8213/ad21e2

Contini, E., Jeon, S., Rhee, J., Han, S., & Yi, S. K. 2023, ApJ, 958, 72, doi: 10.3847/1538-4357/acfd25

Contini, E., Rhee, J., Han, S., Jeon, S., & Yi, S. K. 2024b, AJ, 167, 7, doi: 10.3847/1538-3881/ad0894

Contini, E., Yi, S. K., & Jeon, S. 2024c, arXiv e-prints, arXiv:2404.01560, doi: 10.48550/arXiv.2404.01560

Contini, E., Yi, S. K., Jeon, S., & Rhee, J. 2024d, ApJS, 274, 41, doi: 10.3847/1538-4365/ad70ac

Contini, E., Yi, S. K., & Kang, X. 2018, MNRAS, 479, 932, doi: 10.1093/mnras/sty1518
—. 2019, ApJ, 871, 24, doi: 10.3847/1538-4357/aaf41f

Contreras-Santos, A., Knebe, A., Pearce, F., et al. 2022,
MNRAS, 511, 2897, doi: 10.1093/mnras/stac275

Contreras-Santos, A., Knebe, A., Cui, W., et al. 2024,
A&A, 683, A59, doi: 10.1051/0004-6361/202348474

Dacunha, T., Mansfield, P., & Wechsler, R. 2025, arXiv e-prints, arXiv:2508.02837, doi: 10.48550/arXiv.2508.02837

De Lucia, G., & Blaizot, J. 2007, MNRAS, 375, 2, doi: 10.1111/j.1365-2966.2006.11287.x

Diego, J. M., Pascale, M., Frye, B., et al. 2023, A&A, 679, A159, doi: 10.1051/0004-6361/202345868

Dolag, K., Murante, G., & Borgani, S. 2010, MNRAS, 405, 1544, doi: 10.1111/j.1365-2966.2010.16583.x

Dolag, K., Remus, R.-S., Valenzuela, L. M., et al. 2025, arXiv e-prints, arXiv:2504.01061, doi: 10.48550/arXiv.2504.01061

Donzelli, C. J., Muriel, H., & Madrid, J. P. 2011, ApJS, 195, 15, doi: 10.1088/0067-0049/195/2/15

Dubois, Y., Volonteri, M., & Silk, J. 2014a, MNRAS, 440, 1590, doi: 10.1093/mnras/stu373

Dubois, Y., Pichon, C., Welker, C., et al. 2014b, MNRAS, 444, 1453, doi: 10.1093/mnras/stu1227

Dubois, Y., Beckmann, R., Bournaud, F., et al. 2021, A&A, 651, A109, doi: 10.1051/0004-6361/202039429

Elahi, P. J., Poulton, R. J. J., Tobar, R. J., et al. 2019, PASA, 36, e028, doi: 10.1017/pasa.2019.18

Ellien, A., Slezak, E., Martinet, N., et al. 2021, A&A, 649, A38, doi: 10.1051/0004-6361/202038419

Ellien, A., Montes, M., Ahad, S. L., et al. 2025, A&A, 698, A134, doi: 10.1051/0004-6361/202554460

Englert, A. M., Dell'Antonio, I., & Montes, M. 2025, ApJL, 989, L2, doi: 10.3847/2041-8213/ade8f1

Euclid Collaboration, Bellhouse, C., Golden-Marx, J. B., et al. 2025, A&A, 698, A14, doi: 10.1051/0004-6361/202553887

Fujita, Y. 2004, PASJ, 56, 29, doi: 10.1093/pasj/56.1.29

Genel, S., Vogelsberger, M., Nelson, D., et al. 2013, MNRAS, 435, 1426, doi: 10.1093/mnras/stt1383

Ghigna, S., Moore, B., Governato, F., et al. 1998, MNRAS, 300, 146, doi: 10.1046/j.1365-8711.1998.01918.x

Golden-Marx, J. B., Zhang, Y., Ogando, R. L. C., et al. 2025, MNRAS, 538, 622, doi: 10.1093/mnras/staf277

Gonzalez, A. H., George, T., Connor, T., et al. 2021, MNRAS, 507, 963, doi: 10.1093/mnras/stab2117

Gonzalez, A. H., Zabludoff, A. I., & Zaritsky, D. 2005, ApJ, 618, 195, doi: 10.1086/425896

Gu, M., Conroy, C., Law, D., et al. 2020, ApJ, 894, 32, doi: 10.3847/1538-4357/ab845c

Guo, Q., White, S., Boylan-Kolchin, M., et al. 2011, MNRAS, 413, 101, doi: 10.1111/j.1365-2966.2010.18114.x

Han, S., Dubois, Y., Lee, J., et al. 2025a, ApJ, 978, 96, doi: 10.3847/1538-4357/ad98f4

Han, S., Smith, R., Choi, H., et al. 2018, ApJ, 866, 78, doi: 10.3847/1538-4357/aadfe2

Han, S., Yi, S. K., Dubois, Y., et al. 2025b, arXiv e-prints, arXiv:2507.06301, doi: 10.48550/arXiv.2507.06301

Helmi, A., & White, S. D. M. 1999, MNRAS, 307, 495, doi: 10.1046/j.1365-8711.1999.02616.x

Iwamoto, K., Brachwitz, F., Nomoto, K., et al. 1999, ApJS, 125, 439, doi: 10.1086/313278

Jeon, S. 2025, syj3514/YoungTree: YoungTree: multi-snapshot merger tree builder (v1.0.0), v1.0.0, Zenodo, doi: 10.5281/zenodo.17793360

Jiménez-Teja, Y., & Dupke, R. 2016, ApJ, 820, 49, doi: 10.3847/0004-637X/820/1/49

Jiménez-Teja, Y., Román, J., HyeongHan, K., et al. 2025, A&A, 694, A216, doi: 10.1051/0004-6361/202452384

Joo, H., & Jee, M. J. 2023, Nature, 613, 37, doi: 10.1038/s41586-022-05396-4

Joo, H., Jee, M. J., Kim, J., et al. 2025, ApJ, 990, 96, doi: 10.3847/1538-4357/adf4d0

Jung, S. L., Choi, H., Wong, O. I., et al. 2018, ApJ, 865, 156, doi: 10.3847/1538-4357/aadda2

Kimm, T., & Cen, R. 2014, ApJ, 788, 121, doi: 10.1088/0004-637X/788/2/121

Kimm, T., Katz, H., Haehnelt, M., et al. 2017, MNRAS, 466, 4826, doi: 10.1093/mnras/stx052

Kimmig, L. C., Brough, S., Dolag, K., et al. 2025, A&A, 700, A95, doi: 10.1051/0004-6361/202554777

Kluge, M., Hatch, N. A., Montes, M., et al. 2025, A&A, 697, A13, doi: 10.1051/0004-6361/202450772

Ko, J., & Jee, M. J. 2018, ApJ, 862, 95, doi: 10.3847/1538-4357/aacbda

Kobayashi, C., Umeda, H., Nomoto, K., Tominaga, N., & Ohkubo, T. 2006, ApJ, 653, 1145, doi: 10.1086/508914

Kochanek, C. S. 2006, in Saas-Fee Advanced Course 33: Gravitational Lensing: Strong, Weak and Micro, ed. G. Meylan, P. Jetzer, P. North, P. Schneider, C. S. Kochanek, & J. Wambsganss, 91–268

Komatsu, E., Smith, K. M., Dunkley, J., et al. 2011, ApJS, 192, 18, doi: 10.1088/0067-0049/192/2/18

Kravtsov, A. V., & Borgani, S. 2012, ARA&A, 50, 353, doi: 10.1146/annurev-astro-081811-125502

Kravtsov, A. V., Vikhlinin, A. A., & Meshcheryakov, A. V. 2018, Astronomy Letters, 44, 8, doi: 10.1134/S1063773717120015
Lovell, M. R., Pillepich, A., Engler, C., et al. 2025, arXiv e-prints, arXiv:2509.07078, doi: 10.48550/arXiv.2509.07078

Madau, P., & Dickinson, M. 2014, ARA&A, 52, 415, doi: 10.1146/annurev-astro-081811-125615

Maeder, A., & Meynet, G. 2000, A&A, 361, 159, doi: 10.48550/arXiv.astro-ph/0006405

Manuwal, A., Avila-Reese, V., Montenegro-Taborda, D., Rodriguez-Gomez, V., & Cervantes Sodi, B. 2025, MNRAS, 543, 4020, doi: 10.1093/mnras/staf1717

Martin, G., Kaviraj, S., Devriendt, J. E. G., et al. 2017, MNRAS, 472, L50, doi: 10.1093/mnrasl/slx136

Martin, G., Pearce, F. R., Hatch, N. A., et al. 2024, MNRAS, 535, 2375, doi: 10.1093/mnras/stae2488

Martin, G., Jackson, R. A., Kaviraj, S., et al. 2021, MNRAS, 500, 4937, doi: 10.1093/mnras/staa3443

Martin, G., Bazkiaei, A. E., Spavone, M., et al. 2022, MNRAS, 513, 1459, doi: 10.1093/mnras/stac1003

Martínez-Lombilla, C., Brough, S., Montes, M., et al. 2023, MNRAS, 518, 1195, doi: 10.1093/mnras/stac3119

Martis, N. S., Sarrouh, G. T. E., Willott, C. J., et al. 2024, ApJ, 975, 76, doi: 10.3847/1538-4357/ad7735

Mayes, R. J., Gómez, F. A., & Monachesi, A. 2025, arXiv e-prints, arXiv:2506.16645, doi: 10.48550/arXiv.2506.16645

Mihos, C. 2003, arXiv e-prints, astro, doi: 10.48550/arXiv.astro-ph/0305512

Mihos, J. C. 2016, in IAU Symposium, Vol. 317, The General Assembly of Galaxy Halos: Structure, Origin and Evolution, ed. A. Bragaglia, M. Arnaboldi, M. Rejkuba, & D. Romano, 27-34, doi: 10.1017/S1743921315006857

Mihos, J. C., Harding, P., Feldmeier, J., & Morrison, H. 2005, ApJL, 631, L41, doi: 10.1086/497030

Monaco, P., Murante, G., Borgani, S., & Fontanot, F. 2006, ApJL, 652, L89, doi: 10.1086/510236

Montenegro-Taborda, D., Avila-Reese, V., Rodriguez-Gomez, V., Manuwal, A., & Cervantes-Sodi, B. 2025, MNRAS, 537, 3954, doi: 10.1093/mnras/staf271

Montenegro-Taborda, D., Rodriguez-Gomez, V., Pillepich, A., et al. 2023, MNRAS, 521, 800, doi: 10.1093/mnras/stad586

Montes, M. 2022, Nature Astronomy, 6, 308, doi: 10.1038/s41550-022-01616-z

Montes, M., Brough, S., Owers, M. S., & Santucci, G. 2021, ApJ, 910, 45, doi: 10.3847/1538-4357/abddb6

Montes, M., & Trujillo, I. 2014, ApJ, 794, 137, doi: 10.1088/0004-637X/794/2/137

—. 2018, MNRAS, 474, 917, doi: 10.1093/mnras/stx2847

—. 2019, MNRAS, 482, 2838, doi: 10.1093/mnras/sty2858

—. 2022, ApJL, 940, L51, doi: 10.3847/2041-8213/ac98c5

Mostoghiu, R., Knebe, A., Cui, W., et al. 2019, MNRAS, 483, 3390, doi: 10.1093/mnras/sty3306

Murante, G., Giovalli, M., Gerhard, O., et al. 2007, MNRAS, 377, 2, doi: 10.1111/j.1365-2966.2007.11568.x

Murante, G., Arnaboldi, M., Gerhard, O., et al. 2004, ApJL, 607, L83, doi: 10.1086/421348

Navarro, J. F., Frenk, C. S., & White, S. D. M. 1997, ApJ, 490, 493, doi: 10.1086/304888

Nelson, D., Pillepich, A., Ayromlou, M., et al. 2024, A&A, 686, A157, doi: 10.1051/0004-6361/202348608

Nelson, D., Springel, V., Pillepich, A., et al. 2019, Computational Astrophysics and Cosmology, 6, 2, doi: 10.1186/s40668-019-0028-x

Oh, S., Kim, K., Lee, J. H., et al. 2018, ApJS, 237, 14, doi: 10.3847/1538-4365/aacd47

Park, C., Lee, J., Kim, J., et al. 2022, ApJ, 937, 15, doi: 10.3847/1538-4357/ac85b5

Pearson, W. J., Santos, D. J. D., Goto, T., et al. 2024, A&A, 686, A94, doi: 10.1051/0004-6361/202349034

Pfister, H., Volonteri, M., Dubois, Y., Dotti, M., & Colpi, M. 2019, MNRAS, 486, 101, doi: 10.1093/mnras/stz822

Pillepich, A., Nelson, D., Hernquist, L., et al. 2018, MNRAS, 475, 648, doi: 10.1093/mnras/stx3112

Poole, G. B., Mutch, S. J., Croton, D. J., & Wyithe, S. 2017, MNRAS, 472, 3659, doi: 10.1093/mnras/stx2233

Puchwein, E., Springel, V., Sijacki, D., & Dolag, K. 2010, MNRAS, 406, 936, doi: 10.1111/j.1365-2966.2010.16786.x

Reina-Campos, M., Trujillo-Gomez, S., Pfeffer, J. L., et al. 2022, arXiv e-prints, arXiv:2204.11861, doi: 10.48550/arXiv.2204.11861

Rhee, J., Smith, R., Choi, H., et al. 2017, ApJ, 843, 128, doi: 10.3847/1538-4357/aa6d6c

Rudick, C. S., Mihos, J. C., Frey, L. H., & McBride, C. K. 2009, ApJ, 699, 1518, doi: 10.1088/0004-637X/699/2/1518

Rudick, C. S., Mihos, J. C., & McBride, C. K. 2011, ApJ, 732, 48, doi: 10.1088/0004-637X/732/1/48

Sampaio-Santos, H., Zhang, Y., Ogando, R. L. C., et al. 2021, MNRAS, 501, 1300, doi: 10.1093/mnras/staa3680

Schaller, G., Schaerer, D., Meynet, G., & Maeder, A. 1992, A&AS, 96, 269

Smith, R., Choi, H., Lee, J., et al. 2016, ApJ, 833, 109, doi: 10.3847/1538-4357/833/1/109

Spavone, M., Iodice, E., D’Ago, G., et al. 2022, A&A, 663, A135, doi: 10.1051/0004-6361/202243290

Srisawat, C., Knebe, A., Pearce, F. R., et al. 2013, MNRAS, 436, 150, doi: 10.1093/mnras/stt1545

Taylor, J. E., & Babul, A. 2001, ApJ, 559, 716, doi: 10.1086/322276
Teyssier, R. 2002, A&A, 385, 337,
doi: 10.1051/0004-6361:20011817

Tinsley, B. M. 1979, ApJ, 229, 1046, doi: 10.1086/157039

Tremmel, M., Quinn, T. R., Ricarte, A., et al. 2019,
MNRAS, 483, 3336, doi: 10.1093/mnras/sty3336

Trujillo, I., & Fliri, J. 2016, ApJ, 823, 123,
doi: 10.3847/0004-637X/823/2/123

Tweed, D., Devriendt, J., Blaizot, J., Colombi, S., & Slyz, A. 2009, A&A, 506, 647,
doi: 10.1051/0004-6361/200911787

Willman, B., Governato, F., Wadsley, J., & Quinn, T. 2004,
MNRAS, 355, 159, doi: 10.1111/j.1365-2966.2004.08312.x

Yoo, J., Shin, J., Hwang, H. S., et al. 2025, ApJ, 988, 229,
doi: 10.3847/1538-4357/ade66f

Yoo, J., Ko, J., Sabiu, C. G., et al. 2022, ApJS, 261, 28,
doi: 10.3847/1538-4365/ac7142

Yoo, J., Park, C., Sabiu, C. G., et al. 2024, ApJ, 965, 145,
doi: 10.3847/1538-4357/ad2df8

Zibetti, S., White, S. D. M., Schneider, D. P., &  
Brinkmann, J. 2005, MNRAS, 358, 949,  
doi: 10.1111/j.1365-2966.2005.08817.x

Zwicky, F. 1937, ApJ, 86, 217, doi: 10.1086/143864
—. 1951, PASP, 63, 61, doi: 10.1086/126318