# Brightest Cluster Galaxy Offsets in Cold Dark Matter

Cian Roche,1 Michael McDonald,1 Josh Borrow,1,2 Mark Vogelsberger,1,3 Xuejian
Shen1, Volker Springel4, Lars Hernquist5, Ruediger Pakmor4, Sownak Bose6,
Rahul Kannan7
1Department of Physics and MIT Kavli Institute for Astrophysics and Space
Research, 77 Massachusetts Avenue, Cambridge, MA 02139, USA
2Department of Physics and Astronomy, University of Pennsylvania, 209 South
33rd Street, Philadelphia, PA, USA 19104
3The NSF AI Institute for Artificial Intelligence and Fundamental
Interactions, Massachusetts Institute of Technology, Cambridge, MA 02139, USA
4Max-Planck-Institut für Astrophysik, Karl-Schwarzschild-Str. 1, D-85748,
Garching, Germany
5Harvard-Smithsonian Center for Astrophysics, 60 Garden Street, Cambridge, MA
02138, USA
6Institute for Computational Cosmology, Department of Physics, Durham
University, South Road, Durham DH1 3LE, UK
7Department of Physics and Astronomy, York University, 4700 Keele Street,
Toronto, ON M3J 1P3, Canada<EMAIL_ADDRESS>
(Accepted XXX. Received YYY; in original form ZZZ)

###### Abstract

The distribution of offsets between the brightest cluster galaxies of galaxy
clusters and the centroid of their dark matter distributions is a promising
probe of the underlying dark matter physics. In particular, since this
distribution is sensitive to the shape of the potential in galaxy cluster
cores, it constitutes a test of dark matter self-interaction on the largest
mass scales in the universe. We examine these offsets in three suites of
modern cosmological simulations; IllustrisTNG, MillenniumTNG and BAHAMAS. For
clusters above $10^{14}\,\rm{M_{\odot}}$, we examine the dependence of the
offset distribution on gravitational softening length, the method used to
identify centroids, redshift, mass, baryonic physics, and establish the
stability of our results with respect to various nuisance parameter choices.
We find that offsets are overwhelmingly measured to be smaller than the
minimum converged length scale in each simulation, with a median offset of
$\sim 1\,\rm{kpc}$ in the highest resolution simulation considered, TNG300-1,
which uses a gravitational softening length of $1.48\,\rm{kpc}$. We also find
that centroids identified via source extraction on smoothed dark matter and
stellar particle data are consistent with the potential minimum, but that
observationally relevant methods sensitive to cluster strong gravitational
lensing scales, or those using the the “light traces mass” approach, in this
context meaning gas is used as a tracer for the potential, can overestimate
offsets by factors of $\sim 10$ and $\sim 30$, respectively. This has the
potential to reduce tensions with existing offset measurements which have
served as evidence for a nonzero dark matter self-interaction cross section.

###### keywords:

cosmology: dark matter – galaxies: clusters – gravitational lensing

††pubyear: 2024††pagerange: Brightest Cluster Galaxy Offsets in Cold Dark
Matter–Brightest Cluster Galaxy Offsets in Cold Dark Matter

##  1 Introduction

The first evidence for the existence of a dominant and invisible matter
component in the universe came from observations of the Coma galaxy cluster in
1933 (Zwicky, 2009), for which the kinematics of the galaxies within suggested
an average density 400 times the inferred density from visible matter. The
evidence for this dark matter (DM) has since accumulated to include for
example the rotation curves of galaxies (Rubin & Ford, 1970; Ou et al., 2023),
their escape velocities (Necib & Lin, 2022), the cosmic microwave background
(CMB; Planck Collaboration et al. 2020), and gravitational lens modelling
(Clowe et al., 2006; Sharon et al., 2020). The particle nature of DM remains
one of the primary questions of modern astrophysics, but a cold, collisionless
DM in tandem with a cosmological constant $\Lambda$ constitutes the current
standard cosmological model (Weinberg et al., 2013; Vogelsberger et al.,
2020). The $\Lambda$ cold dark matter (hereafter $\Lambda$CDM) cosmological
paradigm has been remarkably successful in explaining the large scale
structure of the universe via the Lyman-alpha forest (Evslin, 2017; Addison et
al., 2018; Cuceu et al., 2019), the matter power spectrum and galaxy
clustering (Reid et al., 2010; Eisenstein et al., 2011), the accelerated
expansion of the universe (Riess et al., 1998; Perlmutter et al., 1999) and
observations of the cosmic microwave background anisotropies and polarizations
(Page et al., 2003).

Despite its successes, the $\Lambda$CDM model has been subject to increasing
tensions (Bull et al., 2016; Perivolaropoulos & Skara, 2022), many of which
relate to small scales such as the distribution, structure and diversity of
the low-mass DM halos of dwarf galaxies (see Sales et al. 2022 for a
comprehensive review). Certain problems such as the “core-cusp” problem, in
which CDM simulations predict greater dark matter density in the centers of
dwarf galaxies than inferred from observations (Flores & Primack, 1994; Moore,
1994), can be resolved with either baryonic feedback (Navarro et al., 1996a;
Weinberg et al., 2015) or with the introduction of a nonzero dark matter self-
interaction cross section (Spergel & Steinhardt, 2000; Davé et al., 2001;
Vogelsberger et al., 2016).

The properties of DM halos and their baryons in the self-interacting dark
matter (SIDM) paradigm have been well-studied in numerical simulations (e.g.
Vogelsberger et al., 2012; Rocha et al., 2013; Zavala et al., 2013; Elbert et
al., 2015; Robles et al., 2017; Elbert et al., 2018; Brinckmann et al., 2018;
ZuHone et al., 2019; Despali et al., 2019; Fitts et al., 2019; Robertson et
al., 2019; Vogelsberger et al., 2019; Shen et al., 2022; Nadler et al., 2023)
but are generally described by a flattening of the potential in the cores of
halos, particularly in DM-dominated systems like dwarf galaxies. For this
reason most efforts in measuring the shape of density profiles as probes of
the DM self-interaction cross section have focused on small scales, but
significant opportunity also exists on large mass scales in the cores of
galaxy clusters.

Galaxy clusters form at the densest nodes in the cosmic web, with the highest
dark matter densities in their cores. These cores constitute ideal testing
grounds to search for cored dark matter profiles, for example via strong lens
modelling (Zitrin et al., 2012; Limousin et al., 2022), the internal
kinematics of the brightest cluster galaxy (BCG; Kelson et al. 2002), weak
lensing, and combinations of these probes (Newman et al., 2013). Evidence for
both cored and cuspy cluster core profiles has been discovered using these
methods (Biviano et al., 2023), with the source of the measured coring not yet
clearly attributable to baryonic feedback or dark matter particle interactions
alone. This uncertainty is due in part to the fact that baryonic effects have
been demonstrated to both increase DM density slopes in the cores of clusters
via adiabatic contraction (Schaller et al., 2015a), and decrease these slopes
via dynamical friction (Inoue & Saitoh, 2011; Martizzi et al., 2012), galaxy
formation processes (Laporte et al., 2012) and active galactic nucleus
feedback (Ragone-Figueroa et al., 2012; Peirani et al., 2017).

The bulk motion of the BCG also serves as a promising probe of the cluster
core potential. In a simplified picture one may use the BCG, typically found
in the core of a galaxy cluster, as a test particle whose motion reveals the
shape of the underlying potential set by the DM. In a relaxed cluster, the BCG
is expected to rest close to the cluster core and oscillate on small length
scales; the amplitude of oscillation relative to the potential minimum then
reflects the inner slope of the density profile and therefore potentially the
dark matter particle physics. This probe has been referred to as “BCG
wobbling” (see Harvey et al. 2017, 2019). Since in real observations only a
snapshot of each cluster is available (due to dynamical timescales being much
larger than those relevant to observations), a statistical approach is taken
whereby the distribution of observed offsets between BCGs and their respective
cluster potential minima is understood to encode the shape of the potential
dictating the oscillations.

Alternatively, active mergers of galaxy clusters represent a similar but
distinct testing ground for SIDM physics, in particular by considering an
effective drag force on the DM halos of two merging clusters, which would
result in an offset between either the centroid and/or peak of the DM and
assumed collisionless stellar component of each merger member. Whether an
offset is expected to develop between the centroid and/or peak is dependent on
the macroscopic model of the action of the drag force (Markevitch et al.,
2004; Harvey et al., 2014), in addition to aspects of the microscopic
interactions such as the scattering rate and typical scattering angles
(Kahlhoefer et al., 2013). These models are summarized and tested by Kim et
al. (2017), in which the development of BCG-DM offsets are studied in the
context of staged equal-mass mergers in DM-only simulations, with stars
treated as tracer particles. These mergers are examined both in CDM and SIDM,
with varying self-interaction cross section per unit dark matter mass
$\sigma_{\rm{SI}}/m\in[0,1,3,10]\,\rm{cm^{2}\,g^{-1}}$ where zero is the CDM
case. Here it is observed that the BCG-DM offset oscillates with very little
damping post-merger, with typical offset sizes monotonically increasing with
self-interaction cross section. At $1\,\rm{cm^{2}\,g^{-1}}$ offsets are
measured to be on the order of $\sim 20\,\rm{kpc}$, whereas at
$10\,\rm{cm^{2}\,g^{-1}}$ offsets are $\sim 100\,\rm{kpc}$. Note that since
typical dark matter particle velocities increase with cluster mass, and
velocity-dependent dark matter self interaction models are preferred
(Kaplinghat et al., 2016; Sagunski et al., 2021), these studies which consider
velocity-independent dark matter self interaction on the cluster scale should
be understood as probing the high-velocity interaction cross section.

In each case, the expectation for CDM is that there are zero offsets between
bound DM and stellar distributions, with the assumption that they are both
effectively collisionless. A precise quantification of “zero” however has not
yet been fully established, and will be necessary to make the statement that
any observation or set thereof is consistent with SIDM offset distributions at
some self-interaction cross section and simultaneously inconsistent with CDM
predictions. In Kim et al. (2017) the CDM staged mergers exhibit oscillation
amplitudes on the order of $\sim 10\,\rm{kpc}$, which, depending on the
systematics of observations, may be indistinguishable from the $\sim
20\,\rm{kpc}$ offsets of the $\sigma_{\rm{SI}}/m=1\,\rm{cm^{2}\,g^{-1}}$ dark
matter model. Current observational constraints place upper limits on the
self-interaction cross section on the order of $\sim 1\,\rm{cm^{2}\,g^{-1}}$
(Randall et al., 2008) or even $\sigma_{\rm{SI}}/m<0.19\,\rm{cm^{2}\,g^{-1}}$
at 95% confidence (Andrade et al., 2021). For such a cross section, the SIDM
oscillation amplitude is likely to be indistinguishable from that of CDM when
viewed in the context of these staged merger results. This highlights the
necessity for increased precision in the understanding of CDM offsets.

Figure 1: Demonstration of projected dark matter, gas and stellar density
distributions in IllustrisTNG300-1 on scales relevant both to observations and
to measuring faithful offsets between centroid measures of each distribution
in simulations. All panels are centered on the group position (location of
most bound particle) of group 20, which has a mass of
$M_{200,\rm{mean}}=6.4\times 10^{14}\,\rm{M_{\odot}}$ and is shown at redshift
zero. _(Upper left)_ projected dark matter mass density on large scales in a
cluster core, in addition to an isophote of the projected density map and its
associated center (green). The selected isophote is matched in area to a
circle of radius equal to 70kpc, which is the typical size of the Einstein
radius for observations of a cluster on this mass scale. _(Upper right)_ the
large-scale gas which has been used as a tracer of the potential, and its
centroid as determined by SourceExtractor. The solid white ellipses in this
panel and the middle row represent a scaled version of the source’s elliptical
shape, also determined by SourceExtractor. _(Middle row)_ the small-scale dark
matter, stellar and gas projected mass densities, and their centroids as
determined by SourceExtractor. _(Bottom)_ comparison of the various centroid
identification methods and their relationship to the convergence radius, given
by $\sim 2.8\,\epsilon_{\rm{DM,\star}}$ where $\epsilon_{\rm{DM,\star}}$ is
the gravitational softening length of dark matter and star particles in the
simulation (see text). The circle with radius equal to the convergence radius
is centered on the most bound particle, ie. the center of the panel. The
lensing scale center and large scale gas center are greatly offset from this
position and are not visible on this panel.

The relaxed “wobbling” scenario has been studied in Harvey et al. (2017) by
comparing BCG offsets in observations of relaxed clusters to oscillation
amplitudes in the BAHAMAS suite of CDM simulations (McCarthy et al., 2017).
The observed offsets are calculated as the on-sky distance between the
positions of the BCGs and strong lensing centers, used as proxies for the
potential minimum or dark matter center (Medezinski et al., 2013; Markevitch
et al., 2004; Zitrin et al., 2013). Here it was found that offsets measured
via mimicking wobbles analytically and performing mock strong lensing analyses
on the simulation data were too small to be consistent with the observations,
suggesting some degree of coring and thus self-interaction in the
observations. This was developed into a constraint on the self-interaction
cross section of DM in Harvey et al. (2019) using offset distributions in
cosmological SIDM simulations and making ansatzes about the convergence of
median offsets below the gravitational softening length, finding that the
observational data used in Harvey et al. (2017) is consistent with CDM at the
$1.5\sigma$ level and prefers a cross section of
$\sigma_{\rm{SI}}/m<0.39\,\rm{cm^{2}\,g^{-1}}$ at the 95% confidence level.

A hallmark of offset measurements made in relaxed clusters in simulations is
that they are overwhelmingly below the gravitational softening length, and
thus constitute measurements in the unconverged regime. The gravitational
softening introduces an artificial coring independent of and degenerate with
the coring due to self-interactions of the DM, but is necessary to avoid
issues such as shot noise in the mass distribution unphysically affecting
small-scale dynamics (Ludlow et al., 2019). Extrapolating below the softening
length or even operating on scales similar to the softening length is known to
produce untrustworthy results (Power et al., 2003), suggesting that to make
reliable claims about dynamics and thus DM particle physics on scales of $\sim
1-10\,\rm{kpc}$, simulations with sufficiently small softening lengths111Note
that the softening length cannot be arbitrarily decreased, and appropriate
choices for this parameter depend on other factors such as the 2-body
relaxation time set primarily by particle number. See Ludlow et al. (2019) for
a detailed discussion. should be utilized, and the convergence of offsets with
respect to this parameter should be well established.

In this paper, we investigate the offsets between the BCGs of galaxy clusters
and various measures of their underlying oscillation centers across several
cold dark matter cosmological hydrodynamical simulations. We make no attempt
to distinguish between relaxed and unrelaxed clusters, instead aiming to
produce a distribution of offsets over a whole cosmological volume, enabling
comparisons with observations or simulation data which do not depend on
relaxation criteria. The remainder of the paper is structured as follows: In
Sec. 2 we outline the simulations used to calculate offset distributions. In
Sec. 3 we describe the procedure by which we generate mass maps which can be
used to extract centroids for each particle species in these simulations. In
Sec. 4 we define various centroid identification methods to be applied to the
mass maps, including methods relevant to real observations. We discuss the
offset distribution of each simulation considered, and the dependence of this
distribution on centroid identification, redshift, mass, baryonic physics and
various nuisance parameter choices in Sec. 5, and summarize these results in
Sec. 6.

Simulation |  Box Size [cMpc$\,h^{-1}$] | $N_{\rm{part}}$ |  $M_{\rm{b}}$ $[\rm{M_{\odot}}\,$$h^{-1}]$ |  $M_{\rm{DM}}$ $[\rm{M_{\odot}}\,$$h^{-1}]$ |  $\epsilon_{\rm{DM,\star}}$ [kpc$\,h^{-1}$] |  $\epsilon_{\rm{gas,min}}$ [ckpc$\,h^{-1}$] |  $N_{\rm{groups}}$ $>10^{14}\,\rm{M_{\odot}}$ |  $N_{\rm{groups}}$ $>10^{15}\,\rm{M_{\odot}}$
---|---|---|---|---|---|---|---|---
IllustrisTNG300-1 | 205 | $2500^{3}$ | $7.6\times 10^{6}$ | $4.0\times 10^{7}$ | 1 | 0.25 | 460 | 6
IllustrisTNG300-2 | 205 | $1250^{3}$ | $5.9\times 10^{7}$ | $3.2\times 10^{8}$ | 2 | 0.5 | 458 | 6
IllustrisTNG300-3 | 205 | $625^{3}$ | $4.8\times 10^{8}$ | $2.5\times 10^{9}$ | 4 | 1 | 441 | 6
MillenniumTNG740 | 500 | $4320^{3}$ | $2.0\times 10^{7}$ | $1.3\times 10^{8}$ | 2.5 | 0.25 | 6664 | 98
MillenniumTNG185 | 125 | $1080^{3}$ | $2.0\times 10^{7}$ | $1.3\times 10^{8}$ | 2.5 | 0.25 | 98 | 1
BAHAMAS | 400 | $1024^{3}$ | $8.1\times 10^{8}$ | $4.5\times 10^{9}$ | 4 | 4 | 2624 | 30
Table 1: Physical parameters of the simulations considered in this article.
Parameters are (from left to right): the simulation identifier, length of one
side of the cubic simulation box, $N_{\rm{part}}$ is both the initial number
of gas particles (some of which will be converted to stellar particles as time
evolves), and also the number of DM particles (which remains unchanged),
initial baryonic mass for moving-mesh code $M_{\rm{b}}$, dark matter particle
mass $M_{\rm{DM}}$, gravitational softening length of the dark matter and star
particles $\epsilon_{\rm{DM,\star}}$ at $z=0$, minimum gravitational softening
length of the adaptive gas softening $\epsilon_{\rm{gas,min}}$ in comoving
coordinates, number of friends-of-friends identified groups with mass above
$10^{14}\,\rm{M_{\odot}}$ at $z=0$, and number of friends-of-friends
identified groups with mass above $10^{15}\,\rm{M_{\odot}}$ at $z=0$.

##  2 Data

We make use of three modern suites of cosmological simulations; IllustrisTNG
(TNG; Nelson et al. 2021), MilleniumTNG (MTNG; Pakmor et al. 2023), and
BAHAMAS (McCarthy et al., 2017). By analyzing offsets in the context of
varying baryonic physics models and simulation parameters such as mass
resolution, box size, and gravitational softening length, stronger statements
can be made about the offset distribution prediction of CDM in simulations and
the convergence of offset distributions with respect to various simulation
nuisance parameters. A comparison of the simulation parameters for the data
used from each suite can be found in Table 1. In each simulation we utilize
only a minimum mass cut of $M_{200,\rm{mean}}>10^{14}\,\rm{M_{\odot}}$, where
$M_{200,\rm{mean}}$ is the mass enclosed within $R_{200,\rm{mean}}$, defined
as the radius within which the mean mass density is 200 times the mean matter
density of the universe. We use this mass scale to match clusters used to
infer offset measurements in the past, such as the bullet cluster (Randall et
al., 2008), baby bullet (Bradač et al., 2008), Abell 3827 (Williams & Saha,
2011) and El Gordo cluster (Jee et al., 2014), which occupy masses from $\sim
10^{14}-10^{15}\,\rm{M_{\odot}}$. See Ng et al. (2017) for a detailed
comparison of existing observations.

Since the WMAP 9-yr cosmology cosmology used in the BAHAMAS simulations
(Bennett et al. 2013,
$\Omega_{\rm{m}}=\Omega_{\rm{CDM}}+\Omega_{\rm{b}}=0.2793,\>\Omega_{\rm{b}}=0.0463,\>\Omega_{\Lambda}=0.7207,\>\sigma_{8}=0.812,\>n_{s}=0.972,\>h=0.700$)
is inconsistent with that used in both TNG and MTNG (Planck Collaboration et
al. 2016;
$\Omega_{\rm{m}}=\Omega_{\rm{CDM}}+\Omega_{\rm{b}}=0.3089,\>\Omega_{\rm{b}}=0.0486,\>\Omega_{\Lambda}=0.6911,\>\sigma_{8}=0.8159,\>n_{s}=0.9667,\>h=0.6774$),
we choose to analyze the data from each simulation with the corresponding
cosmology, and it should be understood that there is a minor degree of
incompatibility in all following comparisons, but we expect this to have no
bearing on the measurements considered here.

###  2.1 IllustrisTNG

The IllustrisTNG project (TNG; Springel et al. 2017; Pillepich et al. 2017b;
Marinacci et al. 2018; Naiman et al. 2018; Nelson et al. 2021) is a suite of
cosmological hydrodynamical simulations spanning various simulation box sizes,
mass resolutions, and physics models based on the AREPO code (Springel, 2010).
This project succeeds the Illustris project (Vogelsberger et al., 2014a; Genel
et al., 2014; Vogelsberger et al., 2014b; Sijacki et al., 2015), improving
upon the galaxy formation model (Vogelsberger et al., 2013) and expanding the
scope to larger volumes and higher resolutions. The TNG simulation suite
consists of TNG50, TNG100 and TNG300, characterized by their varying
simulation box sizes of $35$, $75$ and $205\,\rm{cMpc}$$\,h^{-1}$
respectively.

We focus on TNG300 in this paper as we are interested in examining the largest
galaxy clusters, and this series of simulations admits the largest mass scales
in TNG. TNG300 contains 3 full-physics runs in a $205\,\rm{\,cMpc}$$\,h^{-1}$
box named TNG300-1, TNG300-2 and TNG300-3, and their DM-only counterparts
TNG300-1-Dark, TNG300-2-Dark, and TNG300-3-Dark, which share the box size and
softening lengths of their corresponding full physics runs (see Table 1) but
at slightly lower mass resolutions of $7\times 10^{7}$, $3.8\times 10^{8}$ and
$3\times 10^{9}\,\rm{M_{\odot}}$$\,h^{-1}$ than the full-physics counterparts
respectively. This collection of mass scales, softening scales and physics
models permit detailed resolution studies and the investigation of the effect
of baryonic physics on the galaxies and galaxy clustering properties, with
these simulations being used extensively in studies of galaxy clusters and
their environments: eg. cool core formation in galaxy clusters (Barnes et al.,
2018, 2019), supermassive black hole feedback (Weinberger et al., 2018), and
gas stripping dense environments in galaxy clusters (Yun et al., 2019).

Groups in IllustrisTNG are identified via a friends-of-friends algorithm with
a linking length $b=0.2$ in units of the mean interparticle separation
(Springel et al., 2001) with substructure identified via subfind (Dolag et
al., 2009). This results in a total of $\sim 450$ groups above
$M_{200,\rm{mean}}=10^{14}\,\rm{M_{\odot}}$ in TNG300 at $z=0$ and 6 above
$10^{15}\,\rm{M_{\odot}}$ which are considered in this article, as in Table 1.
The original Illustris simulations have been used to study DM-BCG offsets and
their dependence on centroid identification methods (Ng et al., 2017), but a
detailed study of BCG offsets has not previously been conducted in
IllustrisTNG.

###  2.2 MillenniumTNG

The MillenniumTNG project (MTNG; Hernández-Aguayo et al. 2023; Pakmor et al.
2023; Barrera et al. 2023) is the successor to the dark matter-only Millennium
simulation (Springel et al., 2005), the first 10 billion particle simulation
in a $500\,\rm{Mpc}$$\,h^{-1}$ box. MillenniumTNG combines the volume of
Millennium with the hydrodynamical modelling222With some small fixes in
addition to changes relative to the TNG model, implemented to save memory in
the much larger box (Pakmor et al., 2023). developed for the IllustrisTNG
project (Weinberger et al., 2016; Pillepich et al., 2017a) and is run using an
improved version of the AREPO code (Pakmor et al., 2015; Weinberger et al.,
2020). The MillenniumTNG simulations involve high time resolution merger trees
and direct light cone output, enabling detailed cosmological inference via
comparisons to large scale surveys (Contreras et al., 2023; Kannan et al.,
2023; Bose et al., 2023) and studies of the impact of baryons and massive
neutrinos on weak lensing observations (Ferlito et al., 2023). The
MillenniumTNG suite consists of numerous dark matter-only runs at various
resolutions, two hydrodynamical runs named MTNG740 and MTNG185 in analogy with
the IllustrisTNG naming scheme, and a collection of dark matter + massive
neutrino runs over a range of resolutions and neutrino masses. In this article
we consider only the hydrodynamical runs.

The significantly larger box size of MTNG enables the study of much more
massive clusters, with $98$ groups above
$M_{200,\rm{mean}}=10^{15}\,\rm{M_{\odot}}{}$ at $z=0$ in comparison to the
$6$ of TNG. In the context of DM-BCG offsets, this permits an analysis of
offset distributions at various mass cuts relevant to strong lensing
observations while maintaining sufficient statistics to study the tail of the
distribution, a feature necessary to quantify in order to perform hypothesis
tests using real observations. The dependence of offset distributions on mass
cuts is discussed in Section 5. The two hydrodynamical simulation runs in MTNG
possess significantly different box sizes but identical gravitational
softening length and mass resolution (see Table 1). The motivation for these
parameter choices is a study of box size effects on the matter power spectrum
and the halo mass function (Hernández-Aguayo et al., 2023). While there is no
expectation that the offset distribution should be meaningfully influenced by
the box size (other than the aforementioned potential dependence of the
distribution on mass), we still make use of this data to perform a consistency
check in Sec. 5.

###  2.3 BAHAMAS

The BAHAMAS suite of cosmological hydrodynamical simulations (McCarthy et al.,
2017; Robertson et al., 2019) has been a valuable asset in testing BCG-DM
offsets, both in the CDM (Harvey et al., 2017) and SIDM paradigms (Harvey et
al., 2019). BAHAMAS (BAryons and HAloes of MAssive Systems) aims to study
observables of large scale structure cosmology, and is based upon a version of
gadget-3 (last described in Springel 2005) which is extended with subgrid
physics including radiative heating/cooling, chemical evolution and both
stelllar and active galactic nucleus feedback. These feedback prescriptions
were calibrated to reproduce the present-day galaxy stellar mass function and
hot gas fractions of clusters and groups, and the simulations have been
successful in reproducing observations of stellar mass fractions in central
galaxies, satellite galaxy dynamics, the local stellar mass autocorrelation
function, the X-ray and SZ scalings of local groups and clusters, and black
hole mass to velocity dispersion and stellar mass relations, among others.

The BAHAMAS simulations use a $400\,\rm{kpc}$$\,h^{-1}$ box following cosmo-
OWLS (Le Brun et al., 2014; McCarthy et al., 2014) and masses for the DM and
baryonic particles a factor of $\sim 100$ larger than the highest resolution
simulation considered here, IllustrisTNG300-1. As in Table 1, the
gravitational softening length for the DM and stars in this simulation is
identical (other than a small disagreement due to chosen cosmology) to that of
IllustrisTNG300-3 at $\epsilon_{\rm{DM,\star}}=4\,\rm{kpc}$$\,h^{-1}$. This
enables a comparison of offset measurements at a given softening length on
similar mass scales, but with a different cosmology, feedback model and
simulation code. For example the hydro method of BAHAMAS is based upon the
smoothed particle hydrodynamics of gadget-3 whereas the hydro method utilized
in IllustrisTNG is the moving voronoi mesh of AREPO. This provides a valuable
opportunity to establish whether the findings of the offset measurements in
this paper depend on some specific choice in the modelling which differs
between TNG and BAHAMAS, or if despite the differences in simulation
methodology the offset distributions are consistent.

##  3 Mass Maps

This paper uses particle data from N-body simulations to produce 2D mass maps
of galaxy clusters. These mass maps will be used to measure the offset between
the centroid of the stellar matter and measures of the oscillation center such
as the potential minimum or centroids of the DM or gas distributions, with the
purpose of illuminating the shape of the potential at the center of galaxy
clusters. In this section, we describe the process by which we transform the
particle position data of a galaxy cluster in an N-body simulation to a mock
pixel image of the mass column density, discretized into pixels.

A mass map is a 2D function on the plane of the sky which represents the
column mass density (units mass per area or per pixel). In order to make a
mass map from simulation data which can be meaningfully compared to similar
measurements in real observations, it is necessary to calculate the projected
mass column density over a grid of pixels. In each pixel the density is
calculated by smoothing each particle close to the line of sight with a
kernel, and adding the contributions from all such kernels at the position of
the pixel. The length of the line of sight is chosen to be $10\,\rm{Mpc}$,
centered on the most bound particle in the cluster: this ensures the entire
cluster is included in the column density without erroneously including
contributions from nearby, unbound matter. To ensure this choice does not bias
our results, we test offset distributions in IllustrisTNG300-1 at various
line-of-sight depths in Sec. 5.

To perform the smoothing, we use the SWIFTsimIO projection module with the
subsampled backend (Borrow & Borrisov, 2020; Borrow & Kelly, 2021), which
smooths particles to a fixed resolution grid using the Wendland-C2 kernel
(Wendland, 1995; Dehnen & Aly, 2012) with $\gamma=1.936$. We choose this
method for smoothing as it ensures well-converged results for the mass maps
even at low pixel resolutions. To calculate the smoothing lengths for the star
and DM particles, we use the SWIFTsimIO generate_smoothing_lengths method and
for the gas we calculate the smoothing length for a Voronoi cell $i$ as
$h_{i}=(m_{i}/V_{i})^{1/3}$ where $m_{i}$ and $V_{i}$ are the total mass and
volume of the Voronoi cell, respectively.

For each galaxy cluster of interest, we project the DM, stars and gas
particles onto a grid with pixels of side length $0.25\,\rm{kpc}$, centered on
the position of the most bound particle. We choose this pixel size to be
smaller than the physical size of Hubble Space Telescope pixels for objects at
redshifts $z\sim 0.5$, which is approximately $0.1^{\prime\prime}\sim
0.6\,\rm{kpc}$. The zoom box is chosen to have a side length of $50\,\rm{kpc}$
for the higher-resolution simulations (IllustrisTNG300-1, IllustrisTNG300-2
and MillenniumTNG), and a side length of $100\,\rm{kpc}$ otherwise. For the
large scale measurements described in Sec. 4 (lensing and large scale gas) we
use a $300\,\rm{kpc}$ zoom box at the same physical pixel size. Example mass
maps for the highest-resolution simulation considered here, IllustrisTNG300-1,
are shown in Fig. 1 along with the ellipses and various centroid
identification methods described in the next section. The choice of pixel
resolution and image field of view have the potential to impact offset
measurements via their impact on centroid identification, and so we test the
dependence of offset measurements on these parameters in Sec. 5.

In Fig. 1 and the remainder of this paper, we use the terminology “convergence
radius”, a parameter which varies in definition (Power et al., 2003; van den
Bosch & Ogiya, 2018) but also potentially from halo to halo in a simulation
once a definition is fixed. Following Ludlow et al. (2019), the convergence
radius may be estimated as $r_{\rm c}=0.055\,l$ where $l$ is the mean comoving
interparticle separation of dark matter, $L_{\rm{box}}/N_{\rm{part}}^{1/3}$,
with $L_{\rm{box}}$ the simulation box size and both parameters available in
Table 1 for the simulations considered here. This considers the two-body
relaxation between dark matter particles as the major mechanism affecting the
convergence of halo central density profiles. Meanwhile, gravitational
softening can also naturally “smooth” the matter distribution at a scale
smaller than a few times $\epsilon_{\rm{DM,\star}}$. For simulations studied
here, $\epsilon_{\rm{DM,\star}}$ values were specifically chosen to be of the
same order as (and slightly smaller than) $r_{\rm c}$. Therefore, for
simplicity in the rest part of this paper, we choose to take a conservative
estimate of $r_{c}=2.8\epsilon_{\rm{DM,\star}}$, motivated by the “spline
softening length" used by the GADGET code as the length scale above which
pairwise forces become exactly Newtonian (Ludlow et al., 2019), but note that
this choice is not impactful for the present analysis.

For the DM distribution, the mass per area is a reasonable metric for
extracting centroids comparable to observations, however in the case of
baryonic matter it is not stellar or gas mass that contributes to observations
directly, but rather the luminosity. Basing these conclusions on mass maps is
equivalent to a “mass traces light” assumption, not to be confused with the
“light traces mass” assumption used to infer the position of DM or trace the
potential using gas. The quality of this assumption can be tested by creating
mock images of the baryonic matter and assessing the difference in centroid
between this mock image and that obtained from a mass map. However for the
galaxy population considered here, that is giant ellipticals in the cores of
clusters, the observed ages and thus mass to light ratios will be highly
uniform, with little scatter between galaxies and within galaxies themselves
(Groenewald & Loubser, 2014). As a result, we identify centroids directly from
the mass maps and leave a detailed treatment of mock images to future work.

##  4 Centroid Identification

To measure offsets between the stellar mass distribution and a measure of the
oscillation center, one must first establish a centroid for the stars and
potentially other species like the DM or gas. In this section, we outline how
we identify the center of the BCG in addition to various measures of the
oscillation center. These are centroids obtained via the position of the most-
bound particle in a cluster, source extraction via thresholding as a means of
centroid fitting, and the centers of isophotes of cluster DM mass maps on
scales relevant to cluster strong lensing.

One could also consider the usage of the cluster center of mass as a proxy for
an oscillation center, however this probe is too sensitive to large-scale
structure which is unrelated to the small-scale oscillation of the BCG in the
core. In Ng et al. (2017), one method of determining the stellar centroid
considered was the position of the “BCG”, defined therein by the peak of the
i-band flux. Other choices considered included a weighted centroid of many
member galaxies, the weighting being the i-band luminosity, and the peaks of
heavily smoothed kernel density estimates (KDEs) of the galaxy luminosities.
The authors observe that the BCG method exhibits by far the least scatter,
therefore in the present work we do not consider the other methods, as when
fitting centroids directly on small scales we are interested in the most
faithful representation of an underlying “true” center.

###  4.1 Most Bound Particle

One method of approximating the oscillation center is using the potential
minimum, identified by the position of the particle in the group of interest
with the lowest gravitational potential energy333This position is also used as
the “group position” field in the AREPO-based simulations (Springel, 2010;
Weinberger et al., 2020) and the center of the mass maps in Sec. 3.. This
centroid identification method has the advantage of being easily interpretable
via the simplified picture of a BCG oscillating in a fixed potential, however
since gravity is improperly resolved below the softening length in a
simulation, this position is uncertain on a scale close to the softening
length. The simplified picture of a fixed potential in which the BCG test mass
oscillates ignores the contribution to the potential of the BCG itself, whose
size ($\sim 100\,\rm{kpc}$) is significantly larger than the length scale of
oscillations ($\sim 10\,\rm{kpc}$). It is possible that this combined
potential minimum oscillates with the BCG but at lower amplitude, perhaps even
with discontinuities due to the most bound particle changing at each time
step, artificially reducing offset measurements. Whether these offsets are in
fact systematically smaller than other centroid identification methods will be
examined in Sec. 5. All references to an “oscillation center” in this paper
should therefore be treated with some nuance, as the underlying potential is
time-dependent.

Figure 2: Dark matter - brightest cluster galaxy centroid offsets in the cores
of galaxy clusters above $10^{14}\,\rm{M_{\odot}}$ in the IllustrisTNG,
MillenniumTNG and BAHAMAS cold dark matter cosmological hydrodynamical
simulations. Offsets are measured by finding the dark matter and stellar
density centroids using SourceExtractor on projected mass density maps
smoothed from the raw particle data using a smoothed particle hydrodynamics
interpolation scheme. (_Inset_) the median and $68\%$ bounds of the offset
distribution for each simulation, including an additional data point for the
MillenniumTNG185 run, which possesses the same gravitational softening length
but a smaller box than its MillenniumTNG740 counterpart. Convergence refers to
the convergence length, calculated as $2.8\epsilon_{\rm{DM}}$ where
$\epsilon_{\rm{DM}}$ is the gravitational softening length (see text).

###  4.2 Centroid Fitting

The offsets between centroids of different particle species, such as DM-stars
or gas-stars offsets, can also be used as a proxy for the shape of the
potential. To locate faithful centroids of these distributions, one can
utilize the mass maps of Sec. 3. A centroid can be extracted from these maps
by various means, for example by assuming a profile shape and fitting the
profile using a minimization or posterior sampling scheme. One can also
extract centroids without assuming a shape for the profile using an algorithm
like shrinking spheres in 3D (eg. Schaller et al. 2015b), shrinking apertures
in 2D, (eg. Ng et al. 2017), or by thresholding, whereby an underlying
statistical variation is identified as a “background”, and any pixel a certain
threshold above this background is classified as belonging to an object (eg
Harvey et al. 2019). In the thresholding approach, the object pixels can then
be divided into different sources, if present. A standard software package
which performs this thresholding procedure is SourceExtractor (SE; Bertin &
Arnouts 1996; Barbary 2016) and we choose to fit centroids of the DM, star and
gas distributions using this method as opposed to profile fitting or
convergent shrinking algorithms. This has the potential advantage of mimicking
some systematics of the procedure applied to astronomical images to identify
centroids of stellar and gas distributions, and we treat the DM in the same
manner for consistency, noting that such a measurement would not be possible
for real observations. This choice is consistent with the treatment of the
stellar maps of Harvey et al. (2019), in which the stellar centroid was
obtained via thresholding with SourceExtractor and the dark matter
centroid/oscillation center was taken as the potential minimum (Sec. 4.1).
Centers of stellar mass maps in this article are always identified using this
thresholding procedure.

After extracting sources in a mass map, we choose the center of the source
with the highest “flux” as the centroid in our image, since the highest flux
is very likely to be associated with the flux of the large-scale source
including the center. The term “flux” is borrowed from the extraction methods
of SourceExtractor, and should be understood here in the context of the mass
traces light assumption. Examples of ellipses describing the center and
orientation of such sources are shown in Fig. 1 in the “DM”, “stars” and “gas”
panels.

The SE object identification algorithm relies on a number of nuisance
parameter choices, for example the kernel size used to determine the
“background” noise level, and the threshold above this level for a detection
to be made. This threshold and other nuisance parameter choices are not
expected to have a significant impact on the determination of the position of
the central source, however we test the stability of the offset measurement to
these choices in Sec. 5. When extracting sources and thus identifying centers,
we choose background extraction parameters as a function of the pixel
resolution of the mass map, denoted $N_{\rm{pix}}\times N_{\rm{pix}}$ as
bw=bh=$\lceil
N_{\rm{pix}}/4\rceil$,fw=$\lceil$bw$/10\rceil$,fh=$\lceil$bh$/10\rceil$ (where
$\lceil x\rceil$ represents the ceiling function) and threshold=$1.5$ which
denote the widths and heights of the background boxes and filters of the SE
detection algorithm, respectively, and the signal to noise threshold above
which an object is considered detected.

###  4.3 Area-Matched Isophote

In observations, the oscillation center has been estimated using the center of
the primary halo in a strong lensing model (Markevitch et al., 2004;
Medezinski et al., 2013; Zitrin et al., 2013). Strong lens modelling is
sensitive to the large-scale potential of the system set by the DM, and not
necessarily to the precise shape of the core. The large-scale information can
however be used to “triangulate” the center under the assumption that the DM
mass distribution possesses at least two axes of symmetry on the sky. An
example of such a distribution would be a normal distribution444The large-
scale halos in strong lens modelling are typically modelled by other profiles
such as generalized forms of the Navarro-Frenk-White profile or the pseudo-
isothermal sphere (Navarro et al., 1996b, 1997). in two dimensions $(x,y)$
with variance along each axis $\sigma_{x}$ and $\sigma_{y}$, not necessarily
equal. This distribution has symmetry axes along $x$ and $y$ and its contours
are ellipses with a shared and well-defined center, such that precise
knowledge of an outer contour yields precise knowledge of the distribution
center. More generally, the center can be identified even if the shape of the
distribution is unknown, provided it possesses this symmetry.

The shape of DM halos in cosmological simulations are typically described by
triaxial ellipsoids (Allgood et al., 2006), and the projection of such a
distribution onto the plane of the sky is again an ellipse, so we expect the
deviation from the above assumption to be small for relaxed systems, and thus
the lensing center to be a good approximation for the DM centroid. However,
for systems undergoing oscillations, the small-scale oscillations may not
propagate to larger scales on sufficiently short timescales, and can lead to
an offset between the center of the small-scale information and the large-
scale. In this case, the large-scale lensing center will be an inappropriate
proxy for the oscillation center and thus lead to spuriously large offset
measurements. This will be exacerbated in unrelaxed systems such as that shown
in Fig. 1, for which the strong lensing scale encompasses more than one large
DM halo and the small scale centroid is significantly offset ($\gtrsim
20\,\rm{kpc}$) from the large scale DM halo center.

To mimic this feature, that is an estimation of the DM centroid on a scale to
which lensing is sensitive, we introduce a new centroid identification method.
We utilize the mass maps of Sec. 3 and extract the mass isophotes555A term
typically used for photometry, but we treat the mass maps as analogous to
images. via photutils (Bradley et al., 2023). Unlike contours, the isophotes
possess a well-defined center and yield an approximation of the center
assuming elliptical symmetry, mimicking the assumptions of the lensing
analysis without the requirement that the isophote ellipses share a center. To
perform the isophote fitting, we use the fit_image method of the
isophote.Ellipse class, with parameters as follows: the minimum and maximum
isophote semi-major axes are set by the mass map pixel resolution
$N_{\rm{pix}}$ as $a_{\rm{min}}=0.1N_{\rm{pix}}$,
$a_{\rm{max}}=0.7N_{\rm{pix}}$. Note that
$N_{\rm{pix}}=(50\,\rm{kpc})/(0.25\,\rm{kpc\,pixel^{-1}})=200$ for the typical
$50\,\rm{kpc}$ mass map size considered here. We also use a minimum of 20 and
maximum of 50 iterations of the fitting algorithm, and use lenient convergence
parameters of maxgerr=$1.0$, fflag=$0.8$ and conver=$0.01$ to obtain estimates
even when isophote identification is challenging due to the presence of
substructure, for example. This lenience is justified by the lack of
uncertainty in the mass map of the simulation, and the fact that dealing with
regions of low signal to noise in the faint outskirts of images, for which
parameters such as maxgerr are intended, is less relevant in this context.

Upon identifying the isophotes, we choose the isophote with the same on-sky
area as a circle formed by the Einstein radius (Narayan & Bartelmann, 1997).
Instead of calculating an Einstein radius by assuming the cluster is a point
mass, or by assuming a shape for the mass distribution such as the NFW profile
or the isothermal sphere (Navarro et al., 1996b, 1997) for which analytical
Einstein radii have been calculated (Narayan & Bartelmann, 1997; Dúmet-Montoya
et al., 2013), we instead choose to estimate an approximate Einstein radius
from observations. For sources at redshift $z_{S}\sim 2$ and 37 galaxy cluster
lenses at $z_{L}\sim 0.5$ in a similar mass range to those clusters considered
in this paper, Sharon et al. (2020) find a median Einstein angle of
$\theta_{E}=10.8^{\prime\prime}$ which corresponds to a lens-plane radius of
approximately $70\,\rm{kpc}$. We use this distance as the Einstein radius when
selecting the area-matched isophote.

This Einstein area-matched isophote center will be subject to the same
potential issues of miscentering between small and large scales as
observations using strong lensing, but the precise nature of this effect will
not be properly modelled here. In future work, we will examine mock lensing
images of clusters in simulations and compare the offset distribution to that
measured on smaller scales, in particular assessing if this effect leads to a
discrepancy between observed offset distributions and those of simulations.
This can be accomplished by assuming statistics of a background distribution
of sources, or by utilizing light cones such as those output by MillenniumTNG
which ensures self-consistency of the background distribution with the lens.

##  5 Results

We now apply the centroid identification techniques of Sec. 4 to the mass maps
of Sec. 3 to identify the distribution of offsets in each simulation suite,
and investigate its dependence on gravitational softening length, centroid
identification method, redshift and the mass range of clusters considered. We
also briefly review the impact of the baryonic model in IllustrisTNG on the
density profiles of the cluster halos, in particular comparing the coring of
the profiles in the full-physics runs to those of the DM-only runs.

###  5.1 Dependence on Softening Length

The gravitational softening length introduces a flattening of the
gravitational potential below a scale $\epsilon_{i}$, which is identical for
the dark matter and stellar particles, though note that the gas uses an
adaptive softening length scheme with a minimum value
$\epsilon_{\rm{gas,min}}$ in IllustrisTNG and MillenniumTNG, and a fixed
softening length in BAHAMAS. This artificial coring could induce oscillations
of the BCG degenerate with the signal of interest, that is the oscillation due
to cored central potentials arising from DM self-interaction. We calculate the
offset distribution in the three resolution levels of IllustrisTNG300, two box
sizes with identical resolution of MillenniumTNG, and the single BAHAMAS
resolution considered, where offsets are measured as the 2D on-sky distance
from the SourceExtractor-identified DM centroid to the SourceExtractor-
identified stellar centroid. The effect of the centroid identification method
is investigated in Sec. 5.2.

In Fig. 2 we present these offset distributions across all simulations
considered, as measured by SE centroids. Here it can be seen that the majority
of the offset distributions are below the softening length scale in each case,
and that the IllustrisTNG 300-3 and BAHAMAS distributions are consistent,
while possessing almost identical softening lengths. There also exists a clear
trend of decreasing offsets with decreasing softening length, which is the
expectation in a physics model whose “true” offset distribution is consistent
with zero, or at least smaller than the smallest softening length considered
here, namely $1.48\,\rm{kpc}$ in IllustrisTNG 300-1. We also represent the
median and 68th percentile offsets as a function of softening length in the
inset panel of Fig. 2. Here it is made clear that with centroids measured by
SourceExtractor, the offset measurements are limited by the gravitational
softening and convergence properties, consistent with the findings of Schaller
et al. (2015b), Ng et al. (2017) and Harvey et al. (2019) in which offsets
measured between the DM centroid (or potential minimum) and BCG are typically
on the order of or below the softening length. The lower resolution
simulations do exhibit some offsets above softening, but these may be due to
resolution and zoom box nuisance parameters, as the relatively poor mass
resolution necessitates a larger zoom box used to create the mass maps. As
discussed in Sec. 5.5 the size of the zoom box positively correlates with the
median offset, as measuring the center on a larger scale has the potential to
deviate due to substructure.

We also note that the two MillenniumTNG runs yield highly consistent offset
distributions, despite the large difference in simulation box size and thus
mass scale. This feature will be further tested with the mass cuts of Sec.
5.3.

###  5.2 Dependence on Centroid Identification

Figure 3: Comparison of offset distributions in IllustrisTNG300-1 as measured
by the four oscillation center identification measures/proxies studied in this
article. Each distribution is measured via the offsets between that method’s
centroid and the SourceExtractor-identified BCG position, applied to all
groups with mass above $M_{200,\rm{mean}}=10^{14}\,\rm{M_{\odot}}{}$ at $z=0$.
Note that “Gas” refers to the large-scale gas measurement, and colors are
roughly matched to Figure 1. The cumulative probability is the probability
that an offset is measured as less than or equal to the corresponding $x$-axis
value.

The offset measurement relies fundamentally on the identification of two
points, one center for the stellar matter (BCG), and another point which
infers an oscillation center. We use SourceExtractor to extract the centroid
of the stellar matter, mimicking the procedure in real observations, but there
are multiple choices for estimating the oscillation center, such as the
potential minimum, a SourceExtractor centroid for the DM, or the potential
isophote method described in Sec. 4. It has also been assumed (eg. Cross et
al. 2023) that “light traces mass” (LTM), which is used as motivation to
estimate the oscillation center using that of the gas in an indirect
measurement. Together, this constitutes a set of four BCG offset measurements
which we discuss in this Section.

In Fig. 3 we represent the offset distribution of IllustrisTNG300-1 as
measured between SourceExtractor stellar centroids (“BCG”) and these four
choices for the oscillation center. Here one can see that the smallest offsets
are measured between the BCG positions and the potential minimum. This is
followed by the BCG-DM offsets, with the majority of the offset distributions
for these two methods measured on sub-softening scales. That the SE-SE offsets
are approximately double the potential minimum-SE offsets on average can be
seen as a twofold measurement of the same noise in source extraction, once in
identifying the BCG and again in identifying the DM centroid. This could also
be understood conversely by considering the potential minimum-SE offsets as a
reduction of the more faithful SE-SE offset measurements. This can be
understood in terms of the discussion of Sec. 4.1. The stability of the source
extraction to the choice of background identification threshold is established
in Fig. 6, demonstrating that in the interpretation of potential min-SE
offsets being more faithful, the source extraction has some intrinsic noise
for a given mass map that is unrelated to nuisance parameter choice. This
noise however is on sub-softening scales, and so does not meaningfully bias
any inference made using these simulation offset distributions.

Figure 4: Variation in dark matter-brightest cluster galaxy offsets in
MillenniumTNG740, as measured by centroid fitting with SourceExtractor, with
the redshift (left) and mass range (right) of clusters considered. Offset
distributions at varying redshift are obtained using offset measurements in
all clusters above $M_{200,\rm{mean}}=3\times 10^{14}\,\rm{M_{\odot}}{}$, the
scale of mass distribution decreasing with increasing redshift. The cumulative
probability is the probability that an offset is measured as less than or
equal to the corresponding $x$-axis value. The number of clusters above this
limit at each redshift is 1197 ($z=0$), 688 ($z=0.25$), and 311 ($z=0.5$). The
mass dependence is assessed at $z=0$, and the number of groups in each mass
bin are 438 ($[1,1.05]\times 10^{14}\,\rm{M_{\odot}}{}$), 600 ($[5,6]\times
10^{14}\,\rm{M_{\odot}}{}$), and 681 ($[1,3]\times
10^{15}\,\rm{M_{\odot}}{}$).

When the position of the DM centroid is measured via the area-matched
isophote, the offset distribution extends to offsets roughly a factor of 10
larger than the potential minimum and SE-SE distributions, with roughly $80\%$
of the lensing scale offset distribution being measured above the softening
length. The miscentering between scales leading to this discrepancy can be
seen in the example of Fig. 1, wherein the Einstein isophote center falls far
outside of the convergence radius.

This hints at significant implications for lensing analyses, in particular
those of a small number of systems. In such a case, it is difficult to make a
statistically significant claim that the observations are inconsistent with
the offset distribution of CDM, given that scale miscentering can increase
offset measurements by an order of magnitude relative to more faithful
measures. The total uncertainties on strong lensing centers relevant to these
analyses are typically on the order of $\pm\sim 1^{\prime\prime}$ (Zitrin et
al., 2015; Harvey et al., 2017), corresponding to a $\pm 5\,\rm{kpc}$ lens-
plane distance. If the lensing center in this approximate picture can be
offset from the underlying faithful DM center by up to $\sim 10\,\rm{kpc}$, as
in Fig. 3, then the relatively small uncertainties on the observed lensing
center may lead to offset measurements falsely interpreted as large
oscillation amplitudes, when in fact it is assumptions (eg. symmetry) which
are breaking down. This highlights again the caution which much be taken when
comparing offset distributions in simulations to those of observations, as not
only are many of these offsets in an unconverged regime, but also systematics
of the centroid measurement method can significantly and artificially inflate
offsets to appear inconsistent with CDM. A full treatment of the offset
distribution obtained via mock lensing analysis is necessary however to
compare this feature directly to observations.

Finally, the large scale gas-BCG offsets are the largest of any centroid
identification method, indicating that the gas is a poor tracer of the
stellar, and therefore dark matter in the core. This can also be seen in Fig.
1, as in the large-scale gas the centroid is offset from the potential minimum
by $\sim 50\,\rm{kpc}$, while the small scale structure of the gas is highly
clumpy, making a small scale measurement of the gas centroid only loosely
correlated with the underlying potential. That the gas is an unreliable tracer
of the potential minimum in simulations is a robust statement as the offsets
measured are above the convergence radius of the underlying dark matter.
Though it is possible that the gas physics is improperly resolved on small
scales, the large scale gas centroid is more relevant to observations and is a
poorer tracer of the potential minimum than even the small scale centroids.
Thus inferences made using the gas as a tracer for the potential center, on
both large and small scales, can lead to misinterpretations of the underlying
dark matter physics.

###  5.3 Dependence on Redshift and Mass Scale

The offset distribution can in principle depend on the redshift of the
simulation snapshot examined, in that the distribution of cluster masses
changes over cosmological time, as does the major merger rate (Fakhouri et
al., 2010; Rodriguez-Gomez et al., 2015). We test the redshift evolution of
the offset distribution using MillenniumTNG740 at redshifts of $z=0,0.25$ and
$0.5$. We choose MillenniumTNG740 for this test as the redshift evolution is
coupled with the mass cut tests, and those tests are more readily conducted
with the larger mass scale and greater statistics of MTNG. We present these
results in Fig. 4, where it can be seen that the offset distribution at each
redshift is entirely consistent, and furthermore the offset distribution in
each case is almost entirely sub-softening. If these methods are primarily
probing noise in the centroid identification method coupled to a lower bound
set by the softening scale, as demonstrated by Fig. 2, then it is not
surprising that the offset distribution is not changed in physical coordinates
with redshift, since the softening length is fixed at these redshifts in
MillenniumTNG.

Also in Fig. 4 are the $z=0$ offset distributions at various mass cuts in
MTNG, with the goal of testing whether the offset distribution depends on the
mass scale of clusters considered. Previous work in BAHAMAS suggests that the
median offset should increase with cluster mass between $10^{14}$ and
$10^{15}\,\rm{M_{\odot}}{}$ (Harvey et al., 2019) in SIDM runs, but that no
mass trends are found in CDM. We test the CDM prediction over a slightly wider
mass range, permitted by the larger box size of MTNG. We find no significant
difference between offset distributions at $M\sim 10^{14}\,\rm{M_{\odot}}{}$
and $M\sim 10^{15}\,\rm{M_{\odot}}{}$, in particular finding the vast majority
of offsets below the softening length, consistent with the results of Fig. 3
in IllustrisTNG300-1. This finding is supported by the comparison between the
MTNG740 and MTNG185 offset distributions in Fig. 2, whereby despite the
significantly different mass statistics, the offset distributions are highly
consistent.

Understanding the oscillation dynamics in the merging and relaxed phases (Kim
et al., 2017) of a galaxy cluster on mass scales relevant to observations
necessitates the study of high time-resolution simulation data of a merger
involving a halo on the $10^{14}-10^{15}\,{\rm{M}}_{\odot}$ mass scale, with a
softening length comparable to the high-resolution simulation suites
considered here ($\lesssim 3\,\rm{kpc}$) and a full treatment of the baryonic
physics. The snapshot output rate of all simulations considered here is
insufficient to successfully resolve merger dynamics in time, as for a length
scale of $100\,\rm{kpc}$ the crossing time of cluster members travelling at
$\sim 10^{3}\,\rm{km\,s^{-1}}$ in the center of a cluster will be on the order
of $\sim 100\,\rm{Myr}$. IllustrisTNG mini-snapshots are stored at $\sim
100-200\,\rm{Myr}$ intervals, meaning that a factor of $10-100$ increase in
snapshot output frequency would enable this analysis. Thus future works may
study the evolution of offsets at this higher time resolution to understand
both the merger and relaxed phases of oscillation, along with the faithfulness
of centering methods in each phase.

Figure 5: Density profiles of example halos in full-physics IllustrisTNG300
runs at three resolution levels, with the density profiles of their matched
dark matter (DM)-only run counterparts. Shown are the total density profiles
and DM density profiles in the full-physics run, and the DM density profiles
in the DM-only run. Halos are also matched across simulations (along columns)
in both mass and radius to below displayed precision. The convergence radius
is calculated as $2.8\epsilon_{\rm{DM}}$ where $\epsilon_{\rm{DM}}$ is the
gravitational softening length (see text).

###  5.4 Dependence on Baryons

We also test the influence of the baryonic physics model on the shape of
cluster core density profiles. If the baryonic physics model contributes to a
significant coring or cusping of the central density profile, this could lead
to offset measurements which are more effective probes of the baryonic physics
than the underlying DM self-interaction properties. In Fig. 5 we examine the
density profiles of IllustrisTNG300 at three resolution levels, and compare
these density profiles to their counterpart halos in the DM-only run, where
the full-physics to DM-only matching is described in Lovell et al. (2018);
Nelson et al. (2021). We also match halos according to their mass and radius
across simulation resolutions to assess if any trend in cluster size exists,
with the matching in both parameters at the $\lesssim$1% level.

In Fig. 5 it can be seen that the effect of the baryonic model on the DM
density profiles (dotted lines) is to increase the density in the core,
consistent with the adiabatic contraction caused by gas condensation and mass
accretion seen in other simulations (eg. Schaller et al., 2015a), but only
below the convergence radius. Since these changes take place primarily below
the convergence radius, robust claims about the impact of this baryonic model
on offset measurements cannot be made. For example offsets are overwhelmingly
measured below the softening scale, which may be in part due to increased
central densities as a result of this contraction, but this cannot be
concluded reliably from this data. It is also the case that no significant
trend appears in cluster size, with the evidence of contraction being
comparable at all mass levels shown.

One factor not explicitly considered here is the efficiency with which AGN
feedback ejects particles on the resulting offset distributions. It may be
beneficial to perform a systematic study of offsets in an isolated system with
variable AGN feedback parameters, or to establish the offset distribution in a
simulation with stronger feedback model such as SIMBA (Davé et al., 2019),
whose cluster mass distribution is similar in scale to that of MTNG185 due to
its $147.1\,\rm{cMpc}$ box. Such a study is however outside the scope of this
paper.

Figure 6: Nuisance parameter tests for BCG offset measurements presented as
cumulative probabilities, all performed in IllustrisTNG300-1 with the fiducial
parameters of Sec. 3 and Sec. 4 unless otherwise specified. The cumulative
probability is the probability that an offset is measured as less than or
equal to the corresponding $x$-axis value. All offsets calculated via
SourceExtractor centroids for stars and dark matter, except for the lensing
scale panel. Note the difference in scale for the rightmost column. (_Upper
left_) depth of the zoom box used to create the mass maps along the line of
sight. (_Upper middle_) pixel size used to generate mass maps. Note that for
this test only groups above $3\times 10^{14}\,\rm{M_{\odot}}$ were considered,
justified by the lack of dependence of the distribution on mass. (_Upper
right_) side length of the zoom box used to create the mass maps. (_Lower
left_) sightline used for creation of mass maps, for example “$z$” means the
mass map was calculated from the perspective of the positive $z$-axis. (_Lower
middle_) threshold used in SourceExtractor detection algorithm, in units of
the standard deviation of the identified background. (_Lower right_) Einstein
radius used as a lensing scale for the area-matched isophote method. Lensing
scale tests are all performed with a $300\,\rm{kpc}$ field of view. Arrows on
a panel represent the value of the parameter for that offset distribution
(matched in color) where appropriate.

###  5.5 Parameter Tests

We also examine the dependence of offsets on the various nuisance parameter
choices used in this analysis, such as the parameters of the box used to
create the mass maps, the parameters of source extraction thresholding
procedure, or Einstein radius used in identifying the area-matched isophote.
We show the results of these tests in Fig. 6, and describe the results below:

  1. 1.

_Depth_ : There is some potential for intervening matter along the line of
sight to increase offset measurements, but we see no evidence that this
meaningfully alters the offset distribution. However, properly quantifying the
effect of intervening matter requires examining the lightcone data for the
group of interest, and we leave this treatment to a future study.

  2. 2.

_Pixel Size_ : We test pixel sizes both significantly smaller and larger than
the choice adopted in this paper, but find no dependence of the offset
distribution on this choice above the softening length. Below the softening
length there is a non-monotonic variance in the distribution, but this regime
is unconverged and thus will not affect inference made with these results. In
the case of $10\,\rm{kpc}$ pixels the source extraction and thresholding are
likely breaking down, with the background window sizes and filters approaching
the size of the mass map.

  3. 3.

_Field of View_ : The field of view is a significant factor in the measurement
of offsets via mass maps, with $\sim 50$% of offsets above the softening
scale, reaching up to $10\,\rm{kpc}$ at a field of view of $150\,\rm{kpc}$.
This trend is generally described by offsets being measured up to $\lesssim
10$% of the field of view, but this may be correlated with choices of the
background identification parameters described in Sec. 4. This correlation
does not significantly alter the results of this paper, as source extraction
is used in part when attempting to obtain faithful representations of the
centroid for a species on small scales (where most offsets are below the
softening scale regardless). It is also utilized on large scales intended to
mimic real observations, for which the offsets measured via the large scale
gas centroid are larger by a factor of $\sim 10$ than the offsets inflated due
to field of view effects in the SE-SE measurement.

  4. 4.

_Sightline_ : We do not expect the offset distribution to depend on the
sightline used throughout all halos considered in the simulation box, and
indeed we observe no evidence of a dependence.

  5. 5.

_Detection Threshold_ : We observe a small scatter in offset distribution as a
function of the source extraction detection threshold nuisance parameter, in
particular with offsets becoming slightly smaller on average as the threshold
becomes more strict. This matches expectations, as this is equivalent to using
only the central portions of detected objects, which likely exhibit smaller
degrees of miscentering, though it should be noted that these offsets are
overwhelmingly in the unconverged regime, so differences on this scale are not
robustly interpretable.

  6. 6.

_Lensing Scale_ : We also test the dependence of the area-matched isophote to
BCG offset distribution as a function of the assumed lensing scale, all tested
in the same $300\,\rm{kpc}$ field of view zoom box. Offsets increase
monotonically with lensing scale, with $100\,\rm{kpc}$ lensing scale offsets
on the order of $\sim 10\,\rm{kpc}$. Though this is a greatly simplified
picture of lensing centers, even on scales of 20 (or 50) kpc the effect of
miscentering pushes 30% (or 60%) of the offset distribution above the
softening scale.

##  6 Conclusions

The distribution of offsets between the BCG of a galaxy cluster and an
appropriate measure of the underlying oscillation center encodes important
information about the shape of the cluster core potential, and therefore the
underlying dark matter particle physics. Given the challenges in identifying
appropriate centers in real observations, and the unknown nature of of the
dark matter underlying those observations, cosmological simulations constitute
a valuable test bed for developing centroid identification methods and testing
both cold dark matter and self-interacting dark matter physics. However, the
usage of this technique to understand dark matter presents a number of
challenges, such as the issue of appropriate centroid identification, relevant
length scales being potentially unconverged and the difficulty of connecting
offset distributions in simulations to sets of observed data.

We test three modern suites of cosmological hydrodynamical simulations with
zero dark matter self-interaction cross section at various resolutions, namely
IllustrisTNG, MillenniumTNG and BAHAMAS. We make the following observations:

  * •

In each simulation, offset measurements are overwhelmingly found below the
gravitational softening length when measured directly via source extraction of
the smoothed particle data or using the position of the most bound particle as
a measure of the cluster center, consistent with and quantifying the
expectation that these cluster-scale BCG offsets should be consistent with
“zero” in CDM. Additionally, we find no evidence for a dependence of the
offset distribution on cluster mass or redshift in simulations.

  * •

When measuring centroids via the “light traces mass” approach, whereby gas is
assumed to trace either the DM centroid or oscillation center, we find the gas
to be significantly offset from the DM in cluster cores in simulations on
scales above the softening length, and thus find this centroid identification
technique to likely be unreliable for identifying evidence of dark matter
self-interaction.

  * •

We introduce a novel centroid identification technique of area-matched
isophotes, intended to probe offsets on scales to which strong lensing
measurements are sensitive, the method typically employed in observations.
With this centroid identification method, we find offsets significantly larger
than the more faithful source extraction or potential minimum offsets,
increasing the median offset by a factor of $\sim 10$. This insight may serve
to weaken evidence of dark matter self-interaction obtained by comparing
observations to CDM simulation offset distributions.

We also find that the baryonic physics model does contribute to some
contraction of DM halos but is not clearly responsible for the overwhelmingly
sub-softening offset measurements, and that the measurements are made with
nuisance parameter choices in robust regimes that do not meaningfully alter
the conclusions of the analysis, with the most impactful nuisance parameter
(other than the scale chosen to assess miscentering relevant to lensing) being
the field of view of the data considered.

To constrain the nature of dark matter, and in particular its self-interaction
cross section, it is necessary to improve the simulation data, analysis
techniques and the catalogue of observations. In simulations, improvements can
be made by examining high time resolution data to study the dynamics in the
cores of clusters at low gravitational softening length. In the analysis of
the data, it will be valuable to conduct extensive mock lensing observations
in simulated clusters to establish the CDM BCG-strong lensing offset
distribution. Finally, the catalogue of observations with high-quality lens
models for which the BCG position is not chosen to coincide with the DM halo
center has great potential to be expanded both in the near term via a
homogeneous re-analysis of _Hubble_ archival data, and in the future via large
lensing surveys conducted by the Euclid and Roman space telescopes.

## Acknowledgements

This work makes use of the following software: Python (Van Rossum & Drake,
2009), numpy (Harris et al., 2020), scipy (Virtanen et al., 2020), astropy
(Astropy Collaboration et al., 2013, 2018), jupyter (Kluyver et al., 2016).
This research also made use of Photutils, an Astropy package for detection and
photometry of astronomical sources (Bradley et al., 2023). CR thanks Meredith
Neyer for numerous helpful discussions.

SB is supported by the UK Research and Innovation (UKRI) Future Leaders
Fellowship [grant number MR/V023381/1].

## References

  * Addison et al. (2018) Addison G. E., Watts D. J., Bennett C. L., Halpern M., Hinshaw G., Weiland J. L., 2018, ApJ, 853, 119
  * Allgood et al. (2006) Allgood B., Flores R. A., Primack J. R., Kravtsov A. V., Wechsler R. H., Faltenbacher A., Bullock J. S., 2006, Monthly Notices of the Royal Astronomical Society, 367, 1781
  * Andrade et al. (2021) Andrade K. E., Fuson J., Gad-Nasr S., Kong D., Minor Q., Roberts M. G., Kaplinghat M., 2021, Monthly Notices of the Royal Astronomical Society, 510, 54
  * Astropy Collaboration et al. (2013) Astropy Collaboration et al., 2013, A&A, 558, A33
  * Astropy Collaboration et al. (2018) Astropy Collaboration et al., 2018, AJ, 156, 123
  * Barbary (2016) Barbary K., 2016, Journal of Open Source Software, 1, 58
  * Barnes et al. (2018) Barnes D. J., et al., 2018, MNRAS, 481, 1809
  * Barnes et al. (2019) Barnes D. J., et al., 2019, MNRAS, 488, 3003
  * Barrera et al. (2023) Barrera M., et al., 2023, Monthly Notices of the Royal Astronomical Society, 525, 6312
  * Bennett et al. (2013) Bennett C. L., et al., 2013, ApJS, 208, 20
  * Bertin & Arnouts (1996) Bertin E., Arnouts S., 1996, A&AS, 117, 393
  * Biviano et al. (2023) Biviano A., et al., 2023, The Astrophysical Journal, 958, 148
  * Borrow & Borrisov (2020) Borrow J., Borrisov A., 2020, Journal of Open Source Software, 5, 2430
  * Borrow & Kelly (2021) Borrow J., Kelly A. J., 2021, preprint (arXiv:2106.05281)
  * Bose et al. (2023) Bose S., et al., 2023, Monthly Notices of the Royal Astronomical Society, 524, 2579–2593
  * Bradač et al. (2008) Bradač M., Allen S. W., Treu T., Ebeling H., Massey R., Morris R. G., von der Linden A., Applegate D., 2008, ApJ, 687, 959
  * Bradley et al. (2023) Bradley L., et al., 2023, astropy/photutils: 1.10.0, doi:10.5281/zenodo.1035865, https://doi.org/10.5281/zenodo.1035865
  * Brinckmann et al. (2018) Brinckmann T., Zavala J., Rapetti D., Hansen S. H., Vogelsberger M., 2018, MNRAS, 474, 746
  * Bull et al. (2016) Bull P., et al., 2016, Physics of the Dark Universe, 12, 56
  * Clowe et al. (2006) Clowe D., Bradač M., Gonzalez A. H., Markevitch M., Randall S. W., Jones C., Zaritsky D., 2006, ApJ, 648, L109
  * Contreras et al. (2023) Contreras S., et al., 2023, Monthly Notices of the Royal Astronomical Society, 524, 2489–2506
  * Cross et al. (2023) Cross D., et al., 2023, Examining the Self-Interaction of Dark Matter through Central Cluster Galaxy Offsets (arXiv:2304.10128)
  * Cuceu et al. (2019) Cuceu A., Farr J., Lemos P., Font-Ribera A., 2019, Journal of Cosmology and Astroparticle Physics, 2019, 044
  * Davé et al. (2001) Davé R., Spergel D. N., Steinhardt P. J., Wandelt B. D., 2001, ApJ, 547, 574
  * Davé et al. (2019) Davé R., Anglés-Alcázar D., Narayanan D., Li Q., Rafieferantsoa M. H., Appleby S., 2019, Monthly Notices of the Royal Astronomical Society, 486, 2827
  * Dehnen & Aly (2012) Dehnen W., Aly H., 2012, MNRAS, 425, 1068
  * Despali et al. (2019) Despali G., Sparre M., Vegetti S., Vogelsberger M., Zavala J., Marinacci F., 2019, MNRAS, 484, 4563
  * Dolag et al. (2009) Dolag K., Borgani S., Murante G., Springel V., 2009, MNRAS, 399, 497
  * Dúmet-Montoya et al. (2013) Dúmet-Montoya H. S., Caminha G. B., Makler M., 2013, A&A, 560, A86
  * Eisenstein et al. (2011) Eisenstein D. J., et al., 2011, AJ, 142, 72
  * Elbert et al. (2015) Elbert O. D., Bullock J. S., Garrison-Kimmel S., Rocha M., Oñorbe J., Peter A. H. G., 2015, MNRAS, 453, 29
  * Elbert et al. (2018) Elbert O. D., Bullock J. S., Kaplinghat M., Garrison-Kimmel S., Graus A. S., Rocha M., 2018, ApJ, 853, 109
  * Evslin (2017) Evslin J., 2017, Journal of Cosmology and Astroparticle Physics, 2017, 024
  * Fakhouri et al. (2010) Fakhouri O., Ma C.-P., Boylan-Kolchin M., 2010, MNRAS, 406, 2267
  * Ferlito et al. (2023) Ferlito F., et al., 2023, Monthly Notices of the Royal Astronomical Society, 524, 5591–5606
  * Fitts et al. (2019) Fitts A., et al., 2019, MNRAS, 490, 962
  * Flores & Primack (1994) Flores R. A., Primack J. R., 1994, ApJ, 427, L1
  * Genel et al. (2014) Genel S., et al., 2014, MNRAS, 445, 175
  * Groenewald & Loubser (2014) Groenewald D. N., Loubser S. I., 2014, MNRAS, 444, 808
  * Harris et al. (2020) Harris C. R., et al., 2020, Nature, 585, 357
  * Harvey et al. (2014) Harvey D., et al., 2014, Monthly Notices of the Royal Astronomical Society, 441, 404
  * Harvey et al. (2017) Harvey D., Courbin F., Kneib J. P., McCarthy I. G., 2017, Monthly Notices of the Royal Astronomical Society, 472, 1972
  * Harvey et al. (2019) Harvey D., Robertson A., Massey R., McCarthy I. G., 2019, Monthly Notices of the Royal Astronomical Society, 488, 1572
  * Hernández-Aguayo et al. (2023) Hernández-Aguayo C., et al., 2023, Monthly Notices of the Royal Astronomical Society, 524, 2556
  * Inoue & Saitoh (2011) Inoue S., Saitoh T. R., 2011, MNRAS, 418, 2527
  * Jee et al. (2014) Jee M. J., Hughes J. P., Menanteau F., Sifón C., Mandelbaum R., Barrientos L. F., Infante L., Ng K. Y., 2014, The Astrophysical Journal, 785, 20
  * Kahlhoefer et al. (2013) Kahlhoefer F., Schmidt-Hoberg K., Frandsen M. T., Sarkar S., 2013, Monthly Notices of the Royal Astronomical Society, 437, 2865
  * Kannan et al. (2023) Kannan R., et al., 2023, Monthly Notices of the Royal Astronomical Society, 524, 2594–2605
  * Kaplinghat et al. (2016) Kaplinghat M., Tulin S., Yu H.-B., 2016, Phys. Rev. Lett., 116, 041302
  * Kelson et al. (2002) Kelson D. D., Zabludoff A. I., Williams K. A., Trager S. C., Mulchaey J. S., Bolte M., 2002, The Astrophysical Journal, 576, 720
  * Kim et al. (2017) Kim S. Y., Peter A. H. G., Wittman D., 2017, Monthly Notices of the Royal Astronomical Society, 469, 1414
  * Kluyver et al. (2016) Kluyver T., et al., 2016, in Loizides F., Schmidt B., eds, Positioning and Power in Academic Publishing: Players, Agents and Agendas. pp 87 – 90
  * Laporte et al. (2012) Laporte C. F. P., White S. D. M., Naab T., Ruszkowski M., Springel V., 2012, Monthly Notices of the Royal Astronomical Society, 424, 747
  * Le Brun et al. (2014) Le Brun A. M. C., McCarthy I. G., Schaye J., Ponman T. J., 2014, Monthly Notices of the Royal Astronomical Society, 441, 1270
  * Limousin et al. (2022) Limousin M., Beauchesne B., Jullo E., 2022, A&A, 664, A90
  * Lovell et al. (2018) Lovell M. R., et al., 2018, Monthly Notices of the Royal Astronomical Society, 481, 1950
  * Ludlow et al. (2019) Ludlow A. D., Schaye J., Bower R., 2019, MNRAS, 488, 3663
  * Marinacci et al. (2018) Marinacci F., et al., 2018, Monthly Notices of the Royal Astronomical Society, 480, 5113
  * Markevitch et al. (2004) Markevitch M., Gonzalez A. H., Clowe D., Vikhlinin A., Forman W., Jones C., Murray S., Tucker W., 2004, The Astrophysical Journal, 606, 819
  * Martizzi et al. (2012) Martizzi D., Teyssier R., Moore B., Wentz T., 2012, Monthly Notices of the Royal Astronomical Society, 422, 3081
  * McCarthy et al. (2014) McCarthy I. G., Le Brun A. M. C., Schaye J., Holder G. P., 2014, Monthly Notices of the Royal Astronomical Society, 440, 3645
  * McCarthy et al. (2017) McCarthy I. G., Schaye J., Bird S., Le Brun A. M. C., 2017, MNRAS, 465, 2936
  * Medezinski et al. (2013) Medezinski E., et al., 2013, The Astrophysical Journal, 777, 43
  * Moore (1994) Moore B., 1994, Nature, 370, 629
  * Nadler et al. (2023) Nadler E. O., Yang D., Yu H.-B., 2023, ApJ, 958, L39
  * Naiman et al. (2018) Naiman J. P., et al., 2018, Monthly Notices of the Royal Astronomical Society, 477, 1206
  * Narayan & Bartelmann (1997) Narayan R., Bartelmann M., 1997, Lectures on Gravitational Lensing (arXiv:astro-ph/9606001)
  * Navarro et al. (1996a) Navarro J. F., Eke V. R., Frenk C. S., 1996a, MNRAS, 283, L72
  * Navarro et al. (1996b) Navarro J. F., Frenk C. S., White S. D. M., 1996b, ApJ, 462, 563
  * Navarro et al. (1997) Navarro J. F., Frenk C. S., White S. D. M., 1997, ApJ, 490, 493
  * Necib & Lin (2022) Necib L., Lin T., 2022, ApJ, 926, 189
  * Nelson et al. (2021) Nelson D., et al., 2021, The IllustrisTNG Simulations: Public Data Release (arXiv:1812.05609)
  * Newman et al. (2013) Newman A. B., Treu T., Ellis R. S., Sand D. J., 2013, The Astrophysical Journal, 765, 25
  * Ng et al. (2017) Ng K. Y., Pillepich A., Wittman D., Dawson W. A., Hernquist L., Nelson D. R., 2017, Offsets between member galaxies and dark matter in clusters: a test with the Illustris simulation (arXiv:1703.00010)
  * Ou et al. (2023) Ou X., Eilers A.-C., Necib L., Frebel A., 2023, arXiv e-prints, p. arXiv:2303.12838
  * Page et al. (2003) Page L., et al., 2003, The Astrophysical Journal Supplement Series, 148, 233
  * Pakmor et al. (2015) Pakmor R., Springel V., Bauer A., Mocz P., Munoz D. J., Ohlmann S. T., Schaal K., Zhu C., 2015, Monthly Notices of the Royal Astronomical Society, 455, 1134
  * Pakmor et al. (2023) Pakmor R., et al., 2023, Monthly Notices of the Royal Astronomical Society, 524, 2539
  * Peirani et al. (2017) Peirani S., et al., 2017, Monthly Notices of the Royal Astronomical Society, 472, 2153
  * Perivolaropoulos & Skara (2022) Perivolaropoulos L., Skara F., 2022, New Astronomy Reviews, 95, 101659
  * Perlmutter et al. (1999) Perlmutter S., et al., 1999, The Astrophysical Journal, 517, 565
  * Pillepich et al. (2017a) Pillepich A., et al., 2017a, Monthly Notices of the Royal Astronomical Society, 473, 4077
  * Pillepich et al. (2017b) Pillepich A., et al., 2017b, Monthly Notices of the Royal Astronomical Society, 475, 648
  * Planck Collaboration et al. (2016) Planck Collaboration et al., 2016, A&A, 594, A13
  * Planck Collaboration et al. (2020) Planck Collaboration et al., 2020, A&A, 641, A6
  * Power et al. (2003) Power C., Navarro J. F., Jenkins A., Frenk C. S., White S. D. M., Springel V., Stadel J., Quinn T., 2003, Monthly Notices of the Royal Astronomical Society, 338, 14
  * Ragone-Figueroa et al. (2012) Ragone-Figueroa C., Granato G. L., Abadi M. G., 2012, Monthly Notices of the Royal Astronomical Society, 423, 3243
  * Randall et al. (2008) Randall S. W., Markevitch M., Clowe D., Gonzalez A. H., Bradač M., 2008, The Astrophysical Journal, 679, 1173
  * Reid et al. (2010) Reid B. A., et al., 2010, MNRAS, 404, 60
  * Riess et al. (1998) Riess A. G., et al., 1998, The Astronomical Journal, 116, 1009
  * Robertson et al. (2019) Robertson A., Harvey D., Massey R., Eke V., McCarthy I. G., Jauzac M., Li B., Schaye J., 2019, MNRAS, 488, 3646
  * Robles et al. (2017) Robles V. H., et al., 2017, MNRAS, 472, 2945
  * Rocha et al. (2013) Rocha M., Peter A. H. G., Bullock J. S., Kaplinghat M., Garrison-Kimmel S., Oñorbe J., Moustakas L. A., 2013, Monthly Notices of the Royal Astronomical Society, 430, 81
  * Rodriguez-Gomez et al. (2015) Rodriguez-Gomez V., et al., 2015, MNRAS, 449, 49
  * Rubin & Ford (1970) Rubin V. C., Ford W. Kent J., 1970, ApJ, 159, 379
  * Sagunski et al. (2021) Sagunski L., Gad-Nasr S., Colquhoun B., Robertson A., Tulin S., 2021, Journal of Cosmology and Astroparticle Physics, 2021, 024
  * Sales et al. (2022) Sales L. V., Wetzel A., Fattahi A., 2022, Nature Astronomy, 6, 897
  * Schaller et al. (2015a) Schaller M., et al., 2015a, Monthly Notices of the Royal Astronomical Society, 452, 343
  * Schaller et al. (2015b) Schaller M., Robertson A., Massey R., Bower R. G., Eke V. R., 2015b, Monthly Notices of the Royal Astronomical Society: Letters, 453, L58–L62
  * Sharon et al. (2020) Sharon K., et al., 2020, The Astrophysical Journal Supplement Series, 247, 12
  * Shen et al. (2022) Shen X., Brinckmann T., Rapetti D., Vogelsberger M., Mantz A., Zavala J., Allen S. W., 2022, MNRAS, 516, 1302
  * Sijacki et al. (2015) Sijacki D., Vogelsberger M., Genel S., Springel V., Torrey P., Snyder G. F., Nelson D., Hernquist L., 2015, MNRAS, 452, 575
  * Spergel & Steinhardt (2000) Spergel D. N., Steinhardt P. J., 2000, Phys. Rev. Lett., 84, 3760
  * Springel (2005) Springel V., 2005, MNRAS, 364, 1105
  * Springel (2010) Springel V., 2010, MNRAS, 401, 791
  * Springel et al. (2001) Springel V., White S. D. M., Tormen G., Kauffmann G., 2001, Monthly Notices of the Royal Astronomical Society, 328, 726
  * Springel et al. (2005) Springel V., et al., 2005, Nature, 435, 629
  * Springel et al. (2017) Springel V., et al., 2017, Monthly Notices of the Royal Astronomical Society, 475, 676
  * Van Rossum & Drake (2009) Van Rossum G., Drake F. L., 2009, Python 3 Reference Manual. CreateSpace, Scotts Valley, CA
  * Virtanen et al. (2020) Virtanen P., et al., 2020, Nature Methods, 17, 261
  * Vogelsberger et al. (2012) Vogelsberger M., Zavala J., Loeb A., 2012, MNRAS, 423, 3740
  * Vogelsberger et al. (2013) Vogelsberger M., Genel S., Sijacki D., Torrey P., Springel V., Hernquist L., 2013, MNRAS, 436, 3031
  * Vogelsberger et al. (2014a) Vogelsberger M., et al., 2014a, Monthly Notices of the Royal Astronomical Society, 444, 1518
  * Vogelsberger et al. (2014b) Vogelsberger M., et al., 2014b, Nature, 509, 177
  * Vogelsberger et al. (2016) Vogelsberger M., Zavala J., Cyr-Racine F.-Y., Pfrommer C., Bringmann T., Sigurdson K., 2016, MNRAS, 460, 1399
  * Vogelsberger et al. (2019) Vogelsberger M., Zavala J., Schutz K., Slatyer T. R., 2019, MNRAS, 484, 5437
  * Vogelsberger et al. (2020) Vogelsberger M., Marinacci F., Torrey P., Puchwein E., 2020, Nature Reviews Physics, 2, 42
  * Weinberg et al. (2013) Weinberg D. H., Mortonson M. J., Eisenstein D. J., Hirata C., Riess A. G., Rozo E., 2013, Phys. Rep., 530, 87
  * Weinberg et al. (2015) Weinberg D. H., Bullock J. S., Governato F., Kuzio de Naray R., Peter A. H. G., 2015, Proceedings of the National Academy of Science, 112, 12249
  * Weinberger et al. (2016) Weinberger R., et al., 2016, Monthly Notices of the Royal Astronomical Society, 465, 3291
  * Weinberger et al. (2018) Weinberger R., et al., 2018, MNRAS, 479, 4056
  * Weinberger et al. (2020) Weinberger R., Springel V., Pakmor R., 2020, The Astrophysical Journal Supplement Series, 248, 32
  * Wendland (1995) Wendland H., 1995, Advances in Computational Mathematics, 4, 389
  * Williams & Saha (2011) Williams L. L. R., Saha P., 2011, Monthly Notices of the Royal Astronomical Society, 415, 448
  * Yun et al. (2019) Yun K., et al., 2019, MNRAS, 483, 1042
  * Zavala et al. (2013) Zavala J., Vogelsberger M., Walker M. G., 2013, MNRAS, 431, L20
  * Zitrin et al. (2012) Zitrin A., et al., 2012, The Astrophysical Journal, 749, 97
  * Zitrin et al. (2013) Zitrin A., Menanteau F., Hughes J. P., Coe D., Barrientos L. F., Infante L., Mandelbaum R., 2013, The Astrophysical Journal Letters, 770, L15
  * Zitrin et al. (2015) Zitrin A., et al., 2015, ApJ, 801, 44
  * ZuHone et al. (2019) ZuHone J. A., Zavala J., Vogelsberger M., 2019, ApJ, 882, 119
  * Zwicky (2009) Zwicky F., 2009, General Relativity and Gravitation, 41, 207
  * van den Bosch & Ogiya (2018) van den Bosch F. C., Ogiya G., 2018, MNRAS, 475, 4066