© The Authors 2025

1. Introduction

Satellite galaxies experience significant changes as they plunge into denser galaxy ecosystems and interact¹ with more massive structures. Many studies have extensively discussed important related subjects, such as the hierarchical assembly and clustering of galaxies (White & Rees 1978; Merritt 1985; Navarro et al. 1996; Mo et al. 2010, and references therein), and the environmental processes that shape satellite galaxies, such as the interactions with the harsh gas environment in galaxy clusters and groups through ram pressure stripping, instabilities, and thermal heating (e.g. Gunn et al. 1972; Mori & Burkert 2000; Boselli & Gavazzi 2006; Peng et al. 2010; Armillotta et al. 2017; Boselli et al. 2022).

Satellite evolution can also be driven by the gravitational field of their host systems and/or by interactions with other satellites (harassment; Toomre & Toomre 1972; Merritt 1983, 1984; Gnedin 2003; Smith et al. 2016). The host’s gravitational tidal field can change the satellite mass distribution through tidal (de-)compressions, which can change the mass densities in satellites (Valluri 1993; Dekel et al. 2003; Masi 2007; Renaud et al. 2009), and through tidal stripping and total satellite disruption (cannibalism; Merritt 1985; Read et al. 2006a,b; Fellhauer et al. 2006a,b, 2008; Penarrubia et al. 2008; Blaña et al. 2015; Ogiya et al. 2022a; Errani & Navarro 2021; Errani et al. 2022, 2023; Stücker et al. 2023; Montero-Dorta et al. 2024). As a consequence, the long-time accreted satellites that continue to survive will only have their central dense regions gravitationally bound, whereas their external, less dense regions, as well as most satellites, are cannibalised by the hosting cluster, especially in the central regions of clusters. First-infall satellites, on the other hand, will have their outer regions stripped only where the host tidal field becomes stronger than the local gravity of the satellite. This establishes a relation between the mass densities of satellite galaxies and their distances from the host (the $\overline{\rho }-r$ relation). Previous studies of the Milky Way (MW) satellite galaxies located within ∼80 kpc found such a relation (Kaplinghat et al. 2019; Robles & Bullock 2021; Pace et al. 2022; Cardona-Barrero et al. 2023; Hammer et al. 2024).

In this article, we investigate this relation further by analysing more massive environments, such as galaxy clusters, while also analysing the mean mass densities of different components (stellar, gas, and dark matter). We also adopt three approaches to study this relation: (1) semi-analytical satellite orbital models with a simplified stripping model, (2)cosmological galaxy simulations, and (3) observations of satellites in galaxy clusters and in the Local Group. The sections of the article are ordered as follows: In Sect. 2 we explain the theoretical framework. In Sect. 3 we describe the methodology used for modelling, analysis, and the data from simulation and observations. In Sect. 4 we present the results, followed by a discussion in Sect. 5. We end in Sect. 6 with our conclusions.

2. Theoretical background

2.1. The relation between the satellite mean mass densities, pericentres, and the host tidal field

When satellite galaxies enter their hosting cluster or group, they are affected by the gravitational tidal field of the main host. The closer the satellite is to the centre of the host, the stronger the tidal field becomes, stripping and accreting the material in the outskirts of the satellite. In consequence, satellites that have long been accreted and still survive the tidal forces have only their dense central regions remaining and gravitationally bound. This is expected to produce the ($\overline{\rho }-r$) relation, which is an (anti-)correlation between satellite mass densities and their orbital pericentres (r_peri). For example, Kaplinghat et al. (2019, see their Fig. 2) and Cardona-Barrero et al. (2023, their Fig. 1) quantified the correlation between the central dynamical mass densities of MW satellites measured at a radius of 150 pc from the satellite centre (ρ₁₅₀), and their orbital pericentres for a few MW satellites within 90 kpc. Similarly, Hayashi et al. (2020, their Fig. 8) compared the values of the MW satellite ρ₁₅₀ with data from cosmological galaxy simulations. Although these studies find a measurable correlation, they used the dynamical mass densities of the satellites measured at a fixed radius of 150 pc, which does not consistently trace environmental tidal effects. In fact, the quantity ρ₁₅₀ was instead proposed to probe the impact of stellar feedback on the central dark-matter densities of dwarf galaxies (Read et al. 2019). Furthermore, Pace et al. (2022, their Fig. 5) do explore the density-pericentre relation ($\overline{\rho }-r$) comparing the satellite mean dynamical mass densities measured within their stellar half-mass radii (${\overline{\rho }}_{h\star }$) with the orbital pericentres and with the MW tidal field, finding that most satellites have densities above the tidal field strength equivalent. Moreover, the analysis of the stellar mass densities of MW star clusters also shows a correlation with their pericentres and the MW tidal field strength, as these systems experience the same environmental tidal forces as satellite galaxies. (Gieles et al. 2011; Carballo-Bello et al. 2012; Webb et al. 2013, and references therein). This has also been demonstrated in the globular cluster system of NGC 1399, the central giant ellitpical galaxy in the Fornax galaxy cluster, where the globular clusters in the circum-nuclear region were found to be systematically smaller and denser as a result of the strong tidal field (Puzia et al. 2014). However, while Gaia has made it possible to estimate the orbital pericentres of the MW satellites in the inner region and in its outskirts, albeit with large observational uncertainties (McConnachie & Venn 2020; McConnachie et al. 2020) (see also Blaña et al. 2020), this becomes technically impossible for distant extra-galactic satellite systems.

In this work, we develop a model to investigate how the gravitational tidal field of a host galaxy cluster or group influences the mean mass densities of its satellite galaxy population, and to identify potentially observable signatures of these effects. To this end, we employ a simplified tidal stripping model presented in Binney & Merrifield (1998, see Eq. (8.91), (8.92), (8.107), (8.108)) (for a more detailed model see Read et al. 2006a). A satellite of total mass m, subject to the tidal field |τ| of its host cluster or group, defines a region around the satellite with a maximum volume of $\frac{4}{3}\pi r_{\mathrm{J}}^3$ with a (Jacobi) radius r_J, in which the material remains bound to the satellite. The combination of this volume and the satellite mass m yields the Jacobi mean mass density (${\overline{\rho }}_{\mathrm{J}}$), which depends on the tidal field as follows:

$$\begin{array}{cc}& {\overline{\rho }}_{\mathrm{J}}=|\tau |.\hfill \end{array}$$(1)

The tidal field (magnitude) of the host system (τ) depends on the gradient of the acceleration field of the host’s mass distribution. For simplicity, we adopt the approximation for a spherically symmetric system, which depends on the mean mass density profile of the host (${\overline{\rho }}^{\mathrm{H}}$) as a function of the distance from the host centre r given as

$$\begin{array}{cc}& |\tau |:=3f^{-3}{\overline{\rho }}^{\mathrm{H}}\left(r\right)=3f^{-3}\frac{M^{\mathrm{H}}(r)}{43\pi r^3};\mathrm{where}f={(1-\frac{1}{3}\frac{\mathrm{d}ln{M}^{\mathrm{H}}}{\mathrm{d}lnr})}^{-1/3},\hfill \end{array}$$(2)

where M^H(r) is the cumulative mass profile of the host and the function f improves the modelling for extended host mass distributions, which simplifies to f = 1 for a point such as mass or the far-field approximation. Considering typical mass profiles for cluster haloes, such as the NFW density profile (Navarro et al. 1996), results in a tidal field profile that decreases with distance (|τ|∝1/r).

The Jacobi density of the satellite (${\overline{\rho }}_{\mathrm{J}}$) is the mean mass density calculated within r_J, i.e.:

$$\begin{array}{cc}& {\overline{\rho }}_{\mathrm{J}}:=\frac{m\left(r_{\mathrm{J}}\right)}{43\pi r_{\mathrm{J}}^3},\mathrm{where}{r}_{\mathrm{J}}={\left(\frac{m}{43\pi |\tau (r)|}\right)}^{1/3}\hfill \end{array}$$(3)

is the instantaneous Jacobi radius, or also called Hills radius or tidal radius.

As a satellite orbits within a galaxy cluster and approaches its centre, the tidal field strength |τ| increases, causing the Jacobi radius to shrink. If r_J becomes smaller than the satellite’s mass truncation radius (r_tr), the satellite begins to lose mass through tidal mass stripping, further reducing r_J. Using Eq. (1) and assuming a mass m(r′) and mean mass density $\overline{\rho }(r\prime )$ profiles for the satellite, where $\overline{\rho }(r\prime )=m(r\prime )/(\frac{4}{3}\pi r^{\prime 3})$ depends on the satellite-centric coordinate r′, we can measure the value of r′ at which the satellite’s mean mass density profile exceeds the tidal field and, thus, remains gravitationally bound to the satellite; namely, where

$$\begin{array}{cc}\hfill {\overline{\rho }}_{\mathrm{tr}}& \ge {\overline{\rho }}_{\mathrm{J}}=|\tau |,\hfill \end{array}$$(4)

with ${\overline{\rho }}_{\mathrm{tr}}=\overline{\rho }(r\prime =r_{\mathrm{tr}})$. With this density we can determine the mass that is still bound and within the truncation radius m_tr​ = ​m(r′​ = ​r_tr), and the truncation radius constrained by r_J according to

$$\begin{array}{cc}\hfill r_{\mathrm{tr}}& \le r_{\mathrm{J}}={\left(\frac{m_{\mathrm{tr}}}{43\pi |\tau (r)|}\right)}^{1/3}.\hfill \end{array}$$(5)

Using Eq. (4) we can predict how the mean mass density of a satellite will evolve as it falls into a cluster. For example, if a satellite has an extended mass distribution with a large initial truncation radius r_tr, init and a low initial mean mass density (${\overline{\rho }}_{\mathrm{tr},\mathrm{init}}$), we can predict that, as soon as the satellite approaches the cluster centre, the tidal field becomes stronger and the outer mass layers of the satellite becomes stripped. This shrinks the truncation radius (r_tr < r_tr, init) and the mean mass density within the new truncation radius increases. This shows that for satellite tidal stripping to occur, neither mass nor size individually matters, but the resulting mean mass density.

Furthermore, as a satellite orbits its host, it can experience different strengths of the tidal field for non-circular orbits, where typically the host tidal field anti-correlates with distance (|τ|∝1/r). In particular, a satellite that approaches its host for the first time with a radial orbit experiences an increasingly stronger tidal field strength as it approaches the host, reaching its maximum at the orbital pericentre (r_peri):

$$\begin{array}{cc}\hfill {|\tau |}_{\max }& =|\tau (r_{\mathrm{peri}})|.\hfill \end{array}$$(6)

This implies that the maximum tidal field experienced by a satellite is inversely proportional to its pericentre distance (|τ|_max ∝ 1/r_peri). Therefore, according to Eq. (3), the maximum mean mass density due to tidal truncation is reached at the pericentre, where the tidal field is maximised, i.e.

$$\begin{array}{cc}\hfill {\overline{\rho }}_{\mathrm{tr}}^{\max }& ={\overline{\rho }}_{\mathrm{J}}^{\max }={|\tau |}_{\max },\hfill \end{array}$$(7)

which is expected for most mass density profiles that we can adopt for the hosting cluster or group (e.g. NFW, Burkert 2000). Consequently, the smallest mass and truncation radius are reached at the pericentre, $m_{\mathrm{tr}}^{\min }=m(r\prime =r_{\mathrm{tr}}^{\min })$ and $r_{\mathrm{tr}}^{\min }$, and the maximum density anti-correlates with the pericentres, ${\overline{\rho }}_{\mathrm{tr}}^{\max }\propto 1/r_{\mathrm{peri}}$, allowing us to obtain the anti-correlated mass-density-distance relation ($\overline{\rho }-r$).

However, since measuring the truncation radius to determine the total mean mass densities can be impractical or difficult in observations or simulations, one alternative is to use the half-mass radius of the remaining bound mass of the satellite within its truncation radius. This is defined as

$$\begin{array}{c}\hfill m(h:=r\prime )=12m_{\mathrm{tr}},\end{array}$$(8)

which depends not only on the shape of the satellite’s mass profile, but also on the truncation radius h​ = ​h(m, r_tr). This allows us to define the mean mass densities within one and twice the half-mass radius, as follows:

$$\begin{array}{cc}\hfill {\overline{\rho }}_h& =\frac{12m_{\mathrm{tr}}}{43\pi h^3},\hfill \end{array}$$(9)

$$\begin{array}{cc}\hfill {\overline{\rho }}_{2h}& =\frac{m_{2h}}{43\pi {(2h)}^3},\hfill \end{array}$$(10)

where m_2h​ := ​m(r′​​ = ​​2h). Moreover, it can be shown that for standard mass profiles (e.g. NFW, Plummer, Exponential), the mean mass densities are higher in the internal regions of a satellite than in its outskirts, i.e. ${\overline{\rho }}_{\mathrm{tr}}\le {\overline{\rho }}_{2h}\le {\overline{\rho }}_h$. This implies that Eq. (4) also imposes

$$\begin{array}{cc}& |\tau |\le {\overline{\rho }}_{\mathrm{tr}}<{\overline{\rho }}_{2h}<{\overline{\rho }}_h.\hfill \end{array}$$(11)

Furthermore, when the satellite reaches its orbital pericentre, the half-mass radius and mean mass densities can also reach their respective minimum ($r_h^{\min }$) and maximum values (${\overline{\rho }}_h^{\max }$ and ${\overline{\rho }}_{2h}^{\max }$).

In the following, we explore whether we can find similar relations for the stellar mass component of satellite galaxies. Eqs. (1) and (4) relate the host’s tidal field to the mean mass density of a satellite within its truncation radius (${\overline{\rho }}_{\mathrm{tr}}$). We can reformulate ${\overline{\rho }}_{\mathrm{tr}}$ as the sum of the mean mass densities of dark matter (χ) and stars (⋆), i.e. ${\overline{\rho }}_{\mathrm{tr}}={\overline{\rho }}_{\mathrm{tr}}^{\star }+{\overline{\rho }}_{\mathrm{tr}}^{\chi }$. For simplicity, we consider a gas-free satellite galaxy, which allow us to re-write Eq. (4) as

$$\begin{array}{cc}& {\overline{\rho }}_{\mathrm{tr}}^{\star }+{\overline{\rho }}_{\mathrm{tr}}^{\chi }\ge |\tau |.\hfill \end{array}$$(12)

Now, in an hypothetical case, we can consider a pristine satellite galaxy that begins its first infall. Using Eq. (12) we can predict that as the satellite approaches the host halo and |τ| increases, the satellite begins to lose its outer layers of material (initially mostly dark matter) and, as a consequence, it shrinks its initial truncation radius r_tr, init according to Eq. (5), resulting in a mean mass density larger than before its infall (i.e. ${\overline{\rho }}_{\mathrm{tr}}^{\chi }>{\overline{\rho }}_{\mathrm{tr},\mathrm{init}}^{\chi }$). Later, as the satellite sinks deeper into the host’s tidal field, the outer layers of the stellar material are increasingly more strongly stripped as well, and both mass components reach the same truncation radius $r_{\mathrm{tr}}^{\star }=r_{\mathrm{tr}}$. As a result, given the typical mass density profiles of satellite galaxies (e.g. exponential or Hernquist profiles), a stripped satellite is expected to have a mean stellar mass density higher than its pre-infall value (i.e. ${\overline{\rho }}_{\mathrm{tr}}^{\star }>{\overline{\rho }}_{\mathrm{tr},\mathrm{init}}^{\star }$), which is due to the shrinking of the truncation radius (r_tr < r_tr, init).

In addition, given that the direct measurement of the stellar mass truncation radius ($r_{\mathrm{tr}}^{\star }$) is also difficult to determine in observations of galaxies with faint substructures, we use instead the stellar half-mass radius (h⋆), and/or multiples of this quantity, as this depends on the mass profile, but also on $r_{\mathrm{tr}}^{\star }$, as follows:

$$\begin{array}{c}\hfill m^{\star }(h\star :=r\prime )=12m_{\mathrm{tr}}^{\star }.\end{array}$$(13)

Similar to the case of the mean dynamical mass densities, h_⋆ allows us to define the mean stellar mass densities within one and two h_⋆ as

$$\begin{array}{cc}\hfill {\overline{\rho }}_{h\star }^{\star }& =\frac{12m_{\mathrm{tr}}^{\star }}{43\pi h_{\star }^3},\hfill \end{array}$$(14)

$$\begin{array}{cc}\hfill {\overline{\rho }}_{2h\star }^{\star }& =\frac{m_{2h\star }^{\star }}{43\pi {(2h_{\star })}^3},\hfill \end{array}$$(15)

where $m_{2h\star }^{\star }:=m^{\star }(r\prime =2h\star )$, with a dependency on $r_{\mathrm{tr}}^{\star }$ through $m_{\mathrm{tr}}^{\star }$. Moreover, following the arguments of Eq. (11) between inner and outer densities for typical stellar mass densities, we can also state that for each satellite galaxy:

$$\begin{array}{c}\hfill {\overline{\rho }}_{\mathrm{tr}}^{\star }<{\overline{\rho }}_{2h\star }^{\star }<{\overline{\rho }}_{h\star }^{\star }.\end{array}$$(16)

Satellite’s stellar masses can be estimated more easily than their dynamical masses, as the former requires adopting stellar-mass-to-light ratios from stellar population models, while the latter requires resolved kinematic observations. Stellar masses can be a useful direct tracer of the radial variation of the mean mass densities of the satellite population.

Modelling the tidal shape of the mass densities of a satellite population, with the equations described above, implies an assumption that satellites behave as rigid profiles and mass stripping removes the outer shells, layer by layer, shrinking r_tr and, therefore, also h and h⋆ to their minimum values when the satellite passes through the pericentre, r_peri. More complex models in the literature also find that the half-mass radius reduces during tidal stripping. This is expected, as the tidal field and stripping can also induce an internal restructuring of the bound material in a satellite galaxy, as its density reacts to its own potential (ϕ_sat) and the host’s potential (ϕ_host) as ϕ_sat, tot = ϕ_sat + ϕ_host. Studies of the tidal tracks of Errani et al. (2022, see their Fig. 2 and 13) and Stücker et al. (2023) reveal that for typical satellite mass distributions (e.g. NFW, Plummer, Exponential), the density profiles are reshaped by the tidal field, preferentially reducing the size of the half-mass radii of the total and stellar mass distributions during the stripping process. Moreover, similar results were obtained for MW star clusters, as expected, because they experience the same environmental tidal forces (Webb et al. 2013). The studies mentioned above suggest that the definitions of h and h_⋆ in Eqs. (8) and (13) capture the overall shrinking of these half-mass radii produced by tidal stripping. Therefore, we use the previous equations to later build in Sect. 3.1 a simple model to understand the effects of tidal stripping on a satellite population in galaxy clusters and groups. Nevertheless, more extreme tidal mechanisms can generate more complex dynamical configurations, such as tidal shocks (Penarrubia et al. 2008; Hammer et al. 2023, 2024), or stages before and after complete tidal disruption that can produce ultradiffuse galaxies (Ogiya 2018) and streams (Smith et al. 2013; Blaña et al. 2015).

2.2. Tracing the distribution of mass densities and additional properties of a satellite galaxy population

Most of the studies in the literature that explore the effects of the host tidal field on the mass densities of its satellite galaxies use their pericentre distances (r_peri; Kaplinghat et al. 2019; Hayashi et al. 2020; Pace et al. 2022; Cardona-Barrero et al. 2023), following Eqs. (6) and (7). However, while the Gaia observatory has allowed us to estimate orbital pericentres of satellite galaxies within the MW, and in its outskirts with larger uncertainties, this becomes impossible with the current technology for extragalactic systems beyond the Local Group. Therefore, the use of observed cluster-centric distances of satellite galaxies (r) is the only alternative.

First infalling satellites have not yet reached their pericentres and, therefore, they are expected to exhibit mass densities that correlate with the local tidal field at their current positions (r) according to Eq. (4). Moreover, satellite orbital distributions in equilibrium can show a correlation between radial satellite distances (r) and their pericentres (r ∝ r_peri), which can have a spread with a size that depends on the shape of the distribution. Furthermore, since we are interested in studying the effects of the tidal field of a galaxy cluster or group, we require selecting systems that have a dominant central potential where |τ|∝1/r. In addition, we require that the host potential evolves slowly compared to the satellite orbital evolution and, therefore, preferentially avoid clusters that are currently undergoing subcluster mergers where complex tidal forces vary on dynamical timescales similar to those of the satellite orbits. Studies find that most low-redshift clusters and groups have their central regions already formed and in equilibrium (Merritt 1985), with the accretion of low-mass satellites and only a few cases of ongoing galaxy cluster major mergers having little effect on the large-scale halo potential (Tempel et al. 2017; Łokas 2023). Therefore, in the most common cases, it is possible to derive a relation between the distributions of the satellites’ cluster-centric distances r and their pericentres, r_peri, showing that, on average, there is a direct correlation between these quantities, such that r ∝ r_peri.

Moreover, given that we are interested in inspecting the collective effect of the host tidal field on its entire satellite population, we use the moving average of the satellite mass densities to characterise their radial transformation. With this we can probe the (anti)correlation between the distribution of satellite mass densities and cluster-centric distances $⟨\overline{\rho }⟩\propto ⟨1/r⟩$. For this, we define the moving average function as follows:

$$\begin{array}{cc}\hfill ⟨Y⟩(r_i)& =K(r_i,Y_i,q)={10}^{\frac{1}{q}{\displaystyle {\sum }_{j=-q/2}^{j=+q/2}}log{Y}_j},\hfill \end{array}$$(17)

which depends on the distance to the host’s centre, r, the number of neighbours, q, and the variable, Y_i. This function is used to calculate a logarithmic moving averaged filter of the same variable ⟨Y_i⟩. We applied this to different satellite variables, Y: dynamical mass densities (${\overline{\rho }}_{\mathrm{tr}}$, ${\overline{\rho }}_{2h}$), stellar mass densities (${\overline{\rho }}_{2h\star }^{\star }$), stellar mass ($M_{2h\star }^{\star }$), pericentres (r_peri), among others. Here, we adopted a logarithmic weighted average filter because: i) the typical mass densities in the satellite populations considered in this study are found to show a log-normal behaviour (see Sect. 4.2) and span two orders of magnitude, and ii) the presence of mass density outliers that could introduce strong variations if the kernel were linear. Therefore, we aim at measuring the bulk of the distribution to determine the distance behaviour of the mass density of the whole satellite population. In addition, we also tested a median filter that yields similar results, leaving our main conclusions unchanged.

3. Methodology

To explore and probe the relation of the mass densities of satellite galaxies and their distances to the host system and the tidal field we used three approaches: (1) a semi-analytical satellite orbital toy model, (2) cosmological galaxy simulations, and (3) observational data of four systems. Each approach allowed us to investigate different aspects. For example, the toy model allowed us to obtain a simple interpretation of the equations in Sect. 2 on the effects of tidal truncation in an environment with controlled variables without including other physical processes that are involved in the evolution of satellites, such as ram-pressure stripping, starvation, etc. Moreover, access to the full orbital history of a satellite allowed us to calculate peri-, apo-centres, and timescales. We could also repeat the calculations while considering different setups and potentials, capabilities that are not possible in cosmological galaxy simulations. However, in this simplified toy model, the satellites were modelled as rigid systems, where the tidal forces strip only outer layers of its dark and baryonic material in an ‘onion’-like interpretation of Eqs. (4) and (7), leaving their central mass distributions unperturbed. Therefore, we also used cosmological galaxy simulations to probe the mass-density distance relation. These simulations included multiple physical processes, such as star formation, dynamical friction, and other processes mentioned above. Lastly, we included observed satellite systems to compare with the models. However, we cannot here access all the information, such as their full dark-matter distributions, 3D spatial position or velocity, or the full orbital histories of a limited sample. We explain the setup of these approaches in the following sections.

3.1. A simplified satellite tidal stripping toy model

This model is an orbital construction of satellite galaxies in a galaxy cluster. We chose a system with the mass and size of the Fornax cluster, as this system is larger in number of satellites and in total mass than the Local Group, with a well-studied and resolved satellite population. The satellites were modelled as test particles that have analytical mass density profiles for the stellar and dark matter components evolving under the influence of the tidal field of the cluster according to the equations of Sect. 2.1. For simplicity, we considered gas-less satellites only. For the toy model, we began by defining the gravitational force of the host, where we set up a Fornax-type cluster with an NFW profile with a virial mass and radius of M_vir = 7 × 10¹³ M_⊙ and R_vir​ = ​1 Mpc (Drinkwater et al. 2001; Schuberth et al. 2010), with a concentration of c = 8.5, and where we used Eq. 2 to calculate its tidal field profile |τ(r)|.

To have an observationally motivated satellite population we used the Schechter Luminosity Function to draw their luminosities for a total of 2 × 10³ satellites, ranging between M_r​ = ​ − 22 and −9 mag. For simplicity, we adopted a constant stellar mass-to-light ratio of ${\mathrm{{\rm Y}}}_r=2{\mathrm{M}}_{\odot }{{\mathrm{L}}_{\odot }}^{-1}$ to obtain the stellar masses, and tested variations between 1 and $4{\mathrm{M}}_{\odot }{{\mathrm{L}}_{\odot }}^{-1}$. To model the stellar mass density profile ρ^⋆(r′), mean stellar mass density ${\overline{\rho }}^{\star }(r\prime )$ and stellar mass m^⋆(r′) profile of each satellite and the impact on the results, we tested three different profiles as a function of distance to the satellite centre (r′): exponential, Plummer, and Hernquist, with initial scale lengths taken to be 10% of the initial satellite size (initial truncation radius r_{tr,  init}). We used the galaxy-halo abundance matching models of Moster et al. (2013) to estimate the maximum dark matter mass that each satellite galaxy contains prior to its infall. Guided by cosmological galaxy simulations and observations of satellite beyond R_vir from their hosts, we sampled the initial total mean mass densities from a log-normal distributions centred at ${\overline{\rho }}_{\mathrm{tr},\mathrm{init}}=2\times {10}^5{\mathrm{M}}_{\odot }{\mathrm{kpc}}^{-3}$ with a half width of 0.1dex, and tested different initial conditions (IC) as well. To model the dark-matter profiles of each satellite, we adopted the NFW model, with concentrations from the c-M halo relation (Correa et al. 2015). The positions and velocities of the satellites were sampled from the NFW equilibrium distribution of the host cluster (test particles) using the Jeans equations. For simplicity, we used an isotropic spherical distribution for the satellites and the cluster potential, and we neglected the increase in the mass of the cluster by satellite accretion since the total satellite mass corresponds to ∼3% of the cluster mass.

Lastly, we performed the orbital calculation with an updated version of the DELOREAN code (Blaña et al. 2020), integrating the orbits for 10 Gyr to phase-mix the initial distribution, to ensure that all satellites passed through their pericentres. From this we calculated the tidal field and the tidal mass stripping for each satellite from Eq. (4) in Sect. 2 and determined their half-mass radii and mean mass densities for their total and stellar mass components (see also Ogiya et al. 2022b). From this we obtained the mean dynamical² mass density of each satellite within its truncation radius r_tr (${\overline{\rho }}_{\mathrm{tr}}$), and within 2h (${\overline{\rho }}_{2h}$), as well as the mean stellar mass density of each satellite within r_tr (${\overline{\rho }}_{\mathit{tr}}^{\star }$) and within 2h⋆ (${\overline{\rho }}_{2h\star }^{\star }$). Moreover, we considered two cases: when 100% of the satellites pass through their orbital pericentres (r_per) obtaining their maximum densities according to Eq. (7), and when 30% of the satellites are on their first infall, and, therefore, their truncated mass densities do not yet reached the maximum values. Then, we explored the resulting distribution of the satellite mass densities with distance using Eq. (17), while comparing with the host’s tidal field. We also used the moving average kernel of Eq. (17) applied to the pericentre distances to compare the pericentre distribution of the satellite population as a function of the cluster-centric distance (i.e. ⟨r_peri⟩∝r). Using the toy model, we show in Sect. 4.1 how the mean mass densities of the satellite population correlate with the cluster-centric distances and the host tidal field.

3.2. Simulated systems

We used data from cosmological galaxy simulations to probe the host’s tidal effects on satellites, which have the advantage over our simple model of having self-gravitating satellite systems that react to the environmental tides self-consistently. For this, we used clusters from ILLUSTRIS TNG50 (Nelson et al. 2019; Pillepich et al. 2019), and the TNG100 and TNG300 simulations (Marinacci et al. 2018; Naiman et al. 2018; Nelson et al. 2018; Pillepich et al. 2018; Springel et al. 2018). Physical and numerical parameters are reported in several publications and on the official webpage³, but we briefly detail the main properties in the following. The periodic box comoving sizes of TNG50, TNG100 and TNG300 are 51.7 cMpc, 106.5 cMpc and 302.6 cMpc, respectively. TNG50 has stellar and dark matter particle masses of 8 × 10⁴ M_⊙ and 4.5 × 10⁵ M_⊙ with a spatial particle resolution of 290 pc and a moving adaptive mesh for the gas, with star-forming cells with a median size of 138 pc, and a recorded minimum value of 8 pc. TNG100 has stellar and dark matter particle masses of 1.4 × 10⁶ M_⊙ and 7.5 × 10⁶ M_⊙, with a resolution of 740 pc and a moving adaptive mesh for the gas, with star-forming cells with a median size of 355 pc, and a recorded minimum value of 14 pc. TNG300 has stellar and dark matter particle masses of 11 × 10⁶ M_⊙ and 59 × 10⁶ M_⊙, with a resolution of 1.48 kpc and a moving adaptive mesh for the gas, with star-forming cells with a median size of 715 pc, and a recorded minimum value of 47 pc.

We analysed different redshifts, but in this work, we present the results for redshifts z​ = ​0 and 1. After requiring subhaloes to be composed of > 100 stellar particles each (van den Bosch et al. 2018; van den Bosch & Ogiya 2018), and each cluster to have at least 60 resolved subhaloes, we obtained from TNG50, TNG100, and TNG300 a total of 235 galaxy clusters (and groups) at z = 0 and 60 at z = 1, which we used for the subsequent analysis. We tested the effects of resolution and sampling by increasing the limit by one order of magnitude or the effects of massive subhalo outliers from the main distribution by setting a maximum stellar mass of 10¹⁰ M_⊙. We found similar results, leaving our main conclusions unchanged. From these simulations, we selected two clusters for a more detailed analysis as our fiducial models, taking a Fornax-type cluster analogue from the TNG50 simulation with a virial mass and radius of 7.6 × 10¹³ M_⊙ and 1115 kpc, respectively, with 492 satellite galaxies, and a Virgo-type cluster from TNG100 with a virial mass and radius of 4.9 × 10¹⁴ M_⊙ and 2079 kpc, respectively. To explore the impact of different subgrid models in the simulations, we also made a comparison with a cluster from the EAGLE simulations (Schaye et al. 2015). The selected cluster and satellite galaxies correspond to a Virgo-type cluster with a virial mass and radius of 3 × 10¹⁴ M_⊙ and 1422 kpc, respectively, containing 271 resolved satellites. The cluster was taken from the L100N1504 run, which consists of a periodic box of 100 cMpc in comoving size, initially containing 1504³ gas particles with an initial mass of 1.81 × 10⁶ M_⊙, and the same number of dark matter particles with a mass of 9.70 × 10⁶ M_⊙. The halo catalogues provided in the public database (McAlpine et al. 2016) were built using a friend-of-friends (FoF) algorithm, which identifies dark-matter overdensities following Davis et al. (1985), considering a linking length of 0.2 times the average inter-particle spacing. Baryonic particles were assigned to the FoF halo of their closest dark-matter particle. Subhalo catalogues were built using the SUBFIND algorithm (Springel et al. 2001; Dolag et al. 2009), which identifies local overdensities using a binding energy criterion for particles within a FoF halo.

The tidal field profiles of the host clusters (|τ|) were calculated from their total mass profiles (gas, stars, and dark matter). We verify that these clusters are the most massive structures within four times their virial radii. Analysis of more and less massive clusters shows the same behaviour as in our fiducial case.

We computed the virial radius R_vir of a cluster from the virial mass M_vir defined in Mo et al. (2010) as an over-dense region with the physical virial radius defined as

$$\begin{array}{cc}\hfill R_{\mathrm{vir}}& ={\left(\frac{M_{\mathrm{vir}}(\mathrm{z})}{4/3\pi {\rho }_{\mathrm{crit}}(\mathrm{z}){\mathrm{\Omega }}_{\mathrm{m}}(\mathrm{z}){\mathrm{\Delta }}_{\mathrm{vir}}(\mathrm{z})}\right)}^{1/3},\hfill \end{array}$$(18)

where the cosmological parameters ρ_crit, Ω_m, Δ_vir are the critical density, the matter to critical density ratio, and the spherical collapse over-density criterion Δ_vir ≈ (8π²+82ϵ−39ϵ²)(ϵ+1)⁻¹ with ϵ = Ω_m − 1.

3.3. Observed systems

We included four observed systems to explore the empirical effect of tidal fields on the mean mass densities of satellites: the galaxy clusters Fornax and Virgo, and the Local Group satellites of the MW and Andromeda (M31) galaxies. We selected these systems for tracing as large a host mass range as possible (10¹²–10^14.6 M_⊙), to ensure deep photometric observations with a large spatial coverage that conveys homogenous and high sample completeness, and because they are publicly available. We modelled the tidal field of observed clusters using mass profiles with parameters based on the literature. In the following, we briefly describe the main characteristics of each dataset:

• Fornax cluster: For the Fornax cluster mass profile we adopted parameters determined by Schuberth et al. (2010) and derived from dynamical models fitted to star cluster kinematics. We selected their NFW model a5 as our fiducial model, which has a virial radius of R_vir​ ≃ ​R₂₀₀​ = ​1075 kpc and mass of M_vir = 7.3 × 10¹³ M_⊙ that agree well with previous estimates (Drinkwater et al. 2001). Our fiducial model also includes the stellar mass component of NGC 1399 modelled according to Schuberth et al. (2010). For the analysis of the Fornax cluster satellites, we used the Next Generation Fornax Survey (NGFS) observations obtained with the wide-field Dark Energy Camera mounted on the 4-m Blanco telescope at the CTIO in Chile (Muñoz et al. 2015; Eigenthaler et al. 2018; Ordenes-Briceño et al. 2018a,b), having the advantage of being among the widest-deepest surveys on a single galaxy cluster, covering ∼50deg² with a wide photometric filter range between NUV and NIR. We used the combined catalogues of Eigenthaler et al. (2018) and Ordenes-Briceño et al. (2018a) with 630 satellite galaxies, which contain photometric information and stellar masses based on colour information. Taking the NGC 1399 galaxy as the centre of the cluster, which is at 20.4 Mpc (Blakeslee et al. 2009), the sample extends to R​ = ​1017 kpc (2.8° ) in projected radius, which includes faint satellites down to M_i′​ = ​ − 8.8 mag (m_i′​ = ​22.73 mag), with objects as extended as R_e​ = ​2.7 kpc (27.9 arcsec) in the effective radius of Sérsic. To probe the relation of Eq. (12), we calculated the stellar mean densities of the Fornax satellites within twice their deprojected effective radii r_e (adopting the conversion r_e​ = ​4/3R_e; see Wolf et al. 2010), obtaining the mean stellar mass densities $({\overline{\rho }}_{2h\star }^{\star })$ for each satellite, where, for simplicity, we approximated h ⋆ ​ ≈ ​r_e. We tested luminosity-weighted densities, adopting a constant stellar mass-to-light ratio Υ_i′ in the i′-band, finding similar results.

• Virgo cluster: To probe the tidal relation in a more massive environment, we selected the Virgo cluster. This cluster has two main substructures, the more massive Virgo A subgroup with a virial mass and radius of $M_{200}^{\mathrm{VA}}=4.2\times {10}^{14}{\mathrm{M}}_{\odot }$ and ${R_{\mathrm{vir}}}^{\mathrm{VA}}=1550\mathrm{kpc}$ respectively, and the Virgo B subgroup with a smaller virial mass and radius of $M_{200}^{\mathrm{VB}}={10}^{14}{\mathrm{M}}_{\odot }$ and ${R_{\mathrm{vir}}}^{\mathrm{VB}}=960\mathrm{kpc}$ (McLaughlin 1999; Ferrarese et al. 2012). For the analysis of the Virgo cluster satellites, we used the publicly available extended catalogue of Kim et al. (2015) with ∼1500 satellites that extend out to R ≃ 3.7R_vir (5730 kpc) in projected radius, containing galaxies brighter than m_r​ ≃ ​18.76 mag, or in absolute magnitude of approximately −12.32 mag, adopting a distance of 16.5 Mpc (Mei et al. 2007). We focused on the Virgo A subgroup and took the central galaxy M87 as the origin, while excluding the subgroup B centred on M49 by masking out satellites within $R={R_{\mathrm{vir}}}^{\mathrm{VB}}/2\simeq 480\mathrm{kpc}$ in projection. We note that by including all satellites only mildly changes the Virgo A satellite distribution, as Virgo B is separated by 4.4deg (1266 kpc) in projection. We determined the mean stellar mass densities $({\overline{\rho }}_{2h\star }^{\star })$ for each satellite using its de-projected stellar half-mass radius h⋆, using the r-band half-light radii and a constant Υ_r. We point out that changing Υ_r would only shift vertically the mass-density-distance relation ($\overline{\rho }-r$) in the density axis, while differences of Υ_r within the satellite population introduce insignificant variations, three orders of magnitude smaller than the spread of the distribution in stellar densities, as shown in Sect. 4.3.

• Local Group systems: To explore the tidal field relation in lower-mass environments, we probed the M31 and MW galaxy systems. The mass models of these hosts include their dark haloes, bulges, stellar and gaseous disc components with parameters adapted from Tamm et al. (2012) and Blaña et al. (2017, 2018) for M31 (M₂₀₀ = 1.04 × 10¹² M_⊙, R₂₀₀​ = ​207 kpc), and from Bland-Hawthorn & Gerhard (2016), Shen et al. (2022) for the MW (M₂₀₀ = 1.08 × 10¹² M_⊙, R₂₀₀ = 216 kpc). We used the catalogue of Putman et al. (2021) with 55 satellites for the MW and 41 for M31 and calculated their mean stellar mass densities (${\overline{\rho }}_{2h\star }^{\star }$) within twice their half-mass radii 2h⋆, and their central mean dynamical mass densities (${\overline{\rho }}_{h\star }$) within h⋆.

4. Results

4.1. Satellite mass densities: Toy model predictions

We start with the results of the toy model in Fig. 1, presenting a snapshot after the 10 Gyr integration. There we show the total and stellar mass mean densities of each satellite measured at their truncation radii (r_tr) and twice their half-mass radii (2h, 2h⋆) and plotted as functions of their distances to the centre of the cluster. We include the moving average of the mean mass densities of the satellite distribution applying the function in Eq. (17). In addition, we show the tidal field of the cluster, which has the largest values in the cluster centre.

Fig. 1. Results of the toy model showing the imprint of the host tidal field on the mean mass density distribution (circles) of nsat​ = ​2 × 103 satellite galaxies shown as a function their distances to the cluster centre. The host has an NFW halo mass and size of a Fornax-type cluster, showing its tidal field (|τ|) in both panels (purple curves). The satellites are taken from a snapshot after a 10 Gyr orbital integration with DELOREAN (Blaña et al. 2020). Left panel: We show the final total (i.e. dynamical) mass mean densities ${\overline{\rho }}_{\mathrm{tr}}$ (large grey circles) and mean stellar mass densities ${\overline{\rho }}_{\mathrm{tr}}^{\star }$ (yellow circles) after 10 Gyr. The moving average profiles of these densities are calculated with Eq. (17) using $q=2\sqrt{n_{\mathrm{sat}}}\simeq 89$ neighbours and are shown with solid thick black and orange curves. The IC of the total and stellar mass densities of the satellites are shown with small points in grey and brown, respectively, showing also their moving average profiles (dashed curves). In addition, we show the effect of including a population of 30% of the first in-falling satellites and the resulting moving average (thin black curve). Right panel: For the same satellites, we show the total and stellar mass mean densities measured in the central regions of each satellite within twice the half-mass radius of the total mass (2h) and the stars (2h⋆), as well as their moving averages (solid curves). The increase in the moving average profiles towards the cluster centre shows that the shrinking of the truncation radius of each satellite (rtr) produces an increasing central mean mass density as a result of the shrinking half-mass radii, h and h⋆. Moreover, the effect is more noticeable in the inner cluster regions (r < 0.5Rvir) where the tidal field is stronger and the satellites have smaller pericentres, making them more susceptible to stronger tidal stripping. The horizontal dot-dashed lines included in both panels are added for reference to illustrate the deviation of the mean stellar mass densities from constant values.

As Fig. 1 reveals, we find that the mean mass densities of the satellite population have systematically higher values towards the centre of the cluster, resulting in a spatial mean mass density segregation. In contrast, we also show the initial satellite mass densities, before being truncated by the tidal field, which show the moving averages centred around the IC. This distribution is produced in this model by the synergy between two effects: (i) the profile of the tidal field of the cluster (|τ|) that imprints its shape according to Eqs. (1) and (7), and (ii) the distribution of the pericentre distances that depend on the satellite orbital distribution. We analysed the distribution of the satellite pericentres applying the moving average function of Eq. (17), finding that these correlate with the distance ⟨r_peri⟩∝⟨r⟩ with a scatter that depends on the orbital distribution. A small fraction of satellites (1.3%) with small pericentres (r​ < ​100 kpc) are at large distances (r​ > ​1 Mpc) in the snapshot taken to produce Fig. 1. These satellites have excentric radial orbits, which explains the presence of some satellites with high mean mass densities in the outskirts of the cluster, since they previously entered the cluster central regions where the tidal forces are stronger. These could correspond to the backsplash satellite population detected in cosmological galaxy simulations and observations (Gill et al. 2005; Teyssier et al. 2012; Garrison-Kimmel et al. 2014; Blaña et al. 2020; Haggar et al. 2020; Diemer 2021). In general, this small satellite sub-population has insignificant effects on the conclusions.

In addition, we explore how the mass density distribution is affected when we include 30% of first-infall satellites, initially with low mass densities. In this case, as shown in Fig. 1, we find that the moving average slightly shifts the overall distribution towards the initial condition (IC) values. Similarly, the minimum distance of the first-infall satellites (r_min) corresponds to their current position (r), yielding a tighter correlation between r_min and r of the entire satellite population (since these distances are the same for the 30%).

We show the effect of the tidal stripping on the mean mass densities measured in the satellite’s central regions, measured at twice the total half-mass radius (2h) and at twice the stellar half-mass radius (2h⋆). As Fig. 1 (right panel) shows, the central mean mass densities increase towards the cluster centre. This is a consequence of the shrinking of h and h⋆ due to the mass stripping of the outer layers of the total and stellar masses that define them (Eqs. (8) and (13)). In addition, as expected from Eqs. (11) and (14), the mean mass density of each satellite measured at its half-mass radius (e.g. h⋆) is higher than at larger radii (e.g. r_tr).

Furthermore, we see in Fig. 1 that the moving average of the total and stellar mass mean density starts to increase within r < 0.5 R_vir, being slightly steeper for the stellar mass densities. This can be seen in the figure for mass densities measured at the truncation radius (left panel), or at their half-mass radii (h or h⋆, right panel). Besides the Hernquist profile used for the satellite mass models shown in Fig. 1, we also tested exponential and Plummer profiles, resulting in a similar shape of the moving average of the satellite mass densities. In the following Sect. 4.2 we show the moving average of stellar mass densities for models that include more realistic physical processes.

4.2. Satellite mass densities: Galaxy simulations

The toy model offers a simplified modelling of the tidal evolution of satellite galaxies, predicting the resulting radial distribution of the mean mass densities of the satellite population, and the strength of the tidal field of the host cluster that drives the evolution. However, the satellites are represented by a rigid model, which has limitations in including additional physical processes, such as tidal shocks, galaxy harassment, star formation, and gas and ram-pressure stripping. Therefore, we performed an analysis of satellite galaxies in cosmological galaxy simulations that include these processes.

We start by presenting the analysis of a Fornax-type analogue cluster in detail, selecting a cluster with a mass of M = 7.6 × 10¹³ M_⊙ from the Illustris TNG50 simulation. We see in Fig. 2 (top panel) that the dynamical mean mass density profile of the satellites $({\overline{\rho }}_{2h})$ follows the tidal relation of Eq. (11), with satellites being tidally shaped and having mass densities that are systematically higher in the central region of the cluster. As the tidal field strength of the host cluster increases towards the cluster centre, it clearly demonstrates its effect on shaping the mass density profile of the satellite population ($⟨{\overline{\rho }}_{2h}⟩$). Moreover, we observe a similar behaviour to the satellite mean mass density profile predicted in the toy model (Fig. 1), where most satellites have ${\overline{\rho }}_{2h}$ values larger than the tidal field at any distance r, because they have already passed their pericentres where |τ| was larger. The presence of some satellites having mean mass densities lower than the tidal field strength in Fig. 2 (top panel) reveals that secondary effects, such as the satellite’s internal processes (star-formation, morphological substructure, etc.), orbital characteristics, or the presence of substructures within the galaxy cluster, may play some role for individual satellites. Furthermore, a careful inspection of the density profile in Fig. 2 (top panel) reveals that the gradient increases within r < 0.5 R_vir towards the cluster centre. It worthy to point out that at these short distances, the tidal field strength of the cluster surpasses the lowest mean mass densities of satellites coming from the outskirts of the cluster i.e. where $|\tau |>{\overline{\rho }}_{2h}\sim {10}^5{\mathrm{M}}_{\odot }{\mathrm{kpc}}^{-3}$.

Fig. 2. Mean mass densities of 492 satellite galaxies (circles) as a function of distance to their hosting TNG50 Fornax-type cluster, with colours indicating their stellar masses (colour bar). The cluster has a virial mass and radius of Mvir = 7.6 × 1013 M⊙ and Rvir​ = ​1075 kpc, respectively, and corresponds to a snapshot at redshift z​ = ​0. Its tidal field radial profile (|τ|) is shown in all panels (purple curve). Shown as a function of cluster-centric distance from the top to bottom panels are: the dynamical mass densities of each satellite ${\overline{\rho }}_{2h}$ (top), measured within 2h; the central mean dynamical mass densities ${\overline{\rho }}_{2h\star }$ within 2h⋆ (middle), and the central stellar mass densities ${\overline{\rho }}_{2h\star }^{\star }$ within 2h⋆ (bottom). The blue curves in each plot show the moving average (Mov. Av.) of different mass densities for all the satellites using Eq. (17). The measurement of the transition radius (ℜ⋆ metric [ii]) is shown with a vertical dot-dashed line (middle and bottom panels). Histograms are shown in the right column sub-panels to show the shape of the mean mass density distributions.

Next, we examine the central mean dynamical-mass densities of the satellites (${\overline{\rho }}_{2h\star }$; see Fig. 2, middle panel). We find that in the cluster outskirts (r > 0.5R_vir) the shape of the moving average of the satellite mass densities $⟨{\overline{\rho }}_{2h\star }⟩$ is consistent with a small positive slope and a roughly constant value beyond r​​ ≈ ​R_vir. Inside r​​ ≈ ​​0.5 R_vir, the central regions within the satellites become increasingly more tidally affected, having their outer layers stripped, which shifts the population’s moving average to higher densities that increase towards the cluster centre. The mass density profile $⟨{\overline{\rho }}_{2h\star }⟩$ has a negative gradient, similar to the behaviour of the mean (total) dynamical-mass density profile $⟨{\overline{\rho }}_{2h}⟩$ (top panel).

Following the analysis of our toy model, we also examine the central mean stellar mass densities of satellites (${\overline{\rho }}_{2h\star }^{\star }$) as a function of cluster-centric distance⁴, shown in Fig. 2 (bottom panel). We find a distribution with more scatter than the central dynamical-mass densities, which is due to the variations of the baryon-to-dark matter fractions between satellites. The moving average of the stellar mass density has a roughly constant value of $⟨{\overline{\rho }}_{2h\star }^{\star }⟩\approx {10}^6{\mathrm{M}}_{\odot }{\mathrm{kpc}}^3$ at large distances r ≳ 1R_vir. It decreases slightly within r ≲ 1 R_vir, reaching a value of $⟨{\overline{\rho }}_{2h\star }^{\star }⟩\approx 5\times {10}^5{\mathrm{M}}_{\odot }{\mathrm{kpc}}^3$. Further in at r ≲ 0.5 R_vir, the average stellar mass density profile starts to increase significantly, following the cluster’s tidal field strength. The correlation between the central dynamical and stellar mass density profiles can be appreciated by comparing the middle and bottom panels in Fig. 2, and in Fig. 3 where we directly compare both density variables, finding that the scatter is less than 1 dex. Overall we find that the moving average of the total, central, and stellar mass density profiles ($⟨{\overline{\rho }}_{2h}⟩$, $⟨{\overline{\rho }}_{2h\star }⟩$, $⟨{\overline{\rho }}_{2h\star }^{\star }⟩$) behave in a similar way within r < 0.5 R_vir, where the tidal field of the cluster dominates.

Fig. 3. Comparison between the mean central dynamical and stellar mass densities for satellites in the TNG50 Fornax-type cluster. The symbol colour encodes stellar mass, and the identity function is shown as the dot-dashed line.

Comparing with the toy model predictions (Fig. 1) we see that the mean stellar mass densities of the toy model are shifted to values lower compared to those in the TNG simulations. However, the shape of both profiles are very similar, especially showing the same transition at r ≃ 0.5 R_vir.

The TNG simulations and our toy model indicate that mass density profiles undergo a transition at r ≈ 0.5 R_vir, within which the cluster’s tidal field strength appears to dominate satellite densities, while further out we observe roughly constant, perhaps slightly increasing satellite densities as a function of cluster-centric distance. In the following, we seek to quantify this distance at r​​ ≈ 0.5 R_vir as the ‘transition radius’ parameter ℜ_⋆. We defined two metrics to measure ℜ_⋆: using the minimum of the moving average of the mean stellar mass profile $⟨{\overline{\rho }}_{2h\star }^{\star }⟩$ within r < R_vir (ℜ_⋆[i]), and using the minimum of the derivative of the moving average profile (ℜ_⋆[ii]). Later, in Sect. 5.1, we provide a complete detailed description of two metrics that we define in Eq. (19) and 20 to measure this feature, as well as estimates of the uncertainties in different simulated and observed systems (Table 1). Here we report that the transition radius measured using $⟨{\overline{\rho }}_{2h\star }^{\star }⟩$ for this TNG cluster is ℜ_⋆​ = ​0.49 R_vir (542 kpc), which is marked with the vertical line in Fig. 2 (bottom panel). In the same figure, we also plot ℜ_⋆ in the middle panel to show that it matches the profile shape of central dynamical mass densities ($⟨{\overline{\rho }}_{2h\star }⟩$).

Furthermore, we also analyse in detail two additional clusters with masses similar to the Virgo cluster (M_vir = 4.2 × 10¹⁴ M_⊙). One cluster is selected from Illustris TNG100, while the other was selected from the EAGLE simulations (Schaye et al. 2015) to explore how different star-formation prescriptions and sub-grid models impact our results. In general, we find that the satellites in these two cluster models exhibit similar behaviour in their mean mass density profiles (${\overline{\rho }}_{2h}$, ${\overline{\rho }}_{2h\star }^{\star }$), having progressively declining values within r ≲ 0.5 R_vir and a transition radius ℜ_⋆ in the mean stellar mass density profiles.

Galaxy clusters and groups are in constant dynamical evolution presenting significant substructures due to massive subgroups of galaxies being accreted, and with occasional major mergers at different redshifts (McGee et al. 2009; Pallero et al. 2019; Kuchner et al. 2022). Therefore, we also analyse a total of 235 galaxy clusters and groups with a range of masses between 9 × 10¹² M_⊙ and 3 × 10¹⁵ M_⊙, and at different redshifts from the Illustris TNG simulations, to explore possible perturbations in the mass density profiles of satellite populations due to infalling substructures in the cluster/group environments. In Fig. 4 we show the dynamical mass densities and the central stellar mass densities of satellites of 155 clusters in the TNG300 simulation at redshifts z​ = ​0 and 31 clusters at z = 1, which have the largest number of satellites.

Fig. 4. Mass densities of satellite galaxies vs their cluster centric distances for 155 host galaxy clusters from a snapshot at redshift z​ = ​0 (left column) and for 31 clusters at z​ = ​1 (right column), taken from TNG300. The tidal field profiles of all the clusters are shown in each panel (|τ| in purple curves). Top panels show the total mass mean densities of the subhaloes (circles, colour-coded by total subhalo masses), and the bottom panels show the central mean stellar mass densities (circles, colour-coded by stellar masses). To avoid overcrowding, we plotted a sub-sample with 10% (left column) and 20% (right column) of the satellites (circles). Circle sizes are proportional to h (top panels) or h⋆ (bottom panels). In each panel the moving average satellite mean mass density profile of each cluster is shown in with an orange curve. In addition, it is shown the global moving average in dashed black curve that is calculated with the moving averages of all clusters. All are moving averages are calculated with Eq. 17. The transition radii (ℜ⋆), measured in each cluster, are located within the vertical grey region (bottom panels), which are determined from the moving average of the mean stellar mass densities $⟨{\overline{\rho }}_{2h\star }^{\star }⟩$, marking where the $⟨{\overline{\rho }}_{2h\star }^{\star }⟩$ start to increase further in. We present a similar analysis of TNG100 and TNG50 in Figs. A.1 and A.2, respectively.

In order to study how the satellite mean mass densities behave in lower mass host environments, we present in Figs. A.1 and A.2 similar results analysing the TNG100 and TNG50 simulations. In general, we find consistent results in the three sets of simulations (TNG50, TNG100, TNG300), where the satellite populations exhibit the same cluster-centric behaviour in their mass density profiles as discussed earlier for the Fornax-type cluster. We measured the transition radius in each of these clusters and groups in TNG300, TNG100, and TNG50, finding typical values in the range of ℜ_⋆ = 0.4R_vir and 0.6R_vir. Moreover, clusters can present substructures due to different processes, such as major or minor mergers of groups (visible as wiggles in the orange curves in Figs. 4, A.1 and A.2). The analysis at redshifts z​ = ​1 and 0 shows a consistent mass density profile behaviour with cluster-centric distance. However, we defer the characterisation of the temporal evolution of the satellite mass density profiles to future works, as this requires the careful assembly of observed satellite samples at corresponding redshifts, which are currently unavailable.

Until now, we have only examined the total and stellar mass densities of satellite galaxies. Although the total (dynamical) mass densities include the gas component, we have not mentioned how the gas properties of satellites may change with distance to the host cluster. As rapid gas-mass loss could be driving internal changes in the stellar and dark matter distribution, we proceed to measure the gas mass content of satellites and the resulting mass densities. We find that within r ≲ 0.5 R_vir only a small percentage of satellites (≲2%) have gas mass densities comparable to the values of stellar mass densities, having gas-to-baryon mass fractions of up to 10%. Later in Sect. 5.3 we discuss the gas mass properties in more detail, but in general we find that the gas-rich satellites are not the dominant population in the central regions of clusters.

4.3. Satellite mass densities: Observations

In this section we explore the ($\overline{\rho }-r$) relation for a selection of satellite galaxies in observed systems. For this we show the mean stellar mass densities of satellites in the Fornax and Virgo galaxy clusters, which have the advantage of having a large observed satellite samples. We also include data from Local Group satellites that have accurate measurements of their internal kinematics and photometry, as well as 3D distances to their MW and M31 host, which allow us to determine their central dynamical and stellar mass densities.

We present the radial mass density distribution profiles for these four observed satellite systems in Fig. 5. Overall, they share a distribution of the stellar mass densities similar to simulations, where the moving average $⟨{\overline{\rho }}_{2h\star }^{\star }⟩$ decreases slightly from the outskirts of the hosts (r​ ≳ ​R_vir) inwards until r ≈ 0.4 − 0.5 R_vir, where a minimum is reached and $⟨{\overline{\rho }}_{2h\star }^{\star }⟩$ starts to increase further inwards. Later in Sect. 5.1 we show in detail how this transition radius (ℜ_⋆) in the moving average profile of the satellite mean stellar mass densities is revealed on a linear scale (see Fig. 6). Furthermore, we find that the location of ℜ_⋆ in the observations (Fornax, Virgo, and M31) is 10% smaller in R_vir units than in the simulations. Such differences could be attributed to total mass overestimation in observations, which could result in larger values of R_vir (Eq. (18)), and/or due to projection effects of the satellite spatial distribution. In fact, in all simulations the parameter ℜ_⋆, measured in projection, has values systematically smaller in R_vir units than measured in three dimensions, as expected for most radial mass distributions (see Sect. 5.4).

Fig. 5. Satellite galaxy mean stellar mass densities (${\overline{\rho }}_{2h\star }^{\star }$ as circles) of four observed systems: the Fornax galaxy cluster (top left panel), the Virgo cluster (top right panel), the MW (bottom left panel), and M31 (bottom right panel), with the satellite stellar masses indicated in colour. In all the panels, we show the tidal field for the host systems (|τ| in purple curves). In each panel we include the moving average profiles of the mean stellar mass densities using Eq. (17) with $q=\sqrt{n_{\mathrm{sat}}}$ (blue curves, with the blue areas casting one standard deviation). The measured transition radii ℜ⋆ are shown for each system (dash-dotted vertical lines), and the expected location of the splashback radii ∼1.4Rvir (vertical dotted lines). For MW and M31, due to the low number of satellites, we add thinner blue curves that show moving averages $⟨{\overline{\rho }}_{2h\star }^{\star }⟩$ where each is calculated with number of neighbours between $q=0.5\sqrt{n_{\mathrm{sat}}}$ and $2\sqrt{n_{\mathrm{sat}}}$. This indicates that, independent of q, all the curves show a low value of $⟨{\overline{\rho }}_{2h\star }^{\star }⟩$ at ℜ⋆. In addition, for MW and M31 we show the central mean dynamical mass densities of each satellite measured within h⋆ (triangles) and the respective moving average profiles (brown curves, bottom panels). This can be compared with other studies in the literature on the $\overline{\rho }-r$ relation for the MW satellites mentioned in Sect. 2, which show that using r instead of the pericentre also reveals an anti-correlation with distance ($\overline{\rho }\propto 1/r$). Moreover, as predicted by the simulations, the central dynamical mass densities show a change in the moving average profile where ℜ⋆ is measured (see Fig. 2, middle panel). Histograms included in the top panels show the overall shape and peak of the distribution of ${\overline{\rho }}_{2h\star }^{\star }$ for the satellites of each cluster.

Fig. 6. Measurement examples of the transition radius ℜ⋆: the Fornax-type cluster (top panel) from TNG50, and the Fornax cluster (bottom panel). From a total of 106 sub-sample realisations and the 106 moving average curves of the satellite mean stellar mass densities, we show here 200 (blue curves). Each curve is calculated with Eq. (17) with a sub-sample of $n_{\mathrm{sub},\mathrm{sat}}=n_{\mathrm{sat}}-\sqrt{n_{\mathrm{sat}}}$ satellites from a total of nsat. The neighbour number is a random number between $q=2\sqrt{n_{\mathrm{sub}}}$ and $4\sqrt{n_{\mathrm{sub}}}$. The first method of measurement (ℜ⋆[i]) uses Eq. (19) and is shown with vertical dot-dashed blue lines and their 1-sigma confidence region. The second method (ℜ⋆[ii]) uses Eq. (20) and is shown with vertical dot dash orange lines and their 1-sigma confidence region. Bootstrapping reveals a small scatter in the moving average profiles, due to the satellite-to-satellite sub-sample variations and the different number of neighbours used. Lastly, we also show the moving average after first randomising (shuffling) the mean mass densities of the satellites, which erases the signature in the satellite population (see Null-hypothesis test in Sect. 5.2).

Towards the cluster centres, the moving average profiles of Fornax and Virgo in Fig. 5 reveal a dip in the mass densities at R ∼ 0.1R_vir, and increasing values towards the centre within. This lower mass density dip may be attributed to projection effects, where distant interlopers with lower mass densities positioned along the line-of-sight in radial filaments decrease the projected moving average of the satellite mean mass densities in the central cluster regions (see Sect. 5.4). Another cause could be the presence of central massive cD galaxies in these two clusters (M87 in Virgo and 1399 in Fornax). The tidal shocks and dynamical friction that these central galaxies produce could be strong enough to quickly cannibalise neighbouring satellites, reducing their masses and densities more efficiently. Future, more detailed comparisons with simulations could reveal whether these deviations are produced by stochastic substructure fluctuations and line-of-sight effects, or whether important physical processes will require further improvement in simulations (e.g. see van den Bosch et al. 2018).

Here we mostly focused on the satellite stellar mass densities, which can be readily determined in large surveys using stellar population synthesis models. The dynamical mass densities, on the other hand, require relatively high-resolution spectroscopy to measure the velocity dispersion of each satellite galaxy, which is observationally challenging to perform for the hundreds or thousands of satellites in galaxy clusters and groups. Nevertheless, we present here two important examples, the MW and M31 systems, the proximity of which allows for more complete kinematic estimates, despite their lower sample statistics containing a few dozen satellites each. An inspection of the bottom panels in Fig. 5 shows that the central mean dynamical-mass densities ${\overline{\rho }}_{h\star }$ of satellites in the MW and M31 exhibit systematically higher average densities towards to the centre of the hosting galaxies. This demonstrates that the ($\overline{\rho }-r$) relation found in simulations is also traceable in the Local Group, as other studies have shown using satellite orbital pericentres (Kaplinghat et al. 2019; Hayashi et al. 2020; Pace et al. 2022; Cardona-Barrero et al. 2023). Moreover, while the dynamical-mass densities have a radial trend similar to their stellar counterparts, the former are shifted to larger values because of their dominant dark matter components, in addition to the measurement being performed within 1h⋆ instead 2h⋆, decreasing the volume by a factor of eight. The $⟨{\overline{\rho }}_{h\star }⟩$ profiles show a similar behaviour as the simulations (Fig. 2), with a weak decrease in the average satellite density from the outskirts towards the transition radius ℜ_⋆. Within this region, the mean satellite densities (${\overline{\rho }}_{h\star }$ and ${\overline{\rho }}_{2h\star }^{\star }$) increase significantly towards the host centre. However, the low number of satellites (55 for the MW and 41 for M31) introduces fluctuations due to stochastic satellite-to-satellite density variations. As example, in the MW the moving average of the satellite stellar and dynamical mass densities decreases in the very central region (r < 0.2 R_vir). This is a consequence of including satellites that are strongly tidally disrupted, such as the dSph Sagittarius and Boötes III, which have large stellar half-mass radii at h ⋆  = 2.6 kpc and h ⋆  = 1.3 kpc, respectively, as well as low luminosities producing low stellar mass densities. However, their central dynamical mass densities are still higher than the local MW tidal field, following the ($\overline{\rho }-r$) relation (see purple curves in Fig. 5).

5. Discussion

Here we discuss in more detail the measurement procedure of the transition radius ℜ_⋆ and possible systematic effects that could influence its estimation. We also discuss physical processes that could impact the stellar mass densities of satellite galaxies, and how these results relate to other works in the literature.

5.1. Measuring the transition radius ℜ_⋆

We show, using the toy model, simulations, and observations, that the stellar and dynamical mean mass densities of satellite galaxies tend to be higher in the central regions of their hosting environments. Furthermore, within a transition radius ℜ_⋆ the moving average profile of the satellite central stellar mass densities $⟨{\overline{\rho }}^{\star }⟩$ increases towards the cluster/group centre.

In Sect. 4.2 we qualitatively describe this transition radius ℜ_⋆ using the profiles of $⟨{\overline{\rho }}^{\star }⟩$, and plot the corresponding ℜ_⋆ values in the figures where we presented the data of simulations and observations (see Figs. 2, 4, and 5). Here we describe in detail how we measure ℜ_⋆ and the possible systematic uncertainties of its measurement. We define two metrics, showing two examples in Fig. 6 with simulations and observations, where we measure ℜ_⋆ according to the following equations:

$$\begin{array}{cc}\hfill {\mathfrak{R}}_{\star }[\mathrm{i}]& :=R\left({⟨{\overline{\rho }}_{2h\star }^{\star }⟩}_{\mathfrak{R}\star }\right)\hfill \end{array}$$(19)

$$\begin{array}{cc}\hfill {\mathfrak{R}}_{\star }[\mathrm{ii}]& :=R({\partial }_r⟨{\overline{\rho }}_{2h\star }^{\star }⟩=0)\hfill \end{array}$$(20)

The first method ([i]) measures a minimum in the moving average of the satellite mean stellar mass density profile in the region 0.2 < R/R_vir < 1, that is, where ${⟨{\overline{\rho }}_{2h\star }^{\star }⟩}_{\mathfrak{R}\star }:=⟨{\overline{\rho }}_{2h\star }^{\star }⟩(\min )$. For this, we show in Fig. 6 data from simulations (Fornax-type cluster, top panel), and from observations (Fornax cluster, bottom panel) showing the application of this metric in a linear scale where the minimum is better revealed.

The second method ([ii]) uses a smoothing function (polynomial of order 7) that is fitted directly to the data to calculate the profile of the derivative. The order of the polynomial is obtained from optimising it to well approximate the moving average profile near its minimum. From this ℜ_⋆ is found where the derivative in the radial range of 0.2 < R/R_vir < 1 becomes zero. Systems that have substructures can produce additional zero roots, as seen in the Fornax cluster (Fig. 6, bottom panel). Therefore, we select the first zero that appears in the range 0.2 < R/R_vir < 1 starting from the lower bound, which is where the polynomial function fits the profile well near the minimum. Even though both methods are arbitrary and empirical they agree well, as shown in Fig. 6, with both capturing the trends of the moving averages. Such methods with derivatives of mass density profiles are also commonly applied to measure changes in density profile slopes to detect and quantify substructures in dark and stellar haloes (O’Neil et al. 2021; Aung et al. 2021).

Given that the first method can be affected by the arbitrary number of neighbours (q) used to calculate the moving average, and that both methods can be affected by the sample number of satellites, substructures, and outliers, we adopted a Monte Carlo technique to measure ℜ_⋆ and its uncertainties. For this we performed 10⁶ sub-sample realisations to measure ℜ_⋆ each time applying a bootstrap (Jackknife) method taking a sub-sample with $n_{\mathrm{sub}\mathrm{sat}}=n_{\mathrm{sat}}-\sqrt{n_{\mathrm{sat}}}$ satellites for each system with n_sat satellites, while taking a range of neighbours between $q=2\sqrt{n_{\mathrm{sub},\mathrm{sat}}}$ and $4\sqrt{n_{\mathrm{sub},\mathrm{sat}}}$. This generates 10⁶ moving average profiles $⟨{\overline{\rho }}_{2h\star }^{\star }⟩$ that can slightly vary due to the data sub-sample. We show in Fig. 6 the moving averages for 200 realisations as an example.

We tested the propagation of photometric errors (20% for the luminosity and 10% for the half-light radius) into the sub-sampling by adding noise to mean stellar mass density of each satellite (observations and simulation) from a normal distribution with a standard deviation of 30% (${\sigma }_{\overline{\rho }\star }$). This results in a slight increase in the spread of $\mathrm{\Delta }⟨{\overline{\rho }}_{2h\star }^{\star }⟩\approx {10}^5{\mathrm{M}}_{\odot }{\mathrm{kpc}}^{-3}$, as shown in Fig. 6. Furthermore, we also tested a modified version of Eq. (17) to calculate the moving averages of the mean stellar mass densities by weighting them with the inverse of the photometric errors ($\mathit{w}=1/{({\sigma }_{\overline{\rho }\star })}^2$), finding similar values of ℜ_⋆. Therefore, in our fiducial procedure, we decide to omit the contribution of photometric errors. The resulting distributions of the measurements of ℜ_⋆ for the different observed and simulated systems and their confidence ranges are presented in Fig. 7 and summarised in Table 1. The distributions of ℜ_⋆ show that the metric that uses the derivative (ℜ_⋆[ii]) is more stable than that that uses the minimum (ℜ_⋆[i]). This is because ℜ_⋆[i] directly measures the minimum in the moving average, which can be affected by outliers during the bootstrap re-sampling, while ℜ_⋆[ii] is obtained from fitting a function of the re-sampled distributions, reducing stochasticity. Furthermore, we also find that ℜ_⋆[i] can be more stable when the spatial completeness of the satellite information is limited, or significant substructures perturb the moving average of the satellite mass densities.

Fig. 7. Distribution of measurements of ℜ⋆, marginalising over the distribution of ${⟨{\overline{\rho }}_{2h\star }^{\star }⟩}_{{\mathfrak{R}}_{\star }}$ (minimum in the moving average of the mean stellar mass densities profiles) using methods [i] (left panel) with Eq. (19) and [ii] (right panel) with Eq. (20). The top four rows of panels show the results for the observed systems (Fornax, Virgo, MW, and M31), and the bottom row panels present the simulated TNG50 Fornax-type cluster, with ℜ⋆ measured in projection (2D) and in 3D.

We find that the transition radius ℜ_⋆ has consistent values in virial radius units (ℜ_⋆​ ≃ ​0.5 R_vir) across the studied environments in simulations and observations. Furthermore, we also measure ℜ_⋆ in the simulations in projection (2D), where for each of the 10⁶ realisations we randomly choose a projection (Table 1). The comparison with ℜ_⋆ measured from the 3D distribution reveals that the projected value is systematically smaller by ∼15%, as generally expected for a projected mass distribution (see Wolf et al. 2010).

5.2. Null-hypothesis test

The signature of ℜ_⋆ is detected at ℜ_⋆​ ≃ ​0.5 R_vir, where we find a minimum and a change of the moving average of the mean mass density profiles of the satellite population as a function of cluster-centric distance. Within this distance, the average density profile of satellites increases towards the centre. Therefore, a coherent global structure where this mass density minimum occurs is encoded by the mean mass densities of satellites. If this coherent structure were produced randomly, statistical realisations would be expected to produce such a minimum, or one of even larger magnitude, very frequently. Therefore, we performed a null hypothesis test to estimate the significance of the ℜ_⋆ measurement in the four analysed observed systems (MW, M31, Fornax, Virgo) and in two clusters from the TNG simulations (Fornax-type and Virgo-type clusters). For this, we randomly shuffled the satellite stellar-mass densities, erasing the profile encoded in the satellite mass density distribution $⟨{\overline{\rho }}_{2h\star }^{\star }⟩$. In Fig. 6 we present two examples, each with 200 moving averages of total 10⁶ realisations of the satellite populations, showing data from simulations (Fornax-type), observations (Fornax) and randomised data. Then we measure ℜ_⋆ and the moving average stellar density $⟨{\overline{\rho }}_{2h\star }^{\star }⟩$ at ℜ_⋆ with methods [i] (Eq. (19)) and [ii] (Eq. (20)) in the original data and the randomised data. With this we calculate the probability that a minimum of this significance is randomly generated, shown in Fig. 8 for the Fornax cluster sample. In this exercise, we find that the ℜ_⋆ in Fornax has a p_NH​ = ​4.0% probability of being randomly generated with a minimum density as low as the observed one. Following the same procedure for the other observed systems, we find probabilities of ℜ_⋆ being randomly generated of 0.1% in Virgo, 17% in the MW, and 26% in M31, with the latter estimates being larger than in Virgo and Fornax due to the limited sample statistics. The compound probability that ℜ_⋆ is generated simultaneously at random in four observed systems is 1.8 ⋅ 10⁻⁶, an unlikely scenario corresponding to ∼5σ. The compound probability would become even smaller if we included the probabilities from all clusters and groups from the cosmological galaxy simulation data.

Fig. 8. Null-hypothesis test: Example with the Fornax galaxy cluster. For each of the 106 measurements of ℜ⋆ ([ii]), the corresponding distribution of the minimum values of the moving average curves of the mean stellar mass densities is shown (purple histogram). The distribution of the measurements of ℜ⋆ for the randomised (shuffled) satellite densities (turquoise histogram) is produced by randomly shuffling the density values between the satellites while keeping their spatial positions within Fornax. A more visual comparison is shown in Fig. 6. Using these distributions we calculate that the probability of the randomised distributions resulting in densities as low as or lower than the most probable value measured in Fornax (i.e. ∼ 5.1 × 105 M⊙ kpc−3) is of pNH​ = ​4.0%. We use the same procedure to calculate these probabilities, pNH, for Virgo, MW, M31 and certain simulation examples.

5.3. Satellite sample completeness, mass range, gas, and other effects on the tracers of ℜ_⋆

We test the impact of satellite completeness on the mass density profiles, and on the measurement of ℜ_⋆. For the observations, we tested the satellite samples of the Fornax cluster, which has satellites as faint as 22.7 mag (i-band), and the Virgo cluster, with a satellite sample as faint as 18.7 mag (r-band). We moved our detection threshold to 2 mag higher luminosities and found similar values of ${\mathfrak{R}}_{\star }[ii]={379}_{-9}^{+14}\mathrm{kpc}$ for Fornax, which is only 3% smaller than for the full dataset. Similarly, for simulations we tested increasing the threshold of the minimum total mass (number of particles) of the subhaloes of the Fornax-type cluster by one order of magnitude larger than our fiducial configuration, finding a similar value of ℜ_⋆ within the scatter.

We also explored how different satellite density tracers may change the measurement of ℜ_⋆. For example, in the case of the observations, we tested the Fornax cluster by using the satellite luminosity mean densities derived from their i-band observations, finding a value of ℜ_⋆ similar to the value estimated with the mean stellar mass densities. In the case of the simulations, we find that the value of ℜ_⋆ varies only within the scatter when using different mass density tracers, such as the central dynamical mass mean density (${\overline{\rho }}_{2h\star }$) shown in Fig. 2 (middle panel), or the central mean stellar mass density within one stellar half-mass radius ${\overline{\rho }}_{h\star }^{\star }$. This co-variation of the central mean stellar density measured at different radii (h⋆ or 2h⋆) and with different tracers (light, stellar, and dynamical mass) is a consequence of the half-mass radius variation (Eq. (8) and (13)) and the mass stripping driven by the cluster tidal field. More importantly, these tests show that, regardless of how the tracer is chosen, as long as it is used consistently for the entire satellite sample, it can trace the mass density profile and the location of the transition radius ℜ_⋆.

We also tested whether the satellite mass densities profiles may be driven by mass segregation, with the most massive galaxies located in the centre of clusters. For this we analysed in detail the Fornax-type cluster from the TNG50 simulation, which shows a correlation between stellar mass densities and stellar masses (Fig. 9). The relation shows a large scatter without any systematic bias of finding a preferred satellite mass at a given cluster-centric radius. We calculate the moving average of the stellar masses as a function of distance using Eq. (17), where Fig. 10 (top panel) shows that within the transition radius ℜ_⋆ the stellar masses start to increase on average. To test whether the stellar masses could be driving the mass density profile, we masked out the massive satellites, leaving only satellites with stellar masses smaller than $M_{2h\star }^{\star }={10}^9{\mathrm{M}}_{\odot }$. We find that the increase in the stellar mass moving average in the centre is negligible when the masking is considered (Fig. 10, middle panel). Furthermore, the stellar mass densities of the satellite sample with the masked massive satellites still shows a significant increment of the moving average (Fig. 10, bottom panel), increasing from a minimum of $⟨{\overline{\rho }}_{2h\star }^{\star }⟩\simeq 3\times {10}^5{\mathrm{M}}_{\odot }{\mathrm{kpc}}^{-3}$ near ℜ_⋆ to ∼ 20 × 10⁵ M_⊙ kpc⁻³ within r < 0.2 R_vir. This shows that low-mass satellites still have systematically higher mass densities in the central region of the cluster. Moreover, as the dynamical friction strength depends directly on the satellite’s mass and is inversely proportional to the satellite orbital velocity (f_DF ∝ m_sat/v_sat²) it can drive massive satellites to the central region of the cluster faster than less massive satellites. However, this only changes the orbital distribution of the satellites, while the mean mass density of each satellite continues to be shaped by the tidal field of the host during its orbit, stripping the outer material of the satellite according to Eq. (12).

Fig. 9. Satellite galaxies in the TNG50 Fornax-type cluster, showing their central stellar mass densities as a function of their stellar masses. The colour-coding here shows their cluster-centric distances, which indicate that massive satellites can be found both in the outskirts and in the internal regions of the cluster.

Fig. 10. Satellite galaxies in the TNG50 Fornax-type cluster. Top and Middle panels: Satellite stellar masses and moving average (calculated with Eq. (17)) shown as a function of their cluster-centric distances and colour-coded by their stellar mass densities. In the middle panel, we masked out satellites with $M_{2h\star }^{\star }>{10}^9{\mathrm{M}}_{\odot }$. Bottom panel: Central mean stellar mass densities of the satellites and the moving average as a function of their distances shown for satellites with stellar masses within $M_{2h\star }^{\star }\le {10}^9{\mathrm{M}}_{\odot }$. The transition radius ℜ⋆ is also plotted (vertical blue line), indicating that within ℜ⋆ the mean stellar mass densities still increase despite the removal of massive satellite. The cluster tidal field is also shown (τ). Histograms are plotted in the right sub-panels to show their distributions.

A similar argument can be stated for galaxy harassment, where the cumulative effects of satellite-satellite collisions within a satellite population must obey the tidal shaping of the overall cluster potential, as shown by simulations.

We also estimated whether the gas component in satellite galaxies is significantly affecting the stellar and dark matter distribution. Satellite gas depletion through ram-pressure stripping, tidal stripping, and stellar feedback, can perturb the satellite’s potential, while star formation converts gas into stars, changing the stellar-mass fraction. Therefore, we inspect the radial distribution of the satellite gas masses and densities in the studied environments. For this, we again used data from cosmological galaxy simulations that provide access to a complete satellite sample with all their gas phases. Overall, our analysis of the gas properties in satellite galaxies shows that within the transition radius (r ≲ ℜ_⋆), most satellites have very low gas-to-baryon mass fractions, with only 1.5% of the satellites having gas-to-baryon fractions above 10%, and increasing to 13% within the splashback radius (r ≲ R_sp). Therefore, we find that the impact of the gas mass densities on the stellar and total mass densities of the analysed environments is minimal in the central region (r ≲ 0.5R_vir), but it matters for larger distances (r ≳ R_vir). Moreover, data from the observed environments show that most satellites located in central regions lack significant amounts of gas, as shown in the MW and M31 (Putman et al. 2021), the Fornax cluster (Loni et al. 2021; Kleiner et al. 2023), and in the Virgo cluster (e.g. Chung et al. 2009; Yoon et al. 2017). Further analysis goes beyond the scope of this work, but there is a vast literature related to the detailed quenching of satellite galaxies at different mass regimes (e.g. Gunn et al. 1972; Mori & Burkert 2000; Faerman et al. 2013; Wetzel et al. 2013; Emerick et al. 2016; Donnari et al. 2018; Lotz et al. 2019; Roberts et al. 2019; Tonnesen 2019; Tonnesen & Bryan 2021; Pallero et al. 2022; Herzog et al. 2022; Rodríguez et al. 2022; Blaña et al. 2024).

5.4. Projection effects and sky observables

We explored with simulations how the satellite stellar mass density distribution is affected by projection effects to understand its impact on the moving average satellite density profiles and the transition radius ℜ_⋆. For this, we analysed the Fornax-type cluster from TNG50, where we randomly project the cluster to measured the projected moving average of the satellite stellar mass density profile $⟨{\overline{\rho }}_{2h\star }^{\star }⟩$. As shown in Fig. 11, the projected profile of $⟨{\overline{\rho }}_{2h\star }^{\star }⟩$ shows the same behaviour as in the case measured with the 3D distances r (Fig. 2, bottom panel) reaching increasingly higher values in the central region of the cluster. However, the transition radius ℜ_⋆ is slightly shifted towards smaller cluster-centric distances at ℜ_⋆​ = ​0.42 R_vir (462 kpc; see Table 1). This smaller projected ℜ_⋆ is likely produced by the projected distribution of the satellites, similar to the projected 3D mass profiles (Wolf et al. 2010). Additionally, in the very centre of the cluster (r < 0.2R_vir) the mean mass densities appear lower because of the distant satellites that appear closer in projection (interlopers). This systematic could partially explain the variations in the shape of the observed satellite mass density profiles in Fig. 5 (see Sect. 4.3). In addition, the random orientation of the cluster along the observer’s line of sight can also affect the profile of $⟨{\overline{\rho }}_{2h\star }^{\star }⟩$, especially in clusters that have massive filaments and merging substructures.

Fig. 11. Mean stellar mass densities (${\overline{\rho }}_{2h\star }^{\star }$) of 492 satellite galaxies (circles) as a function of their projected distances to their hosting TNG50 Fornax-type cluster, with colours indicating their stellar masses. The deprojected profile is shown in Fig. 2 (bottom panel). The host tidal field profile (|τ|) is shown with the purple curve. The blue curve shows the moving average profile using Eq. (17). The measurement of the transition radius (ℜ⋆) is performed using the projected profile Rproj (vertical dot-dash line). The histogram in the right sub-panel shows the shape of the satellite mass density distribution.

5.5. Our results and the literature

As covered in Sect. 2.1, the mean mass densities of satellite galaxies are expected to be a function of the distance to the centre of the host system because of its varying tidal field. Studies focused on MW satellites have identified this correlation in the central dynamical mass densities of satellites and their orbital pericentres ($\overline{\rho }-r_{\mathrm{peri}}$; see Kaplinghat et al. 2019; Hayashi et al. 2020; Pace et al. 2022; Cardona-Barrero et al. 2023). Here, we expand on this further, finding that: i) the mean stellar mass densities (and light) of the satellite population also trace this relation at their current positions, giving rise to the stellar-mass density-distance relation (${\overline{\rho }}^{\star }-r$), and relate this to the host’s tidal field; ii) the behaviour of the ${\overline{\rho }}^{\star }-r$ relation changes at the transition radius at ℜ_⋆​ ≃ ​0.5 R_vir where the host tidal field is weak; and iii) these correlations can be identified in systems with a wide range of host masses (10¹²-10^15.5 M_⊙) considering simulations and observations.

Furthermore, to investigate the pristine scaling relations of satellite galaxies, one must account for the effects of the host tidal field and the observed $\overline{\rho }-r$ relation, which reflects how satellites have been transformed since their infall into denser environments. In particular, tidal evolution impacts the mass-size scaling relations such as the fundamental plane in early-type galaxies (Djorgovski & Davis 1987), the (stellar) mass–size relation in early-type galaxies (Kauffmann et al. 2003a,b) and in late-type galaxies (Kormendy & Ho 2013; Saglia et al. 2016), which can also depend on more detailed morphologies of the galaxies (spheroid versus disc-like) (Dressler 1980; Kuchner et al. 2017). Moreover, observational evidence suggests that the environment affects the galaxy mass-size relation (Kuchner et al. 2017; Demers et al. 2019; Matharu et al. 2019; Venhola et al. 2019; Joshi et al. 2020). For example, Chamba & Hayes (2024) report a similar finding to our analysis of the Fornax cluster satellites. They find that satellites within the cluster are systematically smaller than in the field, which they attribute to rapid gas removal in the cluster environment, truncating satellites (see also Asencio et al. 2022). Instead here we find that, with our toy model and simulations, tidal truncation by the host cluster can explain the relation between the satellite spatial distribution and their stellar and total mean mass densities. The effect on the stellar mass densities strongly depends on cluster-centric distance, starting within the transition radius at ℜ_⋆​ ≃ ​0.5 R_vir, and becoming rapidly stronger towards the centre. However, projection effects can increase the scatter in cluster centres at R_proj < 0.1R_vir.

Moreover, the $\overline{\rho }-r$ relation is intimately connected to the morphology-density relation (Dressler 1980; Fasano et al. 2015); both depend on the cluster-centric distances to the most crowded regions. A growing body of work shows that a satellite’s orbital history leaves a clear imprint on its evolution, with the satellite’s present orbital phase being encoded in its phase-space coordinates and stellar population parameters (e.g. Wetzel 2011; Yoon et al. 2017; Rhee et al. 2017; Smith et al. 2022; Montero-Dorta et al. 2024). In particular, the timing of first infall and pericentric passage strongly influences the star formation, mass, and morphological properties of satellites across different regions of the cluster (Rhee et al. 2017; Lotz et al. 2019; Pallero et al. 2019; Smith et al. 2022), a conclusion reinforced by the strong radial colour gradient of satellites found within 300 kpc of the Fornax cluster centre (Eigenthaler et al. 2018; Spavone et al. 2020). Our results indicate that the interplay between the host’s tidal-field topology and the satellites’ orbital histories shapes the observed radial pattern in their mean mass densities.

6. Summary and conclusions

In this study, we investigated how the mean mass densities of satellite galaxies vary with their cluster-centric distances ($\overline{\rho }-r$) and with the strength of the host’s tidal field. To this end, we combined three complementary approaches: an idealised tidal-stripping toy model (Sect. 4.1), satellite populations drawn from cosmological hydrodynamical simulations (Sect. 4.2), and deep panchromatic observational data (Sect. 4.3). Our main findings are the following:

1. The tidal-stripping toy model and the IllustrisTNG simulations consistently show that satellites become denser towards the cluster centre. We find that the total mean mass densities scale as $⟨{\overline{\rho }}_{2h}⟩\propto 1/r$ and $⟨{\overline{\rho }}_{2h}⟩\propto |\tau |$, where |τ| is the tidal field strength (Figs. 1 and 2). Moreover, the same relation holds for the satellites’ central mean stellar mass densities, $({\overline{\rho }}_{2h\star }^{\star })$ out to a well-defined transition radius at ℜ_⋆​​ ≈ ​0.5R_vir. Beyond ℜ_⋆ the host tidal field weakens below the densities of the least-bound satellites, where stripping is negligible. As a result, ${\overline{\rho }}_{2h\star }^{\star }$ flattens or rises slightly towards the virial and splash-back radii (∼1.4R_vir).

2. We examined 296 clusters in the ILLUSTRIS TNG50, TNG100, and TNG300 cosmological galaxy simulation volumes at redshifts z​ = ​0 and 1 (Figs. 4, A.1, A.2). The ${\overline{\rho }}_{2h\star }^{\star }-r$ transition appears in every run, with ℜ_⋆ ∼ 0.5 R_vir in 3D, and somewhat smaller values once projection effects are included. Detailed comparisons with Fornax-type and Virgo-type clusters in both TNG and EAGLE simulations confirm these trends.

3. Four observed satellite galaxy populations – in the Virgo and Fornax galaxy clusters, as well as in the MW and the M31 galaxies – show the same behaviour in their ${\overline{\rho }}_{2h\star }^{\star }-r$ scaling relations (Fig. 5). The measured ℜ_⋆/R_vir match the simulations, although inside r ≲ 0.2R_vir the observed densities are lower. This offset may reflect projection bias, but could also signal a difference in the actual stripping and/or disruption efficiency, which occurs in cluster cores, warranting further study.

4. Our null-hypothesis statistical tests of the observations (Sect. 5.2, Figs. 6, 8) show that the density transition radius ℜ_⋆ is highly significant. The joint probability of obtaining this feature in all the observed environments purely by chance is approximately ∼2​ ⋅ ​10⁻⁶ (≃5σ).

5. Where dynamical masses are available for the Local Group satellites and the simulations, the central dynamical-mass density correlates with both their central stellar-mass density and tidal field strength (Fig. 2 middle panel, Fig. 5 bottom panels). The resulting ℜ_⋆ values agree with those derived from the stellar mass components, reinforcing the connection between the satellite structure and the host’s tidal field strengths.

We find that the distribution of the satellites’ mean mass densities is governed by three interlocking factors acting simultaneously: (1) the tidal field of the host cluster, (2) the satellites’ orbital histories and their cluster-centric distances, and (3) the satellites’ initial mass densities, with this last point revealing that neither mass nor size on their own dictate how tidal forces reshape satellites.

These findings tighten the link between the satellites’ internal properties and external tidal environment, which drive their evolution. Extending the study of the ($\overline{\rho }-r$) relation to loose groups and rich clusters will reveal how the tidal evolution of satellites is interconnected with large-scale processes, such as the baryon-to-dark matter fraction, star formation, and ram pressure, among others. This in turn will refine classic scaling laws, including the morphology-density and the galaxy mass-size relations.

Despite these insights, our analysis still carries limitations and opens several promising avenues for further investigation. Preliminary snapshots from the IllustrisTNG simulations at z​ = ​1 suggest that the radial density profile and the transition radius are in place early. Ongoing work (Blaña et al. in prep.) will pin down when the ${\overline{\rho }}_{2h\star }^{\star }-r$ scaling relation is established. Our current analysis adopts a spherically symmetric host potential. Mapping the full tidal-tensor field will allow us to quantify anisotropic compression, stretching, and shear, and to predict where satellite galaxies and star clusters are most susceptible to stripping.

Because the transition at ℜ_⋆ is encoded in the satellites’ stellar-mass density alone, it can be measured without kinematic data. This provides a novel and practical way to constrain the tidal field and perhaps even the dark-matter distribution in distant systems where spectroscopic observations are unavailable or limited.

Acknowledgments

The authors deeply thank Mary Putman and her collaborators, as well as Suk Kim and collaborators, for making their respective datasets freely available. Matías Blaña also thanks Andreas Burkert, Raúl Angulo, and Claudia Pulsoni for their helpful comments and discussions, and Carolina Agurto, Julio Carballo, Claudio Blaña and Ana Díaz for their support. The authors also thank the anonymous referee for comments that improved the manuscript. The authors thank the TNG collaboration for making their data, software and server capability freely available. M. Blaña thanks ANID for funding through the Fellowship FONDECYT Postdoctorado 2021 No. 3210592. Y. Ordenes-Briceño thanks ANID for the FONDECYT Postdoctorado 2021 No. 3210442 and ESO comité mixto 2024. T. Puzia acknowledges support through FONDECYT Regular No. 1201016, P.B. Tissera acknowledges funding through FONDECYT Regular No. 1200703, and M. Mora acknowledges funding through FONDO CONICYT-GEMINI Posición postdoctoral. E. Johnston acknowledges support from FONDECYT de Iniciación en Investigación 2020 Project No. 11200263. D. Pallero thanks ANID for funding through the Fellowship FONDECYT Postdoctorado 2024 No. 3230379. T. Ziliotto acknowledges funding from the European Union’s Horizon 2020 research and innovation programme under the Marie Skłodowska-Curie Grant Agreement No. 101034319 and from the European Union – NextGenerationEU. We also acknowledge support through ANID CATA-Basal FB210003. The authors deeply thank the citizens of Chile for their tax contributions to the national development of science and this project in these difficult post-pandemic times. The authors extend their gratitude to all the researchers whose studies have been crucial to this work.

References

Appendix A: Complementary material

We present the mean mass density profiles of satellite galaxies as a function of their cluster-centric distances in Fig. A.1, taken from TNG100 for 54 host galaxy clusters from a snapshot at redshift z​ = ​0 (left column) and 17 clusters at z​ = ​1 (right column). In Fig. A.2 we show the same relations for 26 host galaxy clusters/groups from a snapshot at redshift z​ = ​0 (left column) and 12 clusters at z​ = ​1 (right column) taken from TNG50.

As with the TNG300 cluster (Fig. 4 in Sect. 4.2), we see that the mean dynamical-mass densities of the satellites (${\overline{\rho }}_{2h}$) in TNG100 and TNG50 (top panels in Figs. A.1 and A.2) are shaped by the host tidal field, where the low density subhaloes are clearly limited by the tidal field strength according to Eq. 4 (${\overline{\rho }}_{2h}>|\tau |$). Moreover, the mean stellar-mass mean densities also show a correlation with the tidal field strength that extends out to the transition radius at ℜ_⋆ ≃ 0.5R_vir. Some clusters present high-density variations due to infalling substructure around the main trend, as shown by the density moving averages (wiggles in some orange curves in Figs. 4, A.2 and A.1). However, the main trend for the total- and stellar-mass density profiles remains the same across the TNG simulations, as shown by the moving average of the stacked profiles (dashed black curves in all panels in Figs. 4, A.1 and A.2.

Fig. A.1. Mass densities of satellite galaxies vs their cluster centric distances for 54 host galaxy clusters/groups from a snapshot at redshift z​ = ​0 (left column) and for 17 clusters/groups at z​ = ​1 (right column), taken from TNG100 (equivalent analysis for TNG300 in Fig. 4). The tidal field profiles of all the clusters are shown in each panel (|τ| in purple curves). Top panels show the total mass mean densities of the subhaloes (circles, colour-coded by total subhalo masses), and the bottom panels show the central mean stellar mass densities (circles, colour-coded by stellar masses). To avoid overcrowding, we plotted a sub-sample with 10% (left column) and 20% (right column) of the satellites (circles). Circle sizes are proportional to h (top panels) or h⋆ (bottom panels). In each panel the moving average satellite mean mass density profile of each cluster is shown in with an orange curve. In addition, it is shown the global moving average in dashed black curve that is calculated with the moving averages of all clusters. All are moving averages are calculated with Eq. 17. The transition radii (ℜ⋆), measured in each cluster, are located within the vertical grey region (bottom panels), which are determined from the moving average of the mean stellar mass densities $⟨{\overline{\rho }}_{2h\star }^{\star }⟩$, marking where the $⟨{\overline{\rho }}_{2h\star }^{\star }⟩$ start to increase further in. We present a similar analysis of TNG50 in Figs. A.2.

Fig. A.2. Mass densities of satellite galaxies vs their cluster centric distances for 26 host galaxy clusters/groups from a snapshot at redshift z​ = ​0 (left column) and for 12 clusters/groups at z​ = ​1 (right column), taken from TNG50 (equivalent analysis for TNG300 in Fig. 4). The tidal field profiles of all the clusters are shown in each panel (|τ| in purple curves). Top panels show the total mass mean densities of the subhaloes (circles, colour-coded by total subhalo masses), and the bottom panels show the central mean stellar mass densities (circles, colour-coded by stellar masses). To avoid overcrowding, we plotted a sub-sample with 10% (left column) and 20% (right column) of the satellites (circles). Circle sizes are proportional to h (top panels) or h⋆ (bottom panels). In each panel the moving average satellite mean mass density profile of each cluster is shown in with an orange curve. In addition, it is shown the global moving average in dashed black curve that is calculated with the moving averages of all clusters. All are moving averages are calculated with Eq. 17. The transition radii (ℜ⋆), measured in each cluster, are located within the vertical grey region (bottom panels), which are determined from the moving average of the mean stellar mass densities $⟨{\overline{\rho }}_{2h\star }^{\star }⟩$, marking where the $⟨{\overline{\rho }}_{2h\star }^{\star }⟩$ start to increase further in. We present a similar analysis of TNG100 in Fig. A.1.

All Tables

All Figures

Fig. 1. Results of the toy model showing the imprint of the host tidal field on the mean mass density distribution (circles) of nsat​ = ​2 × 103 satellite galaxies shown as a function their distances to the cluster centre. The host has an NFW halo mass and size of a Fornax-type cluster, showing its tidal field (|τ|) in both panels (purple curves). The satellites are taken from a snapshot after a 10 Gyr orbital integration with DELOREAN (Blaña et al. 2020). Left panel: We show the final total (i.e. dynamical) mass mean densities ${\overline{\rho }}_{\mathrm{tr}}$ (large grey circles) and mean stellar mass densities ${\overline{\rho }}_{\mathrm{tr}}^{\star }$ (yellow circles) after 10 Gyr. The moving average profiles of these densities are calculated with Eq. (17) using $q=2\sqrt{n_{\mathrm{sat}}}\simeq 89$ neighbours and are shown with solid thick black and orange curves. The IC of the total and stellar mass densities of the satellites are shown with small points in grey and brown, respectively, showing also their moving average profiles (dashed curves). In addition, we show the effect of including a population of 30% of the first in-falling satellites and the resulting moving average (thin black curve). Right panel: For the same satellites, we show the total and stellar mass mean densities measured in the central regions of each satellite within twice the half-mass radius of the total mass (2h) and the stars (2h⋆), as well as their moving averages (solid curves). The increase in the moving average profiles towards the cluster centre shows that the shrinking of the truncation radius of each satellite (rtr) produces an increasing central mean mass density as a result of the shrinking half-mass radii, h and h⋆. Moreover, the effect is more noticeable in the inner cluster regions (r < 0.5Rvir) where the tidal field is stronger and the satellites have smaller pericentres, making them more susceptible to stronger tidal stripping. The horizontal dot-dashed lines included in both panels are added for reference to illustrate the deviation of the mean stellar mass densities from constant values.

In the text

Fig. 2. Mean mass densities of 492 satellite galaxies (circles) as a function of distance to their hosting TNG50 Fornax-type cluster, with colours indicating their stellar masses (colour bar). The cluster has a virial mass and radius of Mvir = 7.6 × 1013 M⊙ and Rvir​ = ​1075 kpc, respectively, and corresponds to a snapshot at redshift z​ = ​0. Its tidal field radial profile (|τ|) is shown in all panels (purple curve). Shown as a function of cluster-centric distance from the top to bottom panels are: the dynamical mass densities of each satellite ${\overline{\rho }}_{2h}$ (top), measured within 2h; the central mean dynamical mass densities ${\overline{\rho }}_{2h\star }$ within 2h⋆ (middle), and the central stellar mass densities ${\overline{\rho }}_{2h\star }^{\star }$ within 2h⋆ (bottom). The blue curves in each plot show the moving average (Mov. Av.) of different mass densities for all the satellites using Eq. (17). The measurement of the transition radius (ℜ⋆ metric [ii]) is shown with a vertical dot-dashed line (middle and bottom panels). Histograms are shown in the right column sub-panels to show the shape of the mean mass density distributions.

In the text

	Fig. 3. Comparison between the mean central dynamical and stellar mass densities for satellites in the TNG50 Fornax-type cluster. The symbol colour encodes stellar mass, and the identity function is shown as the dot-dashed line.
In the text

Fig. 4. Mass densities of satellite galaxies vs their cluster centric distances for 155 host galaxy clusters from a snapshot at redshift z​ = ​0 (left column) and for 31 clusters at z​ = ​1 (right column), taken from TNG300. The tidal field profiles of all the clusters are shown in each panel (|τ| in purple curves). Top panels show the total mass mean densities of the subhaloes (circles, colour-coded by total subhalo masses), and the bottom panels show the central mean stellar mass densities (circles, colour-coded by stellar masses). To avoid overcrowding, we plotted a sub-sample with 10% (left column) and 20% (right column) of the satellites (circles). Circle sizes are proportional to h (top panels) or h⋆ (bottom panels). In each panel the moving average satellite mean mass density profile of each cluster is shown in with an orange curve. In addition, it is shown the global moving average in dashed black curve that is calculated with the moving averages of all clusters. All are moving averages are calculated with Eq. 17. The transition radii (ℜ⋆), measured in each cluster, are located within the vertical grey region (bottom panels), which are determined from the moving average of the mean stellar mass densities $⟨{\overline{\rho }}_{2h\star }^{\star }⟩$, marking where the $⟨{\overline{\rho }}_{2h\star }^{\star }⟩$ start to increase further in. We present a similar analysis of TNG100 and TNG50 in Figs. A.1 and A.2, respectively.

In the text

Fig. 5. Satellite galaxy mean stellar mass densities (${\overline{\rho }}_{2h\star }^{\star }$ as circles) of four observed systems: the Fornax galaxy cluster (top left panel), the Virgo cluster (top right panel), the MW (bottom left panel), and M31 (bottom right panel), with the satellite stellar masses indicated in colour. In all the panels, we show the tidal field for the host systems (|τ| in purple curves). In each panel we include the moving average profiles of the mean stellar mass densities using Eq. (17) with $q=\sqrt{n_{\mathrm{sat}}}$ (blue curves, with the blue areas casting one standard deviation). The measured transition radii ℜ⋆ are shown for each system (dash-dotted vertical lines), and the expected location of the splashback radii ∼1.4Rvir (vertical dotted lines). For MW and M31, due to the low number of satellites, we add thinner blue curves that show moving averages $⟨{\overline{\rho }}_{2h\star }^{\star }⟩$ where each is calculated with number of neighbours between $q=0.5\sqrt{n_{\mathrm{sat}}}$ and $2\sqrt{n_{\mathrm{sat}}}$. This indicates that, independent of q, all the curves show a low value of $⟨{\overline{\rho }}_{2h\star }^{\star }⟩$ at ℜ⋆. In addition, for MW and M31 we show the central mean dynamical mass densities of each satellite measured within h⋆ (triangles) and the respective moving average profiles (brown curves, bottom panels). This can be compared with other studies in the literature on the $\overline{\rho }-r$ relation for the MW satellites mentioned in Sect. 2, which show that using r instead of the pericentre also reveals an anti-correlation with distance ($\overline{\rho }\propto 1/r$). Moreover, as predicted by the simulations, the central dynamical mass densities show a change in the moving average profile where ℜ⋆ is measured (see Fig. 2, middle panel). Histograms included in the top panels show the overall shape and peak of the distribution of ${\overline{\rho }}_{2h\star }^{\star }$ for the satellites of each cluster.

In the text

Fig. 6. Measurement examples of the transition radius ℜ⋆: the Fornax-type cluster (top panel) from TNG50, and the Fornax cluster (bottom panel). From a total of 106 sub-sample realisations and the 106 moving average curves of the satellite mean stellar mass densities, we show here 200 (blue curves). Each curve is calculated with Eq. (17) with a sub-sample of $n_{\mathrm{sub},\mathrm{sat}}=n_{\mathrm{sat}}-\sqrt{n_{\mathrm{sat}}}$ satellites from a total of nsat. The neighbour number is a random number between $q=2\sqrt{n_{\mathrm{sub}}}$ and $4\sqrt{n_{\mathrm{sub}}}$. The first method of measurement (ℜ⋆[i]) uses Eq. (19) and is shown with vertical dot-dashed blue lines and their 1-sigma confidence region. The second method (ℜ⋆[ii]) uses Eq. (20) and is shown with vertical dot dash orange lines and their 1-sigma confidence region. Bootstrapping reveals a small scatter in the moving average profiles, due to the satellite-to-satellite sub-sample variations and the different number of neighbours used. Lastly, we also show the moving average after first randomising (shuffling) the mean mass densities of the satellites, which erases the signature in the satellite population (see Null-hypothesis test in Sect. 5.2).

In the text

Fig. 7. Distribution of measurements of ℜ⋆, marginalising over the distribution of ${⟨{\overline{\rho }}_{2h\star }^{\star }⟩}_{{\mathfrak{R}}_{\star }}$ (minimum in the moving average of the mean stellar mass densities profiles) using methods [i] (left panel) with Eq. (19) and [ii] (right panel) with Eq. (20). The top four rows of panels show the results for the observed systems (Fornax, Virgo, MW, and M31), and the bottom row panels present the simulated TNG50 Fornax-type cluster, with ℜ⋆ measured in projection (2D) and in 3D.

In the text

Fig. 8. Null-hypothesis test: Example with the Fornax galaxy cluster. For each of the 106 measurements of ℜ⋆ ([ii]), the corresponding distribution of the minimum values of the moving average curves of the mean stellar mass densities is shown (purple histogram). The distribution of the measurements of ℜ⋆ for the randomised (shuffled) satellite densities (turquoise histogram) is produced by randomly shuffling the density values between the satellites while keeping their spatial positions within Fornax. A more visual comparison is shown in Fig. 6. Using these distributions we calculate that the probability of the randomised distributions resulting in densities as low as or lower than the most probable value measured in Fornax (i.e. ∼ 5.1 × 105 M⊙ kpc−3) is of pNH​ = ​4.0%. We use the same procedure to calculate these probabilities, pNH, for Virgo, MW, M31 and certain simulation examples.

In the text

	Fig. 9. Satellite galaxies in the TNG50 Fornax-type cluster, showing their central stellar mass densities as a function of their stellar masses. The colour-coding here shows their cluster-centric distances, which indicate that massive satellites can be found both in the outskirts and in the internal regions of the cluster.
In the text

Fig. 10. Satellite galaxies in the TNG50 Fornax-type cluster. Top and Middle panels: Satellite stellar masses and moving average (calculated with Eq. (17)) shown as a function of their cluster-centric distances and colour-coded by their stellar mass densities. In the middle panel, we masked out satellites with $M_{2h\star }^{\star }>{10}^9{\mathrm{M}}_{\odot }$. Bottom panel: Central mean stellar mass densities of the satellites and the moving average as a function of their distances shown for satellites with stellar masses within $M_{2h\star }^{\star }\le {10}^9{\mathrm{M}}_{\odot }$. The transition radius ℜ⋆ is also plotted (vertical blue line), indicating that within ℜ⋆ the mean stellar mass densities still increase despite the removal of massive satellite. The cluster tidal field is also shown (τ). Histograms are plotted in the right sub-panels to show their distributions.

In the text

Fig. 11. Mean stellar mass densities (${\overline{\rho }}_{2h\star }^{\star }$) of 492 satellite galaxies (circles) as a function of their projected distances to their hosting TNG50 Fornax-type cluster, with colours indicating their stellar masses. The deprojected profile is shown in Fig. 2 (bottom panel). The host tidal field profile (|τ|) is shown with the purple curve. The blue curve shows the moving average profile using Eq. (17). The measurement of the transition radius (ℜ⋆) is performed using the projected profile Rproj (vertical dot-dash line). The histogram in the right sub-panel shows the shape of the satellite mass density distribution.

In the text

Fig. A.1. Mass densities of satellite galaxies vs their cluster centric distances for 54 host galaxy clusters/groups from a snapshot at redshift z​ = ​0 (left column) and for 17 clusters/groups at z​ = ​1 (right column), taken from TNG100 (equivalent analysis for TNG300 in Fig. 4). The tidal field profiles of all the clusters are shown in each panel (|τ| in purple curves). Top panels show the total mass mean densities of the subhaloes (circles, colour-coded by total subhalo masses), and the bottom panels show the central mean stellar mass densities (circles, colour-coded by stellar masses). To avoid overcrowding, we plotted a sub-sample with 10% (left column) and 20% (right column) of the satellites (circles). Circle sizes are proportional to h (top panels) or h⋆ (bottom panels). In each panel the moving average satellite mean mass density profile of each cluster is shown in with an orange curve. In addition, it is shown the global moving average in dashed black curve that is calculated with the moving averages of all clusters. All are moving averages are calculated with Eq. 17. The transition radii (ℜ⋆), measured in each cluster, are located within the vertical grey region (bottom panels), which are determined from the moving average of the mean stellar mass densities $⟨{\overline{\rho }}_{2h\star }^{\star }⟩$, marking where the $⟨{\overline{\rho }}_{2h\star }^{\star }⟩$ start to increase further in. We present a similar analysis of TNG50 in Figs. A.2.

In the text

Fig. A.2. Mass densities of satellite galaxies vs their cluster centric distances for 26 host galaxy clusters/groups from a snapshot at redshift z​ = ​0 (left column) and for 12 clusters/groups at z​ = ​1 (right column), taken from TNG50 (equivalent analysis for TNG300 in Fig. 4). The tidal field profiles of all the clusters are shown in each panel (|τ| in purple curves). Top panels show the total mass mean densities of the subhaloes (circles, colour-coded by total subhalo masses), and the bottom panels show the central mean stellar mass densities (circles, colour-coded by stellar masses). To avoid overcrowding, we plotted a sub-sample with 10% (left column) and 20% (right column) of the satellites (circles). Circle sizes are proportional to h (top panels) or h⋆ (bottom panels). In each panel the moving average satellite mean mass density profile of each cluster is shown in with an orange curve. In addition, it is shown the global moving average in dashed black curve that is calculated with the moving averages of all clusters. All are moving averages are calculated with Eq. 17. The transition radii (ℜ⋆), measured in each cluster, are located within the vertical grey region (bottom panels), which are determined from the moving average of the mean stellar mass densities $⟨{\overline{\rho }}_{2h\star }^{\star }⟩$, marking where the $⟨{\overline{\rho }}_{2h\star }^{\star }⟩$ start to increase further in. We present a similar analysis of TNG100 in Fig. A.1.

In the text