Velocity dispersion of brightest cluster galaxies in cosmological simulations

The preliminary step to classify substructures in cluster halos is to run, in each high-resolution Lagrangian region, the Friends-of-Friends (FoF) algorithm to arrange particles in groups with a linking length of 0.16 the mean interparticle separation. Substructures are then identified by applying the Subfind algorithm (Dolag et al., 2009) on each FoF group. We consider a substructure to be resolved if it includes a minimum of 50 (DM or stellar) particles and if it has a minimum DM mass of $5\times 10^{9}M_{\odot}$ in the Hydro-10x run ($5\times 10^{10}M_{\odot}$ in the Hydro-1x case). Identified substructures that fulfill these criteria are considered as "bona fide" cluster galaxies.

In addition to substructures, Subfind also identifies the central halo and all the stars belonging to it. In its original formulation, Subfind does not distinguish between stars belonging to a well-identified BCG and stars belonging to the ICL.

In order to overcome this limitation of Subfind, we implement two different algorithms to determine the dynamical distinction between ICL and BCG in the stellar envelope of the main halo of each group.

The first one is a modified version of Subfind, as presented by Dolag et al. (2010) that relies on the assumption that the two velocity distributions of stars belonging to the ICL and the BCG have each a Maxwellian shape so that the overall velocity distribution of stars is described by a double Maxwellian,

Figure 1 shows the particle velocity distributions of the stars in the main halo of one cluster ($M_{200}$ $=1.65\times 10^{15}$$M_{\odot}$) at $z=0$. We observe that the single Maxwellian (best-fit $\sigma=1164$ km s^-1), represented by the dashed black line, does not provide a good fit to the particle distribution. On the other hand, when the fitting procedure is attempted with a double Maxwellian (purple dashed line) the agreement is much more evident, and the residuals can be considered negligible. The diffuse ICL is associated with the Maxwellian with the larger velocity dispersion ($\sigma^{\star}_{\textrm{ICL}}=1417$ km s^-1), in contrast, the BCG, having colder dynamics, populates the distribution at lower dispersion ($\sigma_{\textrm{BCG}}^{\star}=685$ km s^-1). Furthermore, we tested that a triple Maxwellian does not improve the results. To assign each stellar particle to either one of the two dynamical components, the algorithm follows an unbinding procedure, by iteratively computing the gravitational potential given by all particles within a certain radius. This identifies two stellar populations, classified as "bound" and "unbound", which are separately fitted by a single Maxwellian distribution. If the best-fit parameters of the initial double Maxwellian coincide with those obtained in the single Maxwellians, then the procedure is completed, otherwise, the radius of the sphere is changed accordingly and the computation is remade. Hence, by measuring differences in the dynamics we spatially separate the diffuse ICL and the stars bound to the BCG. This first technique is applied to the Hydro-1x set of simulations.

The second algorithm is a novel one and exploits machine learning methods to identify the two dynamical components. This method is described in detail in a forthcoming paper (Marini et al. 2021, in preparation). A brief description is provided in Appendix A, where we illustrate a few key concepts on which it is based. The algorithm is a supervised Random Forest classifier (Breiman, 2001) trained on the results provided by Subfind when applied to the clusters in the Hydro-1x set to separate BCG and ICL stellar populations. For each stellar particle in the main halo, the classifier takes as input vector the cluster mass $M_{200}$, the information on the radial distance from the center of the cluster and the particles velocity, while it gives as output the corresponding label (ICL or BCG). The main reasoning behind using these specific features can be found in the technique employed by Subfind, since it is mainly based on an unbinding procedure supported by both particles’ spatial and velocity distribution. To estimate the accuracy of the algorithm, we define the accuracy score between the predicted $\hat{y}$ and true $y$ label as

where N is the number of objects classified and $\delta$ is the Dirac delta function. Employing this definition, the algorithm is estimated to be accurate at $\sim$ 80-90 percent. We apply this technique to the star particles in the main halos of the cluster in the Hydro-10x.

Distinguishing the BCG from the stellar diffuse envelope within clusters is notoriously non-trivial. Although BCGs are usually the brightest (and most massive) galaxies in their cluster, positioned in a privileged position at the bottom of the gravitational potential, many factors contribute to complicating the analysis, especially in observations. In some cases, the central galaxy may not coincide with the brightest galaxy in the cluster, or its position may be shifted with respect to the geometrical center of the cluster halo (De Propris et al., 2021, e.g.). This adds up to the difficulties of defining physical boundaries to the galaxy, given that its stellar population is wrapped in a diffuse stellar component, the so-called ICL, which is gravitationally bound to the cluster potential. A complete discussion on the methods to discern the two is beyond the scope of this paper, we refer the interested reader to Rudick et al. (2011) for a comprehensive discussion.

To account for possible shifts between the BCG center of mass and the center of the halo (identified by Subfind as the position of the particle with the minimum in the gravitational potential), we iteratively compute the stellar center of mass within $0.1R_{500}$, centering it at each iteration on these coordinates and recomputing it within $0.1R_{500}$ from this new center. Avoiding this correction causes the measurement of (non-physically) large stellar velocity dispersions in several clusters. Furthermore, given that our simulations provide us with the full phase-space information on the distribution of particles, we employ for each simulated cluster the projections along the three orthogonal directions.

In the analysis of observed clusters, all the properties of the BCG are measured within a circular aperture centered on the brightest galaxy member in a given band (Sohn et al., 2020, e.g) or on the galaxy closest to the X-ray peak (Loubser et al., 2018). The choice of the aperture radius is somewhat arbitrary in the literature, still masses measured within a fixed physical radius allow us to neatly compare predictions from simulations with observations. We choose to identify all BCG observables as the properties yielded by the star particles in a cylinder long $R_{200}$ centered on the BCG and within a 50 kpc physical radius (we note that by reducing the aperture to 30 kpc we do not expect major differences in the BCG properties out to $z\simeq 1.5$; Ragone-Figueroa et al. (2018)). We select only stars in the main halo, excluding those which Subfind identifies as bound to other substructures, and we measure the physical properties within this circular area. We point out that with this method we include not only the dynamical information on the BCG but also the contribution from the ICL. In the following, we will refer to this definition as “BCG+ICL”.

As anticipated in Section 2.1, we include a second definition of the BCG, based on the dynamical distinction between BCG and ICL in the stellar population of the main halo hosted in galaxy clusters. Although this technique is applicable only in simulations, it provides us with a substantial understanding of the assumptions made when we take the BCG in a fixed aperture instead of retrieving it in a more physically motivated way. For consistency, we compute all physical properties within a cylinder long $R_{200}$, centered on the BCG, and within a 50 kpc aperture, this time including only stars which are dynamically associated with the BCG. We will refer to this second definition as “BCG” or “BCG-only” interchangeably.

Therefore, we define the line-of-sight (LOS) stellar velocity dispersion $\sigma_{\mathrm{BCG}}^{\star}$ as the r.m.s. of the velocity distribution of the star particles (identified as BCG+ICL or BCG-only) within the cylindrical projection. We find the r.m.s. estimator to obtain stellar velocity dispersions to closely agree with the observational data. All the observational datasets employed in this analysis comprise velocity dispersion measurements recovered from the galaxy spectra with the Gauss-Hermite series(Cappellari & Emsellem, 2004).

We note that using the biweight estimator (which, on the other hand, is used to compute the total cluster velocity dispersion) does not guarantee the same level of agreement. Additionally, velocity dispersions are measured within the three orthogonal projections. The BCG stellar mass $M_{\mathrm{BCG}}^{\star}$ is computed as the median of the sum of the selected star particles within the same bi-dimensional aperture.

Observed cluster properties are usually retrieved from the dynamics of galaxy members. In our simulated clusters, bona fide galaxies correspond to gravitationally bound substructures, which we identify through the SubFind algorithm (see Section 2.1). However, given that substructures in simulations sometimes fail in reproducing the phase-space structure of the real galaxy population in clusters (Hirschmann et al., 2016; Marini et al., 2021, e.g.), it is advisable to also employ DM particles to characterize the cluster dynamical properties. These particles have been selected from the main halo of the central group, excluding the contribution from those bound to substructures. In this regard, the dynamics followed by DM particles should be a robust tracer of the cluster potential, unlike substructures that are dependent on the Halo Finder, and possibly also suffering from the effect of dynamical friction. In Section 4 we further comment on the impact of choosing DM particles with respect to galaxies.

Following the observational approach, we select DM particles within a cylinder long $R_{200}$ whose base of radius $R_{200}$ is centered on the cluster center, corresponding to the position of the particle with the minimum value of the gravitational potential. We determine the DM particles kinematics (i.e. the cluster velocity dispersion $\sigma_{200}$) with the biweight estimator. Given that galaxies are most likely affected by fore-/back-ground interlopers in the observational analysis, observers often choose to use the biweight because it underweights the tails of the velocity distribution which are expected to be populated by outliers (Beers et al., 1990). Analogously, we define the cluster mass $M_{\mathrm{200}}$ as the total particle mass within the same cylinder.

Figure 2 illustrates (from left to right panel) $\sigma_{200}$vs $M_{200}$, $M_{\mathrm{BCG}}^{\star}$ vs $M_{200}$ and $\sigma_{\mathrm{BCG}}^{\star}$ vs $M_{\mathrm{BCG}}^{\star}$ from the cluster sample in our simulations at $z=0$ compared to the observed population provided by Sohn et al. (2020). In this regard, we decided to employ the simulations at $z=0$, given that in any case, our results do not significantly change at late times. Legend is common to all panels: dark green and golden points indicate respectively the Hydro-1x and Hydro-10x samples using all the stars in the main halo enclosed in a bi-dimensional aperture of 50 kpc (BCG+ICL sample). We remind that the BCG+ICL sample is the sample selected with an approach closest to the observational one and shall be regarded as the right sample for comparison. The gray circles stand for the BCG measurements studied by Sohn et al. (2020). More in detail, in the left panel, we examine the relationship between the cluster LOS velocity dispersion $\sigma_{200}$and cluster mass $M_{200}$, to understand the capability of our simulations to recover this well-constrained correlation (Evrard et al., 2008; Saro et al., 2013, e.g.). Indeed, once calibrated, the cluster velocity dispersion is an observable extremely sensitive to the cluster mass. We find that DM particles trace a relation fully consistent with both the theoretical expectation, a power-law with slope 1/3 (black dashed line), and the observational data. Given this tight correlation, in the following, we will often refer to “high-mass” or “low-mass” cluster samples as those with high or low-velocity dispersion values respectively. We point out that the rationale of choosing the DM particles as tracers of the internal dynamics of simulated clusters, instead of using the galaxies identified as substructures, lies in the fact that substructures suffer from a velocity bias possibly larger than in real clusters. In this regard, Figure 3 shows the dependence of the velocity bias with halo mass $M_{200}$ in the same sample of clusters. In our simulations, we measure a consistent velocity bias of the substructures with respect to the DM particles of $\sim 10$ percent at all selected masses which is in line with previous measurements (Lau et al., 2009; Munari et al., 2013, e.g.)

As for the central and right panels, we verify whether the BCG properties in simulations have realistic values compared to the observational data. In the central panel, we examine the relation between $M_{\mathrm{BCG}}^{\star}$ and cluster mass $M_{200}$, a relation generally known to be difficult to reproduce in simulations, as BCG stellar masses are sensitive to the details of the feedback processes. We note that Bassini et al. (2020) and Ragone-Figueroa et al. (2018) previously discussed at length these same results on both resolutions, finding that the Hydro-10x set tends to overestimate the stellar mass of the BCG with respect to the observational data. On the contrary, the Hydro-1x clusters have BCG stellar masses consistent with the observed ones. As discussed by Bassini et al. (2020), rather than a resolution effect, the difference in the BCG stellar masses between the Hydro-10x and Hydro-1x Dianoga clusters lies in the different implementation of AGN feedback. Furthermore, in Ragone-Figueroa et al. (2018) the $M_{200}$ vs $M_{\mathrm{BCG}}^{\star}$ relation is compared also to other simulation results, presented by other groups (Bahé et al., 2017; Hahn et al., 2017; Pillepich et al., 2018), yielding comparable results. We note that this discrepancy is larger at the lower masses.

On the other hand, the BCG stellar mass in the Hydro-1x is compatible with the measurements by Sohn et al. (2020). The right panel shows a tension between the $M_{\textrm{BCG}}^{\star}$– $\sigma_{\textrm{BCG}}^{\star}$ scaling relation for observed and simulated clusters, at least for the Hydro-10x set. Among other factors, an overestimate of the BCG stellar mass contributes to a shift towards the right for the scaling relation predicted by our simulations. In the same plot, we note the presence of three isolated points derived from a cluster in a Hydro-10x simulation with $\sigma_{\mathrm{BCG}}^{\star}$$>$ 600 km s^-1. These points correspond to the three orthogonal projections of the LOS velocity dispersion of the same cluster. To understand the origin of this outlier, we analyze its recent evolution. These extreme features can be traced back to a sequence of major mergers at late times ($z\leq 0.5$). By studying the velocity distribution of the star particles associated with the main halo we note that this structure did not yet reach dynamical equilibrium by $z=0$, with the distribution showing two separate peaks: a central one, corresponding to the most massive subhalo, and a smaller one, corresponding to the merging structure. The reason why this is particularly evident in the Hydro-10x case, and not in the Hydro-1x, is twofold. Firstly, it is due to how Subfind assigns particles to one halo or the other at two resolutions. Secondly, as the resolution changes the details of the timing of the merging also slightly change. We point out that this feature had already emerged in Marini et al. (2021), where we computed the scaling relation with the integrated pseudo-entropy and compare it to the one with the velocity dispersion.

We note that the stellar masses of the BCGs and the total cluster masses are not directly observable quantities, instead, they are inferred from other observable quantities. For instance, scaling relations applied on different observables, such as luminosity (Zhang et al., 2011, e.g.), temperature in the X-ray (Vikhlinin et al., 2006) and velocity dispersion (Sohn et al., 2020) allow to obtain a measure of cluster masses. As for stellar masses, they are often derived by comparing synthetic spectral energy distribution models for single galaxies or through spectroscopy. For this reason, it is useful to also provide a comparison between the stellar velocity dispersion of the BCG and the global cluster velocity dispersion, which are both quantities directly measured from spectroscopic observations.

In Figure 4, we show the LOS velocity dispersion of the BCG as a function of the cluster LOS velocity dispersion. Our simulations span the full range of velocity dispersion probed by observations, extending to even smaller $\sigma_{200}$($<$ 200 km s^-1) or, equivalently smaller masses. We expect these low–$\sigma_{200}$galaxies to be isolated early-type galaxies. On the other hand, in the mass range probed by observations ($\sigma_{200}$$>$ 200 km s^-1), at least for the case of the Hydro-10x, simulated structures show a rather good agreement with the observed stellar dynamics with a steepening of the relation on the very high-mass end.

Quite interestingly, the agreement that we find in Figure 4 between simulations and observations is not in line with the conclusions reached by Sohn et al. (2020) and Sohn et al. (2021), who pointed out a significant tension with the results from the hydrodynamical simulations by Dolag et al. (2010). In the latter (and more recently, including the galaxy group mass scale, in Remus et al. 2017), the stellar BCG velocity dispersions are derived from the Maxwellian fits of the velocity distribution of BCG and ICL separately. The authors find a self-similar scaling between both BCG and ICL velocity dispersions with cluster mass which suggests the existence of self-similarity between the two components over the entire mass range probed. However, such a self-similar scaling is not retrieved in the observed cluster set. In the following, we will extensively discuss the main differences between our analysis of simulated galaxy clusters and that presented by Dolag et al. (2010) to understand where these differences specifically lie, in the effort to clarify the origin of the previously reported disagreement between observations and simulations.

Since the comparison in Sohn et al. (2020) is performed between observed velocity dispersions (thus, including ICL) and three-dimensional BCG-only velocity dispersions, as a first step we test whether accounting for projection effects and the inclusion of the ICL contribution affect the results, but these effects turn out to be negligible in our analysis. Therefore, we test whether our cluster sets differ in the virial scaling relations. For this purpose we show the stellar velocity dispersion of BCG (in blue) and ICL (red) as a function of the cluster mass $M_{200}$ in Figure 5. The scaling between velocity dispersion and cluster mass is modeled as

where $\alpha=3$ is value expected from viral equilibrium and it is kept fixed in the analysis by Dolag et al. (2010). In our analysis, we decided to treat $\alpha$ as a fitting parameter, given that a simple comparison by eye between the virial expectation and our dataset, unlike that from Dolag et al. (2010), showed a consistent discrepancy in the slope. We distinguish between results from the Hydro-1x (simple circles) and Hydro-10x (crossed circles) simulations to clarify the dependence on resolution. Additionally, we report the best-fit relations as provided by Dolag et al. (2010) (dashed lines) and as obtained from the analysis of our set of simulated clusters (solid lines). The best-fit curves from the two distinct fits in the scaling $\sigma^{\star}_{\mathrm{ICL}}$– $M_{200}$ are comparable to each other, and our fitting value is very close to the virial expectation, with $1/\alpha\simeq 0.339$. The difference in the normalization is minimal and still within the internal scatter. As for the $\sigma_{\mathrm{BCG}}^{\star}$– $M_{200}$ relation, our BCG set has a significantly different behavior with respect to the best-fit line from Dolag et al. (2010), especially when expanding the dynamical range down to the mass range of galaxy groups. It is precisely this difference that explains the difference with respect to Dolag et al. (2010) and Remus et al. (2017). A higher BCG velocity dispersion at fixed cluster mass (or global velocity dispersion) in the mass range of galaxy groups brings our measurements in agreement with observational results. This difference is most relevant in the low-mass end, $\sim 5\times 10^{13}$ $M_{\odot}$. The slope of our best-fit relation is shallower than the virial value ($1/\alpha\simeq 0.289$) and indicates substantial differences in the mechanisms operating in the BCG accretion with respect to those in the ICL, which traces more the overall virial dynamics of the host halo.

These two effects are consistent with the different dynamical evolution of two stellar components. The ICL is associated with the stellar envelope bound to the cluster potential, and thus we expect this component to dynamically co-evolve with the DM halo, following the self-similar scaling described by the virial theorem. On the other hand, the BCG is expected to form from the dissipative collapse of the gas cooling and fueling star formation within the proto-BCG building blocks, which later assemble into the BCG dry mergers.

Furthermore, we analyzed the correlation $\sigma_{\mathrm{BCG}}^{\star}$– $\sigma^{\star}_{\mathrm{ICL}}$. In their analysis, Dolag et al. (2010) found these two velocity dispersions to be proportional to each other, according to

We expect this relation not to hold for our set of simulated clusters and groups, given that our $\sigma_{\mathrm{BCG}}^{\star}$distribution is different in the low-mass range. Figure 6 presents our datapoints (circles mark the clusters in the Hydro-1x run, while circled crosses are for the Hydro-10x) with the corresponding best-fit curve (dashed line) for a power-law with varying normalization and slope and the best-fit in Dolag et al. (2010) (solid line). In the light of previous results, it is not surprising that the best-fit slope for our set of simulated clusters is different than 1, as implied by Eq. 4, and the power-law is in fact shallower. Interestingly, we color-coded the data points with the ICL fraction over the total stellar mass associated with the central galaxy, given that Remus et al. (2017) found that outliers in their cluster distribution were due to clusters with below 5 percent BCG mass fractions. We find a color gradient along the best-fit curve: data points laying on this line have larger ICL mass fractions as we move towards higher $\sigma_{\mathrm{ICL}}^{\star}$.

In conclusion, we highlighted the origin of the different conclusions reached by Sohn et al. (2020); Sohn et al. (2021) and by our analysis on the comparison between observed and simulated results for the BCG velocity dispersion. Since we expect the dynamics of the BCG to be different from what was observed in Dolag et al. (2010), the measurements of the stellar velocity dispersions within the innermost regions will also be different, as we have seen in Figure 4. This still holds even if the stellar mass of the BCG in our simulations is overestimated when compared to observational results. Hence, while it is clear that the modeling of BCG stellar mass accretion is extremely sensitive to several different factors, yet the "cold" stellar dynamics of the stars in the central regions is well-reproduced at $z=0$ in our "Hydro-10x" simulations.

As a final test, we investigate the effects of excluding the contribution of the ICL to the stellar velocity dispersion of the BCG, according to what is discussed in Section 2.2. This test is based on the BCG sample which, contrary to the BCG+ICL sample, includes only the stars bound to the BCG in a cylinder long $R_{200}$, centered on the BCG, and within a 50 kpc aperture. We stress that this test has the only purpose of displaying the impact of removing the ICL with a dynamically motivated procedure, rather than for a comparison with observational results, which cannot separate the effect of the diffuse stellar component to the velocity dispersion measured in projection. In Figure 7, we show the comparison of this dataset with the BCG+ICL sample from the Hydro-10x simulation only, since we checked that results do not change with the lower resolution settings. Legend is common to all panels: brown crosses mark the BCG sample, while golden points label the BCG+ICL sample. From left to right we plot the $\sigma_{\mathrm{BCG}}^{\star}$– $\sigma_{200}$, the $M_{\textrm{BCG}}^{\star}$– $M_{200}$ and $\sigma_{\mathrm{BCG}}^{\star}$– $\sigma_{200}$ relations respectively.

The plot in the left panel is complementary to Figure 4 and it shows that the velocity dispersions of the "BCG" and "BCG+ICL" cases are different albeit not significantly (within 50 kpc from the center). Similarly, the central panel examines the amount of ICL mass that is still accounted for in the 50 kpc bi-dimensional aperture centered on the central galaxy. Unsurprisingly stellar BCG masses in the BCG-only sample are lower than the BCG+ICL set, given that the former does not take into account the mass contribution from the stellar component assumed to be part of the ICL envelope. This mass difference accounts for a factor $\sim$3. As for the right panel, the stellar velocity dispersion as a function of the stellar mass, we note that the reduction of the stellar mass is only partially compensated by the reduced velocity dispersion, once the ICL contribution to such two quantities is removed. As a consequence, the scaling relation between velocity dispersion and stellar mass has a higher normalization for the "BCG" case. It is striking how the predictions from the virial theorem are still conserved in both cases (a power-law with slope 1/3 as reported by the dashed colored lines).

In conclusion, we find that excluding the ICL component from the 50 kpc aperture on the central galaxy implies lower velocity dispersions, compatible with a cut in dynamically hot stellar population which makes up the ICL, and lower stellar masses.

Having established that the observed $\sigma_{\mathrm{BCG}}^{\star}$– $\sigma_{200}$ relation is correctly predicted by our simulations, we investigate its redshift evolution, similar to what is done for the evolution of the stellar mass of BCGs. We select five different redshifts ($z=2$, 1, 0.5, 0.2, 0.0) from the Hydro-1x and Hydro-10x simulations for a total of 113, 80, 84, 80, and 93 clusters at each redshift. For all clusters, the BCG stellar velocity dispersion is computed in projection by including the ICL contribution (BCG+ICL case) with a 50 kpc aperture. At each redshift, we group clusters in bins of $\sigma_{200}$, with a minimum of 10 clusters per bin. The value of $\sigma_{\mathrm{BCG}}^{\star}$assigned to each bin corresponds to the median value among the clusters belonging to that bin. We report the results of this analysis in Figure 8. Additionally, we plot the standard deviation at $z=0$ as a function of the cluster velocity dispersion with the grey shaded band, noting that the standard deviation is similar at all considered redshifts.

We observe that the $\sigma_{\mathrm{BCG}}^{\star}$– $\sigma_{200}$ relation is almost redshift-independent with no significant variation within its intrinsic scatter. At $z=2$ there is a mild indication for a slightly higher normalization in low-$\sigma_{200}$halos. This result is possibly consistent with the picture of a rapid BCG formation through dry mergers (i.e. gas-poor and negligible star formation), whose assembly anticipates that of the host halo.

As a second test of evolution, we select the most massive FoF group for the eight regions from the Hydro-10x simulation set, and for each of them, we follow the velocity dispersion as a function of redshift up to $z=2$. Figure 9 shows the evolution of $\sigma_{\mathrm{BCG}}^{\star}$vs $\sigma_{200}$ both normalized at $z=0$, with the color-coding indicating the redshift. On average, the cluster velocity dispersion doubles from $z=2$ to $z=0$. The BCG velocity dispersion also increases with time, with most of the evolution occurring at $z>1$. The evolution of the cluster velocity dispersion is rather smooth, although episodes of merger may temporarily cause a sudden increase of $\sigma_{200}$, followed by relaxation. For instance, this is the case for cluster 5 at $z=0.2$, whose $\sigma_{200}$ is higher than at $z=0$. The same can also happen for $\sigma_{\mathrm{BCG}}^{\star}$, especially if these merging events affect the inner regions. Much like the BCG stellar mass accretion, the stellar velocity dispersion has very little change at late times thus confirming that, once formed through dry mergers, it remains almost independent of the cluster mass growth. In other words, we expect the late-time mass accretion in the central regions not to significantly affect the dynamics of the BCG.

Recently Sartoris et al. (2020) determined the full dynamical analysis of Abell S1063 ($R_{200}=$ 2.36 Mpc, as inferred from weak lensing measurements) combining two different tracers: the kinematics of the cluster galaxy members and the stellar velocity dispersion profile of the BCG. For this cluster at $z=0.3$, the observations consisted of an intensive spectroscopic campaign using the VIMOS and MUSE spectrographs at the VLT. With our analysis, we want to verify whether our simulations of clusters of comparable size reproduce the cluster dynamics as traced by the galaxy population and, in the innermost region, by the BCG. For this comparison, we extend the aperture radius to include all the stars in the main halo out to $R_{200}$. This has the advantage to provide us information on the dynamics traced by both the stellar component of the main cluster halo and by the cluster galaxies at large radii.

Figure 10 shows the comparison between observational results from Sartoris et al. (2020) and one of the simulated clusters from the Hydro-10x runs at redshift $z=0$. The latter was selected for its closeness in size to the observed cluster ($R_{200}$$=$ 2.34 Mpc). Solid lines are for the sample obtained from simulations, while the filled and open circles correspond to the stars and galaxies of Abell S1063, respectively. Error bars on observational data points refer to the reported 68 percent uncertainties. The accuracy of our simulated cluster in reproducing the observed data points is quite remarkable, especially for the velocity dispersion profile of stars. In the outer regions, we see that DM and galaxies profiles from the simulations both capture the observed negative gradients, albeit with a slightly lower normalization for $r<1$ Mpc. This difference can be due to several factors. Indeed, we stress that the cluster selected from simulations has been chosen for its size (the most massive in our cluster sample and the closest to the observed cluster mass), thus differences in the dynamics of its components can be ascribed to distinct dynamical states or specific cluster formation histories.

A further step can be taken by extending the comparison to include the eight clusters from the Hydro-10x set of simulations to infer the general characteristics of the distinct tracers when combining all the simulated clusters. The results are shown in Figure 11, with shaded areas showing the 16^thand 84^thpercentile given by the entire cluster sample.

Except for particularly disturbed systems, as it has been discussed in Section 4, BCGs present minimal dynamical differences from cluster to cluster. On the contrary, as we move towards the outskirts, the contribution from the ICL increases. This component feels more the effects of the overall cluster potential, which is more sensitive to the global dynamical state of the cluster and the accretion pattern from the surrounding large-scale structure (Evrard et al., 2008; Saro et al., 2013, e.g.), thus producing a larger scatter in the velocity dispersion profiles. Given that we are especially interested in examining the inner regions, we decided to stack the stellar profiles without any suitably rescaling of their amplitudes. Conversely, DM particles and galaxies velocity dispersion profiles require to be normalized to the cluster radius before stacking (provided that $\sigma(r)\propto r^{-1}$, as discussed in details in Appendix A of Marini et al. 2021). Once we have the stacked median profiles of simulated clusters, we multiply the resulting profile by the virial radius of Abell S1063, $R_{200}$ = 2.36 Mpc, to properly compare to observational results. We see that in the internal regions, DM particles have a larger velocity dispersion than the stellar component, a consequence of the dissipative collapse which forms the BCG stars and determines their "colder" dynamics. As for substructures, since they undergo tidal effects during the merger process and they are disrupted when reaching the cluster center, the profile abruptly falls for $r/R_{200}<0.1$ . Moving away from the cluster center, the three components converge to similar profiles. At large radii, the median profile spans values of the velocity dispersion which are lower than those in the observed cluster, similarly to what already encountered in Figure 10. Differences at this level can be again ascribed to several factors which may include distinct dynamical state and/or projection effects. At these scales, we find the substructures velocity dispersion profile to be consistently higher than that from the DM case. Analogously to what discussed in Section 4, we observe a velocity bias between the two tracers due to the effect of tidal stripping which is more effective in substructures with relatively low orbital velocity. The resulting mechanism is the selective removal of lower-velocity substructures, and thus an increase in the total velocity dispersion profile as traced by substructures.

Velocity dispersion profiles of the BCGs relate to the dynamical mass profiles of clusters and are expected also to be determined by the formation history of such extremely massive galaxies. Therefore, it is important to investigate and deepen the study on the velocity dispersion profile and its gradient. Up to now, studies of velocity dispersion profiles of BCGs exhibit a large variety of slopes, with a significantly larger fraction of positive slopes for BCGs, with respect to other early-type galaxies and brightest group galaxies (BGGs) (Von Der Linden et al., 2007; Bernardi, 2009; Huertas-Company et al., 2013, e.g.).

We present here the analysis of the slopes of the velocity dispersion profiles of BCGs in our set of "Hydro-10x" simulated clusters and compare with observational results obtained by Loubser et al. (2018). Following these authors, we adopt the following model to describe the projected velocity dispersion profiles of BCGs:

Here the central velocity dispersion $\sigma_{0}$ is measured within an aperture $r_{0}$ from the BCG center. Loubser et al. (2018) (and Newman et al. 2013 therein) obtained the BCG velocity dispersion within a radius $r_{lim}=15$ kpc and apply a central aperture $r_{0}=$ 5 kpc. This choice is motivated by previous studies (Graham et al., 1996) which measured a typical half-light radius $R_{e}$ of 16.7 kpc in 119 Abell clusters.

As for the analysis of simulations, we decided to follow two different approaches. In one case, for the comparison with observational results, we select a bi-dimensional radial aperture $r_{lim}=$ 15 kpc and compute the profiles on the Hydro-10x runs only. We exclude from this analysis the Hydro-1x clusters since their force resolution does not allow to adequately resolve the scales involved in the observational analyses. In the second case, we let $r_{lim}$ to vary with the size of cluster by taking $r_{lim}=0.1R_{500}$: this value proved to include the BCG self-consistently with cluster size (Ragone-Figueroa et al., 2018, and references therein). In this case, the central aperture is also chosen to scale with the cluster radius, while being significantly larger than the softening length of both simulations with $r_{0}=0.05R_{500}$. While this analysis can not be directly compared with observational results, it allows us to quantify the role of the ICL in determining the profiles over fixed fractions of the characteristic cluster scale radius.

The velocity dispersion slope $\eta$ as a function of the central velocity dispersion $\sigma_{0}$ (or cluster velocity dispersion $\sigma_{200}$) is reported in the top left (right) panel of Figure 12. The legend is common to all panels: in golden the Hydro-10x, in dark green the BCG+ICL sample of the Hydro-1x and the BCG sample with light green crosses; maroon triangles and turquoise circles are the Newman et al. (2013) and Loubser et al. (2018) observational sets respectively. The top panels show the results of the analysis with $r_{lim}=$ 15 kpc, while the bottom panels refer to the analysis based on the different cut on the radial aperture $r_{lim}$. In the former, simulations predict slightly positive slopes, a result generally consistent with observational results from Loubser et al. (2018) and Newman et al. (2013). We also note that our set of simulated BCGs shows few cases of negative slopes that often correspond to distinct projections of the same clusters having a three-dimensional negative radial gradient. These clusters do not seem to yield any distinct features with respect to the other clusters in the simulated set, except for the different gradients.

In the bottom panels, we report the results from all runs (i.e. Hydro-1x and Hydro-10x, for both BCG+ICL and BCG cases), given that they should not be contaminated by resolution effects for the chosen radial range over which profiles are analyzed. On the left, we show the slope $\eta$ as a function of the central velocity dispersion $\sigma_{0}$: allowing $r_{0}$ to scale with $R_{500}$ yields a different sampling of $\sigma_{0}$ that extends to larger values than for the observational analysis shown in the upper panels. Furthermore, we observe that the BCG-only sample generally yields lower values of the gradients with respect to their counterpart BCG+ICL. Similar results are obtained in the right bottom panel which analyzes the relationship between $\eta$ the total cluster velocity dispersion $\sigma_{200}$. The generally larger values of $\eta$ at fixed $\sigma_{0}$ (or $\sigma_{200}$) when ICL stars are included is due to the contribution of this dynamically warmer stellar component that increases the velocity dispersion when measured along the line-of-sight. At the same time, removing the ICL contribution still allows finding both positive and negative $\eta$, thus indicating that positive gradients are not necessarily due to the larger velocity dispersion of the ICL, but also reflect intrinsic dynamical properties of the BCGs. Additionally, here both Hydro-1x and Hydro-10x distributions hint at a positive correlation between $\eta$ and $\sigma_{200}$, suggesting that more massive halos have increasing BCG velocity dispersion slopes. This correlation is confirmed by the Spearman coefficient $R_{s}$ computed on 100 bootstrapped samples of the simulated data which gives $R_{s}\simeq 0.5$ and p-value $\sim 0.02$. To help the eye, we report the best-fit lines for each data sample separately, following the color legend used with the points. A steeper slope in larger systems can be interpreted in light of Figure 8, which shows that BCG velocity dispersion increases less than linearly with the overall cluster velocity dispersion. This implies that in more massive clusters the relative contribution of the ICL to the observed $\sigma_{\textrm{BCG}}^{\star}$ increases. As a consequence, we expect this to translate into a larger contribution of the ICL to the projected velocity dispersion as we move toward larger projected radii, thus producing steeper slopes of $\sigma_{\textrm{BCG}}^{\star}(r)$. A similar argument is proposed by Loubser et al. (2018) when including a sample of brightest group galaxies in their analysis. These smaller galaxies – measured within an aperture of 10 kpc and $r_{0}=$1 kpc which cannot be thoroughly compared with our simulations due to resolution limits – exhibit a similar gradient. This is also in line with other observational studies, such as Veale et al. (2017) who investigated the correlation between velocity dispersion profile slopes with galaxy environment and luminosity.