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

# Intracluster light is a biased tracer of the dark matter distribution in clusters

J Butler, G Martin, N.A Hatch, F Pearce, S Brough, Y Dubois

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

## To cite this version:

J Butler, G Martin, N.A Hatch, F Pearce, S Brough, et al.. Intracluster light is a biased tracer of the dark matter distribution in clusters. Monthly Notices of the Royal Astronomical Society, 2025, 539 (3), pp.2279-2291. 10.1093/mnras/staf615. hal-05057646

HAL Id: hal-05057646

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

Submitted on 27 May 2025

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_3.png" /></figure>
<figure><img src="image_4.png" /></figure>

# Intracluster light is a biased tracer of the dark matter distribution in clusters

J. Butler<sup>1,★</sup> G. Martin<sup>1</sup> N. A. Hatch<sup>1</sup> F. Pearce,<sup>1</sup> S. Brough<sup>1,2</sup> and Y. Dubois<sup>3</sup>

Accepted 2025 April 7. Received 2025 April 4; in original form 2025 February 20

## ABSTRACT

The diffuse stellar component of galaxy clusters known as intracluster light (ICL) has been proposed as an observable tracer of the cluster's dark matter (DM) halo. Assessing its reliability as a DM tracer requires understanding how the intracluster stars are energetically linked to the underlying DM distribution, which we investigate at $z \approx 0$ in 12 galaxy clusters with $M_{178} = 1.18 - 3.71 \times 10^{14} M_{\odot}$ from the HORIZON-AGN simulation. We quantify the orbital energies of these components by their mean specific energies $\langle\epsilon\rangle$, and find that this quantity is $\approx 25$ per cent lower for the intracluster stars than the DM, whilst the energetics of the satellite galaxies (a standard DM tracer) are only marginally ($\approx 5$ per cent) higher than the DM. Importantly, the lower $\langle\epsilon\rangle$ of the intracluster stars compared to the DM is robust against the precise separation between the brightest cluster galaxy and the ICL. The specific energy distribution of ICL stars is concentrated towards lower energies and poorly samples the higher energies, where much of the DM resides. Consequently, the intracluster stars have velocity distributions with lower typical speeds and a more centrally concentrated density profile than the DM. We also find that intracluster stars have more radially biased orbits than the DM, indicating these components have distinct orbital distributions. This study demonstrates that although the morphology of the ICL may match the DM halo, the ICL is a biased tracer of DM, and these biases must be understood in order to infer properties of the DM from the ICL.

**Key words:** methods: numerical – galaxies: clusters: general – galaxies: interactions – galaxies: kinematics and dynamics.

## 1 INTRODUCTION

Galaxy clusters are the most massive collapsed structures in the Universe, with their mass distribution dominated by their dark matter (DM) halo. These haloes initially form via the gravitational collapse of the highest peaks of primordial density fluctuations (Zel'dovich 1970), and then grow hierarchically via accretion of less massive DM haloes and the direct accretion of DM from less dense regions. Independent measures of cluster masses via gravitational lensing (e.g. Grossman & Narayan 1989), observations of X-rays emitted by the hot intracluster medium (e.g. Eyles et al. 1991), the Sunyaev-Zel'dovich effect (Sunyaev & Zel'dovich 1972) and galaxy velocity dispersions (Zwicky 1933) imply that DM must account for 80–90 per cent of the mass of these structures (e.g. Gonzalez, Zaritsky & Zabludoff 2007; Gonzalez et al. 2013).

Beyond their integrated properties, the shape of DM haloes offers a probe of the nature of DM. Different models of DM (see Feng 2010 for a review) such as self-interacting DM (e.g. Spergel & Steinhardt 2000; Elbert et al. 2015; Tulin & Yu 2018) and warm DM (e.g. Lovell et al. 2012, 2014) provide possible non-baryonic solutions to small-scale challenges to standard cold DM, such as the core-cusp problem (Moore 1994; Blok 2010). Accurate measurements of the DM mass distribution are therefore needed to constrain these solutions. Since

DM is non-luminous by definition, measurements of the structure, shape, and morphology of cluster haloes must utilize visible tracers of the gravitational potential it generates. One approach is measuring the gravitational lensing of background galaxies, which in turn can be used to reconstruct the mass distribution of the galaxy cluster (see for reviews: Kneib & Natarajan 2011; Hoekstra et al. 2013). Alternatively, visible (baryonic) cluster components can be employed as tracers of the gravitational potential; typical tracers include the dynamics of satellite galaxies (e.g. Gifford, Miller & Kern 2013) and the hot intracluster medium (e.g. Borgani & Guzzo 2001; Ettori et al. 2013), the latter of which can be detected in X-rays emitted via thermal bremsstrahlung. However, all of these methods have limitations. Gravitational lensing maps can be limited by systematics such as foreground and background contamination, misidentified multiple objects or a need for stacking to enhance weak signals. Dynamical mass measurements are limited by the number of satellite galaxies in an individual cluster, resulting in a noisy signal or else requiring multiple clusters to be stacked. Finally, the intracluster medium is collisional, meaning its distribution can differ greatly from that of the DM halo (Smith et al. 2016a), especially in unrelaxed systems (e.g. the Bullet Cluster; Clowe, Gonzalez & Markevitch 2004; Markevitch et al. 2004).

In recent years, the intracluster light (ICL) has been postulated as another possible tracer of the DM halo in galaxy clusters. This light is emitted from a diffuse stellar component consisting of stars gravitationally bound to the cluster but not to any particular galaxy

<footnote n="1"><sup>1</sup>School of Physics and Astronomy, University of Nottingham, University Park, Nottingham NG7 2RD, UK</footnote>

<footnote n="2"><sup>2</sup>School of Physics, University of New South Wales, NSW 2052, Australia</footnote>

<footnote n="3"><sup>3</sup>Institut d'Astrophysique de Paris, CNRS, 98 bis blvd Arago, F-75014 Paris, France</footnote>

<footnote>* E-mail: joseph.butler@nottingham.ac.uk</footnote>
(see Montes 2022 for a review). While the existence of ICL has been
known for a long time (Zwicky 1951), its low-surface-brightness
nature (e.g. $\mu_r(3\sigma, 10'' \times 10'') > 26$ mag arcsec$^{-2}$, Montes et al.
2021; Martínez-Lombilla et al. 2023; Brough et al. 2024) has
prevented large statistical studies of the resolved properties of the
ICL. This has begun to change in the past decade, with the Dark
Energy Survey (DES; The Dark Energy Survey Collaboration 2005)
e.g. Golden-Marx et al. (2023, 2025) and the Hyper-Suprime Cam
Subaru Strategic Program (HSC-SSP; Aihara et al. 2018, 2022).
The advent of the Vera C. Rubin Observatory's Legacy Survey of
Space and Time (LSST;IVEzić et al. 2019) and the European Space
Agency's *Euclid* Wide Survey (EWS; Scaramella et al. 2022) will
deliver systematic deep and wide imaging of the ICL across a range
of redshifts, cluster masses, and dynamical states.

The primary formation mechanisms of the ICL are mergers
between satellite galaxies and the brightest cluster galaxy (BCG;
Willman et al. 2004; Conroy, Wechsler & Kravtsov 2007; Murante
et al. 2007), and the tidal stripping of satellite galaxies (Rudick
et al. 2009; Contini et al. 2014; Brown et al. 2024; Golden-Marx
et al. 2025). For very massive clusters ($M_{178} > 10^{15}$ M$_\odot$), pre-
processing is thought to be an important secondary mechanism,
where the intracluster or intragroup stars of infalling haloes are
accreted (Mihos 2004; Rudick, Mihos & McBride 2006; Contini
et al. 2014; Cañas et al. 2020). The build-up of the ICL is therefore
linked to the hierarchical assembly of galaxy clusters (Golden-Marx
et al. 2025). Since much of the DM accreted onto clusters is also
accreted from infalling galaxy and group haloes, the orbital properties
of the intracluster stars should be related to those of the DM. At
the stellar densities that produce ICL, both intracluster stars and
DM are collisionless particles, i.e. they do not share energy via
collisions between themselves or other cluster components, and thus
their subsequent dynamics are entirely determined by the evolving
potential of the cluster. Since DM is the dominant matter component
(with the exception of within the central regions of galaxies), the
gravitational potential almost entirely reflects the spatial distribution
of DM. Therefore, the ICL (emitted by intracluster stars governed
by this potential) has promise as a luminous tracer of the DM halo.

Observational studies have found a strong connection between
ICL and the host DM halo. Kluge et al. (2021) reported positive
correlations between BCG + ICL brightness and cluster properties
such as mass, radius, and velocity dispersion, indicating that both the
BCG and the ICL correlate with the host cluster properties, which
provides further motivation to investigate the ICL as a possible
DM tracer. Furthermore, the morphologies of the ICL and DM in
clusters has been compared by Montes & Trujillo (2019). They
showed that the shape of the ICL contours are well-matched to the
total mass contours, in the central regions of the Hubble Frontier
Field clusters, better than that of X-ray contours. This finding
has also been replicated in hydrodynamical simulations (Alonso
Asensio et al. 2020; Yoo et al. 2024). Although the morphologies
of the ICL and DM appear to match, many simulations have shown
that the radial density profiles do not. The ICL tends to be more
centrally concentrated, such that its radial density profile is much
steeper than the DM (Alonso Asensio et al. 2020; Contreras-Santos
et al. 2024).

The state of the cluster is described by its 6D phase space, but can
be more conveniently described in terms of orbits by transforming
to orbital energy and angular momentum. If the distributions of the
orbital energies and orbital angular momenta of the ICL stars match
the DM well, then the ICL could be used as an unbiased tracer of
the DM. However, if these distributions are significantly different,
then these differences must be quantified and their physical origin

understood in order to safely infer properties of the DM from the
ICL.

In this paper, we investigate the orbital energetics and orbital
angular momenta – as inferred by the orbital anisotropy – of the
intracluster stars and DM, using 12 clusters with total masses 1.18 –
3.71 × 10¹⁴ M$_\odot$ at $z = 0.06$ from the HORIZON-AGN cosmological
hydrodynamic simulation.

This paper is organized as follows. In Section 2, we outline the
details of the simulation and selection of the cluster sample. In
Section 3, we define the BCG, ICL, and satellite galaxies, describe the
calculation of liberation time for the intracluster stars, and define the
specific energies and orbital anisotropies of the cluster components.
In Section 4, we present our results, showing that the intracluster
stars generally have lower specific energies and a more radially
biased orbital anisotropy than both the satellite galaxies and DM.
In Section 5, we discuss potential reasons why the orbital properties
of the intracluster stars are different to the satellite galaxies and
the DM, and the observable implications of these differences. We
conclude in Section 6 by summarizing our main results.

## 2 SIMULATION

In this work, we utilize the cosmological hydrodynamic simulation, HORIZON-AGN (Dubois et al. 2014). In brief, a standard ΛCDM cosmology is adopted with the following parameters: $\Omega_m = 0.272$, $\Omega_\Lambda = 0.728$, $\Omega_{baryon} = 0.045$, $\sigma_8 = 0.81$, $H_0 = 70.4$ km s$^{-1}$ Mpc$^{-1}$, and $n_s = 0.967$, which are compatible with the WMAP7 cosmology (Komatsu et al. 2011), and also with Planck Collaboration (2014) within a 10 per cent relative variation. The simulation box with side length $100\,h^{-1}$ comoving Mpc contains $1024^3$ DM particles with a mass resolution of $M_{\text{DM, res}} = 8 \times 10^7$ M$_\odot$, and an initial gas resolution of $M_{\text{gas, res}} = 1 \times 10^7$ M$_\odot$. The adaptive mesh refinement code RAMSES (Teyssier 2002) is used as the N-body and hydrodynamical solver, which has up to 7 levels of refinement according to a quasi-Lagrangian criterion down to a spatial resolution of $\Delta x = 1$ kpc. An extra level of refinement is added every time the expansion scale factor $a$ doubles in order to keep the minimum cell size approximately constant in proper units. Gas heating from a uniform UV background takes place from $z=10$ (following Haardt & Madau 1996) and gas cooling occurs via H, He, and metals (following Sutherland & Dopita 1993). Details of the implementation of star formation, stellar feedback, and active galactic nuclei (AGN) feedback can be found in Kaviraj et al. (2017).

Cosmological simulations must make a compromise between box
size (thus sample size) and mass resolution. Given the significant im-
pact of resolution on tidal stripping (Martin et al. 2024), we prioritize
resolution in our choice of simulation to study the resolved properties
of the ICL. HORIZON-AGN has high-stellar mass ($M_{\text{star, res}} \approx 2 \times 10^6$ M$_\odot)$ and spatial resolution, comparable to other cosmological
simulations (e.g. TNG100, Pillepich et al. 2018) and cluster zoom-in
simulations (e.g. Hydrangea, Bahé et al. 2017) with similar sample
sizes in our chosen cluster mass regime. Additionally, many studies
have used HORIZON-AGN to investigate the ICL (Cañas et al. 2020;
Brough et al. 2024; Brown et al. 2024; Kimmig et al. 2025).

Apart from choosing black hole (BH) feedback parameters to
match observed BH mass relations at $z=0$, HORIZON-AGN is not
otherwise calibrated to the local Universe. Kaviraj et al. (2017) find
that the aggregate star formation histories of the model galaxies
in HORIZON-AGN broadly reproduce that of galaxies in the real
Universe, such that important bulk properties (luminosity functions,
stellar mass functions, the star formation main sequence, the cos-
mic star formation history, and rest-frame ultraviolet-optical-near
infrared colours) are in good agreement with observational data in the redshift range $0 < z < 6$.

The simulation begins at $z \approx 38$, with the initial conditions generated using MPGRAFIC (Prunet et al. 2008). The particle data is saved with a time resolution of $\approx 20$ Myr, resulting in 782 fine snapshots up to $z = 0.018$. In each snapshot, DM haloes are identified using the ADAPTAHOP structure finder (Aubert, Pichon & Colombi 2004) applied to the DM particles. Specifically, the updated version (Tweed et al. 2009) is used which ensures substructures always have lower mass than their host (sub)haloes. First, the local density associated with each particle is computed using the 20 nearest neighbours. Local maxima are identified by walking through the density field. Particles with a density above the threshold of 80 times the total matter density $\rho_M$ of the entire simulation box are then linked with their closest local maxima. Groups of particles containing at least 100 particles are then defined as (sub)structures, resulting in a minimum detectable DM (sub)halo mass of $8 \times 10^9 M_\odot$. Galaxies are identified in a similar way independently from the DM haloes, with the original ADAPTAHOP structure finder applied to the star particles. However, the density threshold is instead $178 \rho_M$, and the minimum number of particles in a (sub)structure is 50. This corresponds to a minimum detectable galaxy mass of $\sim 10^8 M_\odot$. Galaxies are then associated with DM haloes by linking the most massive galaxy within $0.1 R_{178}$ of each DM halo, where $R_{178}$ is the radius of a sphere around the halo centre with average DM density $\bar{\rho}_{DM} = 178 \rho_M$.

We generate merger trees from these fine snapshots using the TREEMAKER algorithm (originally developed by Hatton et al. 2003). This scheme identifies descendants of a structure by the maximum mass overlap between the structure in a given snapshot and structures in the subsequent snapshot, and similarly identifies main progenitors of a structure by the maximum mass overlap between the structure in a given snapshot and structures in the previous snapshot. The updated version (Tweed et al. 2009) expands this scheme to include substructures in the merger trees. From these 782 fine snapshots, we select 12 at intervals of $\sim 1$ Gyr from $z = 2.871$ to $z = 0.056$; for the remainder of this paper, the term ‘snapshots’ will refer to these 12 coarse snapshots. We chose the 14 most massive DM haloes from this final snapshot ($z = 0.06$) as cluster haloes, with a $M_{178}$ range of $1.18 - 8.30 \times 10^{14} M_\odot$, where $M_{178}$ is the total (inclusive) mass within $R_{178}$. The two most massive clusters are undergoing major mergers, and thus their DM and ICL distribution have multiple spatially distinct peaks. Much of the later analysis will involve evolution with BCG-centric radial distance and assumes both a level of spherical symmetry and that the system is approximately virialized – these assumptions are clearly broken for these two clusters, and thus they are removed from the sample.

The remaining 12 clusters, with a $M_{178}$ range of $1.18 - 3.71 \times 10^{14} M_\odot$, are chosen as our cluster sample. The properties of these clusters are displayed in Table A1. The galaxy considered to be the BCG of each cluster is the most massive galaxy within $0.1 R_{178}$ of the halo centre in the final snapshot, and the main progenitor of this galaxy is considered to be the BCG in prior snapshots. For each cluster we extract cubes with side lengths of 8 Mpc centred on the BCG for all 12 snapshots.

# 3 METHOD

## 3.1 Defining the BCG, ICL, and satellite galaxies

We investigate the orbital energies and orbital angular momenta of the intracluster stars and dark matter. There is no straightforward way to separate the ICL from the BCG, both in simulations and

observations (Cui et al. 2014; Pillepich et al. 2014; Jiménez-Teja & Dupke 2016; Cañas et al. 2019; Montes 2022; Canepa et al. 2025; see Brough et al. 2024 for a recent comparison). Thus we begin by extracting the combined BCG and ICL component, along with the main DM halo for each cluster.

To achieve this, substructures (galaxies in the stellar component and subhaloes in the DM component) must be isolated. While ADAPTAHOP does identify substructures, it does not use an unbinding procedure to determine which particles are bound to the substructures, and therefore may attribute some loosely bound star particles on the outskirts of galaxies to the BCG + ICL component. Moreover, even structure finders that incorporate unbinding procedures or leverage full 6D phase-space information can yield differing results (Knebe et al. 2013).

Since the orbits of those star particles just outside of the ADAPTAHOP structure boundary may still follow the local potential well generated by the associated satellite galaxy instead of the global cluster potential well, which reflects the global DM distribution, we impose a more stringent cut on the star particles that make up the BCG + ICL. We do this by imposing an additional spherical cut around the centre of each satellite of radius $4 r_{sat}$, and attribute all star particles in this region to the satellite. For $M_{stellar} < 10^{11} M_\odot$, $r_{sat}$ is a constant $r_{const}$, set as the 84th percentile of the satellite effective radii, $r_{eff}$, in this $M_{stellar}$ range within all clusters in the given snapshot, where $r_{eff}$ is the half-stellar mass radius. For $M_{stellar} \ge 10^{11} M_\odot$, $r_{sat} = r_{eff}$ if $r_{eff} > r_{const}$, otherwise $r_{sat} = r_{const}$.

These choices are informed by the following factors. First, there are orphan galaxies within the simulation. This justifies using a quantity related to the stellar $r_{eff}$ rather than the subhalo $R_{178}$. Secondly, we find that the $r_{eff}$ distribution at a given redshift is approximately constant<fnref n="1" /> with stellar mass up to $M_{stellar} \approx 10^{11} M_\odot$, but with some significant outliers. Thirdly, the substructure mass recovery fraction of most structure finders decreases towards the centre of clusters (Muldrew, Pearce & Power 2011), so ADAPTAHOP may underestimate $r_{eff}$ for satellites near the cluster centre (although this is less of an issue for stars than DM). These two prior reasons, along with the aim to minimize contamination of the ICL from galactic outskirts, justify using a constant value $r_{const}$ that is related to the upper end of the $r_{eff}$ distribution. Finally, although the factor of 4 in $4 r_{sat}$ is arbitrary, Ahvazi et al. (2024a) also impose a similar cleaning procedure on the satellites galaxies using a factor of 4 times $r_{eff}$. We test factors of 6, 8, and 10 and find our results involving the ICL are robust to this choice. Further details of these tests are given in Section 4.1.

The BCG + ICL for each cluster is defined as any star particle within $R_{178}$ from the centre of the BCG that is not attributed to another galaxy. The main DM halo is defined in a similar way, with all DM particles attributed to subhaloes or beyond $R_{178}$ removed, but centred on the BCG rather than the centre of the DM halo itself. An example can be seen in the projected surface density map of the cluster Hrz001 at $z \approx 0$ (Fig. 1), where the middle panels show the BCG + ICL (top) and main DM halo (bottom), and the rightmost panels show the removed galaxies (top) and subhaloes (bottom).

To separate the intracluster stars from those that belong to the BCG at $z \approx 0$, a spherical aperture of radius 100 kpc is placed around the centre of the BCG, and all star particles outside this radius are considered intracluster stars. We define the spherical shell with inner radius 100 kpc and outer radius $R_{178}$ as the ‘ICL region’. This aperture

<footnote n="1">¹This is primarily a numerical effect due to gravitational force softening at the spatial resolution limit of the simulation.</footnote>
<figure><img src="image_5.png" /><figcaption>**Figure 1.** Projected surface density map of the stars (top, green) and DM (bottom, purple) along two arbitrary spatial axes centred on the BCG for Hrz001 at $z \approx 0$. In all panels, darker colours represent higher projected surface densities on a logarithmic scale. The panels show: all star particles in a cube of volume ($R_{178}$)$^3$ (top left); the BCG + ICL (top middle); the satellite galaxies (top right); all DM particles in a cube of volume ($R_{178}$)$^3$ (bottom left); the main halo (bottom middle); and the subhaloes (bottom right). The inner and outer black circles indicate 100 kpc and $R_{178}$ around the centre, respectively.</figcaption></figure>

is a relatively conservative radial cut for the ICL – for example, Brough et al. (2024) use 30, 50, and 100 kpc cuts between the BCG and ICL. The transition radius, where the density profile of the ICL becomes dominant over the BCG, is typically identified to be between 30 and 100 kpc (Zhang et al. 2019; Montes et al. 2021; Kluge et al. 2024; Zhang et al. 2024; Golden-Marx et al. 2025). Proctor et al. (2023) identify a transition radius for $10^{14} M_{\odot}$ haloes of $\sim 100$ kpc, matching well with our BCG-ICL cut.

We discuss how different BCG-ICL radial separations affect our results in Section 4.1.

## 3.2 Liberation time of the intracluster stars

This work will not be concerned with the progenitor objects of the ICL (see Brown et al. 2024 for this analysis on a similar sample of HORIZON-AGN clusters); however, the time at which a star particle becomes part of the ICL, the so-called liberation time, is relevant. Since the ICL is assembled from many individual stripping and merging events, the orbital energy and angular momentum of the progenitor galaxies and the cluster itself at the time of liberation are imprinted on to the phase space distribution of each of these ICL subpopulations. Although these individual ICL subpopulations will evolve with time, the particles are collisionless so only phase mixing will occur. Once a star particle is removed from a galaxy, its orbital energy should only be impacted by the large-scale potential, so trends in the $z \approx 0$ ICL properties with liberation time can reveal

how the galaxy or cluster properties at the time of liberation impact the present ICL properties.

The $z \approx 0$ intracluster stellar particles are assigned a liberation time based on the following criteria:

(i) The particle's first appearance in the cutout cube must be in a satellite.

(ii) The particle must undergo a transition from a satellite to the BCG + ICL (a 'liberation event').

(iii) If a particle undergoes multiple liberation events, the first instance is selected as the liberation time.

In extreme scenarios the main progenitor branch of the merger trees may break (Srisawat et al. 2013), which we find to occur for clusters Hrz078, Hrz137, and Hrz157 in our sample, so we remove them for the liberation time analysis.

For this part of the method, the definition of a satellite in every snapshot is extended to encompass a sphere of radius $8 r_{\text{sat}}$. A factor of 8 is used here instead of the factor of 4 used to define galaxies in Section 3.1, to remove spurious liberation events due to particles with highly radial orbits continually passing between the inside and outside of the galaxy border. It is emphasized that the term 'liberation' is used here not to imply a binding energy criterion has been applied, but instead as a neutral term that describes particles that have been removed from their satellites without attempting to identify the process by which this happened (e.g. tidal stripping or satellite-BCG mergers). A particle may therefore appear to undergo
multiple liberation events if it overlaps spatially with another satellite in a subsequent snapshot, and so the first instance is selected as the liberation time.

The fraction of $z \approx 0$ ICL particles assigned a liberation time averaged between the 9 clusters is $f_{\text{liberation}} = 0.83$ with an uncertainty between clusters of $\pm 0.04$. Not all ICL particles are assigned a liberation time because the identification requires a star particle to be a satellite member in one snapshot and a BCG + ICL member in the subsequent snapshot, and so naturally does not include stars of infalling groups that have been liberated from a satellite before it was first seen inside the cutout box. This material would typically be referred to as pre-processed ICL which makes up the intragroup light of the infalling group. Note that pre-processed stars always account for a minority of the ICL in our sample (see Brown et al. 2024). There may also be contributions from star particles that were formed in a satellite but liberated before the next snapshot, and star particles formed in the BCG that later became part of the ICL, however, these channels account for less than 5 per cent (Brown et al. 2024) of the ICL at $z \approx 0$. Expanding the satellite definition from $4 r_{\text{sat}}$ to $8 r_{\text{sat}}$ will also reduce $f_{\text{liberation}}$. However, we prioritize purity over completeness, so a conservative identification resulting in $f_{\text{liberation}} < 1$ is not a concern.

### 3.3 Defining the orbital properties

The orbits of particles can be described in terms of kinetic and potential energies. These quantities vary throughout a particle's orbit, but their sum (the orbital energy) is conserved. A typical DM particle will have a higher orbital energy than a typical star particle simply due to its mass ($M_{\text{DM}}/M_{\text{star}} \approx 40$), so to allow comparisons we divide the energies by particle mass, and thus use specific kinetic energy ($\epsilon_k$) and specific potential energy ($\epsilon_p$). The specific total energy ($\epsilon$) is defined as

$$ \epsilon = \epsilon_k + \epsilon_p = \frac{v^2}{2} + (\phi - \phi_{\text{min}}), \qquad (1) $$

where $v$ is the particle speed, $\phi$ is the potential at the position of the particle, and $\phi_{\text{min}}$ is the potential at the position of the most bound particle in the cluster. The potential is computed following the Barnes & Hut (1986) method.<fnref n="2" /> The normalization to $\phi_{\text{min}}$ means the specific orbital energy will always be positive, thus the usual implication that a positive total energy corresponds to an unbound particle does not apply. In a static potential, $\epsilon$ would be constant throughout the orbit for a given particle.

We also compute $\epsilon_k$ and $\epsilon_p$ for the satellite galaxies, treating them as point masses. The radial position, $r$, and speed, $v$, of each galaxy is calculated using the mean position and velocity of the star particles within $r_{\text{sat}}$ of the galaxy centre. Calculating $\epsilon_k$ from $v$ is straightforward, but for calculating $\epsilon_p$, we must approximate the potential of the cluster at the position of the galaxy's centre, since the presence of the galaxy itself will lower the potential at that point. We do this by taking the mean of $\epsilon_p$ of the DM particles on a thin shell at $r$, excluding particles that are within $8 r_{\text{sat}}$ of any satellite.

The orbital properties of the intracluster stars and DM are not fully defined by their orbital energy; we also need to examine their angular momenta, which we explore through the orbital anisotropy, given by

$$ \beta(r) = 1 - \frac{\sigma_{\theta}(r)^2 + \sigma_{\phi}(r)^2}{2 \sigma_{r}(r)^2}, \qquad (2) $$

where $\sigma_{\theta}$, $\sigma_{\phi}$, and $\sigma_r$ are the velocity dispersions of the sample of particles in spherical coordinates (Binney & Tremaine 2008). $\beta(r) = 0$ corresponds to a population of isotropic orbits, $\beta(r) = 1$ corresponds to a population on purely radial orbits (with $\beta > 0$ generally said to be a orbits with a radial bias), whereas $\beta(r) \to -\infty$ corresponds to circular orbits (with $\beta < 0$ describing orbits with a tangential bias).

## 4 RESULTS

### 4.1 Energetics

We begin by examining the orbital energetics of a single cluster, Hrz001, at $z \approx 0$. In Fig. 2, we show the energetics of the stellar particles of the BCG and ICL components (left panel) and the main halo DM (right panel), with the galaxies overlaid in magenta. We define the phase space as the specific potential energy versus the specific kinetic energy. The galaxies preferentially sample almost the same regions of phase space as the DM. Both populations preferentially occupy the higher potential energy and lower kinetic energy regions of the phase space as any non-circular orbit will spend more time near its apocentre, where potential energy is at maximum, than its pericentre.

Conversely, the distribution of the BCG+ICL stellar particles is fundamentally different to the DM in this phase space: BCG + ICL particles have, on average, even higher potential energies relative to their kinetic energy than the DM particles, indicating a larger radial bias. The majority of the BCG + ICL particles also occupy overall lower specific energies than the DM particles. To demonstrate this mismatch, the 90th percentile of the BCG + ICL specific orbital energy distribution is added to both panels as a diagonal, black-dotted line. The significant fraction of DM particles above this line (48 per cent) clearly demonstrates that the ICL poorly samples the full phase space occupied by the DM energy distribution. We note that this mismatch in energetics is not just dominated by the stars associated with the BCG: the vertical, dashed line indicates the approximate boundary between the BCG and ICL, which we define to be at 100 kpc, and there is a large difference in the phase space distribution between the ICL and the DM to the right of that line (in the ICL region).

Clear structure is visible in both the BCG + ICL and DM distributions, much of which is aligned along isoenergy surfaces (diagonal lines in Fig. 2). Some of these will be identifiable coherent structures in 6D phase space such as stellar streams (Rudick et al. 2009), whereas others will have undergone significant phase mixing. Further investigation of these structures is however beyond the scope of this paper.

We plot the specific total energy distribution of the BCG stars, ICL stars, DM, and galaxies from cluster Hrz001 in Fig. 3 to show that the difference in the energetics of the DM and ICL stars is not a simple translation in energy states. We focus on the energetics of the star particles belonging to the ICL and compare them to the DM and galaxies that are located in the ICL region (i.e. at radii greater than 100 kpc from the cluster centre, as defined in Section 3.1).

The specific energies of the galaxy population are in reasonable agreement with the DM, except for a deficiency of galaxies at the lowest specific energies ($< 1.5 \times 10^6 \text{ km}^2 \text{ s}^{-2}$). The sharp decline in the distributions of both DM and galaxy components beyond $\sim 2.1 \times 10^6 \text{ km}^2 \text{ s}^{-2}$ results from truncating the cluster boundary at $R_{178}$. However, the ICL stars display a very different distribution to the DM, with most stars having low specific energies ($< 1.5 \times$

<footnote n="2">²Our implementation can be found at: https://github.com/garrethmartin/bh_ potential.</footnote>
<figure><img src="image_6.png" /><figcaption>**Figure 2.** Mass-weighted 2D histogram of specific kinetic energy ($\epsilon_k$) against specific potential energy ($\epsilon_p$) for the BCG + ICL (left, green) and main halo DM (right, purple) for an individual cluster, Hrz001. Darker colours represent higher densities on a logarithmic scale. The individual galaxies are plotted as magenta circles, the dashed black line indicates the mean $\epsilon_p$ at the boundary between the BCG and ICL, and the dotted black line indicates the 90th percentile of the specific (total) energy distribution of the BCG + ICL.</figcaption></figure>

<figure><img src="image_7.png" /><figcaption>**Figure 3.** Specific energy ($\epsilon$) distribution for DM in the ICL region (purple), BCG (orange), ICL (green), and galaxies in the ICL region (magenta) for an individual cluster, Hrz001.</figcaption></figure>

$10^6 \text{ km}^2 \text{ s}^{-2}$). The sharp decline at $\sim 0.7 \times 10^6 \text{ km}^2 \text{ s}^{-2}$ is due to the separation between the BCG and ICL that occurs at 100 kpc, which results in the ICL stars having a minimum specific total energy of

the $\epsilon_p$ at 100 kpc. We emphasize that the precise shapes of these distributions therefore depend on these cuts at 100 kpc and $R_{178}$, and also varies between clusters. For example, the apparent bimodality seen in the distribution of the ICL is not a consistent feature. However, the broad result that the galaxy distribution is similar to the DM distribution, and the ICL distribution mostly populates lower energies, holds.

We next explore the orbital energies in all 12 of the simulated clusters at $z \approx 0$. The total energy for each cluster varies due to the different cluster masses, therefore, to normalize the orbital information of all the clusters we first compute the specific energy ratio $\langle\epsilon_{\text{component/DM}}\rangle$, defined as the average specific energy $\langle\epsilon\rangle$ for each component (BCG + ICL and galaxies) divided by that of the DM for each cluster. Finally, we calculate the average of these ratios for all 12 clusters. In this calculation a mass-weighted average is taken for stellar and DM components, while galaxies are treated as point tracers, with no mass weighting. We calculate the specific energy ratios for all cluster components within $R_{178}$, but also calculate them while excluding the inner 100, 200, and 300 kpc to explore the impact of including the star particles that are part of the BCG. These specific energy ratios are shown in Fig. 4.

Consistent with Figs 2 and 3, the energetics of the galaxies and the DM halo are similar, with $\langle\epsilon_{\text{galaxies/DM}}\rangle = 1.05 \pm 0.01$ (as seen in Fig. 4). The quoted uncertainties of specific energy ratios are the 1$\sigma$ standard error on the mean, but we also note that the dispersion between clusters (the shaded regions in Fig. 4) is consistently $\le 0.05$. The slightly higher energies of the galaxies compared to the DM halo can be explained by galaxies with low orbital energy preferentially merging with the BCG and thus being effectively removed from the galaxy sample (galactic cannibalism; Ostriker & Tremaine 1975).
<figure><img src="image_8.png" /><figcaption>**Figure 4.** Specific energy ratio of galaxies ($\langle\epsilon_{galaxies/DM}\rangle$, magenta) and ICL ($\langle\epsilon_{ICL/DM}\rangle$, green) at different BCG–ICL separations, such that only the populations beyond this separation are considered for all components (galaxies, ICL, and DM). The points represent the means, the error bars are the 1$\sigma$ standard error on the mean, the shaded regions are the 1$\sigma$ dispersion between clusters, and the purple dashed line indicates a specific energy ratio of 1. The extreme BCG–ICL separation of 300 kpc still resulting in $\langle\epsilon_{ICL/DM}\rangle < 1$ demonstrates that the BCG is not the primary driver of this difference in average specific energies.</figcaption></figure>

Indeed, if only the DM and galaxies beyond 100 kpc are considered, then this ratio drops to $\langle\epsilon_{galaxies/DM}\rangle = 1.01 \pm 0.01$, and remains consistent with 1 when this cut is increased to 200 kpc and 300 kpc. This demonstrates that the satellite galaxies sample the DM halo phase space well and are thus an effective, mostly unbiased tracer of the underlying DM.

If the intracluster stars were unbiased tracer particles of the DM halo, we would expect $\langle\epsilon_{ICL/DM}\rangle$ to be approximately unity. However, taking the ICL region to be that beyond 100 kpc of the BCG, we find that $\langle\epsilon_{ICL/DM}\rangle = 0.75 \pm 0.01$, i.e. the mean specific energy of the ICL is only $\approx 75$ per cent of the DM.

The stars associated with the BCG are expected to have lower energies than those associated with the ICL due to their location near the bottom of the potential well. Therefore, when we include all the BCG stars with the ICL component, we measure $\langle\epsilon_{BCG+ICL/DM}\rangle = 0.56 \pm 0.01$. As previously mentioned, Proctor et al. (2023) identified the BCG-ICL transition region as $\sim 100$ kpc for $10^{14} M_{\odot}$ haloes, so excluding the central 100 kpc region may still leave a contribution of star particles that belong to the BCG outskirts in this definition of the ICL. In Fig. 4, we show that $\langle\epsilon_{ICL/DM}\rangle$ increases as we progressively exclude more of the core of the cluster. Importantly, however, even in the extreme case of excluding out to a radius of 300 kpc (at which point any BCG contribution can safely be assumed negligible), we find that $\langle\epsilon_{ICL/DM}\rangle$ still only reaches $0.86 \pm 0.01$, meaning the contamination from the BCG is not the primary driver of our results. Overall, these show that although the separation between the BCG and ICL will affect the precise $\langle\epsilon_{ICL/DM}\rangle$, the general result that $\langle\epsilon_{ICL/DM}\rangle < 1$ holds.

<figure><img src="image_9.png" /><figcaption>**Figure 5.** Anisotropy $\beta$ against radius for dark matter (purple), ICL (green), and galaxies (magenta). The points represent the mean between clusters, and the error bars are the 1$\sigma$ standard error on the mean. The midpoints of the radial bins for the galaxies are plotted slightly offset for visual clarity, and truncated below $\approx 0.15 R_{178}$ due to low counts.</figcaption></figure>

We also test the robustness of these results to different definitions of satellite galaxies. The factor of 4 in the spherical cut of radius $4 r_{\text{sat}}$ around the centre of each galaxy is relatively arbitrary, so we test factors of 6, 8, and 10. The effect on $\langle\epsilon_{ICL/DM}\rangle$ is a negligible increase, only increasing from $0.75 \pm 0.01$ for a factor of 4 to $0.77 \pm 0.01$ for a factor of 10. However, increasing this factor results in significant decreases of BCG + ICL mass as more stars are attributed to galaxies. Thus, we choose to continue using a factor of 4.

## 4.2 Anisotropy

Now that we have established differences between the orbital energetics of the DM, ICL stars and galaxies, we turn our attention to the orbital anisotropy. Fig. 5 shows the anisotropy profiles of the DM, ICL stars and galaxies averaged between clusters.

We see that the DM anisotropy increases (becomes more radial) with clustercentric radius, increasing from $\beta \approx 0.15$ to $\beta \approx 0.35$, in agreement with previous studies (Willman et al. 2004; Wojtak et al. 2005; Ascasibar & Gottlöber 2008; Lemze et al. 2012; He et al. 2024). The anisotropy profile of the galaxies (magenta points) is truncated below $\approx 0.15 R_{178}$ due to low counts, and has a larger dispersion between clusters, such that a positive trend with radius is less certain. However, the average galaxy anisotropy profile from all 12 clusters is similar to the DM profile, in agreement with figure 10 from Willman et al. (2004). Our average galaxy anisotropy profile is in qualitative agreement with observational studies (Biviano et al. 2013; Mamon et al. 2019) that find $\beta \approx 0$ in the cluster core to $\beta \approx 0.5$ at $R_{200}$.

Other studies investigate the anisotropy profile of the subhaloes (Diemand, Moore & Stadel 2004; He et al. 2024) and find the subhalo anisotropy is typically lower than in our work. For example, He et al. (2024) find the anisotropy of subhaloes increases from $\beta \approx 0$ to $\beta \approx 0.2$ in the halo mass range $\log_{10}(M_{200}/M_{\odot}) = 13.8 - 14.3$. However, it is important to note that the anisotropy profiles of the
galaxies and subhaloes are not directly comparable since orphan galaxies exist in some simulations.

An increase in anisotropy towards more radial orbits can also be seen for the ICL stars, increasing from $\beta \approx 0.4$ to $\beta \approx 0.55$. This is in good agreement with He et al. (2024; see their fig. 6), who also find the same trend and similar values for accreted halo stars in the TNG300 simulation. Willman et al. (2004) also find that $\beta$ increases with radius for unbound stars, albeit encompassing a much larger range, increasing from $\beta \approx -0.2$ to $\beta \approx 0.8$. Importantly, we find that at all radii, the anisotropy of the ICL stars is more radially biased than that of the DM and galaxies.

# 5 DISCUSSION

Our work demonstrates that the orbital energies and anisotropy of the ICL stars inherently differ from both the galaxies and the DM. This is despite the dynamics of each of these components being governed by the same potential, and the ICL stars primarily originating from galaxies. We discuss potential reasons why the orbital properties of the ICL stars are different to the satellite galaxies and the DM, and the observable implications of these differences.

## 5.1 Why do the orbits differ?

### 5.1.1 Evolving cluster potential

We first consider the possibility that the energies of the intracluster stars could be lower than the DM because the intracluster stars' energy reflects the evolving properties of the cluster potential at the time of stripping. Galaxies that entered the cluster at earlier times, when it was less massive, would travel on lower energy orbits compared to galaxies infalling at later times, and thus the orbits of the stars liberated through tidal stripping would also have these lower energies. If these orbital energies at liberation time were approximately maintained until $z = 0$, one might expect that $\langle\epsilon_{ICL/DM}\rangle$ split up by liberation time would steadily increase towards $z = 0$, with the most recently liberated material matching the energetics of the $z = 0$ DM halo. However, in Fig. 6 we show there is no clear dependence of $\langle\epsilon_{ICL/DM}\rangle$ on liberation time until a lookback time of ~6 Gyr, after which there is a slight increase in $\langle\epsilon_{ICL/DM}\rangle$, such that more recently liberated stars tend to have higher energies. However, this trend is not consistent between different definitions of satellite galaxies and BGC-ICL separations, so we do not explore this further. More importantly, even the most recently liberated stars have $\langle\epsilon_{ICL/DM}\rangle = 0.80 \pm 0.02$, and thus do not reach the average energies of the DM. This mechanism, therefore, cannot fully explain the lower energetics of the ICL stars compared with the DM.

### 5.1.2 Differential stellar and DM stripping

An alternative possibility is that a majority of intracluster stars are stripped from galaxies on orbits with lower energies and/or higher angular momenta than the DM. Many studies have found that the ICL is not built up by uniform contributions from all galaxies (e.g. Contini et al. 2014; Contini, Yi & Kang 2019; Chun et al. 2023, 2024; Ahvazi et al. 2024b; Brown et al. 2024) and this also applies to their orbits. As the DM halo of an orbiting satellite is, on average, more weakly bound than its stars, significant amounts of DM can be tidally stripped on orbits that result in very little stellar stripping (e.g. Villalobos et al. 2012; Smith et al. 2016b; Joshi et al. 2019; Haggar et al. 2021; Martin et al. 2024). This would therefore bias

<figure><img src="image_10.png" /><figcaption>**Figure 6.** Specific energy ratio $\langle\epsilon_{ICL/DM}\rangle$ split up by liberation time. This is the $\langle\epsilon_{ICL}\rangle$ at $z \approx 0$ of intracluster stars with a given liberation time, divided by $\langle\epsilon_{DM}\rangle$ at $z \approx 0$. The points represent the mean between clusters, the vertical error bars, the $1\sigma$ standard error on the mean, and the horizontal error bars the width of the individual age bins. This plot demonstrates that even the most recently liberated stars have $\langle\epsilon_{ICL/DM}\rangle < 1$.</figcaption></figure>

the orbital energies and anisotropy of the intracluster stars to those of galaxies on lower energy and more radial orbits, and thus explain why the intracluster stars have lower specific energies (Figs 2 and 3) and appear more radially biased (Fig. 5) compared with both the DM and the galaxies.

Another possibility is that loss of orbital energy over time due to dynamical friction (Chandrasekhar 1943) will result in higher levels of stellar stripping and therefore bias the orbital energies of the intracluster stars to reflect those of more massive galaxies that have fallen into lower orbits. At the same time, galaxies experiencing high levels of dynamical friction will eventually merge with the BCG (Tormen 1997), biasing the satellite population to modestly higher orbital energies (Ostriker & Tremaine 1975; Dolag, Murante & Borgani 2010). However, this mechanism alone cannot account for the more radial orbits observed in the intracluster stars, as dynamical friction typically results in either no change or in the circularization of satellite orbits, dependent on the density profile of the host cluster (van den Bosch et al. 1999; Arena & Bertin 2007).

To summarize, our result that the intracluster stars have lower orbital energies and a more radially biased anisotropy compared to the galaxies and DM particles could be explained by biases in orbits and galaxies that progenate the bulk of the ICL stars. $\langle\epsilon_{ICL/DM}\rangle$ may also depend on how the ICL formed in a particular cluster. For example, if a large fraction of ICL was formed from pre-processed material, which is loosely bound to the infalling subhalo, then the orbital energy difference in such a cluster may be smaller compared to clusters where the ICL formed primarily by mergers with the BCG. $\langle\epsilon_{ICL/DM}\rangle$ and the radially biased $\beta$ likely depend on the specific ICL formation process of the cluster, and therefore there may be a scatter in these parameters within the cluster population.
<figure><img src="image_11.png" /><figcaption>**Figure 7.** Top panel: Density profiles of DM (purple) and ICL (green). The DM density profile is divided by an arbitrary factor of 10 to make visual comparison of the DM and ICL profile shapes clearer. Bottom panel: Density ratio $\rho_{\text{ICL}}/\rho_{\text{DM}}$ as function of radius. Grey lines are individual clusters and the black line is the power-law fit $y = a \cdot x^b$ with $a = (2.16 \pm 0.12) \times 10^{-3}$ and $b = -1.19 \pm 0.04$.</figcaption></figure>

## 5.2 Observable implications

### 5.2.1 Density profile

The nature of DM can be probed by analysing the radial mass distribution of the halo (e.g. Spergel & Steinhardt 2000), and previous studies (e.g. Alonso Asensio et al. 2020; Contreras-Santos et al. 2024) have suggested that the ICL radial profile could be used to infer the DM radial profile. However, a notable consequence of the differing energetics and orbital anisotropy is that the density profile of the ICL should differ from the DM profile. The Jeans equation tells us that there are two competing effects: the greater anisotropy of the intracluster stars means the ICL density profile should be flatter than the DM, but the lower orbital energetics of the intracluster stars means the ICL profile should be more concentrated. In Fig. 7, we show the individual density profiles of both the DM and ICL for each cluster (top panel) and their ratio (bottom panel). We see that the ratio of the density profiles between the ICL and DM can be fit by a power law with a negative exponent, such that the density of the ICL decreases faster than the DM. This demonstrates that

the lower energetics of the ICL are the dominant factor over the anisotropy in determining the radial profile. The gradient of the density ratio is similar for all clusters, but the normalization differs. We test if this normalization is related to $\langle\epsilon_{\text{ICL/DM}}\rangle$ and find no clear correlation.

A comparison of the ICL profile and the DM profile, as measured by gravitational lensing, has been performed for SMACS0723 (Diego et al. 2023) and MACS0416 (Diego et al. 2024). They find that the ICL profile is steeper than the DM profile, in qualitative agreement with our findings. Similarly, Chen et al. 2022 find the surface mass density profile of the BCG + ICL is more centrally concentrated than the DM.

It should be emphasized that the difference in the shapes (density profiles) of the ICL and DM does *not* imply that their morphologies should be different. Since the dynamics of both components are governed by the same cluster potential, it is perhaps unsurprising that both observations (Montes & Trujillo 2019) and simulations (Alonso Asensio et al. 2020) agree that the morphologies match well in the inner region, but crucially this does not guarantee that the two components sample this potential in a similar way.

Our work confirms the findings of Alonso Asensio et al. (2020) and Contreras-Santos et al. (2024) who demonstrated that the ratio of the density profiles of the intracluster stars and DM approximately follows a power-law relationship. However, the exact mechanisms driving this relationship remain unclear. It is uncertain whether this relationship holds universally across different cluster properties such as dynamical states, mass, and redshift. Furthermore, the observed differences between the density profiles of DM and ICL may be sensitive to the dominant ICL formation process as well as the properties of the galaxy populations that progenate the ICL, whose size and stellar mass distribution are themselves sensitive to the implemented galaxy evolution models. Finally, the presence of non-phase-mixed material may cause temporal variations in the density profiles. We therefore caution against using this relation to estimate the DM profile from the observed ICL profile, until the reasons for the differing energetics and orbital anisotropies are better understood.

### 5.2.2 Kinematics

Another potentially observable implication of the differing energetics of the ICL and DM is that their velocity distributions will also differ. The key information about the velocity distributions for each cluster in order of increasing $M_{178}$ are summarized in Fig. 8. For each cluster, the arrays of DM and ICL particle speeds are bootstrapped 100 times and splines fit to each. The peak of each distribution (represented by triangles for DM and stars for ICL) is then computed from the mean of the peaks of each spline, and the full width at half maximum (FWHM, represented by solid bars) from where the median curve through these splines equals half the peak height. The large number of particles in both components ensures the error on the mean speed of their spline peaks are negligible. The top-left inset illustrates this for one of the clusters.

For every cluster, the peak speed is higher for the DM than the ICL stars. This is in agreement with Contreras-Santos et al. (2024), who use a similar aperture-based method of BCG-ICL separation with a 50 kpc radius sphere, and find the velocity dispersion ratio between the ICL and DM to be ≈ 0.75 within $R_{500}$. Notably, our work also agrees with Dolag et al. (2010), who employ a kinematic separation of the BCG and ICL. This method naturally selects the higher speed stars to be associated with the ICL instead of the BCG, and yet the
<figure><img src="image_12.png" /><figcaption>**Figure 8.** Velocity distributions of the DM in the ICL region (purple) and ICL (green) for each cluster, ordered by $z = 0$ mass. The peaks of the velocity distributions are represented by triangles for the DM and stars for the ICL, and the error bars represent the FWHM of the velocity distributions. The inset indicates these quantities for Hrz001. The difference in the peaks of the velocity distributions for each cluster is given in Table A1.</figcaption></figure>

velocity distribution of the ICL still peaks at a lower speed than the DM.

Observationally, measuring the kinematics of the ICL is extremely difficult, due to both its faint nature and uncertainty in the transition radius. Integrated-light absorption spectra has been measured for many BCGs (e.g. Carter et al. 1981, 1985; Kelson et al. 2002; Bender et al. 2015; Boardman et al. 2017; Barbosa et al. 2018; Loubser et al. 2022), but typically the kinematics are limited to within 3 – 4 half-light radii. Some studies (e.g. Bender et al. 2015; Edwards et al. 2020) have found that the stellar velocity dispersion reaches that of the cluster in the BCG outskirts, which may disagree with our result that the velocity dispersion of the ICL should never reach that of the DM. However, the difference in velocity dispersions between the two components may not be pronounced, and therefore be within uncertainties. For example, Dolag et al. (2010) find an ~500 km s⁻¹ difference in the peak speeds of the DM and the diffuse stellar component, but this translates to a difference in velocity dispersion of only ~20 km s⁻¹ (see their fig. 5).

To probe the ICL at larger distances, most studies attempt to identify and measure the kinematics of bright tracers, usually intracluster planetary nebulae (PNe; Arnaboldi et al. 1998; Longobardi et al. 2013; Hartke et al. 2017; Pulsoni et al. 2018) and/or globular clusters (GCs; Schuberth et al. 2010; Strader et al. 2011; Pota et al. 2018). There have only been a handful of successful observations to date, which include the Virgo, Fornax, Hydra I, and Coma clusters – see Arnaboldi & Gerhard (2022) and references therein for full detail on these. However, different ICL tracers may have different energetics and anisotropy than the global intracluster stars. Simulations (Reina-Campos et al. 2022) and observations (Alabi et al. 2016, 2017; Hudson & Robison 2018) show that intracluster GCs are good tracers of dark matter haloes, and so may have energetics and anisotropy more similar to the DM rather than the ICL. On the other hand, Ramos-Almendares et al. (2020) simulate intracluster GCs by associating GCs with DM particles selected to match the

radial distribution of GCs in observed galaxies. They show that the GC kinematics are 25 per cent lower than the DM, which aligns with our results for the ICL.

## 6 CONCLUSIONS

In this study we have investigated the energetics of the ICL and DM at $z \approx 0$, using 12 clusters with masses $1.18 – 3.71 \times 10^{14} M_{\odot}$ from the HORIZON-AGN cosmological hydrodynamic simulation. We summarize our results below:

(i) We calculate the mean specific energies $\langle \epsilon \rangle$ of the DM halo, ICL, and satellite galaxies. We find that the ratio between that of the ICL and DM is $\langle \epsilon_{ICL/DM} \rangle = 0.75 \pm 0.01$, such that the average orbital energy of the intracluster stars is $\approx 75$ per cent of that the DM. On the other hand, galaxies (treated as point masses) have marginally higher energies than the DM with $\langle \epsilon_{galaxies/DM} \rangle = 1.05 \pm 0.01$. Whilst the specific energy distributions of the DM and galaxies appear to match well (Fig. 3), the ICL has a very different distribution, concentrated towards lower energies and poorly sampling the higher energies at which much of the DM and galaxies reside (Fig. 2).

(ii) These results are tested against different BCG–ICL separations (Fig. 4). The energetics of the galaxy population becomes consistent with the DM halo when the innermost 100 kpc of both components are excluded ($\langle \epsilon_{galaxies/DM} \rangle = 1.01 \pm 0.01$), showing that there is a deficit in the low-energy galaxy population in the cluster core, but the galaxies are otherwise a good tracer of DM. As the BCG–ICL separation is increased to an extreme of 300 kpc, where the BCG contribution can safely be assumed to be negligible, $\langle \epsilon_{ICL/DM} \rangle$ still only reaches $0.86 \pm 0.01$. This demonstrates that the lower energetics of the ICL is not simply due to BCG contribution, and thus this qualitative result is robust to the precise definition of the ICL. We also test the impact of expanding the extent of satellite galaxies, which assigns additional stars to the satellites
galaxies instead of to the ICL, and find only small quantitative
differences.

(iii) The angular momentum of the cluster components is quantified by the orbital anisotropy. The ICL is found to be more radially biased (larger $\beta$) than the DM at all clustercentric radii (Fig. 5). This, in combination with the differing energetics, highlights that the orbital configurations of the ICL and DM are inherently different.

(iv) These differences in the energetics and anisotropy of the ICL and DM manifest themselves in their velocity distributions and density profiles. We find that the ICL density profile is more centrally concentrated than that of the DM (Fig. 7), such that the ratio of their densities as a function of radius has a negative gradient. We also find that the peak speed is lower and FWHM is smaller for the ICL velocity distribution than that of the DM for all clusters (Fig. 8). Both of these results qualitatively agree with other simulation work.

While ICL represents a promising tool for exploring the properties of DM, we have shown that it is not an unbiased proxy for the underlying DM distribution. ICL stars have, on average, lower orbital energies and more radial orbits than DM, meaning their density profiles and velocity distributions differ. The relationship between ICL and DM likely also depends on other factors, such as cluster mass, dynamical state, and the interplay of different mechanisms driving ICL production. Given the limitations of traditional methods like weak lensing, incorporating ICL as a complementary tracer could significantly improve our ability to map DM, particularly with upcoming surveys such as LSST and *Euclid*, which will provide deep imaging for large numbers of clusters. A comprehensive understanding of the processes that generate ICL – and the origins of its discrepancies with DM – will therefore be essential for calibrating its use effectively.

## ACKNOWLEDGEMENTS

We thank the anonymous referee for their constructive comments which have improved this work. We also thank members of the NottICL group for helpful discussions and comments, in particular Yannick Bahé, Jesse B. Golden-Marx, Tutku Kolcu, and Harley J. Brown.

JB and NAH gratefully acknowledge support from the Leverhulme Trust through a Research Leadership Award. GM, FRP, and NAH acknowledge support from the UK STFC under grant ST/X000982/1. SB acknowledges funding support from the Australian Research Council through a Discovery Project DP190101943. This research was supported by the International Space Science Institute (ISSI) in Bern, through ISSI International Team project #23-577.

This work was granted access to the HPC resources of CINES under allocations 2013047012, 2014047012, and 2015047012 made by GENCI. This work has made use of the Infinity cluster, hosted by the Institut d'Astrophysique de Paris. We warmly thank S. Rouberol for running it smoothly.

## DATA AVAILABILITY

The raw data products of the HORIZON-AGN simulation are available upon reasonable request through the collaboration's website: https://www.horizon-simulation.org. The data products generated from this work are available upon reasonable request from the corresponding author.

## REFERENCES

Ahvazi N., Sales L. V., Navarro J. F., Benson A., Boselli A., D'Souza R., 2024a, Open J. Astrophys., 7, 111

Ahvazi N., Sales L. V., Doppel J. E., Benson A., D'Souza R., Rodriguez-Gomez V., 2024b, MNRAS, 529, 4666

Aihara H. et al., 2018, PASJ, 70, S4

Aihara H. et al., 2022, PASJ, 74, 247

Alabi A. B. et al., 2016, MNRAS, 460, 3838

Alabi A. B. et al., 2017, MNRAS, 468, 3949

Alonso Asensio I., Dalla Vecchia C., Bahé Y. M., Barnes D. J., Kay S. T., 2020, MNRAS, 494, 1859

Arena S. E., Bertin G., 2007, A&A, 463, 921

Arnaboldi M., Freeman K. C., Gerhard O., Matthias M., Kudritzki R. P., Méndez R. H., Capaccioli M., Ford H., 1998, ApJ, 507, 759

Arnaboldi M., Gerhard O. E., 2022, Front. Astron. Space Sci. 9, 403

Ascasibar Y., Gottlöber S., 2008, MNRAS, 386, 2022

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

Bahé Y. M. et al., 2017, MNRAS, 470, 4186

Barbosa C. E., Arnaboldi M., Coccato L., Gerhard O., Mendes De Oliveira C., Hilker M., Richtler T., 2018, A&A, 609, A78

Barnes J., Hut P., 1986, Nature, 324, 446

Bender R., Kormendy J., Cornell M. E., Fisher D. B., 2015, ApJ, 807, 56

Binney J., Tremaine S., 2008, Galactic Dynamics: Second Edition. Princeton Univ. Press, Princeton

Biviano A. et al., 2013, A&A, 558, A1

Blok W. J. G. d., 2010, Adv. Astron., 2010, 789293

Boardman N. F. et al., 2017, MNRAS, 471, 4005

Borgani S., Guzzo L., 2001, Nature, 409, 39

Brough S. et al., 2024, MNRAS, 528, 771

Brown H. J., Martin G., Pearce F. R., Hatch N. A., Bahé Y. M., Dubois Y., 2024, MNRAS, 534, 431

Cañas R., Elahi P. J., Welker C., Lagos C. D. P., Power C., Dubois Y., Pichon C., 2019, MNRAS, 482, 2039

Cañas R., Lagos C. D. P., Elahi P. J., Power C., Welker C., Dubois Y., Pichon C., 2020, MNRAS, 494, 4314

Canepa L., Brough S., Lanusse F., Montes M., Hatch N., 2025, ApJ, 980, 245

Carter D., Efstathiou G., Ellis R. S., Inglis I., Godwin J., 1981, MNRAS, 195, 15

Carter D., Inglis I., Ellis R. S., Efstathiou G., Godwin J. G., 1985, MNRAS, 212, 471

Chandrasekhar S., 1943, ApJ, 97, 255

Chen X., Zu Y., Shao Z., Shan H., 2022, MNRAS, 514, 2692

Chun K., Shin J., Smith R., Ko J., Yoo J., 2023, ApJ, 943, 148

Chun K., Shin J., Ko J., Smith R., Yoo J., 2024, ApJ, 969, 142

Clowe D., Gonzalez A., Markevitch M., 2004, ApJ, 604, 596

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

Contini E., De Lucia G., Villalobos A., Borgani S., 2014, MNRAS, 437, 3787

Contini E., Yi S. K., Kang X., 2019, ApJ, 871, 24

Contreras-Santos A. et al., 2024, A&A, 683, A59

Cui W. et al., 2014, MNRAS, 437, 816

Diego J. M. et al., 2023, A&A, 679, A159

Diego J. M. et al., 2024, A&A, 690, A114

Diemand J., Moore B., Stadel J., 2004, MNRAS, 352, 535

Dolag K., Murante G., Borgani S., 2010, MNRAS, 405, 1544

Dubois Y. et al., 2014, MNRAS, 444, 1453

Edwards L. O. V. et al., 2020, MNRAS, 491, 2617

Elbert O. D., Bullock J. S., Garrison-Kimmel S., Rocha M., Oñorbe J., Peter A. H. G., 2015, MNRAS, 453, 29

Ettori S., Donnarumma A., Pointecouteau E., Reiprich T. H., Giodini S., Lovisari L., Schmidt R. W., 2013, Space Sci. Rev., 177, 119

Eyles C. J., Watt M. P., Bertram D., Church M. J., Ponman T. J., Skinner G. K., Willmore A. P., 1991, ApJ, 376, 23

Feng J. L., 2010, ARA&A, 48, 495

Gifford D., Miller C., Kern N., 2013, ApJ, 773, 116

Golden-Marx J. B. et al., 2023, MNRAS, 521, 478

Golden-Marx J. B. et al., 2025, MNRAS, 538, 622
Gonzalez A. H., Zaritsky D., Zabludoff A. I., 2007, ApJ, 666, 147

Gonzalez A. H., Sivanandam S., Zabludoff A. I., Zaritsky D., 2013, ApJ, 778, 14

Grossman S. A., Narayan R., 1989, ApJ, 344, 637

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

Haggar R., Pearce F. R., Gray M. E., Knebe A., Yepes G., 2021, MNRAS, 502, 1191

Hartke J., Arnaboldi M., Longobardi A., Gerhard O., Freeman K. C., Okamura S., Nakata F., 2017, A&A, 603, A104

Hatton S., Devriendt J. E. G., Ninin S., Bouchet F. R., Guiderdoni B., Vibert D., 2003, MNRAS, 343, 75

He J. et al., 2024, ApJ, 976, 187

Hoekstra H., Bartelmann M., Dahle H., Israel H., Limousin M., Meneghetti M., 2013, Space Sci. Rev., 177, 75

Hudson M. J., Robison B., 2018, MNRAS, 477, 3869

Ivezić Ž. et al., 2019, ApJ, 873, 111

Jiménez-Teja Y., Dupke R., 2016, ApJ, 820, 49

Joshi G. D., Parker L. C., Wadsley J., Keller B. W., 2019, MNRAS, 483, 235

Kaviraj S. et al., 2017, MNRAS, 467, 4739

Kelson D. D., Zabludoff A. I., Williams K. A., Trager S. C., Mulchaey J. S., Bolte M., 2002, ApJ, 576, 720

Kimmig L. C. et al., 2025, preprint (arXiv:2053.20857)

Kluge M., Bender R., Riffeser A., Goessl C., Hopp U., Schmidt M., Ries C., 2021, ApJS, 252, 27

Kluge M. et al., 2024, preprint (arXiv:2405.13503)

Knebe A. et al., 2013, MNRAS, 435, 1618

Kneib J.-P., Natarajan P., 2011, A&AR, 19, 47

Komatsu E. et al., 2011, ApJS, 192, 18

Lemze D. et al., 2012, ApJ, 752, 141

Longobardi A., Arnaboldi M., Gerhard O., Coccato L., Okamura S., Freeman K. C., 2013, A&A, 558, A42

Loubser S. I., Lagos P., Babul A., O'Sullivan E., Jung S. L., Olivares V., Kolokythas K., 2022, MNRAS, 515, 1104

Lovell M. R. et al., 2012, MNRAS, 420, 2318

Lovell M. R., Frenk C. S., Eke V. R., Jenkins A., Gao L., Theuns T., 2014, MNRAS, 439, 300

Mamon G. A., Cava A., Biviano A., Moretti A., Poggianti B., Bettoni D., 2019, A&A, 631, A131

Markevitch M., Gonzalez A. H., Clowe D., Vikhlinin A., Forman W., Jones C., Murray S., Tucker W., 2004, ApJ, 606, 819

Martin G., Pearce F. R., Hatch N. A., Contreras-Santos A., Knebe A., Cui W., 2024, MNRAS, 535, 2375

Martínez-Lombilla C. et al. 2023, MNRAS, 518, 1195

Mihos J. C., 2004, in Mulchaey J. S., Dressler A., Oemler A. eds, Clusters of Galaxies: Probes of Cosmological Structure and Galaxy Evolution. Cambridge University Press, Cambridge, p.277

Montes M., 2022, Nat. Astron., 6, 308

Montes M., Brough S., Owers M. S., Santucci G., 2021, ApJ, 910, 45

Montes M., Trujillo I., 2019, MNRAS, 482, 2838

Moore B., 1994, Nature, 370, 629

Muldrew S. I., Pearce F. R., Power C., 2011, MNRAS, 410, 2617

Murante G., Giovalli M., Gerhard O., Arnaboldi M., Borgani S., Dolag K., 2007, MNRAS, 377, 2

Ostriker J. P., Tremaine S. D., 1975, ApJ, 202, L113

Pillepich A. et al., 2014, MNRAS, 444, 237

Pillepich A. et al., 2018, MNRAS, 475, 648

Planck Collaboration XIV, 2014, A&A, 571, A16

Pota V. et al., 2018, MNRAS, 481, 1744

Proctor K. L., Lagos C. d. P., Ludlow A. D., Robotham A. S. G., 2023, MNRAS, 527, 2624

Prunet S., Pichon C., Aubert D., Pogosyan D., Teyssier R., Gottlieber S., 2008, ApJS, 178, 179

Pulsoni C. et al., 2018, A&A, 618, A94

Ramos-Almendares F., Sales L. V., Abadi M. G., Doppel J. E., Muriel H., Peng E. W., 2020, MNRAS, 493, 5357

Reina-Campos M., Trujillo-Gomez S., Deason A. J., Kruijssen J. M. D., Pfeffer J. L., Crain R. A., Bastian N., Hughes M. E., 2022, MNRAS, 513, 3925

Rudick C. S., Mihos J. C., McBride C., 2006, ApJ, 648, 936

Rudick C. S., Christopher Mihos J., Frey L. H., McBride C. K., 2009, ApJ, 699, 1518

Scaramella R. et al., 2022, A&A, 662, A112

Schuberth Y., Richtler T., Hilker M., Dirsch B., Bassino L. P., Romanowsky A. J., Infante L., 2010, A&A, 513, A52

Smith G. P. et al., 2016a, MNRAS, 456, L74

Smith R., Choi H., Lee J., Rhee J., Sanchez-Janssen R., Yi S. K., 2016b, ApJ, 833, 109

Spergel D. N., Steinhardt P. J., 2000, Phys. Rev. Lett., 84, 3760

Srisawat C. et al., 2013, MNRAS, 436, 150

Strader J. et al., 2011, ApJS, 197, 33

Sunyaev R. A., Zeldovich Y. B., 1972, Comment. Astrophys. Space Phys., 4, 173

Sutherland R. S., Dopita M. A., 1993, ApJS, 88, 253

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

The Dark Energy Survey Collaboration, 2005, preprint (arXiv:astro-ph/0510346)

Tormen G., 1997, MNRAS, 290, 411

Tulin S., Yu H.-B., 2018, Phys. Rep., 730, 1

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

van den Bosch F. C., Lewis G. F., Lake G., Stadel J., 1999, ApJ, 515, 50

Villalobos Á., De Lucia G., Borgani S., Murante G., 2012, MNRAS, 424, 2401

Willman B., Governato F., Wadsley J., Quinn T., 2004, MNRAS, 355, 159

Wojtak R. a., Łokas E. L., Gottlöber S., Mamon G. A., 2005, MNRAS, 361, L1

Yoo J. et al., 2024, ApJ, 965, 145

Zel'dovich Y. B., 1970, A&A, 5, 84

Zhang Y. et al., 2019, ApJ, 874, 165

Zhang Y. et al., 2024, MNRAS, 531, 510

Zwicky F., 1933, Helv. Phys. Acta, 6, 110

Zwicky F., 1951, PASP, 63, 61

## APPENDIX A: TABLE OF CLUSTER PROPERTIES

We present key properties of the 12 HORIZON-AGN clusters in our sample in Table A1. All clusters are used in our $z \approx 0$ analysis, but, due to breaks in their merger tree main progenitor branch, clusters Hrz078, Hrz137 and Hrz157 are excluded from the liberation time analysis presented in Fig. 6. For the same reason, no $z_{50}$ is available for cluster Hrz137.
<table><caption>**Table A1.** Properties of the 12 HORIZON-AGN clusters. $R_{178}$ is the radius of a sphere around the halo centre with average DM density $\bar{\rho}_{DM} = 178 \rho_M$, $M_{178}$ is the total mass within $R_{178}$, and $M_{178,*}$ is the total stellar mass within $R_{178}$. $z_{50}$ is the redshift at which half of the final mass is assembled into the cluster halo. No $z_{50}$ is available for cluster Hrz137 due to issues constructing a reliable merger tree. $f_{BCG}$, $f_{ICL}$, and $f_{sat}$ are the fractions of total cluster stellar mass in the BCG, ICL, and satellite galaxies, respectively, with a 100 kpc spherical aperture used to separate the BCG and ICL. $\Delta v$ is the difference between the peaks of the DM and ICL velocity distributions, as described in Section 5.2.2.</caption><thead><tr><th>M<sub>178</sub> (×10<sup>14</sup> M<sub>☉</sub>)</th><th>M<sub>178,*</sub> (×10<sup>12</sup> M<sub>☉</sub>)</th><th>R<sub>178</sub> (Mpc)</th><th>z<sub>50</sub></th><th>f<sub>BCG</sub></th><th>f<sub>ICL</sub></th><th>f<sub>sat</sub></th><th>Δv (km s<sup>-1</sup>)</th><th>Cluster Name</th></tr></thead><tbody><tr><td>1.18</td><td>2.84</td><td>1.13</td><td>1.24</td><td>0.35</td><td>0.29</td><td>0.36</td><td>270</td><td>Hrz009</td></tr><tr><td>1.22</td><td>3.58</td><td>1.15</td><td>0.33</td><td>0.24</td><td>0.18</td><td>0.58</td><td>112</td><td>Hrz183</td></tr><tr><td>1.23</td><td>3.28</td><td>1.15</td><td>0.75</td><td>0.24</td><td>0.24</td><td>0.52</td><td>172</td><td>Hrz048</td></tr><tr><td>1.40</td><td>3.53</td><td>1.19</td><td>0.57</td><td>0.24</td><td>0.19</td><td>0.56</td><td>179</td><td>Hrz046</td></tr><tr><td>1.46</td><td>4.09</td><td>1.21</td><td>0.27</td><td>0.24</td><td>0.17</td><td>0.59</td><td>445</td><td>Hrz132</td></tr><tr><td>1.60</td><td>4.78</td><td>1.24</td><td>0.37</td><td>0.17</td><td>0.20</td><td>0.63</td><td>283</td><td>Hrz071</td></tr><tr><td>1.70</td><td>4.73</td><td>1.29</td><td>0.42</td><td>0.12</td><td>0.22</td><td>0.65</td><td>227</td><td>Hrz001</td></tr><tr><td>1.80</td><td>5.12</td><td>1.31</td><td>0.79</td><td>0.28</td><td>0.29</td><td>0.43</td><td>188</td><td>Hrz157</td></tr><tr><td>2.32</td><td>6.74</td><td>1.42</td><td>0.47</td><td>0.20</td><td>0.28</td><td>0.52</td><td>216</td><td>Hrz078</td></tr><tr><td>2.43</td><td>6.15</td><td>1.44</td><td>0.23</td><td>0.10</td><td>0.24</td><td>0.66</td><td>353</td><td>Hrz174</td></tr><tr><td>2.66</td><td>6.43</td><td>1.46</td><td>N/A</td><td>0.12</td><td>0.21</td><td>0.67</td><td>376</td><td>Hrz137</td></tr><tr><td>3.71</td><td>8.66</td><td>1.64</td><td>1.07</td><td>0.17</td><td>0.29</td><td>0.54</td><td>173</td><td>Hrz049</td></tr></tbody></table>

This paper has been typeset from a TeX/LATEX file prepared by the author.