# Dark against luminous matter around isolated central galaxies: a comparative
study between modern surveys and IllustrisTNG

Pedro Alonso  Department of Astronomy, Shanghai Jiao Tong University, Shanghai
200240, China  Shanghai Key Laboratory for Particle Physics and Cosmology,
Shanghai 200240, China  Wenting Wang  Department of Astronomy, Shanghai Jiao
Tong University, Shanghai 200240, China  Shanghai Key Laboratory for Particle
Physics and Cosmology, Shanghai 200240, China  Jun Zhang  Department of
Astronomy, Shanghai Jiao Tong University, Shanghai 200240, China  Shanghai Key
Laboratory for Particle Physics and Cosmology, Shanghai 200240, China  Hekun
Li  Department of Astronomy, Shanghai Jiao Tong University, Shanghai 200240,
China  Shanghai Key Laboratory for Particle Physics and Cosmology, Shanghai
200240, China  Shi Shao  Key Laboratory for Computational Astrophysics,
National Astronomical Observatories, Chinese Academy of Sciences, Beijing
100101, China  Qi Guo  Key Laboratory for Computational Astrophysics, National
Astronomical Observatories, Chinese Academy of Sciences, Beijing 100101, China
Yanqin He  Tianjin Astrophysics Center, Tianjin Normal University, Tianjin
300387, China  Cai-Na Hao  Tianjin Astrophysics Center, Tianjin Normal
University, Tianjin 300387, China  Rui Shi  Department of Astronomy, Shanghai
Jiao Tong University, Shanghai 200240, China  Shanghai Key Laboratory for
Particle Physics and Cosmology, Shanghai 200240, China Wenting Wang
<EMAIL_ADDRESS>
###### Abstract

Based on independent shear measurements using the DECaLS/DR8 imaging data, we
measure the weak lensing signals around isolated central galaxies (ICGs) from
SDSS/DR7 at $z\sim 0.1$. The projected stellar mass density profiles of
satellite galaxies are further deduced, using photometric sources from the
Hyper Suprime-Cam (HSC) survey (pDR3). The signals of ICGs $+$ their extended
stellar halos are taken from Wang et al. (2021a). All measurements are
compared with predictions by the IllustrisTNG300-1 simulation. We find,
overall, a good agreement between observation and TNG300. In particular, a
correction to the stellar mass of massive observed ICGs is applied based on
the calibration of He et al. (2013), which brings a much better agreement with
TNG300 predicted lensing signals at
$\log_{10}M_{\ast}/\mathrm{M_{\odot}}>11.1$. In real observation, red ICGs are
hosted by more massive dark matter halos, have more satellites and more
extended stellar halos than blue ICGs at fixed stellar mass. However, in
TNG300 there are more satellites around blue ICGs at fixed stellar mass, and
the outer stellar halos of red and blue ICGs are similar. The stellar halos of
TNG galaxies are more extended compared with real observed galaxies,
especially for blue ICGs with $\log_{10}M_{\ast}/\mathrm{M_{\odot}}>10.8$. We
find the same trend for TNG100 galaxies and for true halo central galaxies.
The tensions between TNG and real galaxies indicate that satellite disruptions
are stronger in TNG. In both TNG300 and observation, satellites approximately
trace the underlying dark matter distribution beyond $0.1R_{200}$, but the
fraction of total stellar mass in TNG300 does not show the same radial
distribution as real galaxies.

hydrodynamical simulation — galaxy stellar halos — sky surveys — dark matter —
galaxy formation — gravitational lensing shear

##  1 Introduction

Most of the matter in our Universe is in the form of invisible dark matter.
Dark matter can collapse to form cosmic sheets, filaments and clumpy halos.
The growth of large and intermediate scales of cosmic structures can be
successfully modeled by the linear perturbation theory, while on small scales,
galaxies form through the gas cooling and condensation within dark matter
halos (e.g. White & Rees, 1978). Smaller halos and galaxies can merge with
larger halos, becoming the so-called substructures and satellite galaxies,
orbiting around the central dominant galaxy in the host halo. These satellite
galaxies would eventually be disrupted by tidal forces, becoming part of the
central galaxy. The stripped material from these satellites mainly builds the
outskirts of the central galaxy, forming the extended diffuse stellar halos.
Understanding the distribution of galaxies on small scales and their
connections to the host dark matter halos is one of the key topics in
astrophysics, but it is not trivial, since it is beyond the capability of the
linear perturbation theory and the baryonic physics complicates the
understanding.

A few empirical approaches have been developed to model the formation and
evolution of galaxies under the standard framework of cosmic structure
formation, including halo occupation distribution (HOD) models (e.g. Peacock &
Smith, 2000; Seljak, 2000; Scoccimarro et al., 2001; Berlind & Weinberg, 2002;
Zheng et al., 2005), semi-analytical models (SAMs) (e.g. White & Frenk, 1991;
Kauffmann et al., 1993; Springel et al., 2001a; De Lucia & Blaizot, 2007; Guo
et al., 2011b) and hydrodynamical simulations such as Illustris (e.g.
Vogelsberger et al., 2014a, b; Genel et al., 2014; Sijacki et al., 2015),
EAGLE (e.g. Schaye et al., 2015; Crain et al., 2015), and IllustrisTNG (e.g.
Weinberger et al., 2017; Pillepich et al., 2018a, b; Marinacci et al., 2018;
Naiman et al., 2018; Springel et al., 2018; Nelson et al., 2018).

HOD is mathematical, which models the probability of finding an $N$ number of
galaxies in a halo of mass $M$, through the probability distribution $P(N|M)$,
with the free parameters tuned to match the clustering of observed galaxies.
HOD widely assumes that the spatial distribution of satellites trace the dark
matter distribution. While more complicated HOD models, such as those
considering how the occupation number depends on the colors of the central
galaxies, have been developed (e.g. Skibba & Sheth, 2009), physical mechanisms
behind galaxy formation and evolution are not directly considered. SAM solves
the non-linearity of galaxy formation by analytical approximations and by
tracing the merger histories of galaxies, with the merger trees either
generated through Monte-Carlo methods or from dark matter only simulations.
Freedoms in adopted physical processes are tuned to match real observations at
redshift $z=0$ or higher redshifts. SAMs are computationally cheap, but they
have the drawback of depending on assumptions about the physics behind galaxy
formation.

Compared with HOD and SAM, numerical simulations are more sophisticated
approaches to study the galaxy-halo connection, although they are more
computationally expensive. With the inclusion of baryonic physics, dark matter
only simulations have evolved into more complex hydrodynamical simulations.
Further developments in the formation and evolution of galaxies include
detailed treatments of various physical processes, such as cooling, star
formation and evolution, chemical enrichment and gas recycling. The spatial
distribution of star particles and gas cells can be well resolved. However,
many physical processes adopted in hydrodynamical simulations are still below
the resolution limit and have to rely on assumptions of subgrid physics.

Given the uncertainties of the adopted physical processes, it is essential to
perform studies that directly compare and verify that these simulations are as
close as possible to the real Universe. Such comparisons are useful in terms
of guiding to which direction the models are supposed to be improved, whether
the included physics are reasonable and whether there are missing ingredients
(e.g. Guo et al., 2011a, 2012; Wang & White, 2012; Wang et al., 2014; Han et
al., 2015; Henriques et al., 2017; Merritt et al., 2020; Ardila et al., 2021;
Zhang et al., 2022b).

The IllustrisTNG set of simulations are one of the most up-to-date
hydrodynamical simulations of galaxy formation and evolution. It has been
shown in previous studies that TNG predictions show reasonable agreement with
real data, in terms of the global stellar mass and luminosity functions,
global color distribution, galaxy clustering and satellite abundance etc.
(e.g. Nelson et al., 2018; Springel et al., 2018; Pillepich et al., 2018a).
However, in terms ,of a few detailed and stringent comparisons, they still
deviate from real data, especially when the comparison is made against
morphology, the outer stellar halos of massive galaxies and when in
combination with weak lensing measurements. For example, Merritt et al. (2020)
compared a few galaxies from the Dragonfly Nearby Galaxies Survey with TNG
mass-matched counterparts, and found that real galaxies have less mass or
light at large radii, which they denote as a so-called “missing outskirts
problem". More recently and in combination with weak lensing measurements,
Ardila et al. (2021) reported that at the virial mass of $\sim
10^{14}\mathrm{M_{\odot}}$, TNG shows excess in the outer stellar mass.
Moreover, Renneby et al. (2020) claimed that TNG300 tends to predict $\sim$50%
higher lensing signals for red galaxies with
$10.2<\log_{10}M_{\ast}/\mathrm{M_{\odot}}<11.2$.

In a previous study, Wang et al. (2021a) measured the projected radial density
profiles of stellar mass in isolated central galaxies, in their extended
stellar halos, and in surviving bound satellite galaxies. The fractions of
total stellar mass versus total mass as measured from weak lensing signals are
further calculated based on the HSC internal S19 shear measurements, though
the weak lensing signals are noisy due to the small footprint. In this study,
we first improve the weak lensing measurements using a much larger shear
catalog based on galaxy images from the Dark Energy Camera Legacy Survey
(DECaLS; Dey et al., 2019) and the Fourier_Quad Method. We then measure the
projected stellar mass density profiles for satellites following the method of
Wang et al. (2021a), with a new threshold in the stellar mass of satellites to
ensure fair comparisons with TNG predictions. The projected stellar mass
density profiles for centrals and their diffuse stellar halos are directly
taken from Wang et al. (2021a). We compare these profiles as well as their
connections to each other with the corresponding predictions by IllustrisTNG.

In particular, the Fourier_Quad Method (Zhang et al., 2015, 2017, 2019) is a
shear measurement pipeline that calculates the shear in Fourier space. An
advantage of this method is that it does not require any assumption about the
galaxy morphology. In Zhang et al. (2017), a novel shear measurement approach
was introduced, the so called PDF symmetrization method, which recovers the
shear signal by symmetrizing the PDF of the shear estimators. Zhang et al.
(2017) proved that with this new approach, Fourier_Quad was allowed to reach
the Cramer-Rao Bound, which is the lower bound to the statistical variance of
an estimator.

For observational calculations throughout this paper, we adopt as our fiducial
cosmological model the first-year Planck cosmology (Planck Collaboration et
al., 2014), with $H_{0}=67.3\mathrm{kms^{-1}/Mpc}$,
$\Omega_{\mathrm{m}}=0.315$ and $\Omega_{\Lambda}=0.685$, to be consistent
with Wang et al. (2021a).

The structure of this paper is as follows. We introduce our observational and
simulation data in section 2. In section 3, we describe the Fourier_Quad
Method and introduce our methods of calculating the lensing signals, the
projected stellar mass density profiles for satellite galaxies, and extended
stellar halos centered on isolated central galaxies. We present our results in
section 4, and include discussions in section 5. We conclude in section 6.

##  2 Data

In this paper, we will present measurements of the weak lensing signals and
the projected stellar mass density profiles for photometric satellite galaxies
and for central galaxies $+$ their extended stellar halos. Our sample of
central galaxies are selected as those galaxies which are the brightest within
a certain volume, which we call isolated central galaxies, and all the
measurements are centered on these isolated central galaxies binned according
to stellar mass or color. Except for the projected stellar mass profiles for
isolated central galaxies and their stellar halos, which were measured by
stacking galaxy images in Wang et al. (2021a), and we directly take the
measurements from the previous study, all the other measurements are performed
independently in this study. All observational measurements will be compared
with the predictions by IllustrisTNG300 in detail. In the following, we first
introduce the IllustrisTNG suite of simulations, we then move on to describe
our selection of isolated central galaxies in both observation and TNG300,
followed by descriptions about the weak lensing data and the photometric
galaxies from the Hyper Suprime-Cam Imaging Survey used to calculate the
projected density profiles for satellite galaxies.

###  2.1 IllustrisTNG

The IllustrisTNG project111https://www.tng-project.org (Nelson et al., 2018)
is composed by a series of magneto-hydrodynamical simulations, run with the
moving-mesh code AREPO (Springel et al., 2019). These simulations are the
descendants of the original Illustris simulations, with improved galaxy
formation models (Marinacci et al., 2018; Naiman et al., 2018; Nelson et al.,
2018; Pillepich et al., 2018a; Springel et al., 2018; Nelson et al., 2019),
which describe the coupled evolution of gas, dark matter, stars, and black
holes. It includes extensive treatments of various galaxy formation and
evolution processes, including metal line cooling, star formation and
evolution, stellar feedback, chemical enrichment, gas recycling and magneto-
hydrodynamics. The simulations were carried out using the Planck 2015
$\Lambda$CDM cosmological model with parameters $\Omega_{\mathrm{m}}=0.3089$,
$\Omega_{\Lambda}=0.6911$, $\Omega_{\mathrm{b}}=0.0486$, $\sigma_{8}=0.8159$,
$n_{\mathrm{s}}=0.9667$, and $h=0.6774$ (Planck Collaboration et al., 2016a).
Dark matter halos are identified with the friends-of-friends (FoF) algorithm
(Davis et al., 1985). In each FoF group, substructures including galaxies are
identified with the SUBFIND algorithm (Springel et al., 2001b).

IllustrisTNG includes a total of 18 simulations, which can be grouped into
three main categories with different resolutions and box sizes ($\sim$50, 100,
and 300 Mpc). Each hydrodynamical simulation is run together with two lower-
resolution versions of it, and a dark-matter only version accompanying each of
them. The simulation outputs are in the form of 100 snapshots with redshifts
from $z=127$ to $z=0$. Each snapshot provides information about dark matter
particles, star particles, gas cells and black holes. Halo, subhalo (galaxy)
catalogs, and merger trees are also provided for each simulation.

In this paper, we mainly use the redshift $z=0$ snapshot of the TNG300-1
simulation (TNG300 hereafter), which has a box size of approximately 300 Mpc
in comoving units, and it has the highest resolution among the 300 Mpc
category. TNG300 has a dark matter particle mass of $3.98\times
10^{7}\mathrm{M_{\odot}}$/h, an average gas cell mass of $7.44\times
10^{6}\mathrm{M_{\odot}}$/h. The average mass of star particles is $7.8\times
10^{6}\mathrm{M_{\odot}}$/h. The resolution of galaxies, calculated as
approximately 100 times the average star particle mass, is $\sim
10^{9}\mathrm{M_{\odot}}$. Compared with TNG300, the resolution of TNG100 is
approximately an order of magnitude higher. The large box size of TNG300 helps
us to achieve better signals, especially for massive galaxies. However, as it
has been discussed in previous studies, properties of simulated galaxies may
not be necessarily fully converged in TNG300, though it explores larger volume
(e.g. Pillepich et al., 2018b; Weinberger et al., 2017; Pillepich et al.,
2018a). Thus in order to make sure that our results are not limited by
resolution, we also compare the predictions by TNG100-1 (TNG100 hereafter). We
do not show TNG100 results repeatedly. This is because due to the smaller
volume, TNG100 results are noisier. However, we do have double checked that
the conclusions in this paper do not show prominent differences across TNG300
and TNG100.

###  2.2 Isolated central galaxies in observation and TNG

In real observation, it is difficult to directly identify galaxies which are
true central galaxies of dark matter halos. In order to select a sample with a
high fraction of true halo central galaxies, we follow the method of Wang et
al. (2021a) to select those so-called isolated central galaxies, which are the
brightest within a given volume. This is because the actual central galaxies
of dark matter halos are often brighter than other companion galaxies around
them. The sample of isolated central galaxies (hereafter ICGs or ICG) are
selected from the spectroscopic Main galaxies of the seventh data release of
the Sloan Digital Sky Survey (SDSS/DR7; Abazajian et al., 2009). The parent
sample used for selection is from the New York University Value Added Galaxy
(NYU-VAGC; Blanton et al., 2005) catalog, which is flux limited down to $r\sim
17.7$.

An ICG is defined as the brightest galaxy in $r$-band magnitude, within a
cylinder in redshift space, formed by a circular region in projected sky with
a radius of virial radius and the height of the cylinder being $\pm$3 times
the virial velocity along the line of sight. In addition, the ICG should not
be within such a volume of another more massive galaxy in stellar mass. Here
the virial radius and velocity are estimated through the abundance matching
formula of Guo et al. (2010). Throughout this paper, our measurements of the
lensing signals, the projected stellar mass density profiles for these ICGs
$+$ their extended stellar halos, and for surrounding satellites will all be
centered on this sample of ICGs.

To ensure that no ICG candidates have brighter companions without
spectroscopic redshifts, we have used the photoz probability distribution
catalog of Cunha et al. (2009) to look for additional companions missed in the
spectroscopic sample. We further eliminate any candidate with a companion in
this catalog of equal or brighter $r$-band magnitude and projected within the
virial radius, unless the photometric redshift distribution of the potential
companion is inconsistent222The photoz probability distribution of the
companion gives a less than 10% of probability that it shares the same
redshift as the central galaxy. with the spectroscopic redshift of the
candidate. We refer to Wang et al. (2021a) for further details on the galaxy
selection and completeness of this sample.

To select ICGs in TNG, we follow exactly the same isolation criteria as for
real observation. Here we choose the $z$-axis as the line of sight, and the
$z$ coordinate of each galaxy is displaced according to its velocity along the
$z$-axis, to mimic the redshift space distortion of real galaxies. The virial
radius and velocity are estimated in the same way as for real observed ICGs,
instead of directly using the values from the simulation. Notably, the
magnitudes used to select simulated ICGs are those which involved galaxy dust
attenuation and solar neighborhood extinctions (Nelson et al., 2018). Nelson
et al. (2018) claims significant improvements with respect to the original
Illustris simulations after they apply this new dust model (model C in their
paper). With the sample of ICGs selected from TNG300, we can directly check
the fraction of true halo central galaxies. We find the fraction is as high as
89.1%, with 10.9% of satellite contamination.

In our analysis, we will separate our sample of observed and simulated ICGs
into red and blue populations. The separation in real observation is drawn
according to the $g-r$ color-mass diagram of SDSS spectroscopic Main galaxies.
For TNG, the separation is determined according to the $g-r$ color-mass
diagram of TNG galaxies, i.e., different from observation. For TNG galaxies,
the colors have included the dust model (Nelson et al., 2018), which directly
affects the color separation. The explicit color cuts are
$g-r=0.065\log_{10}M_{\ast}/\mathrm{M_{\odot}}+0.1$ for observation and
$g-r=0.03\log_{10}M_{\ast}/\mathrm{M_{\odot}}+0.4$ for TNG300. Table 1 shows
the number of red and blue ICGs in real observation and in TNG300, divided
into different stellar mass bins.

In this paper, when we try to validate our lensing signal measurements, we
will compare our measured signals with the previous measurements in Mandelbaum
et al. (2016) and Wang et al. (2016). Mandelbaum et al. (2016) performed their
measurements by using a sample of so-called locally brightest galaxies (LBGs).
The sample of LBGs are selected in a similar way as our sample of ICGs, but
with different selection criteria. Explicitly, LBGs are the brightest in
$r$-band within 1 Mpc in projection and $\pm$1000 km/s along the line of
sight. We refer to Wang et al. (2016) for details. This sample of LBGs has
been used for stacking weak lensing, X-ray and SZ signals (e.g. Planck
Collaboration et al., 2013; Anderson et al., 2015; Hernández-Monteagudo et
al., 2015; Planck Collaboration et al., 2016b; Mandelbaum et al., 2016) over a
wide range in stellar mass, and historically, they were referred to as LBGs,
so we maintain this nomenclature to distinguish them from our sample of ICGs.
In this paper, the sample of LBGs will only be used when we compare our
lensing signals with previous studies and when we test the effect of stellar
mass corrections (Section 3.4).

Stellar mass | $R_{200}$ | SDSS | TNG300
---|---|---|---
[$\log_{10}M_{\ast}/\mathrm{M_{\odot}}$] | [kpc] | red | blue | red | blue
11.4-11.7 | 758.00 | 1723 | 22,840 | 79 | 852
11.1-11.4 | 459.08 | 9056 | 42,006 | 172 | 2127
10.8-11.1 | 288.16 | 33,010 | 78,461 | 970 | 4741
10.5-10.8 | 214.80 | 50,499 | 61,813 | 5100 | 7984
10.2-10.5 | 173.18 | 43,011 | 31,183 | 13,137 | 3673
Table 1: The virial radius, $R_{200}$, of ICGs in different stellar mass bins,
and the number of red and blue ICGs in observation and TNG300. Here $R_{200}$
is defined as the radius within which the matter density is 200 times the
critical density of the universe. $R_{200}$ values are taken from Wang et al.
(2021a) and are based on ICGs selected from a mock galaxy sample of the Munich
semi-analytical model. Note that we use the same $R_{200}$ for red and blue
ICGs in a given stellar mass bin.

###  2.3 Weak Lensing Data

The Dark Energy Spectroscopic Instrument (DESI; DESI Collaboration et al.,
2016) is a spectroscopic project that aims to obtain optical spectra for
millions of galaxies, quasars and Milky Way stars over a 5 year period of
time. The DESI Legacy Imaging Surveys (Dey et al., 2019) are the precedented
programs to take source images and prepare input photometric target sources
for the DESI spectroscopic observations. It includes three imaging projects
(the Dark Energy Camera Legacy Survey, the Beijing–Arizona Sky Survey, and the
Mayall $z$-band Legacy Survey), and cover a total of $\approx 14,000$ square
degrees of the sky in three optical bands ($g$, $r$, and $z$). Among the three
imaging surveys, the Dark Energy Camera Legacy Survey covers the North
Galactic Cap region of the sky at $\mathrm{Dec}\leq 32^{\circ}$ and the South
Galactic Cap region of the sky at $\mathrm{Dec}\leq 34^{\circ}$, reaching
$5\sigma$ depths of 24.0, 23.4, and 22.5 in magnitudes for the $g$, $r$ and
$z$ bands, respectively.

In this study, our weak lensing shear data is constructed from the Dark Energy
Camera Legacy Survey (DECaLS) DR8 data using Fourier_Quad (see Section 3.1.1
for details). The shear catalog itself has passed the field distortion
test333The field distortion is an optical aberration found in all optical
systems that produces a distortion in the galaxy shapes similar to the cosmic
shear. (Zhang et al., 2019). According to Zhang et al. (2022a), the
multiplicative biases in DECaLS $r$ and $z$-bands are all in reasonable
ranges. However, in DECaLS DR8, the $g$-band data shows a larger bias compared
to $r$ and $z$ bands in the initial field distortion test. Therefore, we only
use the $r$ and $z$-bands data in this work.

6Our galaxy catalog for the above shear measurements is provided by Zou et al.
(2019). Zou et al. (2019) presents a catalog of approximately 0.19 billion
galaxies with $r<23$ mag in DECaLS, in three optical bands ($g$, $r$, and
$z$), and two infrared ones. We use the photometric redshifts provided by Zhou
et al. (2021), which adopts a machine learning approach to compute photometric
redshifts. Concretely, they applied a random forest regression model, a
popular supervised learning algorithm composed by several decision trees. For
the training process, they combined the data from 10 spectroscopic surveys to
form a training set of what they called "truth" redshift values.

###  2.4 The Hyper Suprime-Cam Survey imaging products and photometric sources

The Hyper Suprime-Cam Subaru Strategic Program (HSC-SSP) project (Aihara et
al., 2018a; Miyazaki et al., 2018, 2012; Komiyama et al., 2018; Furusawa et
al., 2018) is a wide-field multi-band (g,r,i,z,y) imaging survey that uses the
Hyper Suprime-Cam (HSC) on the Subaru $8.2m$ telescope. The project works on
three levels of depth: Wide, Deep, and Ultradeep, which aims to cover a total
of $1,400$, $27$, and $3.5$ square degrees of footprints, and with $r=26$,
$r=27$, and $r=28$ depths in magnitudes, respectively. In this paper, we use
photometric galaxies from the wide area of the third public data release
(pDR3) to measure the projected stellar mass density profiles of satellite
galaxies around our sample of ICGs. pDR3 (Aihara et al., 2022) covers 670
square degrees in all five bands, in the Wide layer. Besides, the HSC coadded
imaging product was used by Wang et al. (2019) and Wang et al. (2021a) to
measure the surface brightness and projected stellar mass density profiles of
ICGs and their stellar halos, and the measurements of Wang et al. (2021a) will
be compared with TNG predictions in this study. In the following, we briefly
introduce the HSC imaging products and photometric sources.

HSC-SSP data is processed with the HSC pipeline, which is a specialized
version of the LSST (Axelrod et al., 2010; Jurić et al., 2017) pipeline code.
HSC “visited” the same sky for multiple times. The HSC pipeline involves four
main steps: (1) processing of single exposure (one visit) image (2) joint
astrometric and photometric calibration (3) image coaddition and (4) coadd
measurement. In the first step, bias, flat field and dark flow are corrected
for. Bad pixels, pixels hit by cosmic rays and saturated pixels are masked and
interpolated. The pipeline first performs an initial background subtraction
before source detection. Detected sources are matched to external reference
catalogs in order to calibrate the zero point and a gnomonic world coordinate
system for each CCD. After galaxies and blended objects are filtered out, the
sky background is estimated and subtracted again. A secure sample of stars are
used to construct the PSF model. In the second joint calibration step, the
astrometric and photometric calibrations are refined by requiring that the
same source appearing on different locations of the focal plane during
different visits should give consistent positions and fluxes (Bosch et al.,
2018). In the third step, the HSC pipeline resamples images to the predefined
output skymap. It involves resampling both the single exposure images and the
PSF model (Jee & Tyson, 2011). Resampled images of different visits are then
combined together (coaddition). Images produced through coaddition are called
coadded images.

The last step is the most relevant for the measurement of the projected
stellar mass density profiles of photometric satellites. In this step,
photometric sources are detected, deblended and measured from the coadded
images. A maximum-likelihood detection algorithm is run independently for each
band at first. Background is estimated and subtracted once again. The detected
footprints and peaks of sources are merged across different bands to maintain
detections that are consistent over bands and to remove spurious peaks. These
detected peaks are deblended and a full suite of source measurement algorithm
is run on all objects, yielding independent measurements of source positions
and other properties in each band. A reference band of detection is then
defined for each object based on both the signal-to-noise ratio (SNR) and the
purpose of maximizing the number of objects in the reference band. Finally,
the measurement of sources are run again with the position and shape
parameters fixed to the values in the reference band to achieve the “forced”
measurement. The forced measurement brings consistency across bands.

For analysis in this paper, we select those photometric sources444The sources
are used to calculate projected stellar mass density profiles of photometric
satellites, so different from background photometric sources used to calculate
lensing signals from DECaLS. which are classified as extended, i.e., galaxies.
We only use galaxies which are brighter than 25 in $r$-band. As has been
discussed by Aihara et al. (2018b) and Wang et al. (2021b), the completeness
of photometric sources in HSC is very close to 1 at $r\sim 25$. We also
exclude sources with any of the following flags set as true in $g$, $r$ and
$i$-bands: bad, crcenter, saturated, edge, interpolatedcenter or
suspectcenter. Note, in principle, we can also use the photometric sources
from DECaLS, i.e., the imaging survey adopted for our weak lensing shear
measurements, to measure the projected stellar mass density profiles of
photometric satellites. However, as having been shown by Wang et al. (2021b),
DECaLS and HSC photometric sources give fully consistent results, and although
HSC is much smaller in footprint, the signals measured from HSC are comparable
to or sometimes slightly better than DECaLS, due to its much deeper flux
limit. Besides, the projected stellar mass density profiles of ICGs and their
stellar halos, as measured by Wang et al. (2021a) and will be investigated in
this study, are also based on HSC galaxy images. Thus we choose to use HSC
photometric sources when we investigate satellite galaxies in this study.

##  3 Methods

###  3.1 The Fourier_Quad Method and Lensing Signal Measurements

####  3.1.1 Fourier_Quad

Fourier_Quad (Zhang et al., 2015, 2017, 2019) is a unique shear measurement
pipeline that measures the shear signal in Fourier space, based on the
weighted moments of the power spectrum of the galaxy images. Here we briefly
describe Fourier_Quad and the readers can refer to Zhang et al. (2019) for
more details.

Fourier_Quad calculates the 2D power spectrum of the sources. The shear
estimators are defined as:

$G_{1}=-\frac{1}{2}\int{d^{2}\vec{k}(k^{2}_{x}-k^{2}_{y})T(\vec{k})M(\vec{k})},$
(1) $G_{2}=-\int{d^{2}\vec{k}k_{x}k_{y}T(\vec{k})M(\vec{k})},$ (2)
$N=\int{d^{2}\vec{k}\left[k^{2}-\frac{\beta^{2}}{2}k^{4}\right]T(\vec{k})M(\vec{k})},$
(3)

where $\vec{k}$ is the wave vector, and

$T(\vec{k})=\left|\tilde{W}_{\beta}(\vec{k})\right|^{2}/\left|\tilde{W}_{PSF}(\vec{k})\right|^{2},$
(4)
$M(\vec{k})=\left|\tilde{f}^{S}(\vec{k})\right|^{2}-\tilde{F}^{S}-\left|\tilde{f}^{B}(\vec{k})\right|^{2}+\tilde{F}^{B}.$
(5)

$T(\vec{k})$ represents the ratio between the power spectrum of an isotropic
Gaussian function, $\tilde{W}_{\beta}(\vec{k})$, and the power spectrum of the
original PSF. $\tilde{W}_{\beta}(\vec{k})$ is defined as:

$\tilde{W}_{\beta}(\vec{x})=\frac{1}{2\pi\beta^{2}}\mathrm{exp}\left(-\frac{\left|\vec{x}\right|^{2}}{2\beta^{2}}\right),$
(6)

where $\beta$ is the scale radius of the Gaussian function. $T(\vec{k})$
transforms the form of the original PSF to the desired Gaussian form in order
to correct for the PSF effect. Zhang et al. (2015) claims the ideal value of
$\beta$ to be slightly larger ($20-50\%$) than the scale radius of the
original PSF, to avoid singularities in the conversion. Additionally, the
value of $\beta$ should not be too large, in order to avoid information loss.
$M(\vec{k})$ corrects the power spectrum of the source by subtracting the
contributions of the background and the Poisson noise.
$\tilde{f}^{S}(\vec{k})$ and $\tilde{f}^{B}(\vec{k})$ are the Fourier
transformations of the galaxy and the background noise, respectively.
$\tilde{F}^{\mathrm{S}}$ and $\tilde{F}^{\mathrm{B}}$ are estimates of the
power spectrum of the Poisson noise for the galaxy image and the background,
respectively.

It can be shown that the ensemble averages of the shear estimators so defined
satisfy the following relations (Zhang & Komatsu, 2011) :

$\frac{\left<G_{1}\right>}{\left<N\right>}=g_{1}+O(g_{1,2}^{3}),\frac{\left<G_{2}\right>}{\left<N\right>}=g_{2}+O(g_{1,2}^{3}),$
(7)

where $g_{1}$ and $g_{2}$ are the two shear components. In a more recent work,
Zhang et al. (2017) proposed a new approach for the shear measurement, called
the PDF Symmetrization Method (PDF-SYM). PDF-SYM symmetrizes the probability
distribution function (PDF) of $G_{1}-\hat{g}_{1}\left(N+U\right)$ and
$G_{2}-\hat{g}_{2}\left(N-U\right)$ in order to determine the shear values
$\hat{g}_{1}$ and $\hat{g}_{2}$. This introduces two additional terms ($U$ and
$V$) to the three initially defined ($G_{1}$, $G_{2}$, and $N$). The new terms
are defined as:

$U=-\frac{\beta^{2}}{2}\int{d^{2}\vec{k}(k_{x}^{4}-6k_{x}^{2}k_{y}^{2}+k_{y}^{4})T(\vec{k})M(\vec{k})},$
(8)
$V=-2\beta^{2}\int{d^{2}\vec{k}(k_{x}^{3}k_{y}-k_{x}k_{y}^{3})T(\vec{k})M(\vec{k})}.$
(9)

$U$ involves the parity properties of the shear estimators, whereas $V$ is
needed to transform $U$ when a coordinate rotation appears in the shear
measurement. Zhang et al. (2017) proved that PDF-SYM allowed the shear
estimation to reach the lowest theoretical statistical limit (Cramer-Rao
limit).

Fourier_Quad has been validated by the field distortion test (see Section 2.3
above), and the readers can find more details in Zhang et al. (2019) and Zhang
et al. (2022a).

####  3.1.2 The observed lensing signals

In this section we describe the method used for the lensing signal calculation
in real observation. In weak lensing, the tangential component of the shear
signal $\gamma_{\mathrm{t}}$ is related to the density profile of the
foreground lens via:

$\Delta\Sigma(r_{\mathrm{p}})=\gamma_{\mathrm{t}}\Sigma_{\mathrm{cr}},$ (10)

in which $\Delta\Sigma$ is the excess surface density defined as:

$\Delta\Sigma(r_{\mathrm{p}})=\bar{\Sigma}(<r_{\mathrm{p}})-\Sigma(r_{\mathrm{p}}),$
(11)

where $\bar{\Sigma}(<r_{\mathrm{p}})$ is the averaged projected mass density
within radius $r_{\mathrm{p}}$, and $\Sigma(r_{\mathrm{p}})$ is the projected
mass density or surface density at radius $r_{\mathrm{p}}$.
$\Sigma_{\mathrm{cr}}$ is the critical surface density, defined as

$\Sigma_{\mathrm{cr}}=\frac{c^{2}}{4\pi
G}\frac{d_{\mathrm{s}}}{d_{\mathrm{l}}d_{\mathrm{ls}}},$ (12)

in physical units. Here $d_{\mathrm{s}}$, $d_{\mathrm{l}}$, and
$d_{\mathrm{ls}}$ are the angular diameter distances to the source, to the
lens, and between source and lens, respectively.

Conventionally, the lensing signal is calculated as the weight sum of all
source-lens pairs. In this project, we follow a different approach, based on
the PDF-SYM method. As described in section 3.1.1, Fourier_Quad derives the
shear as the value that better symmetrizes the PDF of the shear estimators.
For the lensing signals, this means that we need to symmetrize the following
PDF,

$G_{\mathrm{t}}-\frac{\widehat{\Delta\Sigma}}{\Sigma_{\mathrm{cr}}(z_{\mathrm{l}},z_{\mathrm{s}})}\left(N+U_{t}\right),$
(13)

in which $G_{\mathrm{t}}$ is the tangential shear estimator, which is derived
from $G_{1}$ and $G_{2}$ defined in Equation 1 and 2 through a coordinate
rotation. Similarly, $U_{t}$ is from rotating the the $U$ and $V$ components
defined in Equation 8 and 9. $\widehat{\Delta\Sigma}$ is the excess surface
density that one seeks to symmetrize the PDF of the tangential shear estimator
defined in Equation 13. We set a photometric redshift cut of $\Delta z=0.2$
between source and lens in order to mitigate the dilution effect of foreground
galaxies.

A few recent works have used this shear measurement method and have calculated
the lensing signals with it. For example, Fong et al. (2022) performed the
first measurements of the so called characteristic depletion radius using weak
lensing data. Wang et al. (2022) made use of the same weak lensing data to
study halo properties of different lens samples.

In Appendix A, we show a comparison between our calculated lensing signals
based on Fourier_Quad shear estimator and the previous measurement of
Mandelbaum et al. (2016) and Wang et al. (2016) based on shear measurements of
SDSS. The signals are centered on the same sample of LBGs. The agreement is
very good, which together with the excellent results on the field distortion
test (Zhang et al., 2019, 2022a), indicating that we can safely use the
lensing results based on the Fourier_Quad method for our analysis throughout
this paper.

####  3.1.3 Lensing signals in TNG

For TNG, the excess surface density can be directly calculated based on the
spatial distribution of all particles, following Equation 11. We include dark
matter particles, stars, wind particles, gas cells, and black holes. Radial
bins are defined in projected distance, $r_{\mathrm{p}}$.
${\bar{\Sigma}(<r_{\mathrm{p}})}$ includes all particles within radius
$r_{\mathrm{p}}$, and $\Sigma(r_{\mathrm{p}})$ includes all particles between
$r_{\mathrm{p}}-\Delta r_{\mathrm{p}}$ and $r_{\mathrm{p}}+\Delta
r_{\mathrm{p}}$, where $\Delta r_{\mathrm{p}}<<r_{\mathrm{p}}$.

ICGs in real observation are flux limited, whereas ICGs in TNG are volume
limited. In order to consider the flux limit in our observational data and
ensure fair comparison between observation and TNG300, we weight each ICG in
TNG by a maximum volume when calculating the excess surface density profiles.
The maximum volume is calculated from the limiting redshift or distance when
one can still observe the galaxy given the SDSS flux limit, by placing
simulated galaxies from the $z=0$ snapshot at a redshift where the apparent
magnitudes reach $r=17.7$ (flux limit of our ICGs from the SDSS Main galaxy
sample). Similar weights are assigned when we compute the projected stellar
mass density profiles for ICGs $+$ their stellar halos and for satellites (see
subsections below). However, we find including or not including this weight
show very small differences in our results, because ICGs are divided into a
few narrow stellar mass bins in our results, and the flux limit barely affects
the median or mean stellar mass and halo mass given the narrow bin width.

###  3.2 Projected Stellar Mass Density Profile of satellites

####  3.2.1 Observed satellite profiles

In real observation, we follow the approach of Wang et al. (2021a) to measure
the projected stellar mass density profiles of satellites (hereafter satellite
profiles). We adopt a satellite mass threshold of
$M_{\ast}>10^{9}\mathrm{M_{\odot}}$, that is, only satellites more massive
than this threshold can contribute to the final signals. The threshold is
chosen to match the resolution limit for TNG to ensure fair comparisons. The
readers can refer to Wang et al. (2021a) and Wang & White (2012) for details,
and here we briefly introduce the method.

The observed satellite profiles are measured by counting photometric galaxy
companions around spectroscopically identified ICGs. The photometric galaxies
are from HSC pDR3 (see Section 2.4). In our method, we first bin ICGs
according to their stellar masses or colors, and for a given bin, we loop
through all ICGs. For each ICG, we first assume that all its companions are at
the same redshift, and we use the spectroscopic redshift of this ICG to
calculate the intrinsic luminosity and rest frame colors of the companions.
Here we adopt the $K$-correction formulas of Westra et al. (2010). The stellar
mass of each companion is further calculated through an empirical relation
between the stellar-mass-to-light ratio ($M/L$) and the rest-frame color (see
Wang & White (2012) for details). The cumulated stellar masses of satellites
as a function of projected radius for ICGs in the same bin are added together
and averaged in the end.

We then correct the contamination by foreground and background sources, we
repeat exactly the same steps around a sample of random points within the HSC
footprint. The random points are assigned the same redshift and stellar mass
distributions as real ICGs. The averaged companion profiles around the random
points in different stellar mass and color bins are statistical estimates of
the foreground $+$ background source contamination, and are subtracted off
from the profiles around real ICGs.

The photometric sources of HSC are approximately fluxed limited down to 26 in
$r$-band apparent magnitude, and we adopt a safe flux limit of $r=25$. To
ensure the completeness of our satellite counts at
$M_{\ast}>10^{9}\mathrm{M_{\odot}}$, the flux limit is transferred to a lower
limit in stellar mass, again based on the redshift of the central ICG and the
$M/L$ estimated from the reddest color allowed at the corresponding redshift.
Here the reddest color is calculated through the red end boundary of all
spectroscopic SDSS galaxies at this redshift, which corresponds to the largest
allowed value of $M/L$, and thus the safest lower limit in stellar mass. Only
when this mass limit is smaller than $M_{\ast}=10^{9}\mathrm{M_{\odot}}$,
companions around this ICG will be cumulated and averaged in the end.

Notably, Wang et al. (2021a) showed that observed satellite profiles can be
severely affected by deblending mistakes of the central dominant ICG within
$0.1R_{200}$, where $R_{200}$ is the virial radius of the host dark matter
halo within which the enclosed matter density is 200 times the critical
density of the Universe555In this paper, our definition of the virial mass,
$M_{200}$, is the total mass enclosed within $R_{200}$.. Thus throughout this
paper, when discussing satellite profiles, we should only focus on radial
ranges beyond $0.1R_{200}$.

The above methods have been robustly developed and carefully tested in Wang et
al. (2011); Wang & White (2012); Wang et al. (2014). In fact, similar
methodology has been independently developed and applied to many different
modern surveys and varying types of galaxies at different redshifts in the
last decade or so (e.g. Lares et al., 2011; Guo et al., 2011a; Jiang et al.,
2012; Guo et al., 2012; Kawinwanichakij et al., 2014; Lan et al., 2016; Xu et
al., 2022c; Xu & Jing, 2022; Xu et al., 2022a, b). The very early idea can be
traced back to Lorrimer et al. (1994).

####  3.2.2 Satellite profiles in TNG

For TNG, in principle, we directly know from the simulation the true
population of satellite galaxies bound to each ICG. However, in order to
ensure fair comparisons with observation, we calculate satellite profiles in a
way similar to real observation. We choose the $z$-axis of the simulation box
as the line-of-sight direction, and count the number of companions around our
selected sample of TNG ICGs in projected radial bins perpendicular to the
artificial line of sight, i.e., the $x-y$ plane. We count all the companions
projected around the ICG, i.e., not only those true satellites, but also those
foreground and background galaxies. When cumulating the stellar mass as a
function of projected radius to the central ICG, the stellar mass of each
companion is the total stellar mass of all bound star particles belonging to
it. Only companions more massive than $M_{\ast}=10^{9}\mathrm{M_{\odot}}$ are
used, which is approximately the resolution limit of satellite galaxies in
TNG300-1. Similar to how we obtain the lensing signals, we weight each ICG by
the maximum volume to which one can still observe this galaxy given the flux
limit of SDSS.

We then subtract the foreground and background contamination, which is
calculated by randomly generating $10^{4}$ points within the simulation box
and calculating their companion profiles in the same way as for ICGs.
Moreover, since the simulation box is regular, we can also simply estimate the
level of background contamination from the average volume density of galaxies
in the simulation box and the volume of each projected radial bin. This gives
nearly identical results as when the background contamination is estimated
from companions around random points in the simulation box.

Rigorously, in order to fully mimic real observation, we have to generate
light-cone mock observations for simulated galaxies, so that we have realistic
flux limits, background contamination and galaxy evolution. However, most of
our ICGs are at $z<0.25$, with the most massive ICGs extending to at most
$z\sim 0.3$. The evolution of galaxies at $z<0.3$ is weak. Besides, as having
been shown by Wang & White (2012), directly projecting the simulation box and
light-cone mocks lead to very similar results in the final stacked/averaged
surface brightness, with the noise in the former suppressed lower, and thus we
choose to directly project the simulation box.

Figure 1: Lensing signals for ICGs in real observation (green points) and
TNG300 (gray lines). Each panel corresponds to ICGs in a different stellar
mass range, as indicated by the text, which is the stellar mass range in unit
of $\log_{10}\mathrm{M_{\odot}}$. Observed lensing signals are measured using
Fourier_Quad, with errorbars calculated as the 1-$\sigma$ scatter of 200
jackknife subsamples. The subsampling is applied to the lens catalog. Errors
for TNG300 are calculated as the 1-$\sigma$ scatters of 100 bootstrap
subsamples, which are compared to the line width. Stellar mass corrections
have been included to account for the missing stellar mass in outskirts of
massive galaxies in observation.

###  3.3 Projected stellar mass density profiles for ICGs and their stellar
halos in observation and TNG

####  3.3.1 Observed central profiles

In real observation, the projected stellar mass density profiles of ICGs and
their stellar halos (hereafter central profiles) are directly taken from Wang
et al. (2021a), and here we only briefly introduce the method.

The ICG surface brightness profiles are obtained by stacking galaxy images
from the HSC coadded images (see Section 2.4 for details), which are then
converted to projected stellar mass density profiles based on the PSF
corrected $g-r$ color profiles (Wang et al., 2021a). Explicitly, image cutouts
are at first created for each ICG, with bad pixels masked, such as those hit
by cosmic rays and saturated pixels. For each cutout, companion sources,
including true satellite galaxies, fore/background sources are detected by
Sextractor and masked. Wang et al. (2021a) have tried different source
detection and masking thresholds, to ensure safe masking of companion sources.
Each cutout is resampled to exactly the same world coordinate system and the
same output grids, which we call as input images. For each pixel in the output
image plane, we at first clip corresponding pixels from all input images
(masked pixels are not included). Say all input images, if not masked at this
pixel, contribute $N$ reads in surface brightness. We discard 10% pixels near
the two distribution tails of all input pixel values, in order to avoid
extreme pixel values biasing our results. After the clipping we take the
average of all remaining pixel values. The choice of 10% follows the value
used by D’Souza et al. (2014) and is empirical. Wang et al. (2021a) tested
that choices of 1% to 10% lead to very similar results. In the end, we repeat
exactly the same steps for the images around a sample of random points within
the HSC footprint. The random stacks are estimates of the residual sky
backgrounds, which are subtracted off from the stacks for real ICGs.

####  3.3.2 Central profiles in TNG

For TNG, central profiles are directly calculated by using all bound star
particles belonging to each ICG and those diffuse star particles not belonging
to any other galaxies but are in the same friends-of-friends (FoF) group.
Again we choose the $z$-axis as the line-of-sight direction. Projected radial
bins are arranged in the same way as for satellite profiles.

###  3.4 Stellar mass correction

The stellar masses we use to divide ICGs into different stellar mass bins are
initially from SDSS, which were measured by fitting the stellar population
synthesis model (Blanton & Roweis, 2007) assuming a Chabrier (2003) initial
mass function to the SDSS photometry. However, as have been pointed out in
many previous studies (e.g. Graham et al., 2005; Von Der Linden et al., 2007;
Aihara et al., 2011; Bernardi et al., 2013; He et al., 2013; D’Souza et al.,
2015), the total stellar mass of massive central galaxies could be
underestimated, due to the fact that their outskirts contribute a non-
negligible amount fraction of the total stellar mass but are often
oversubtracted as backgrounds or are under the detection threshold, especially
for shallow surveys. This leads to disagreement with predictions by numerical
simulations, in terms of the stellar mass functions at the massive end (e.g.
He et al., 2013; D’Souza et al., 2015). The sample of ICGs (also LBGs) used in
our study are selected from the SDSS spectroscopic main galaxies, and thus we
need to correct for the missing part in their total stellar mass. This
correction is im,portant to ensure fair comparisons with TNG300, because for
TNG300 galaxies, we use their total stellar mass from the database, without
any aperture cuts.

To achieve the corrections, we use the data provided by the authors of He et
al. (2013). Explicitly, He et al. (2013) reanalyzed the images of a large
sample of red elliptical galaxies from SDSS (e.g. Heymans et al., 2012), and
the photometry of these galaxies down to different sky limits is compared with
the original SDSS photometry, which can be used as the reference for our
stellar mass correction. We use the difference between their remeasured
magnitudes and the original SDSS Petrosian magnitudes ($\Delta m_{r}$) for the
correction. We convert $\Delta m_{r}$ to the corresponding correction in log
stellar mass ($\Delta\log_{10}M_{\ast}$), simply through Equation 14.

$\log_{10}M_{\mathrm{new}}=\log_{10}M_{\mathrm{old}}+\frac{\Delta
m_{r}}{2.5}.$ (14)

For each massive ICG or LBG in our analysis, we match it to the closest one in
$r$-band absolute magnitude from the data of He et al. (2013), and obtain the
newly corrected stellar mass according to Equation 14. We only apply the
correction for galaxies of $r$-band absolute magnitude brighter than -22.5,
since smaller galaxies are expected to suffer less from such underestimates in
their stellar mass and He et al. (2013) only considers galaxies in this
magnitude range.

A comparison of the lensing signals before and after the correction in stellar
mass, is provided in Appendix B. For the two most massive bins in stellar
mass, the difference before and after the correction is prominent. Without the
correction, the values of stellar mass are smaller than the true stellar mass,
and thus given a stellar mass bin, we are actually selecting more massive
galaxies to fit into this bin. Moreover, at the massive end, a small increase
in stellar mass would correspond to a more significant increase in host halo
mass. As a result, the lensing signals are significantly higher in amplitudes
in the two most massive bins. After the stellar mass correction, the
amplitudes are significantly lower, bringing in good agreement with the
prediction by IllustrisTNG300. Our results thus show that proper corrections
to the stellar mass of massive galaxies in observation are crucial for fair
comparisons with theoretical predictions.

In our analysis, we apply the stellar mass correction above when calculating
the lensing signals and the satellite profiles. The central profiles are
directly taken from Wang et al. (2021a) without including such corrections in
stellar mass. More explicitly, the lensing and satellite profiles are binned
by stellar masses after the correction. The central profiles are binned
according to the original stellar mass without such corrections. However, we
will show that our conclusions based on the central profiles in this paper are
not affected by the stellar mass correction.

Figure 2: Lensing signals for observed blue (blue dots) and red (red dots)
ICGs, and TNG300 blue (cyan curves) and red (magenta curves) ICGs. Each panel
corresponds to a different bin in stellar mass, with the text in each panel
showing the stellar mass range in unit of $\log_{10}\mathrm{M_{\odot}}$.
Errorbars for real observation are calculated as the 1-$\sigma$ scatters of
200 jackknife subsamples applied on the lens catalog. Errors for TNG300 are
calculated as the 1-$\sigma$ scatters of 100 bootstrap subsamples, which are
mostly comparable to the line width. Stellar mass corrections have been
included to account for the missing stellar mass in outskirts of massive
galaxies in observation.

##  4 Results

In this section, we present our measured lensing signals, satellite and
central profiles around our sample of ICGs. The measurements are shown in
parallel for both the real observation and predictions by TNG300, for direct
comparisons. We have also confirmed that the trends are similar though noisier
in TNG100. Part of the TNG100 predictions will be shown in Appendix C. In the
end, we investigate the fraction of total stellar mass and the fraction of
stellar mass in surviving satellites as a function of projected radius.

###  4.1 Lensing signals

Figure 1 shows the lensing signals for real observation and TNG300, centered
on ICGs grouped in five different stellar mass ranges. We scale the projected
distance, $r_{\mathrm{p}}$, by the virial radius, $R_{\mathrm{200}}$, with the
values of $R_{\mathrm{200}}$ for ICGs in different stellar mass bins taken
from Wang et al. (2021a), which are based on ICGs selected from a mock ICG
sample of the Munich semi-analytical model666As we have checked, the
$R_{\mathrm{200}}$ values for different stellar mass bins based on the Guo et
al. (2011b) model and based on the TNG simulations can be slightly different
from each other, which may also differ from the actual $R_{\mathrm{200}}$
values of real galaxies, though in observation $R_{\mathrm{200}}$ cannot be
measured directly. However, this does not affect our comparisons, because we
scale observational results and TNG predictions by the same
$R_{\mathrm{200}}$. 777In Section 2.2, we use the virial radius predicted by
the abundance matching formula of Guo et al. (2010) as the criterion to select
ICGs. For the selection, each galaxy can have different virial radius based on
their stellar mass. From now on, the virial radius, $R_{200}$, used to scale
$r_{\mathrm{p}}$ is from mock ICGs of the Guo et al. (2011b) model, so not
exactly the same as those calculated from the abundance matching formula of
Guo et al. (2010). We make such choices in order to be consistent with Wang et
al. (2021b). For all ICGs in the same stellar mass bin, the $R_{200}$ value
adopted to scale $r_{p}$ is the same, i.e., the mean $R_{200}$ for all mock
ICGs in the bin. Also note although the explicit choice of $R_{200}$ can
affect our ICG sample selection, this does not affect our results, because we
select ICGs in TNG and real observation in exactly the same way, i.e., both
are based on the same Guo et al. (2010) abundance matching $R_{200}$.(Guo et
al., 2011b) (see Wang & White (2012) for details). We scale $r_{\mathrm{p}}$
by $R_{200}$ in order to investigate the density profiles over halo scales for
different mass bins, and this will be the convention when we report our
measurements in this paper. The $R_{200}$ values used to scale
$r_{\mathrm{p}}$ are provided in Table 1 for each stellar mass bin.

We find good agreements between observed and TNG lensing signals. For the
observational lensing signal calculation, a correction in stellar mass has
been applied to massive galaxies (See Section 3.4 and Appendix B for details).
This correction significantly affects the signals in the two most massive
bins, which improves the agreement between the real observation and TNG300 at
the massive end. Without this correction in stellar mass, the observed lensing
profiles are significantly higher than the simulation predictions in the two
most massive panels (See figure 10).

In a previous study, Renneby et al. (2020) calculated the lensing signals
based on the Munich semi-analytical model and TNG300. Good agreements were
reported between the real signals and simulation predictions. Wang et al.
(2016) reached similar results when comparing predictions by different Munich
semi-analytical models (e.g. Guo et al., 2011b; Henriques et al., 2015).
However, the simulation predictions are lower than the observed lensing
signals for the most massive bin. We think this is mainly due to the missing
stellar mass in outskirts of our ICGs selected from SDSS, which were not
properly corrected in Wang et al. (2016).

In Figure 2, we present our lensing signals for red and blue ICGs separately,
and for both real observation and TNG300. Note as having been mentioned in
Section 2.2, we determine red and blue galaxies separately for real
observation and TNG, based on the $g-r$ color-mass diagrams. The division
slightly depends on stellar mass. In general, we find an overall good
agreement between both lensing signals and for red and blue ICGs. However, in
the outskirts of the three most massive panels, TNG predicted that signals
around blue ICGs are higher in amplitudes than those around real observed blue
ICGs.

In previous studies based on numerical simulations or based on direct weak
lensing measurements, it was reported that red galaxies are hosted by more
massive dark matter halos than blue ones at fixed stellar mass (e.g. Peng et
al., 2012; Wang & White, 2012; Mandelbaum et al., 2016; Man et al., 2019).
This is also seen in our observed lensing signals. For example, in the
$10.8<\log_{10}<11.1$ panel, the red dots are above the blue ones by $\sim
3\sigma$ beyond $0.3R_{200}$. In the other two more massive panels there are
also indications that the red dots are above blue ones at
$r_{\mathrm{p}}>0.3R_{200}$, but with a lower significance with respect to the
errors ($\sim 1.8$ and $2\sigma$ between $0.3R_{200}$ and $R_{200}$ for the
most and second massive panels).

For TNG, we do not see prominent differences between the lensing signals
around red and blue ICGs. The magenta curves seem to be only slightly above
the cyan curves at $10.2<\log_{10}M_{\ast}/\mathrm{M_{\odot}}<11.1$. The trend
remains the same if we use true halo central galaxies in TNG300. Besides, we
have explicitly checked that TNG100 predicts very similar results. The
differences between the magenta and cyan curves for TNG results are mostly not
significant. Compared with the errors, TNG red and blue ICGs are hosted by
dark matter halos with more similar masses.

Theoretically, it is expected that red galaxies are hosted by more massive
halos than blue galaxies with the same stellar mass (e.g. Peng et al., 2012;
Wang & White, 2012). This is because the star formation activities are
quenched earlier for red galaxies, but their host dark matter halos keep
growing by accretion. Thus at fixed host halo mass, blue galaxies are more
massive in stellar mass. In other words, if the stellar mass is the same, red
galaxies are hosted by more massive dark matter halos. This is seen in our
observed lensing signals, but TNG does not show prominent difference between
the lensing signals around red and blue ICGs.

In the following subsections, we move on to compare the satellite and central
profiles between real observation and TNG predictions. We will provide more
discussions on the difference in host halo mass of red and blue ICGs in
Section 5.2, by looking at the halo and stellar mass growth histories.

Figure 3: Satellite profiles for observed blue (blue dots) and red (red dots)
ICGs, and TNG300 blue (cyan curves) and red (magenta curves) ICGs. Only
satellites more massive than $10^{9}\mathrm{M_{\odot}}$ are included. Each
panel corresponds to a different stellar mass bin, with the text in each panel
showing the stellar mass range in unit of $\log_{10}\mathrm{M_{\odot}}$.
Shaded areas indicate the errors for TNG300. All errors are calculated as
1-$\sigma$ scatters of 100 bootstrap subsamples. Note that the observed
satellite profiles within 0.1$R_{200}$ are significantly affected by source
deblending mistakes with the central dominant galaxy, and hence are not
robust. This is marked by the black dashed vertical lines. Beyond
0.1$R_{200}$, there are more satellites around red ICGs in the three most
massive panels, but the trend is opposite in TNG300. Stellar mass corrections
have been included to account for the missing stellar mass in outskirts of
massive galaxies in observation. Figure 4: Central profiles for SDSS blue
(blue dots) and red (red dots) ICGs, and TNG300 blue (cyan curves) and red
(magenta curves) ICGs. Each panel represents the results for a different
stellar mass range, and it is indicated by the text in each panel, in unit of
$\log_{10}\mathrm{M_{\odot}}$. All errors are calculated as 1-$\sigma$
scatters of 100 bootstrap subsamples, but they are mostly comparable to the
symbol size and line width, hence not prominently shown.

###  4.2 Satellite profiles

Figure 3 shows the satellite profiles (projected stellar mass density profiles
of satellites), for red and blue ICGs separately and for both real observation
and TNG300. Although there is a good overall agreement between observation and
TNG300, there is a significant disagreement when we look at the profiles
around red and blue ICGs separately. In real observation, the satellite
profiles are higher in amplitude beyond $0.1R_{200}$ around red ICGs than
around the blue ones888Within $0.1R_{200}$, the observational satellite
profiles are strongly affected by deblending mistakes and are thus not robust
(Wang et al., 2021a). In fact, part of the star-forming regions and spiral
arms of blue ICGs are mistakenly deblended as satellites within 0.1$R_{200}$,
making the signals artificially higher around blue ICGs at such scales. for
the three most massive panels. This is in good agreement with the previous
studies of Wang & White (2012) and Wang et al. (2021a). However, TNG300 shows
a different trend, in terms that blue ICGs have satellite profiles higher in
amplitudes than the red ones, except for the lowest mass panel.

As having been reported by a few previous studies, the total luminosity or
total stellar mass in satellites are good proxies to the host halo mass (e.g.
Wang & White, 2012; Wang et al., 2021a; Tinker et al., 2021), and thus we
expect the total stellar mass in satellites around red and blue ICGs show
similar trends as lensing signals. As we have shown in the previous section,
there are indications that red ICGs are hosted by more massive dark matter
halos than blue ICGs in real observation, and thus the observed trend in
satellite profiles is roughly consistent with what we see in the observed
lensing profiles.

However, there are more stellar mass in satellites around blue ICGs in TNG300,
i.e., cyan curves are above magenta ones in Figure 3, whereas we did not see
the cyan curves being prominently above the magenta ones in Figure 2. Note, to
ensure the completeness of satellites in TNG300, the satellite profiles are
based on objects which are more massive than $10^{9}\mathrm{M_{\odot}}$, which
are relatively massive objects, but this does not affect the comparison with
observation because we have adopted the same satellite mass threshold for both
observation and TNG300.

We explicitly repeat our calculation for TNG100, and the results are shown in
Appendix C. TNG100 produces results which are noisier due to the smaller
simulation box size, but the trend remains very similar, i.e., there are more
stellar mass in satellites around blue ICGs. Besides, if using only true halo
central galaxies or true satellites bound to central galaxies to calculate
their satellite profiles, the trend stays very similar as well. So the
inconsistency with observation cannot be trivially explained by the resolution
limit or by satellite contamination.

We have also tried different $g-r$ color divisions when separating ICGs into
red and blue populations, such as shifting the division slightly, moving up
and down, finding the division according to the dip in the $g-r$ color
distribution of each stellar mass bin, and we end up with very similar trends
in TNG satellite profiles.

We postpone more detailed discussions on possible explanations to the
discrepancies in Section 5. To gain more evidences, we move on to investigate
the central profiles in the next subsection. Moreover, we will investigate in
Section 4.4 the fraction of total stellar mass versus total mass, as a
function of the projected radius.

###  4.3 Central profiles

Figure 4 presents the central profiles (projected stellar mass density
profiles of ICGs and their extended stellar halos), for both real observation
and TNG300. The observed central profiles are directly taken from the
measurements of Wang et al. (2021a). TNG300 central profiles exhibit a good
overall agreement with real observations. Nevertheless, we observe the TNG300
profiles are more extended, which is particularly relevant for the radial
range of $r_{\mathrm{p}}\approx 0.1-0.3R_{200}$ and for blue ICGs more massive
than $10^{10.8}\mathrm{M_{\odot}}$. Similar trend also exists in the
$10.5<\log_{10}M_{\ast}/\mathrm{M_{\odot}}<10.8$ panel, but not as
significant. Note that for the observed central profiles, corrections to the
stellar mass (see Section 3.4) are not applied. However, this cannot explain
the discrepancy. In fact, as we have mentioned, without the correction to the
observed stellar mass, we end up selecting more massive galaxies with more
extended stellar halos in each bin and in observation. Hence if including
corrections to the stellar mass, we would select less massive galaxies in
observation, and since less massive galaxies have less extended stellar halos,
the discrepancy between real observation and TNG300 predictions would be even
larger. We have also explicitly repeated our measurements with TNG100, and
have found similar disagreement with real observation.

The more extended central profiles of TNG300 galaxies might be partly related
to the fact that the sizes of TNG galaxies are larger than real observed
galaxies. It was shown by Pillepich et al. (2018a) that the average sizes of
TNG galaxies are larger than those of real galaxies. The mass-size relation
predicted by TNG cannot even touch the 1-$\sigma$ error boundaries of SDSS
galaxies. Part of the large discrepancy could be due to the fact that the
outskirts of massive galaxies might fall below the detection limit of the
shallower SDSS survey. Moreover, the SDSS mass-size relation measured by Shen
et al. (2003) is based on the circularized radius, which is on average
$\sim$1.4 times smaller than the actual major axis length (Furlong et al.,
2017). Nevertheless, the fact that hydrodynamical simulations such as TNG and
EAGLE tend to predict larger galaxy sizes than observation seems to be a
universal problem (e.g. Snyder et al., 2019; Yang et al., 2021; De Graaff et
al., 2022; Valluri et al., 2022).

In addition to the connection to the global mass-size relation, our results
show that the discrepancy of the TNG predicted central profiles from real
observed galaxies is mainly in outskirts, where the stellar material is
dominated by accreted stars from satellite galaxies (e.g. Rodriguez-Gomez et
al., 2016). This is in good agreement with Ardila et al. (2021), in which a
comparison of the most massive galaxies were performed between TNG100 and the
measurements based on the Hyper Suprime-Cam (HSC) imaging data at $z\sim 0.4$.
It was reported that the outer stellar halos of observed massive galaxies
hosted by halos with virial masses of $\sim 10^{14}\mathrm{M_{\odot}}$ are
less extended than TNG100 predictions.

The story is similar for very nearby Milky-Way-mass galaxies. Merritt et al.
(2020) reported a discrepancy between the surface brightness profiles of
Milky-Way-like disk galaxies from the Dragonfly Nearby Galaxies Survey (DNGS)
and TNG100 counterparts. The TNG100 galaxies were matched in stellar mass to
those of DNGS, and surface brightness calculations showed that the DNGS
galaxies have "missing" stellar mass in the outskirts. DNGS galaxies have
significantly lower surface brightness profiles beyond $20\mathrm{kpc}$.
However, as having been pointed out by Merritt et al. (2020), their sample
size of nearby disk galaxies is relatively small, so it is likely that the
results could have been affected by statistical fluctuations. Nevertheless, it
seems that at $z=0$, $z\sim 0.1$, and $z\sim 0.4$, evidences all exist to
prove that TNG tends to predict more extended outer stellar halos than real
galaxies. Compared with the smaller sample of galaxies in the local Universe,
the sample of ICGs is large enough to ensure good statistics, and Wang et al.
(2021a) have performed intensive tests on background subtraction and companion
source masking. In particular, for blue ICGs in the three massive panels, the
disagreement between TNG300 prediction and the Wang et al. (2021a) measurement
is significant compared with the small errors in Figure 4.

In addition to the tension we discussed above, we can observe in real data
that the outer stellar halos are more extended around red ICGs, but TNG300
does not show such a trend. There is no significant difference between the
outer central profiles around red and blue ICGs in TNG. We have checked the
same is true in TNG100.

###  4.4 Total stellar mass versus total mass

Figure 5: Ratios between the projected density profiles of total stellar mass
(central $+$ satellites) versus total mass. Red and blue dots with errors are
centered on observed red and blue ICGs, while magenta and cyan curves are
around red and blue ICGs in TNG, respectively. Shaded areas indicate the
errors for TNG300. The stellar mass range is indicated in each panel, in unit
of $\log_{10}\mathrm{M_{\odot}}$. Errors are propagated from the corresponding
central/satellite profiles (Figures 3 and 4) and lensing signals (Figure 2).
For TNG300 predictions, the projected density profiles of total mass are
estimated from the radial distribution of all types of particles.  Figure 6:
Similar to Figure 5, but showing the ratios between the projected stellar mass
density profiles of satellites versus total mass. Red and blue dots with
errors are centered on observed red and blue ICGs, while magenta and cyan
curves are around red and blue ICGs in TNG, respectively. Shaded areas
indicate the errors for TNG300. Errors are propagated from the corresponding
satellite profiles (Figure 3) and lensing signals (Figure 2). The black dashed
lines mark 0.1$R_{200}$ within which the deblending errors are significant,
and therefore the observed satellite profiles are not robust.

In this section we investigate how the fraction of total stellar mass versus
total mass changes as a function of the projected distance to the galaxy
center, for the total stellar mass in central ICGs $+$ their stellar halos and
satellites. We also check whether the radial distribution of satellite
galaxies trace the underlying dark matter.

In order to obtain the projected density profiles for the total mass
distribution from the observed lensing signals in Figure 2, we fit the
following projected NFW profile (Wright & Brainerd, 2000) to the observed
lensing signals for red and blue ICGs

$X_{ij}=\left\\{\begin{aligned}
&\frac{r_{\mathrm{s}}\delta_{\mathrm{c}}\rho_{\mathrm{c}}}{\Sigma_{\mathrm{c}}}g_{<}(x)&(x<1)\\\
&\frac{r_{\mathrm{s}}\delta_{\mathrm{c}}\rho_{\mathrm{c}}}{\Sigma_{\mathrm{c}}}\left[\frac{10}{3}+4\mathrm{ln}\left(\frac{1}{2}\right)\right]&(x=1)\\\
&\frac{r_{\mathrm{s}}\delta_{\mathrm{c}}\rho_{\mathrm{c}}}{\Sigma_{\mathrm{c}}}g_{>}(x)&(x>1)\end{aligned}\right.$
(15)

,

where $r_{\mathrm{s}}=R_{200}/c$ is the scale radius, and
$x=r_{\mathrm{p}}/r_{\mathrm{s}}$ is the projected radius to galaxy center
divided by the scale radius. The functions of $g_{<}(x)$ and $g_{>}(x)$ take
the following form

$\displaystyle g_{<}(x)=$
$\displaystyle\frac{8\mathrm{arctanh}\sqrt{\frac{1-x}{1+x}}}{x^{2}\sqrt{1-x^{2}}}+\frac{4}{x^{2}}\mathrm{ln}\left(\frac{x}{2}\right)$
(16)
$\displaystyle-\frac{2}{(x^{2}-1)}+\frac{4\mathrm{arctanh}\sqrt{\frac{1-x}{1+x}}}{(x^{2}-1)(1-x^{2})^{1/2}},$
$\displaystyle g_{>}(x)=$
$\displaystyle\frac{8\mathrm{arctan}\sqrt{\frac{x-1}{1+x}}}{x^{2}\sqrt{x^{2}-1}}+\frac{4}{x^{2}}\mathrm{ln}\left(\frac{x}{2}\right)$
(17)
$\displaystyle-\frac{2}{(x^{2}-1)}+\frac{4\mathrm{arctan}\sqrt{\frac{x-1}{1+x}}}{(x^{2}-1)^{3/2}}.$

For ICGs in TNG300, the projected density profiles for the total mass
distribution are directly calculated from the density profiles of all
different types of particles, including star and dark matter particles, gas
cells, and black holes.

We show in Figure 5 the fractions of total stellar mass out of total mass, for
both real observation and TNG300 predictions and as a function of the
projected radius scaled by $R_{200}$. We stop at an inner radius of
$0.03R_{200}$ of the $x$-axis, because for real observation we do not have
weak lensing signals measured on such small scales in Figure 2, due to the
lack of enough lens-source pairs at very small radii. The total stellar mass
is estimated from the contribution of both satellite and central profiles, and
thus it includes the contribution from surviving satellites, central ICGs and
their extended stellar halos formed by disrupted satellites. Note although the
satellite profiles are affected by deblending mistakes within $0.1R_{200}$,
the total stellar mass is dominated by central ICGs within $0.1R_{200}$, and
thus is not affected by deblending mistakes. Besides, for satellites, we only
calculate the stellar mass in those which are more massive than
$10^{9}\mathrm{M_{\odot}}$, so in principle we should take our fractions of
total stellar mass as lower limits. However, as long as the power-law index of
the satellite stellar mass function does not rise steeper than $-2$ at the
low-mass end, we expect the total stellar mass in satellites to be dominated
by the most massive/luminous satellites (e.g. Lan et al., 2016; Wang et al.,
2021b), so we do not expect our fraction to be too different from that of
using all satellites.

In the three most massive panels of Figure 5, the TNG300 curves are higher
than real observations around $r_{\mathrm{p}}\sim 0.15R_{200}$, indicating
there are slightly more total stellar mass in TNG300 than real galaxies at
$r_{\mathrm{p}}\sim 0.15R_{200}$ around ICGs more massive than
$\log_{10}M_{\ast}/\mathrm{M_{\odot}}\sim 10.8$. This discrepancy is mainly
due to the more extended central profiles in TNG compared with real observed
galaxies, as the readers can see from Figure 4 that the TNG central profiles
show the most prominent differences from observation at $\sim 0.15R_{200}$ and
in the three most massive panels.

Besides, we can see in Figure 5 that cyan curves are above magenta ones at
almost every radii and in every panels. In outskirts close to the virial
radius boundary, the total stellar mass is contributed more by those from
surviving satellites, and thus the trend is related to the fact that we see
more stellar mass in satellites around blue ICGs in TNG (Figure 3). In more
central regions, this is mainly because there are slightly more total mass
around red ICGs than blue ones in TNG (see the TNG300 lensing signals in
Figure 2), which increases the denominator and brings down the fraction in
total stellar mass. For the most massive panel, the central profile around
blue ICGs are slightly higher, though this is not very prominently seen in
Figure 4. However, we do not see the same trend in real observation. Blue dots
are above red dots in only central regions, and in fact red dots are mostly
above blue ones at $r_{\mathrm{p}}>0.1R_{200}$, except for the lowest mass
panel.

In outskirts of Figure 5, the results are dominated by surviving satellites,
and the fractions already drop below 1%. It seems observed satellites have
slightly more extended radial distributions, in terms that the outer most red
and blue dots are above the magenta and cyan curves in the three most massive
panels. In the least massive panel, the three red/blue dots at large radius
are significantly above the magenta/cyan curves, but this could be due to the
reason that the TNG300 predicted satellite profiles are too noisy in this
least massive panel (see Figure 3).

To more clearly investigate and separate the contribution from satellites, we
now show in Figure 6 the ratios between the stellar mass in surviving
satellites versus the total mass as a function of projected radius to the
centers of ICGs. Despite the discrepancy in the few outer most points between
observation and TNG, the radial distribution of satellites mildly trace the
underlying dark matter in both real observation and simulation. Note within
$0.1R_{200}$ (black dashed vertical lines), the observed satellite profiles
are significantly affected by deblending mistakes, hence showing some
artificial increases.

Theoretically, we expect the total mass in accreted subhalos and the total
accreted dark matter contributing to the growth of host dark matter halos
(including the smoothly accreted part not in any clumpy structures) to be
proportional to each other. However, subhalos undergo significant tidal
disruptions after infall, so the total mass in surviving subhalos is not
really expected to trace well the underlying dark matter. Instead, the amount
of stellar mass bound to surviving satellite galaxies in the very center of
infalling subhalos are expected to be less significantly stripped, except for
those satellites fell in very early and are close to the peri-centers. It is
thus expected that the total stellar mass in satellites should approximately
trace the spatial distribution of dark matter (e.g. Han et al., 2016; Wang et
al., 2018; McDonough & Brainerd, 2022; Dong et al., 2022). This is proved in
our measurements of real data, and also approximately reproduced by TNG.

However, due to the fact that in TNG300 there are more stellar mass in
satellites around blue ICGs (Figure 3), the cyan curves in Figure 6 are
significantly higher than magenta curves. This is not seen in real
observation, because red ICGs are hosted by more massive dark matter halos
while at the same time have more stellar mass in surviving satellites, and
thus the trend is not obviously monotonic between red and blue.

So far, we have raised a few detailed tensions between TNG300 and observation,
when looking at host dark matter halos, population of satellites and outer
stellar halos around red and blue ICGs separately. Similar tensions also exist
between TNG100 and observation. In the next section, we move on to discuss
possible explanations to these tensions.

##  5 Discussions

In the previous section, we have shown evidences that in real observation red
ICGs are hosted by more massive dark matter halos, have more stellar mass in
their satellites and have more extended outer stellar halos than blue ICGs.
TNG300 predicts in general a good agreement with real observation, but
delicate tensions exist when we look at the red and blue ICGs separately. As
we have checked, similar tensions still exist with TNG100.

First of all, we see that red and blue ICGs in TNG are hosted by more similar
dark matter halos. Secondly, blue ICGs have more stellar mass in their
satellites than red ICGs in TNG. Moreover, the outer stellar halos are similar
between red and blue ICGs in TNG, but overall the outer stellar halos are more
extended than real galaxies. This is particularly true for blue ICGs more
massive than $10^{10.8}\mathrm{M_{\odot}}$. Due to the more extended outer
stellar halos for TNG galaxies, the total stellar mass versus total mass
fractions do not show exactly the same radial distributions as real galaxies.

###  5.1 The tension between observed and TNG predicted central/satellite
profiles

Based on the fact that TNG galaxies have more extended outer stellar halos, we
suggest that satellite galaxies in TNG300 and TNG100 might be disrupted more
efficiently than real observed galaxies, which forms the heavier outer stellar
halos. This is because the outer stellar halos are dominated by stripped stars
from infalling satellites (e.g. Rodriguez-Gomez et al., 2016).

Moreover, faster satellite disruption in TNG may also weaken the correlation
between the total stellar mass in satellites and the host halo mass on a
secondary order. Red galaxies are in relatively more over-dense regions than
blue ones. If in TNG the satellites are accreted on average earlier around red
ICGs than blue ones at fixed stellar mass, the satellites around red ICGs
could have been disrupted more due to their earlier infall. Besides, we have
seen in Figure 2 that red ICGs are hosted by slightly larger dark matter halos
at $10.2<\log_{10}M_{\ast}/\mathrm{M_{\odot}}<11.1$ than blue ICGs in TNG, so
the tidal stripping can be slightly stronger for satellites around these red
ICGs in TNG, due to the more massive host halos. If the tidal disruption in
TNG happens more efficiently than real Universe, maybe this would cause less
stellar mass remained in bound surviving satellites around red ICGs than their
blue counterparts, resulting in an opposite trend compared with real
observation.

Figure 7: The fraction of ex-situ formed stars from disrupted satellites
versus the total accreted stellar mass (disrupted $+$ surviving satellites),
for red and blue ICGs in five stellar mass bins.

We first checked the infall time distribution of satellites around red and
blue ICGs, and indeed identified earlier infall of satellites around red ICGs
(see Appendix D). Moreover, the amount of stripped stellar mass from disrupted
satellites can be directly known for red and blue ICGs in TNG. Here we use the
stellar assembly catalog by (Rodriguez-Gomez et al., 2016, 2017), which
provides the information of ex-situ formed stellar mass for galaxies in TNG.
Figure 7 shows the results for red and blue ICGs in our five stellar mass
bins. For the quantity of the $y$-axis, the nominator is the ex-situ formed
stellar mass of red and blue ICGs, which includes not only those stripped
stars from surviving satellites at $z=0$, but also those entirely disrupted
satellites999In principle, we should only include the stellar mass stripped
from $z=0$ surviving satellites, but Figure 7 here already demonstrates a
useful trend.. The denominator is the ex-situ formed stellar mass plus the
total stellar mass in all101010We only count satellites which are more massive
than $10^{9}\mathrm{M_{\odot}}$ in the satellite profiles, but since massive
satellites dominate over stellar mass, if we exclude the stellar mass smaller
satellites, the trend of Figure 7 remains very similar. $z=0$ surviving
satellites within $R_{200}$. It is prominent that the fraction of stripped
stars is higher for satellites around red ICGs, except for the least massive
bin, so indeed satellites around red ICGs are disrupted faster than those
around blue ICGs in TNG for the four massive bins. For the least massive bin
the cyan triangle is above the magenta one, while we did not see the cyan
curve to be significantly above the magenta one in the corresponding panel of
Figure 3 either. We think this is because tidal effects are weaker around less
massive ICGs, so satellites around red ICGs are no longer strongly stripped in
the least massive bin. Besides, the satellite profile is already quite noisy
in the least massive panel of Figure 3. At least, Figure 7 and Figure 3 are
fully consistent. Thus our explanation that faster satellite disruption in TNG
may turn red ICGs to have less stellar mass in their surviving satellites
seems plausible.

Note, however, in addition to the fast disruption of satellites around red
ICGs, cyan curves in Figure 3 are higher in amplitudes than blue dots in the
two most massive bins, indicating satellites around blue ICGs in TNG are NOT
stripped enough. However, given the fact that the TNG central profiles are
more extended in ourskirts, especially around blue ICG, it seems unlikely that
satellites around blue ICGs are not stripped enough, because this should lead
to less extended outer stellar halos around blue ICGs. In fact, Merritt et al.
(2020) have shown in their results that artificially delaying the disruption
of satellite galaxies and reducing the spatial extent of in-situ stellar
populations result in improved matches between the outer surface brightness
profile shapes and stellar halo masses. We think the more efficient satellite
disruption around red ICGs in TNG is a self-consistent explanation.

The faster satellite disruption in TNG is at least partially related to the
resolution limit. For example, Green et al. (2021) investigated how subhalo
mass functions, number density profiles, and substructure mass fractions are
affected by mass resolution limits. It was shown that the resolution limit
plays a primary role, and lower resolution simulations have more significant
disruptions than those in higher resolution simulations, though the
conclusions are based on N-body simulations. Moreover, to correct for possible
biases introduced by the resolution limit, a rescaling method (Pillepich et
al., 2018b) has been adopted to recalibrate the stellar mass in TNG300, which
is based on the ratio of stellar mass between two higher resolution TNG100
runs and at fixed halo mass. The stellar mass versus halo mass relations and
the stellar mass functions in TNG300 show large differences before and after
the recalibration.

Another possible explanation to why satellite profiles are higher in
amplitudes around blue ICGs might be related to the galactic conformity
phenomenon. In both observation and simulation, people have long been noticed
the so-called galactic conformity phenomenon, that the color distribution of
satellites correlate with the central galaxies, i.e., satellites around blue
centrals are bluer and have stronger star formation rates than satellites
around red centrals (e.g. Weinmann et al., 2006; Yang et al., 2006; Wang &
White, 2012; Kauffmann et al., 2013; Henriques et al., 2017). This might
explain why there are more stellar mass in surviving satellites around blue
ICGs, because more stellar mass could have formed in satellites around blue
ICGs. Nevertheless, if we want to explain this using the galactic conformity
phenomenon, the conformity signal has to be prominently stronger in TNG than
real observation. Based on a previous study of Ayromlou et al. (2022), the
TNG300 conformity signal is in a general good agreement with SDSS and DESI
galaxies. We have repeated the calculation of our conformity signal, based on
the sample of ICGs in TNG and observation, and we did not find the conformity
signal to be significantly stronger in TNG than that of real galaxies. Thus
galactic conformity is unlikely to explain the tension in satellite profiles
between TNG and real observation.

###  5.2 Indications from the mass growth histories of ICGs in TNG

Figure 8: The growth history in stellar mass and $M_{200}$ for TNG300 ICGs.
Red and blue curves show the $M_{200}$ growth histories in bins of
$\log_{10}M_{\ast}$, while magenta and cyan curves show the stellar mass
growth histories in bins of $\log_{10}M_{200}$, for red and blue ICGs,
respectively. Each panel represents the results for a different stellar mass
or $M_{200}$ bin, indicated by the text at the bottom of each panel, which is
in unit of $\log_{10}\mathrm{M_{\odot}}$. Note the chosen range of $M_{200}$
bins approximately corresponds to the common range in $M_{200}$ for red and
blue ICGs in the corresponding stellar mass bin.

As shown in Figure 2, the lensing signals around red and blue ICGs already do
not show exactly the same trend between observation and TNG. As we have
mentioned above, the difference between the lensing signals around red and
blue ICGs is related to the stellar and halo mass growth histories. Thus in
this section we investigate the stellar mass and dark matter growth histories
for red and blue ICGs.

Figure 8 shows the average growth history of $M_{200}$ for red and blue ICGs
in our five stellar mass bins (red and blue curves) and in TNG300. In
addition, we also choose five bins in $M_{200}$, which roughly corresponds to
the common range in $M_{200}$ for red and blue ICGs in the five stellar mass
bins. After selecting red and blue ICGs according to $M_{200}$, the magenta
and cyan curves in each panel of Figure 8 also shows the growth history in
stellar mass for these ICGs in $M_{200}$ bins. The ranges in both log stellar
mass and in log host halo mass are indicated by the text in each panel.

At fixed $M_{200}$, red ICGs indeed stop their growth in stellar mass earlier,
with the magenta curves starting to become flat111111We define the redshift
when the curve starts to become flat, as the time when the change rates in
stellar mass start to be smaller than 25% per Gyr. at around $z\sim 1.6$ for
all stellar mass bins, whereas the cyan curves start to be flat at $z<1.6$,
$z<1.07$, $z<0.58$, $z<0.38$ and $z<0.38$, from the most to the least massive
bins. The stop in the stellar mass growth for blue ICGs happens later than red
ICGs by approximately 1.54, 4.01, 5.51, and 5.51 Gyr on average, from second
most to least massive bins. For the most massive bin, the number of blue ICGs
is not enough, but the difference is small between magenta and cyan curves in
that bin anyway. On the other hand, the growth histories of $M_{200}$ at fixed
stellar mass are more similar between red and blue ICGs. The red curves are at
most slightly above the blue ones at $z<0.8$ and
$10.2<\log_{10}M_{\ast}/\mathrm{M_{\odot}}<11.1$, with a bit larger
differences at $z>1$.

In real observation, Figure 2 above has shown that red ICGs are hosted by more
massive dark matter halos in the
$10.8<\log_{10}M_{\ast}/\mathrm{M_{\odot}}<11.1$ panel, and there are similar
trends in a few other panels, but with low statistical significances. On the
other hand, TNG shows no prominent differences between the lensing signals
around red and blue ICGs. We propose that if red TNG ICGs can stop their
growth in stellar mass earlier, their $z=0$ stellar mass would become smaller,
and we would select more massive red ICGs at $z=0$. As a result, the contrast
in the host halo mass for red and blue ICGs would become larger at fixed
stellar mass. This would bring better agreement with the observed lensing
signals around red and blue ICGs, especially in the
$10.8<\log_{10}M_{\ast}/\mathrm{M_{\odot}}<11.1$ panel.

Based on the mass growth histories for ICGs in TNG300, red ICGs stop the
growth in their stellar mass at $z\sim 1.6$ on average, which correspond to
9.78 Gyrs ago. If our argument above is true, nowadays red central galaxies in
our real Universe with stellar masses greater than
$10^{10.8}\mathrm{M_{\odot}}$ are expected to stop significant growth in their
stellar mass slightly earlier than $\sim$9.78 Gyrs ago, based on the TNG300
model.

##  6 Conclusions

In this study, we measure the weak lensing signals and the projected stellar
mass density profiles of satellite galaxies (satellite profiles) centered on a
sample of isolated central galaxies from SDSS DR7 and at $z\sim 0.1$. We also
take measurements of the projected stellar mass density profiles for these
ICGs $+$ their extended stellar halos (central profiles) from a previous study
of Wang et al. (2021a). The shear measurements used for weak lensing signal
calculations are based on the Fourier_Quad Method (Zhang et al., 2019) and the
eighth data release of the Dark Energy Camera Legacy Survey, while the
satellite and central profiles are based on photometric sources and galaxy
images of the Hyper Suprime-Cam (HSC) survey (pDR3). All the observational
measurements are compared in every detail with predictions by IllustrisTNG300,
with validations from IllustrisTNG100 to confirm that our comparisons are not
affected by the resolution.

Excellent agreement is found between the lensing signals from real observation
and TNG for all ICGs with different stellar mass. In particular, the good
agreement at the massive end relies on a correction to the stellar mass of
observed massive galaxies, to account for the missing stellar mass in
outskirts due to an oversubtraction in the sky background (He et al., 2013).
Otherwise the observed lensing signals are significantly higher than the TNG
predictions.

However, once our sample of ICGs is divided into red and blue populations, we
find more detailed tensions between observed and TNG predicted lensing
signals. For galaxies with $10.8<\log_{10}M_{\ast}/\mathrm{M_{\odot}}<11.1$,
red ICGs are hosted by more massive dark matter halos than their blue
counterparts in real observation, and similar trends exist in other stellar
mass bins, though with low statistical significances. On the other hand, red
and blue ICGs in TNG are hosted by more similar dark matter halos.

Theoretically, the difference in host halo mass for galaxies with the same
stellar mass but different color is due to the early quench of star formation
in red galaxies, but whose host halos keep growing through accretion, i.e.,
the stellar mass of red galaxies are smaller at fixed halo mass. We thus
propose that if TNG galaxies can stop their growth in stellar mass slightly
earlier, the agreement between observed and TNG predicted lensing signals
would be better. If our proposal is correct and based on the stellar mass
growth histories of TNG300 galaxies, we predict that nowadays red central
galaxies in our real Universe with $M_{\ast}>10^{10.8}\mathrm{M_{\odot}}$ are
expected to stop significant growth in their stellar mass slightly earlier
than $\sim$9.78 Gyrs ago.

We find that the satellite profiles around red and blue ICGs do not show a
consistent trend with the lensing signals in TNG300. On the contrary, the
satellite profiles are higher in amplitudes around blue ICGs than red ones at
fixed stellar mass, when real observed satellite profiles are higher in
amplitudes around red ICGs. The same is true in TNG100, and is true for true
satellites bound to true central galaxies.

The outer stellar halo profiles of galaxies in TNG300 and TNG100, which are
dominated by stripped stars from satellites, are more extended than real
observation, especially for blue ICGs at
$\log_{10}M_{\ast}/\mathrm{M_{\odot}}>10.8$. This is in good agreement with
previous studies based on more massive elliptical galaxies at $z=0.4$ (Merritt
et al., 2020) and based on a smaller sample of Milky-Way-like galaxies in the
local volume (Ardila et al., 2021). Moreover, while we see that red ICGs have
more extended outer stellar halos than blue ones in observation, we do not see
the same trend in TNG, with red and blue ICGs showing similar shapes in their
outer stellar halos.

The more extended outer stellar halo profiles of galaxies in TNG300 and TNG100
than those of real galaxies can be explained by a more efficient satellite
disruption compared with that in the real Universe. In addition, we identified
higher fractions of disrupted stellar mass around red ICGs than around their
blue counterparts. If satellite disruptions in TNG happen more efficiently
than in real observation, satellites around red ICGs would have been disrupted
more than those around real galaxies, which can explain why there are less
stellar mass in surviving satellites around red ICGs, in contrary to real
observation.

The more efficient satellite disruption in TNG and the disagreement with real
observation are at least partially due to the resolution limit, in terms that
lower resolution simulations have more significant satellite disruptions than
those in higher resolution simulations (e.g. Pillepich et al., 2018b; Green et
al., 2021).

In both real observation and TNG, satellites approximately trace the
underlying radial distribution of dark matter beyond $0.1R_{200}$. However,
the fractions of total stellar mass versus total mass do not show exactly the
same radial distribution as real galaxies. The fraction of total stellar mass
is slightly higher around blue ICGs than red ones in TNG at every radius, but
we do not see the same trend in real observation.

##  7 Acknowledgments

This work is supported by NSFC (12022307, 12273021, 11621303, 11890691,
12073017), National Key Basic Research and Development Program of China (No.
2018YFA0404504), 111 project (No. B20019), Shanghai Natural Science Foundation
(No. 19ZR1466800) and the science research grants from the China Manned Space
Project (No. CMS-CSST-2021-A02, No. CMS-CSST-2021-B03, No. CMS-CSST-2021-A01).
We gratefully acknowledge the support of the Key Laboratory for Particle
Physics, Astrophysics and Cosmology, Ministry of Education. We thank the
sponsorship from Yangyang Development Fund. SS acknowledges support from NSFC
grants (Nos.11988101, 12273053), CAS Project for Young Scientists in Basic
Research Grant (No. YSBR-062), the K. C. Wong Foundation, and the science
research grants from the China Manned Space Project (No. CMS-CSST-2021-B03).
QG is supported by NSFC (12033008, 11622325), and the science research grants
from the China Manned Space Project (No. CMS-CSST-2021-A03 and CMS-
CSST-2021-A07). CNH acknowledges support from the NSFC grant (No. 11733002).

The computations in this paper were run on the $\pi$ 2.0 cluster supported by
the Center of High Performance Computing at Shanghai Jiaotong University, and
the Gravity supercomputer of the Astronomy Department, Shanghai Jiaotong
University.

## References

  * Abazajian et al. (2009) Abazajian, K. N., Adelman-McCarthy, J. K., Agüeros, M. A., et al. 2009, ApJS, 182, 543, doi: 10.1088/0067-0049/182/2/543
  * Aihara et al. (2018a) Aihara, H., Arimoto, N., Armstrong, R., et al. 2018a, PASJ, 70, S4, doi: 10.1093/pasj/psx066
  * Aihara et al. (2018b) Aihara, H., Armstrong, R., Bickerton, S., et al. 2018b, PASJ, 70, S8, doi: 10.1093/pasj/psx081
  * Aihara et al. (2011) Aihara, H., Allende Prieto, C., An, D., et al. 2011, ApJS, 193, 29, doi: 10.1088/0067-0049/193/2/29
  * Aihara et al. (2022) Aihara, H., AlSayyad, Y., Ando, M., et al. 2022, PASJ, 74, 247, doi: 10.1093/pasj/psab122
  * Anderson et al. (2015) Anderson, M. E., Gaspari, M., White, S. D. M., Wang, W., & Dai, X. 2015, MNRAS, 449, 3806, doi: 10.1093/mnras/stv437
  * Ardila et al. (2021) Ardila, F., Huang, S., Leauthaud, A., et al. 2021, MNRAS, 500, 432, doi: 10.1093/mnras/staa3215
  * Axelrod et al. (2010) Axelrod, T., Kantor, J., Lupton, R. H., & Pierfederici, F. 2010, in Proc. SPIE, Vol. 7740, Software and Cyberinfrastructure for Astronomy, 774015, doi: 10.1117/12.857297
  * Ayromlou et al. (2022) Ayromlou, M., Kauffmann, G., Anand, A., & White, S. D. M. 2022, arXiv e-prints, arXiv:2207.02218.  https://arxiv.org/abs/2207.02218
  * Berlind & Weinberg (2002) Berlind, A. A., & Weinberg, D. H. 2002, ApJ, 575, 587, doi: 10.1086/341469
  * Bernardi et al. (2013) Bernardi, M., Meert, A., Sheth, R. K., et al. 2013, MNRAS, 436, 697, doi: 10.1093/mnras/stt1607
  * Bertin & Arnouts (1996) Bertin, E., & Arnouts, S. 1996, A&AS, 117, 393, doi: 10.1051/aas:1996164
  * Blanton & Roweis (2007) Blanton, M. R., & Roweis, S. 2007, AJ, 133, 734, doi: 10.1086/510127
  * Blanton et al. (2005) Blanton, M. R., Schlegel, D. J., Strauss, M. A., et al. 2005, AJ, 129, 2562, doi: 10.1086/429803
  * Bosch et al. (2018) Bosch, J., Armstrong, R., Bickerton, S., et al. 2018, PASJ, 70, S5, doi: 10.1093/pasj/psx080
  * Crain et al. (2015) Crain, R. A., Schaye, J., Bower, R. G., et al. 2015, MNRAS, 450, 1937, doi: 10.1093/mnras/stv725
  * Cunha et al. (2009) Cunha, C. E., Lima, M., Oyaizu, H., Frieman, J., & Lin, H. 2009, MNRAS, 396, 2379, doi: 10.1111/j.1365-2966.2009.14908.x
  * Davis et al. (1985) Davis, M., Efstathiou, G., Frenk, C. S., & White, S. D. M. 1985, ApJ, 292, 371, doi: 10.1086/163168
  * De Graaff et al. (2022) De Graaff, A., Trayford, J., Franx, M., et al. 2022, MNRAS, 511, 2544, doi: 10.1093/mnras/stab3510
  * De Lucia & Blaizot (2007) De Lucia, G., & Blaizot, J. 2007, MNRAS, 375, 2, doi: 10.1111/j.1365-2966.2006.11287.x
  * DESI Collaboration et al. (2016) DESI Collaboration, Aghamousa, A., Aguilar, J., et al. 2016, arXiv e-prints, arXiv:1611.00036.  https://arxiv.org/abs/1611.00036
  * Dey et al. (2019) Dey, A., Schlegel, D. J., Lang, D., et al. 2019, AJ, 157, 168, doi: 10.3847/1538-3881/ab089d
  * Dong et al. (2022) Dong, F., Zhao, D., Han, J., et al. 2022, ApJ, 929, 120, doi: 10.3847/1538-4357/ac5aaa
  * D’Souza et al. (2014) D’Souza, R., Kauffman, G., Wang, J., & Vegetti, S. 2014, MNRAS, 443, 1433, doi: 10.1093/mnras/stu1194
  * D’Souza et al. (2015) D’Souza, R., Vegetti, S., & Kauffmann, G. 2015, MNRAS, 454, 4027, doi: 10.1093/mnras/stv2234
  * Fong et al. (2022) Fong, M., Han, J., Zhang, J., et al. 2022, MNRAS, 513, 4754, doi: 10.1093/mnras/stac1263
  * Furlong et al. (2017) Furlong, M., Bower, R. G., Crain, R. A., et al. 2017, MNRAS, 465, 722, doi: 10.1093/mnras/stw2740
  * Furusawa et al. (2018) Furusawa, H., Koike, M., Takata, T., et al. 2018, PASJ, 70, S3, doi: 10.1093/pasj/psx079
  * Genel et al. (2014) Genel, S., Vogelsberger, M., Springel, V., et al. 2014, MNRAS, 445, 175, doi: 10.1093/mnras/stu1654
  * Graham et al. (2005) Graham, A. W., Driver, S. P., Petrosian, V., et al. 2005, AJ, 130, 1535, doi: 10.1086/444475
  * Green et al. (2021) Green, S. B., van den Bosch, F. C., & Jiang, F. 2021, MNRAS, 503, 4075, doi: 10.1093/mnras/stab696
  * Guo et al. (2011a) Guo, Q., Cole, S., Eke, V., & Frenk, C. 2011a, MNRAS, 417, 370, doi: 10.1111/j.1365-2966.2011.19270.x
  * Guo et al. (2012) —. 2012, MNRAS, 427, 428, doi: 10.1111/j.1365-2966.2012.21882.x
  * Guo et al. (2010) Guo, Q., White, S., Li, C., & Boylan-Kolchin, M. 2010, MNRAS, 404, 1111, doi: 10.1111/j.1365-2966.2010.16341.x
  * Guo et al. (2011b) Guo, Q., White, S., Boylan-Kolchin, M., et al. 2011b, MNRAS, 413, 101, doi: 10.1111/j.1365-2966.2010.18114.x
  * Han et al. (2016) Han, J., Cole, S., Frenk, C. S., & Jing, Y. 2016, MNRAS, 457, 1208, doi: 10.1093/mnras/stv2900
  * Han et al. (2015) Han, J., Eke, V. R., Frenk, C. S., et al. 2015, MNRAS, 446, 1356, doi: 10.1093/mnras/stu217810.48550/arXiv.1404.6828
  * He et al. (2013) He, Y. Q., Xia, X. Y., Hao, C. N., et al. 2013, ApJ, 773, 37, doi: 10.1088/0004-637X/773/1/37
  * Henriques et al. (2015) Henriques, B. M. B., White, S. D. M., Thomas, P. A., et al. 2015, MNRAS, 451, 2663, doi: 10.1093/mnras/stv705
  * Henriques et al. (2017) —. 2017, MNRAS, 469, 2626, doi: 10.1093/mnras/stx1010
  * Hernández-Monteagudo et al. (2015) Hernández-Monteagudo, C., Ma, Y.-Z., Kitaura, F. S., et al. 2015, Phys. Rev. Lett., 115, 191301, doi: 10.1103/PhysRevLett.115.191301
  * Heymans et al. (2012) Heymans, C., Van Waerbeke, L., Miller, L., et al. 2012, MNRAS, 427, 146, doi: 10.1111/j.1365-2966.2012.21952.x
  * Jee & Tyson (2011) Jee, M. J., & Tyson, J. A. 2011, PASP, 123, 596, doi: 10.1086/660137
  * Jiang et al. (2012) Jiang, C. Y., Jing, Y. P., & Li, C. 2012, ApJ, 760, 16, doi: 10.1088/0004-637X/760/1/16
  * Jurić et al. (2017) Jurić, M., Kantor, J., Lim, K.-T., et al. 2017, in Astronomical Society of the Pacific Conference Series, Vol. 512, Astronomical Data Analysis Software and Systems XXV, ed. N. P. F. Lorente, K. Shortridge, & R. Wayth, 279.  https://arxiv.org/abs/1512.07914
  * Kauffmann et al. (2013) Kauffmann, G., Li, C., Zhang, W., & Weinmann, S. 2013, MNRAS, 430, 1447, doi: 10.1093/mnras/stt007
  * Kauffmann et al. (1993) Kauffmann, G., White, S. D. M., & Guiderdoni, B. 1993, MNRAS, 264, 201, doi: 10.1093/mnras/264.1.201
  * Kawinwanichakij et al. (2014) Kawinwanichakij, L., Papovich, C., Quadri, R. F., et al. 2014, ApJ, 792, 103, doi: 10.1088/0004-637X/792/2/103
  * Komiyama et al. (2018) Komiyama, Y., Obuchi, Y., Nakaya, H., et al. 2018, PASJ, 70, S2, doi: 10.1093/pasj/psx069
  * Lan et al. (2016) Lan, T.-W., Ménard, B., & Mo, H. 2016, MNRAS, 459, 3998, doi: 10.1093/mnras/stw898
  * Lares et al. (2011) Lares, M., Lambas, D. G., & Domínguez, M. J. 2011, AJ, 142, 13, doi: 10.1088/0004-6256/142/1/13
  * Lorrimer et al. (1994) Lorrimer, S. J., Frenk, C. S., Smith, R. M., White, S. D. M., & Zaritsky, D. 1994, MNRAS, 269, 696, doi: 10.1093/mnras/269.3.696
  * Man et al. (2019) Man, Z.-Y., Peng, Y.-J., Shi, J.-J., et al. 2019, ApJ, 881, 74, doi: 10.3847/1538-4357/ab2ece
  * Mandelbaum et al. (2016) Mandelbaum, R., Wang, W., Zu, Y., et al. 2016, MNRAS, 457, 3200, doi: 10.1093/mnras/stw188
  * Marinacci et al. (2018) Marinacci, F., Vogelsberger, M., Pakmor, R., et al. 2018, MNRAS, 480, 5113, doi: 10.1093/mnras/sty2206
  * McDonough & Brainerd (2022) McDonough, B., & Brainerd, T. G. 2022, ApJ, 933, 161, doi: 10.3847/1538-4357/ac752d
  * Merritt et al. (2020) Merritt, A., Pillepich, A., van Dokkum, P., et al. 2020, MNRAS, 495, 4570, doi: 10.1093/mnras/staa1164
  * Miyazaki et al. (2018) Miyazaki, S., Komiyama, Y., Kawanomoto, S., et al. 2018, PASJ, 70, S1, doi: 10.1093/pasj/psx063
  * Miyazaki et al. (2012) Miyazaki, S., Komiyama, Y., Nakaya, H., et al. 2012, in Society of Photo-Optical Instrumentation Engineers (SPIE) Conference Series, Vol. 8446, Ground-based and Airborne Instrumentation for Astronomy IV, ed. I. S. McLean, S. K. Ramsay, & H. Takami, 84460Z, doi: 10.1117/12.926844
  * Naiman et al. (2018) Naiman, J. P., Pillepich, A., Springel, V., et al. 2018, MNRAS, 477, 1206, doi: 10.1093/mnras/sty618
  * Nelson et al. (2018) Nelson, D., Pillepich, A., Springel, V., et al. 2018, MNRAS, 475, 624, doi: 10.1093/mnras/stx3040
  * Nelson et al. (2019) Nelson, D., Springel, V., Pillepich, A., et al. 2019, Computational Astrophysics and Cosmology, 6, 2, doi: 10.1186/s40668-019-0028-x
  * Peacock & Smith (2000) Peacock, J. A., & Smith, R. E. 2000, MNRAS, 318, 1144, doi: 10.1046/j.1365-8711.2000.03779.x
  * Peng et al. (2012) Peng, Y.-j., Lilly, S. J., Renzini, A., & Carollo, M. 2012, ApJ, 757, 4, doi: 10.1088/0004-637X/757/1/4
  * Pillepich et al. (2018a) Pillepich, A., Springel, V., Nelson, D., et al. 2018a, MNRAS, 473, 4077, doi: 10.1093/mnras/stx2656
  * Pillepich et al. (2018b) Pillepich, A., Nelson, D., Hernquist, L., et al. 2018b, MNRAS, 475, 648, doi: 10.1093/mnras/stx3112
  * Planck Collaboration et al. (2016a) Planck Collaboration, Ade, P. A. R., Aghanim, N., et al. 2016a, A&A, 594, A13, doi: 10.1051/0004-6361/201525830
  * Planck Collaboration et al. (2013) —. 2013, A&A, 557, A52, doi: 10.1051/0004-6361/201220941
  * Planck Collaboration et al. (2014) —. 2014, A&A, 571, A16, doi: 10.1051/0004-6361/201321591
  * Planck Collaboration et al. (2016b) —. 2016b, A&A, 586, A140, doi: 10.1051/0004-6361/201526328
  * Renneby et al. (2020) Renneby, M., Henriques, B. M. B., Hilbert, S., et al. 2020, MNRAS, 498, 5804, doi: 10.1093/mnras/staa2675
  * Rodriguez-Gomez et al. (2016) Rodriguez-Gomez, V., Pillepich, A., Sales, L. V., et al. 2016, MNRAS, 458, 2371, doi: 10.1093/mnras/stw456
  * Rodriguez-Gomez et al. (2017) Rodriguez-Gomez, V., Sales, L. V., Genel, S., et al. 2017, MNRAS, 467, 3083, doi: 10.1093/mnras/stx305
  * Schaye et al. (2015) Schaye, J., Crain, R. A., Bower, R. G., et al. 2015, MNRAS, 446, 521, doi: 10.1093/mnras/stu2058
  * Scoccimarro et al. (2001) Scoccimarro, R., Sheth, R. K., Hui, L., & Jain, B. 2001, ApJ, 546, 20, doi: 10.1086/318261
  * Seljak (2000) Seljak, U. 2000, MNRAS, 318, 203, doi: 10.1046/j.1365-8711.2000.03715.x
  * Sérsic (1963) Sérsic, J. L. 1963, Boletin de la Asociacion Argentina de Astronomia La Plata Argentina, 6, 41
  * Shen et al. (2003) Shen, S., Mo, H. J., White, S. D. M., et al. 2003, MNRAS, 343, 978, doi: 10.1046/j.1365-8711.2003.06740.x
  * Sijacki et al. (2015) Sijacki, D., Vogelsberger, M., Genel, S., et al. 2015, MNRAS, 452, 575, doi: 10.1093/mnras/stv1340
  * Skibba & Sheth (2009) Skibba, R. A., & Sheth, R. K. 2009, MNRAS, 392, 1080, doi: 10.1111/j.1365-2966.2008.14007.x
  * Snyder et al. (2019) Snyder, G. F., Rodriguez-Gomez, V., Lotz, J. M., et al. 2019, MNRAS, 486, 3702, doi: 10.1093/mnras/stz1059
  * Springel et al. (2019) Springel, V., Pakmor, R., & Weinberger, R. 2019, AREPO: Cosmological magnetohydrodynamical moving-mesh simulation code.  http://ascl.net/1909.010
  * Springel et al. (2001a) Springel, V., White, S. D. M., Tormen, G., & Kauffmann, G. 2001a, MNRAS, 328, 726, doi: 10.1046/j.1365-8711.2001.04912.x
  * Springel et al. (2001b) —. 2001b, MNRAS, 328, 726, doi: 10.1046/j.1365-8711.2001.04912.x
  * Springel et al. (2018) Springel, V., Pakmor, R., Pillepich, A., et al. 2018, MNRAS, 475, 676, doi: 10.1093/mnras/stx3304
  * Tinker et al. (2021) Tinker, J. L., Cao, J., Alpaslan, M., et al. 2021, MNRAS, 505, 5370, doi: 10.1093/mnras/stab1576
  * Valluri et al. (2022) Valluri, M., Chabanier, S., Irsic, V., et al. 2022, arXiv e-prints, arXiv:2203.07491.  https://arxiv.org/abs/2203.07491
  * Vogelsberger et al. (2014a) Vogelsberger, M., Genel, S., Springel, V., et al. 2014a, Nature, 509, 177, doi: 10.1038/nature13316
  * Vogelsberger et al. (2014b) —. 2014b, MNRAS, 444, 1518, doi: 10.1093/mnras/stu1536
  * Von Der Linden et al. (2007) Von Der Linden, A., Best, P. N., Kauffmann, G., & White, S. D. M. 2007, MNRAS, 379, 867, doi: 10.1111/j.1365-2966.2007.11940.x
  * Wang et al. (2018) Wang, C., Li, R., Gao, L., et al. 2018, MNRAS, 475, 4020, doi: 10.1093/mnras/sty073
  * Wang et al. (2022) Wang, J., Yang, X., Zhang, J., et al. 2022, ApJ, 936, 161, doi: 10.3847/1538-4357/ac8986
  * Wang et al. (2011) Wang, W., Jing, Y. P., Li, C., Okumura, T., & Han, J. 2011, ApJ, 734, 88, doi: 10.1088/0004-637X/734/2/88
  * Wang et al. (2014) Wang, W., Sales, L. V., Henriques, B. M. B., & White, S. D. M. 2014, MNRAS, 442, 1363, doi: 10.1093/mnras/stu988
  * Wang & White (2012) Wang, W., & White, S. D. M. 2012, MNRAS, 424, 2574, doi: 10.1111/j.1365-2966.2012.21256.x
  * Wang et al. (2016) Wang, W., White, S. D. M., Mandelbaum, R., et al. 2016, MNRAS, 456, 2301, doi: 10.1093/mnras/stv2809
  * Wang et al. (2019) Wang, W., Han, J., Sonnenfeld, A., et al. 2019, MNRAS, 487, 1580, doi: 10.1093/mnras/stz1339
  * Wang et al. (2021a) Wang, W., Li, X., Shi, J., et al. 2021a, ApJ, 919, 25, doi: 10.3847/1538-4357/ac0e38
  * Wang et al. (2021b) Wang, W., Takada, M., Li, X., et al. 2021b, MNRAS, 500, 3776, doi: 10.1093/mnras/staa3495
  * Weinberger et al. (2017) Weinberger, R., Springel, V., Hernquist, L., et al. 2017, MNRAS, 465, 3291, doi: 10.1093/mnras/stw2944
  * Weinmann et al. (2006) Weinmann, S. M., van den Bosch, F. C., Yang, X., & Mo, H. J. 2006, MNRAS, 366, 2, doi: 10.1111/j.1365-2966.2005.09865.x
  * Westra et al. (2010) Westra, E., Geller, M. J., Kurtz, M. J., Fabricant, D. G., & Dell’Antonio, I. 2010, PASP, 122, 1258, doi: 10.1086/657452
  * White & Frenk (1991) White, S. D. M., & Frenk, C. S. 1991, ApJ, 379, 52, doi: 10.1086/170483
  * White & Rees (1978) White, S. D. M., & Rees, M. J. 1978, MNRAS, 183, 341, doi: 10.1093/mnras/183.3.341
  * Wright & Brainerd (2000) Wright, C. O., & Brainerd, T. G. 2000, ApJ, 534, 34, doi: 10.1086/308744
  * Xu & Jing (2022) Xu, K., & Jing, Y. 2022, ApJ, 926, 130, doi: 10.3847/1538-4357/ac4707
  * Xu et al. (2022a) Xu, K., Jing, Y. P., & Gao, H. 2022a, ApJ, 939, 104, doi: 10.3847/1538-4357/ac8f47
  * Xu et al. (2022b) Xu, K., Jing, Y. P., Zheng, Y., & Gao, H. 2022b, arXiv e-prints, arXiv:2211.02665.  https://arxiv.org/abs/2211.02665
  * Xu et al. (2022c) Xu, K., Zheng, Y., & Jing, Y. 2022c, ApJ, 925, 31, doi: 10.3847/1538-4357/ac38a2
  * Yang et al. (2021) Yang, H., Gao, L., Frenk, C. S., et al. 2021, arXiv e-prints, arXiv:2110.04434.  https://arxiv.org/abs/2110.04434
  * Yang et al. (2006) Yang, X., van den Bosch, F. C., Mo, H. J., et al. 2006, MNRAS, 369, 1293, doi: 10.1111/j.1365-2966.2006.10373.x
  * Zhang & Komatsu (2011) Zhang, J., & Komatsu, E. 2011, MNRAS, 414, 1047, doi: 10.1111/j.1365-2966.2011.18436.x
  * Zhang et al. (2022a) Zhang, J., Liu, C., Vaquero, P. A., et al. 2022a, AJ, 164, 128, doi: 10.3847/1538-3881/ac84d8
  * Zhang et al. (2015) Zhang, J., Luo, W., & Foucaud, S. 2015, J. Cosmology Astropart. Phys, 2015, 024, doi: 10.1088/1475-7516/2015/01/024
  * Zhang et al. (2017) Zhang, J., Zhang, P., & Luo, W. 2017, ApJ, 834, 8, doi: 10.3847/1538-4357/834/1/8
  * Zhang et al. (2019) Zhang, J., Dong, F., Li, H., et al. 2019, ApJ, 875, 48, doi: 10.3847/1538-4357/ab1080
  * Zhang et al. (2022b) Zhang, Y., Yang, X., & Guo, H. 2022b, MNRAS, 517, 3579, doi: 10.1093/mnras/stac293410.48550/arXiv.2210.05215
  * Zheng et al. (2005) Zheng, Z., Berlind, A. A., Weinberg, D. H., et al. 2005, ApJ, 633, 791, doi: 10.1086/466510
  * Zhou et al. (2021) Zhou, R., Newman, J. A., Mao, Y.-Y., et al. 2021, MNRAS, 501, 3309, doi: 10.1093/mnras/staa3764
  * Zou et al. (2019) Zou, H., Gao, J., Zhou, X., & Kong, X. 2019, ApJS, 242, 8, doi: 10.3847/1538-4365/ab1847

##  Appendix A Consistency of the lensing results

Figure 9: Lensing profiles for locally brightest galaxies (LBGs) selected from
SDSS DR7, grouped in different stellar mass bins. Black points show the
results obtained by using the Fourier_Quad pipeline (FQ), using DECaLS data,
and green points represent the results from Mandelbaum et al. (2016) and Wang
et al. (2016) (Wang 16). The ranges in stellar mass are shown by the text in
each panel, in unit of $\log_{10}M_{\odot}$. Note since Mandelbaum et al.
(2016) and Wang et al. (2016) did not include corrections to the stellar mass
of massive galaxies (see Appendix B below), no stellar mass corrections have
been made to LBGs in this plot, to ensure fair comparisons.

Figure 9 shows the comparison between the lensing signals calculated using
Fourier_Quad and the results from Mandelbaum et al. (2016) or Wang et al.
(2016). Here we use exactly the same sample of locally brightest galaxies
(LBGs) from SDSS DR7 as Wang et al. (2016), and the only difference is the
shear catalog. We used photometric sources from DECaLS, whereas Wang et al.
(2016) used photometric sources from SDSS. Our Fourier_Quad method is also
different from that of the previous studies. Encouragingly, we find a very
good agreement between our results and Wang et al. (2016), which allows us to
confidently use the Fourier_Quad method to calculate the lensing signals in
this paper. Notably, Wang et al. (2016) or Mandelbaum et al. (2016) did not
include any correction to the stellar mass of LBGs, and thus to ensure fair
comparisons with previous studies, stellar mass corrections are not included
in our results of Figure 9.

##  Appendix B Stellar Mass Correction

Figure 10: Lensing signals around locally brightest galaxies (LBGs) in six
different stellar mass bins. The ranges in stellar mass are shown by the text
in each panel, in unit of $\log_{10}M_{\odot}$. Purple points represent the
observational results without stellar mass correction. Orange points represent
the observational signals after the correction in stellar mass has been
applied. The black lines show the TNG300 predictions. Without the stellar mass
corrections, the observed lensing signals are significantly overestimated in
the two most massive panels.

He et al. (2013) performed a photometric analysis of early-type galaxies
(ETGs) from SDSS DR7, with $r$-band magnitude brighter than -22.5, stellar
mass greater than $10^{11}\mathrm{M_{\odot}}$ and redshifts between 0.05 and
0.15. They argue that SDSS underestimated the luminosity due to two possible
reasons. Firstly, the PHOTO algorithm of SDSS overestimates the sky background
subtraction, which leads to an underestimate of the luminosity. This
overestimate of the sky background is mainly due to wrong masking of the
galaxy neighbors and extended stellar halos. Aihara et al. (2011) raised this
problem for SDSS DR7, and introduced an improved sky subtraction algorithm for
SDSS DR8, leading to a smaller but still significant underestimate in
magnitude. Secondly, the surface brightness distribution of galaxies does not
fit well to a Sersic profile (Sérsic (1963)), and this can also cause
underestimates in galaxy luminosity or stellar mass.

He et al. (2013) achieved a more accurate photometry for SDSS ETGs by
improving the sky background subtraction. In their method, they first used
Sextractor Bertin & Arnouts (1996) for object masking. They found that, for
bright galaxies of $M_{\mathrm{r}}<-22.5$, a significant number of companion
objects were not properly masked, mainly due to the galaxies being present in
crowded fields, or having an extended stellar halo, or being near foreground
bright stars. They then correct the masks by hand, one by one. This correction
leads to a much more accurate photometry, and they measured an underestimate
in magnitude by $\sim$0.2-0.3 mag for SDSS DR7.

Following this line of thought, we perform a stellar mass correction (see
Section 3.4 for details). We then re-calculate the lensing signals with the
new stellar mass binning. Figure 10 shows the observational results before and
after the correction in stellar mass, together with the simulated results from
TNG300. Here we use the same sample of locally brightest galaxies (LBGs) from
SDSS DR7 as Wang et al. (2016). We find that the stellar mass correction
solves the discrepancy found in the two most massive bins. Without the
correction in stellar mass, the observed lensing signals are significantly
higher than TNG in the two most massive bins and in outskirts.

##  Appendix C Satellite profiles in TNG100

Figure 11: Satellite profiles around observed blue (blue points) and red (red
points) ICGs, and TNG100 blue (cyan shaded curves) and red (magenta shaded
curves) ICGs, divided in different stellar mass bins, at redshift $z=0$.
Satellites more massive than $10^{8}\mathrm{M_{\odot}}$ are included to
calculate the profiles. The ranges in stellar mass are shown by the text in
each panel, in unit of $\log_{10}M_{\odot}$. Errorbars are the 1-$\sigma$
scatters of 100 bootstrap subsamples.

Figure 11 shows the satellite profiles for TNG100 ICGs, compared to real
observations. The resolution limit of TNG100 is about one order of magnitude
better than TNG300, so we include satellites which are more massive than
$10^{8}\mathrm{M_{\odot}}$ in Figure 11. Similar to the TNG300 results, red
and blue ICGs in TNG100 show different trends from those of real observations,
though noisier. This discards the resolution limit as a possible explanation
to these discrepancies.

##  Appendix D Infall time distribution of satellite galaxies

Figure 12: Cumulative distribution of the infall time for satellites around
red (red curves) and blue (blue curves) ICGs in TNG300. Each panel represents
the results for a different stellar mass bin, indicated by the text at the
bottom of each panel, in unit of $\log_{10}\mathrm{M_{\odot}}$.

Figure 12 shows the cumulative infall time distribution121212We defined the
infall time of a $z=0$ satellite as the highest redshift when it was still a
satellite galaxy. of $z=0$ surviving satellites in TNG. We identify some
differences between the infall time distribution of satellites around red and
blue ICGs. The red curves have slightly earlier infall time distributions than
the blue curves in the three middle panels. This indicates that for ICGs with
$10.5<\log_{10}M_{\ast}/\mathrm{M_{\odot}}<11.4$, earlier infall of satellites
around red ICGs can partially explain the stronger tidal stripping. However,
we also note that the infall time distribution of satellites around red and
blue ICGs looks similar in the most massive panel of Figure 12. And in the
second most massive panel, the difference between red and blue curves is
minor. The stronger stripping of satellites around red galaxies in TNG might
depend on other factors in addition to the infall time distribution, such as
the difference in host halo masses between red and blue ICGs at higher
redshifts, the orbits of infalling satellites and how long can satellites
maintain their star formation after infall and hence stay more self bound etc.