© ESO 2022

1. Introduction

The concordance, or Λ cold dark matter (CDM), cosmological model, dominated by CDM and with a cosmological constant Λ, allows for specific predictions about the physical properties of dark matter (DM) haloes at different scales. Cold dark matter drives the hierarchical growth of structures in the Universe: the most massive DM haloes form from mergers of smaller objects and by accretion (e.g. Tormen 1997; Moore et al. 1999). Many of the less massive haloes are therefore found orbiting around the largest structures: they are usually referred to as sub-haloes or the substructures of a main halo. N-body cosmological simulations can be used to predict the number of sub-haloes around a halo of a given mass. They also show that DM haloes of any mass have a similar mass density profile, with a central cusp, usually described using the Navarro-Frenk-White (Navarro et al. 1996, 1997) or Einasto (Einasto 1965) profiles (e.g. Wang et al. 2020). Any significant discrepancy between these predictions and the observations may imply that the formation of structures does not proceed as forecast by the ΛCDM model.

Galaxy clusters are especially well suited for testing the predictions of cosmological simulations. Around 85 − 90% of their mass is made up of DM located in cluster-scale haloes and galaxy-scale haloes or sub-haloes. The remaining 10 − 15% is composed of baryons: hot X-ray-emitting plasma and stars, which represent only approximately 1% of the total mass (Annunziatella et al. 2017; Sartoris et al. 2020). The entire baryonic mass budget of clusters can be measured from photometric and spectroscopic observations in the near-IR, visible, and X-ray bands. Thus, if one determines the total mass distribution, one can disentangle that of the DM haloes and compare it with simulations (e.g. Diemand & Moore 2011).

In the past few years, gravitational lensing has become one of the most effective techniques for studying the total mass distribution of clusters (e.g. Natarajan & Kneib 1997) as it does not discriminate between baryons and DM. Strong lensing (SL) is especially effective in the central regions, or cores, of massive clusters (i.e. within approximately 1 arcmin, or approximately 350 kpc at z = 0.5), where several multiple images of background sources are usually found. In these regions, two important discrepancies between the predictions of N-body simulations and lensing models have emerged. First, several SL models suggest that the main DM haloes of clusters might have central flat cores (Sand et al. 2004), extended up to around 100 kpc (Grillo et al. 2015; Annunziatella et al. 2017; Bergamini et al. 2019, hereafter B19). This result seems to be confirmed on a homogenous and unbiased sample of massive lens galaxy clusters (Postman et al. 2012). A second discrepancy appears when comparing the amount and distribution of DM substructures (or sub-haloes). As found in Grillo et al. (2015), Munari et al. (2016), and Bonamigo et al. (2018), N-body simulations underestimate the number of massive substructures near the cores of massive clusters. Furthermore, the number of predicted galaxy-scale SL events seems to imply that the mass distributions of these substructures are more concentrated than suggested by simulations (Meneghetti et al. 2020, hereafter M20).

Several Hubble Space Telescope (HST) programmes have recently been dedicated to the observation of SL phenomena caused by massive clusters. For instance, the Cluster Lensing And Supernova survey with Hubble (CLASH; Postman et al. 2012), the Hubble Frontier Fields (HFF; Lotz et al. 2017) campaign, and the Reionization Lensing Cluster Survey (RELICS; Coe et al. 2019) obtained multi-band HST images with an unprecedented depth for 25, 6, and 41 clusters, respectively. Moreover, spectroscopic campaigns such as CLASH-VLT (Dark Matter Mass Distributions of Hubble Treasury Clusters and the Foundations of Λ CDM Structure Formation Models, PI: P. Rosati; see Rosati et al. 2014), complemented with data from the new integral field Multi Unit Spectroscopic Explorer (MUSE; Bacon et al. 2010) on the Very Large Telescope (VLT), have permitted secure identifications of cluster members and precise redshift measurements for hundreds of multiply lensed sources. These data have allowed for the construction of detailed SL models, whose accuracy can be assessed by considering the distance between the observed positions of multiple images and those predicted by a model. The root mean square (rms) of this discrepancy can be as low as 0.3″ (e.g. B19, Caminha et al. 2019); however, this is still higher than the observational uncertainty on the position of the images, which can be less than a tenth of an arcsecond. This might suggest that the total mass modelling of galaxy clusters needs to be further improved (see Meneghetti et al. 2017; Remolina González et al. 2018; Raney et al. 2020).

Strong lensing observations permit a very effective determination of the total mass enclosed within the multiple images of a single background source, but the same amount of mass can be obtained with relatively different mass distribution models (e.g. Limousin et al. 2016; Ghosh et al. 2021). This implies that the parameters that define a model often have a certain degree of degeneracy between them, which corresponds to a possible transfer of mass between the various components. Specifically, galaxy clusters are modelled with extended, cluster-scale mass clumps, to describe the main DM haloes (Natarajan & Kneib 1997) and the intra-cluster medium (ICM), and with smaller haloes, to represent the member galaxies (De Lucia et al. 2004). This last component is generally parametrised in SL models with spherical, truncated, isothermal profiles, defined only by their values of effective velocity dispersion and truncation radius, in order to reduce the number of optimised free parameters. Furthermore, their values are usually linked to the luminosity of each cluster member with two power-law scaling relations (e.g. Richard et al. 2014; Grillo et al. 2015; Monna et al. 2017). This approach is equivalent to using a Faber–Jackson scaling law (Faber & Jackson 1976) to determine the value of the velocity dispersion of an elliptical galaxy.

Information on the kinematics of galaxies can significantly reduce the degeneracy between the parameters of an SL model. Observations of galaxy-scale SL events by the Sloan Lens ACS (SLACS) survey allowed for a precise determination of the slope of the total mass density profile of massive lens galaxies (Koopmans et al. 2006). In recent years, velocity dispersion measurements have also become available for several members in some galaxy clusters. Monna et al. (2015) used spectrography from the Hectospec fibre spectrograph at the MMT and early CLASH data for 21 members of the cluster Abell 383, showing a significant reduction in the degeneracy between the SL model parameters. B19 refined and significantly extended this method to three CLASH clusters: MACS J1206.2−0847, MACS J0416.1−2403, and Abell S1063. The integral field data provided by MUSE permitted velocity dispersion measurements to be obtained for 40−60 galaxies for each cluster. These data have been exploited to estimate the values of the slope and normalisation of the Faber–Jackson relation for those members. This has, in turn, been used as a prior to determine the exponent of the power-law relation that binds the total mass of sub-haloes to their luminosity. The new data and procedure have allowed for a more accurate mass assignment to the diffuse and sub-halo components, reducing the degeneracy between their parameters. This is crucial for both inferring the DM mass density profile in the central regions of a cluster and precisely determining its sub-halo mass function.

This technique has enabled a first step towards a more accurate mass modelling of the cluster members, but a single power-law relation between their total mass and luminosity may be inaccurate on the very large total mass range of member galaxies (10⁹ − 10¹³ M_⊙). In order to improve on the adopted scaling relations for members, in this work we consider the so-called Fundamental Plane (FP) relation (Dressler et al. 1987; Djorgovski & Davis 1987; Bender et al. 1992), which has a lower scatter than the Faber–Jackson law. We focus on one of the three clusters considered by B19, namely Abell S1063 (AS1063). We chose it due to its rather regular morphology and because we have at our disposal line-of-sight velocity dispersion measurements for a sizable subset of cluster members, for which we also have effective radius and magnitude estimates. We can therefore use these galaxies to obtain the best-fit values of the parameters that define the FP.

The paper is organised as follows. In Sect. 2 we give details on our photometric and spectroscopic data and on how they were analysed and reduced for our objectives. In Sect. 3 we describe our SL model built with the FP, while in Sect. 4 we present the analysis and discussion of the results of the optimisation and the comparison with the previous model from B19. In Sect. 5 we show the cumulative, projected mass profiles for all the baryonic and DM components of our model. In Sect. 6 we compare our results with the predictions of stellar-to-halo mass relations derived from N-body simulations. In Sect. 7 we compare the stellar-to-total mass fraction and the compactness of the sub-haloes in our model with those suggested by recent high-resolution hydrodynamical simulations. Finally, in Sect. 8 we summarise our most important conclusions.

In this work, we use a flat ΛCDM cosmology with Ω_m = 0.3 and H₀ = 70 km s⁻¹ Mpc⁻¹, in which 1″ corresponds to a scale of 4.92 kpc at z = 0.348, the redshift of AS1063. All magnitudes are expressed in the AB system.

2. Abell S1063: Data and analysis

Abell S1063, first identified in Abell et al. (1989), is a massive galaxy cluster at redshift z = 0.348. A virial mass value of (2.17 ± 0.06)×10¹⁵ M_⊙ was estimated by Sartoris et al. (2020). Despite its regular shape, it appears to have recently undergone an off-axis merger, as inferred by Gómez et al. (2012) and Rahaman et al. (2021) from X-ray data, and by Mercurio et al. (2021) from the dynamics of the members. Its high mass and redshift make it a very effective gravitational lens, which allows us to use the observed position of 55 spectroscopically confirmed multiple images from 20 background sources, distributed over a wide range of redshifts (z = 0.73 − 6.11; see Caminha et al. 2016), to build an SL model. This high number of spectrocopically confirmed sources is crucial for an accurate SL model. Indeed, the well-known degeneracy between the total mass of a lens and the redshift of a multiply imaged source can be broken when several multiply imaged sources in a wide redshift range are available and their redshift are measured.

The cluster is also a very bright X-ray source (L_X ≈ 2.5 × 10⁴⁵ erg s⁻¹ in the 0.1 − 2.4 keV band), due to the emission of the hot (T ≈ 12 keV) ICM (Rahaman et al. 2021): it was identified as RXJ2248.7−4431 in the ROSAT All-Sky Survey (De Grandi et al. 1999; Guzzo et al. 1999).

2.1. Photometric data

Abell S1063 has been included in the two main photometric campaigns by the HST that, in the last decade, targeted SL phenomena in galaxy clusters. Firstly, it was one of the 25 clusters included in CLASH. Both the Advanced Camera for Surveys (ACS) and the Wide Field Camera 3 (WFC3) instruments were used to observe the clusters with 16 different broadband filters, from the near-UV to the near-IR. An even deeper view on AS1063 was obtained thanks to the HFF programme, which was assigned 840 HST orbits to target a selected sample of only six clusters, in seven different bands (HST/ACS F435W, F606W, and F814W and HST/WFC3 F105W, F125W, F140W, and F160W). Following Tortorelli et al. (2018), we used the HFF photometric data in the F814W band, which corresponds to a rest-frame R-band for the cluster redshift, as it allows for a good signal-to-noise ratio down to magnitudes of around 28 and it has a lower full width at half maximum of the point spread function when compared to the near-IR bands.

In order to obtain the parameters of the FP and consequently derive the values of the velocity dispersions of the cluster members, we need robust measurements of their magnitudes and effective radii. We adopt the prescriptions developed in Tortorelli et al. (2018) and further enhanced in an upcoming work from the same authors to analyse the 222 galaxies included in our model. The galaxy surface brightness profile fitting consists of iterative steps, where we consecutively increase the size of the fitted images. We started by measuring the photometric properties with SExtractor (Bertin & Arnouts 1996). These properties were used as initial guesses for the more accurate galaxy light profile analysis with GALFIT (Peng et al. 2010). We cut stamps around each galaxy, fitting its light profile and those of the closest galaxies. Then, the best-fitting parameter values were used as starting points for the following iterations, where we simultaneously fitted all galaxies falling in the nine regions we separated our image into. Finally, the parameters obtained in the previous step were used as a starting point for the final simultaneous light profile fitting of all galaxies in the whole image.

The surface brightness profiles of the galaxies were fitted with GALFIT using single Sérsic light profiles with seven free parameters (the two centre coordinates, major-to-minor axis ratio, position angle, major-axis effective radius, Sérsic index, and magnitude). For every galaxy, GALFIT was provided with different combinations of point spread functions, background value estimates, and noise images, as described in Tortorelli et al. (2018). This allowed us to average out the contribution of systematic effects in the surface brightness fit.

The spectral energy distribution obtained from the multi-band HST photometry also permits us to infer the stellar mass of the member galaxies, using composite stellar population models. Further details are provided in Mercurio et al. (2021), where the authors used the HST photometry plus the MUSE spectra to fit the spectral energy distributions of the members with a delayed exponential star formation history, solar metallicity, and a Salpeter (Salpeter 1955) stellar initial mass function (IMF). The presence of dust was taken into account following Calzetti et al. (2000).

However, MUSE data are only available for a fraction of our sample. For those galaxies we compare their stellar mass values obtained including information from spectra with those derived without. We notice only a small systematic offset between the two estimates. Therefore, we chose to consider, for all the members, the stellar masses obtained by using only the photometric data so as to have a uniform sample. From this comparison, we also fixed a conservative error of 0.22 on the decimal logarithmic value of the stellar mass.

2.2. Spectroscopic data

Abell S1063 was targeted by several spectroscopic programmes. It was one of the 13 clusters included in the CLASH-VLT survey. A total of 16 masks were used on a field of approximately 25 arcmin across, allowing for the measurement of 3607 reliable redshifts, 1109 of which identify cluster members (Mercurio et al. 2021).

The integral field spectrograph MUSE has complemented the Visible Multi Object Spectrograph (VIMOS) observations in cluster cores. Abell S1063 was observed with two different pointings, centred in the NE and SW regions of the cluster (Karman et al. 2015, 2017; Caminha et al. 2016), with an exposure time of 3.1 and 4.8 h and a seeing of 1.1″ and 0.9″, respectively. MUSE provided a total of 175 additional reliable redshifts, 104 of which refer to cluster members. Following Caminha et al. (2016), we identify as spectroscopic cluster members those galaxies whose redshift falls in the range 0.336 − 0.362, corresponding to a rest-frame velocity range of ±3000 km s⁻¹. Moreover, multi-band CLASH data permit photometric members to be identified with the method presented in Grillo et al. (2015).

The MUSE data are deep enough for the measurement of the line-of-sight stellar velocity dispersion of a considerable fraction of the cluster members (80 out of 222). As in B19, we used the pPXF software (penalised pixel-fitting) by Cappellari & Emsellem (2004) for the measurement, but instead of using a fixed circular aperture of 0.8″ for all members, we chose a wider one, of 2″, and weigh each pixel with its corresponding surface brightness value: this choice has been found to provide velocity dispersion values equivalent to the central ones. The aperture was reduced to 0.5″ for a member whose velocity dispersion measurement was strongly influenced by the effects of the diffuse light of the brightest cluster galaxy (BCG). pPXF estimates the line-of-sight-velocity-distribution (LOSVD) parameters by comparing the observed spectra with a set of stellar templates from the MILES library (Vazdekis et al. 2010) convolved with a LOSVD. The code minimises a χ² function between the measured spectrum and a model.

3. Building the strong lensing model

We built our SL model using that from B19 as a basis and taking advantage of the information provided by the FP relation for a more accurate determination of the total mass of the cluster members. We then compared our model with the B19 model to understand how the new modelling choices regarding the mass of the sub-haloes influence the reconstruction of the diffuse mass component of the cluster.

As in B19, we started from the same 55 spectroscopically confirmed multiple images produced by 20 sources presented in Caminha et al. (2016). Likewise, we modelled the various lens mass components with dual pseudo-isothermal elliptical (dPIE) mass density profiles (Limousin et al. 2005; Elíasdóttir et al. 2007), which are the ellipsoidal generalisation of a truncated isothermal sphere with a central core. As in B19, we used the publicly available SL code LensTool (Kneib et al. 1996; Jullo et al. 2007; Jullo & Kneib 2009) in order to find the best-fit values of the parameters of our model. The dPIE profile is defined by seven free parameters: two for the centre position, the ellipticity e and its position angle θ_e, the central velocity dispersion σ₀, and the core and truncation radii r_c and r_t, respectively. The projected mass density of a spherical dPIE as a function of the two-dimensional distance from the centre R is (Elíasdóttir et al. 2007)

$$\begin{array}{c}\hfill \mathrm{\Sigma }(R)=\frac{{\sigma }_0^2}{2G}\frac{r_{\mathrm{t}}}{r_{\mathrm{t}}-r_{\mathrm{c}}}(\frac{1}{\sqrt{r_{\mathrm{c}}^2+R^2}}-\frac{1}{\sqrt{r_{\mathrm{t}}^2+R^2}})·\end{array}$$(1)

LensTool actually parametrises the dPIE profile with a fiducial velocity dispersion that is related to the central velocity dispersion by ${\sigma }_{\mathrm{LT}}=\sqrt{2/3}{\sigma }_0$ (see Elíasdóttir et al. 2007; Agnello et al. 2014, and Appendix C in B19). As the projected mass of a sub-halo within a given aperture depends on its truncation radius and central velocity dispersion, the same enclosed mass can be obtained with different combinations of values for these two parameters. As shown by B19, the introduction of a kinematic prior significantly reduces this degeneracy, allowing one to obtain total mass values that are consistent with the kinematic data. The diffuse, cluster-scale, DM haloes were represented with an elliptical dPIE, close to the BCG position, and a spherical dPIE with vanishing core radius, whose position was left free to vary on an area of 150″ × 120″ centred in the NE region of the cluster. This second halo is necessary to reproduce some arcs in the NE region and its presence is compatible with the post-merger scenario of the cluster. As in B19, the hot gas mass component is completely fixed by the Chandra X-ray data, reduced in Bonamigo et al. (2018). This component was modelled with three dPIE mass profiles, whose parameters were not optimised. The position of each multiple image and source is identified by two coordinates. As the positions of the sources are not known, 55 multiple images and 20 sources correspond to, respectively, 110 observables and 40 free parameters. The two cluster-scale DM clumps have nine further parameters: the values of the coordinates of the centre and of the velocity dispersions were left free for both haloes, while those of the core radius, ellipticity and orientation angle were only optimised for the main halo. There is therefore a total of 49 free parameters, leading to 61 degrees of freedom.

As anticipated, the main difference with the model obtained in B19 is in the way we assign mass to the cluster members. They were modelled as spherical dPIE profiles with a negligible core radius, so that their total mass is entirely fixed by their values of velocity dispersion and truncation radius. Their three-dimensional total mass density ρ(r) scales as r⁻² inside the truncation radius, and r⁻⁴ outside. Within a sphere with radius equal to the truncation radius, 50% of the total mass is enclosed, while approximately 90% of it is enclosed within a six times larger radius. In order to estimate the velocity dispersion of member galaxies from their measured effective radius and surface brightness, B19 calibrated, for a subset of the members, the Faber–Jackson relation between the observed values of luminosity and line-of-sight velocity dispersion. The resulting relation was used to impose a prior on the exponent of an optimised power-law scaling relation between the two observables, valid for all members. With this approach, a fixed law, with no scatter, is imposed to describe members in a wide mass range. Furthermore, the total mass of members only depends on a single measured quantity, their total luminosity, and might therefore be more sensitive to observational offsets or errors.

The Faber–Jackson law has a scatter of around 0.1 on the logarithmic value of the velocity dispersion (Davies et al. 1983), and the residuals show a significant correlation with galaxy size (e.g. Fraix-Burnet et al. 2010; Nigoche-Netro et al. 2011). This indicates that there is a tighter relation, which includes the effective radius of galaxies, known as the FP (Dressler et al. 1987; Djorgovski & Davis 1987; Bender et al. 1992), of which the Faber–Jackson law is just a projection. The fact that the FP deviates from the expectations for homologous galaxies in virial equilibrium can be explained with the hypothesis of a variation in the total mass-to-light ratio of galaxies with their total mass (Faber et al. 1987), or of a non-homologous structure of the galaxies along the FP (Djorgovski 1995; Hjorth & Madsen 1995; Desmond & Wechsler 2017).

As anticipated, in our model we used the FP relation to completely fix the value of the velocity dispersion of the members from their measured effective radius and total magnitude values. This procedure has several advantages. First, the FP has a lower scatter of 0.07 on the logarithmic value of the line-of-sight velocity dispersion (Jorgensen et al. 1996). The best-fit parameters of the relation were determined with an optimisation performed on a selected sample of members, to avoid the effects of systematics affecting the observations. Furthermore, the value of the velocity dispersion of a member depends on two observed quantities, instead of one. This method is therefore less sensitive to possible measurement errors on a specific quantity. Late-type galaxies were included with early-type galaxies and modelled in the same way. We are aware that this is an approximation, but, as clearly visible in Fig. 1, late-type galaxies are a very small fraction of the members located in the cluster central regions. Moreover, the most massive cluster members, providing the most important contribution to the SL model, are early-type galaxies.

Fig. 1. Colour composite image of the core of AS1063 created from HFF photometry. Blue, green, and red channels are obtained by combining the filters F435W, F606W + F814W, and F105W + F125W + F140W + F160W, respectively.

As for the truncation radii of the members, we chose them to be proportional to their observed half-light radii. This ansatz is suggested by the abundance matching results of Kravtsov (2013), who found that the virial and half-light radii of galaxies are approximately proportional, as previously suggested by some theoretical studies (e.g. Mo et al. 1998). The ratio between the two mostly depends on the angular momentum acquired from tidal torque by the halo during its formation. The relation was observed for galaxies within a very wide range of stellar masses and morphologies, including elliptical galaxies in the Virgo cluster. In our case, rather than the virial radii of members, we are interested in estimating their truncation radii, which are significantly smaller. Galaxies in the cluster core might have a generally self-similar total mass profile, leading to an approximately constant ratio between truncation and virial radii, as a consequence of similar evolution processes. This leads to an approximately proportional relation between effective and truncation radius, which is the ansatz we used in this paper. Compared to previously adopted power-law relations between the truncation radius and the total luminosity of a member, our assumption relates in a more direct way two spatial scales of the members. We calibrated the proportionality relation by testing several values for the truncation-to-effective radius ratio, optimising all the free parameters describing the other mass components of the cluster each time. Finally, we chose the ratio value that allows for the lowest rms value of the discrepancy between the observed and model-predicted positions of the multiple images considered in the optimisation of the model.

3.1. The Fundamental Plane of Abell S1063

The FP is an empirical relation between the logarithm of the effective radius (i.e. the bi-dimensional radius inside which half of the total luminosity is emitted) R_e, the average surface brightness within it, SB_e, and the logarithm of the central stellar velocity dispersion σ₀. The FP can be written as

$$\begin{array}{c}\hfill log{R}_{\mathrm{e}}=\alpha log{\sigma }_0+\beta {\mathrm{SB}}_{\mathrm{e}}+\gamma .\end{array}$$(2)

We consider R_e as measured in kpc, and σ₀ in km s⁻¹, whereas SB_e is defined as (1 mag arcsec⁻²)$(\frac{m}{1\mathrm{mag}}+2.5log\left(\frac{2\pi R_{\mathrm{e}}^2}{1{\mathrm{kpc}}^2}\right))$, where m is the total magnitude measured using GALFIT.

As anticipated, MUSE observations and their subsequent analysis provide us with line-of-sight central stellar velocity dispersions, σ₀, for 80 out of the 222 member galaxies included in the model. As the FP is actually only valid for early-type galaxies, we need to ensure that the values of its parameters are not influenced by galaxies to which the FP relation does not strictly apply. There are different ways of selecting early-type galaxies, based on their morphologies, colours, and stellar populations. In Tortorelli et al. (2018), it has been shown that early-type galaxies selected via a Sérsic index value n > 2.5, elliptical morphology and a passive stellar population constitute a coeval population. Therefore, we selected early-type galaxies as those galaxies with Sérsic index n > 2.5. Furthermore, we abide by the criteria introduced in B19, based on the robustness of the velocity dispersion measurements by pPXF on simulated spectra, thus including only those with a spectrum signal-to-noise ratio greater than 10 and σ₀ > 80 km s⁻¹. We also excluded two further galaxies: one is in the diffuse light of the BCG, and the second is very close to another member, which influences the velocity dispersion measurement. This leaves us with a total of 30 spectroscopically confirmed cluster members, out of the 80 for which we had a measurement of the stellar velocity dispersion, which we used as a basis to estimate the parameters of the FP relation.

To obtain the FP parameter values, we performed an optimisation using the code ltsfit, developed by Cappellari et al. (2013) for robust fits of lines or planes. A linear regression algorithm is used on the data to derive a first value for the parameters. A σ dispersion around the plane is defined as the range in which 68% of the values fall. The points that lie at more than 3σ from the plane determined by the first regression are clipped and the code repeats the procedure iteratively with the remaining points, until no new member is discarded and the final parameter values are obtained. As for the errors on the observables, we used the propagated experimental uncertainty. Tortorelli et al. (2018) find, for the same cluster, an average uncertainty of 0.09 kpc on R_e and of 0.19 mag arcsec⁻² on SB_e for galaxies with m_F814W ≤ 22.5, which are significantly higher than what we obtain through the error propagation for several members. We chose 0.09 kpc and 0.19 mag arcsec⁻² as lower limits for those uncertainties. We list the best-fit values of α, β, and γ obtained with ltsfit in Table 1. The values of α and β, which do not depend on redshift, are compatible with those reported in the literature (e.g. Barr et al. 2006).

As is clear from Fig. 2, there are very few outliers around the FP: their photometric and spectroscopic analyses have been checked and confirmed. They do not however influence the values of the FP parameters, as they are clipped during the optimisation procedure. To evaluate the accuracy of the calibrated FP at predicting cluster member central velocity dispersions, we can compare the measured values of σ₀ with those obtained through the relation for the 30 members used to derive the parameters. We call these values σ_0, FP. The median value of the ratio $\frac{{\sigma }_{0,\mathrm{FP}}}{{\sigma }_0}$ is 0.98, with a standard deviation of 0.20, corresponding to a ${\chi }_{\mathrm{red}}^2=\sum _{i=1}^{30}\frac{{({\sigma }_{0,\mathrm{FP}i}-{\sigma }_{0i})}^2}{30{\sigma }_{0i}^2}=0.05$, denoting a low scatter around the plane, as shown by Fig. 2. Furthermore, we do not notice any systematic offsets between the observed values of σ₀ and those predicted by the FP.

Fig. 2. Distribution around the FP of the 30 members selected to determine its parameters, with the 1σ uncertainties obtained through the error propagation. For the physical quantities, we use the units listed in the main text.

3.2. The best-fit model

From the best-fit FP parameters, we can compute the value of the central stellar velocity dispersion for all members. As anticipated, we take the value of the truncation radius of the members as a multiple of their measured effective radius, testing several proportionality constant values and re-optimising the values of all the parameters describing the diffuse DM components each time.

To determine the values of the parameters of the best-fit model, we used the χ² minimisation implemented in LensTool. The χ² function is defined as

$$\begin{array}{c}\hfill {\chi }^2(\mathit{\theta })={\displaystyle \sum _{j=1}^{N_{\mathrm{fam}}}\sum _{i=1}^{N_{\mathrm{img}}^j}{\left(\frac{||{\mathit{x}}_{\mathrm{obs}i,j}-{\mathit{x}}_{\mathrm{pred}i,j}(\mathit{\theta })||}{{\sigma }_{xi,j}}\right)}^2,}\end{array}$$(3)

where N_fam is the number of multiple image families considered to evaluate the precision of our model and $N_{\text{img}}^j$ is the number of multiple images that compose the jth family. x_obs i, j and σ_x i, j are, respectively, the observed position of the ith image of the jth family and its uncertainty. Finally, x_{pred i, j}(θ) is the position of the same image as predicted by the model defined by the set of parameter values θ. As in Bonamigo et al. (2018) and B19, we assume an error of 0.5″ on the HST multiple images and of 1″ for the multiple images only found in the MUSE data.

We find that the lowest values of the χ² are obtained when the proportionality constant varies between 2.2 and 2.5, as clear from Fig. 3. In particular, we obtain the absolute lowest value of χ² for r_t = 2.3 R_e. In this case, we find that the minimum χ² has a value of 84.92, with 61 degrees of freedom, which corresponds to a rms value of 0.60″.

Fig. 3. Values of the rms difference between the observed and model-predicted positions of the 55 multiple images used in our model as a function of the ratio between the truncation and effective radii of the cluster members.

In Table 2, the best-fit values of the free parameters describing the cluster-scale DM haloes for this model and for that by B19 are listed. In the table, the index 1 refers to the main, central halo, the index 2 to the secondary, spherical one. The positions of the clumps are indicated with respect to the centre of the BCG. The ellipticity of the main halo is defined as $e=\frac{a^2-b^2}{a^2+b^2}$, where a and b are the major and minor semi-axes. The orientation angle is measured counter-clockwise from the positive x axis. For both models, we determine the statistical uncertainties on the parameters using the Monte Carlo Markov chains as implemented by LensTool, based on the BayeSys algorithm (Jullo et al. 2007). To obtain a reliable estimate of the uncertainties on the parameter values from the Monte Carlo sampling, the observational uncertainty on the positions of the multiple images has been increased to get a χ² value close to the number of degrees of freedom.

4. Comparison and discussion

In this section we compare our model with that presented in B19, which is our starting point and reference. Using the same code, LensTool, B19 obtain a minimum value of χ² of 69.90, which corresponds to a value of the rms on the positions of the 55 multiple images of 0.55″, slightly lower than our value of 0.60″. The maximum values of the likelihood for the two models, provided by LensTool, are reported in Table 3. This small increase in the rms value with respect to B19 was expected, as in our model the mass component representing the member galaxies has less freedom. In B19, the truncation radius and the fiducial velocity dispersion of the ith member are obtained through power laws from its luminosity L_i, specifically

$$\begin{array}{c}\hfill \begin{array}{cc}& {\sigma }_{\mathrm{LT},i}={\sigma }_{\mathrm{LT}}^{\mathrm{ref}}{\left(\frac{L_i}{L_0}\right)}^{0.27},\hfill \\ \hfill & r_{\mathrm{t},i}=r_{\mathrm{t}}^{\mathrm{ref}}{\left(\frac{L_i}{L_0}\right)}^{0.66},\hfill \end{array}\end{array}$$(4)

where L₀ is the luminosity of the BCG. The reference values in the relations (${\sigma }_{\mathrm{LT}}^{\mathrm{ref}}$ and $r_{\mathrm{t}}^{\mathrm{ref}}$) are free to vary, increasing the number of free parameters from 49 to 51. This additional freedom corresponds to some extent to a possible mass exchange between the member galaxies and the diffuse cluster-scale DM haloes. In our model, the total mass component of the cluster galaxies is completely fixed by the observations through the FP, and only the two main DM haloes have their mass free to vary.

We first focus on how mass is assigned to the member galaxies. For both models, cluster members are modelled as dPIEs with vanishing core radius and zero ellipticity, so their total mass can be obtained as

$$\begin{array}{c}\hfill M=\frac{\pi {\sigma }_0^2{r}_{\mathrm{t}}}{G}·\end{array}$$(5)

In B19, r_t and σ₀ are derived, for each member, with the two power-law scaling relations reported in Eq. (4), which imply that $M\propto {\sigma }_0^{4.44}$. Such a precise and simple relation between M and σ₀, with no scatter, might be an oversimplification given the large range in total mass value and the great variety in galactic morphology. The use of two measured observables and a more accurate relation (the FP) should allow for a more complex and realistic dependence of M from the other observables. This is shown in Fig. 4: the M(σ₀) relation has a visible scatter, which the model by B19 could not include. The relation also shows a shallower slope (a bi-logarithmic fit on the members with σ_0, FP > 80 km s⁻¹ gives a slope value of around 2.4). The two models are similar in the high-mass regime and agree on the prediction of the total mass of the BCG. Despite being fixed by the FP with the same procedure as for all other cluster members, our model predicts a BCG total mass value of 3.16 × 10¹² M_⊙, very close to that found in B19 (3.35 × 10¹² M_⊙), where it is de facto a free parameter. Furthermore, the new procedure also leads to a different degeneracy between the parameters that describe the total mass of the sub-haloes. While, as anticipated in Sect. 3, the introduction of a kinematic prior in B19 already allows for a significant reduction in the degeneracy between the values of r_t and σ₀ for the various cluster members, the two reference values that appear in Eq. (4) are still clearly degenerate (as shown in Fig. 4 from B19). In this work, instead, the values of σ₀ are fixed a priori from observations, while the values of r_t are obtained at a later stage, testing several values for the proportionality constant with respect to the observed values of the effective radius.

Fig. 4. Total mass of the cluster members as a function of their velocity dispersion as predicted by our model (in blue) compared to values predicted by the model presented in B19 (in red).

As explained in Sect. 2, we measured the stellar mass values of all 222 member galaxies (Mercurio et al. 2021) from HST multi-band photometric data, using a Salpeter IMF. By combining this with the results of Eq. (5), we can estimate the stellar fraction of the cluster members, computed as the ratio between their stellar mass and total mass values. As lensing models are more sensitive to the effects of high-mass members, we focused only on the 63 galaxies with σ_0, FP > 80 km s⁻¹. For this subset of members, the median value of the stellar fraction is 0.47, with a standard deviation of 0.34. From the model presented in B19, using the same stellar mass values, we obtain a value of 0.25. Therefore, our model predicts, on average, slightly higher values for the stellar fraction, as illustrated in Fig. 5.

Fig. 5. Stellar fraction of the cluster members as a function of their central velocity dispersion for the 63 galaxies with σ0, FP > 80 km s−1: comparison between the values obtained with our model (in blue) and those from the model in B19 (in red).

We now compare our stellar fraction values with those of other massive early-type galaxies. Grillo (2010) studied the projected DM fraction within the effective radius for a sample of nearly 1.7 × 10⁵ early-type galaxies from the Sloan Digital Sky Survey (SDSS) Data Release Seven (DR7), selected to reproduce well the physical properties of the lenses in the SLACS survey (Koopmans et al. 2009). He considered the values of the galaxy stellar mass as reported in the Max Planck for Astrophysics and John Hopkins University catalogue, and obtained the values of the galaxy projected total mass within the effective radius, by assuming a one-component isothermal model, parametrised by the value of the measured central stellar velocity dispersion. He found that the best-fit linear relation between the logarithm of the projected total and stellar mass values within the effective radius, M(< R_e) and M_*(< R_e), respectively, is

$$\begin{array}{c}\hfill log(M(<R_{\mathrm{e}}))=-0.58+1.09log(M_{\ast }(<R_{\mathrm{e}})).\end{array}$$(6)

As far as the cluster members of AS1063 are concerned, assuming that light traces stellar mass, M_*(< R_e) is half of their total stellar mass, as R_e is defined as the two-dimensional radius within which half of the total light is emitted. The stellar mass values of our cluster members have been obtained assuming a Salpeter stellar IMF, while those presented by Grillo (2010) were derived adopting a Chabrier stellar IMF (Chabrier 2003). To compare the two samples, we use a constant conversion factor (0.585, from Speagle et al. 2014) between the stellar mass values obtained with a Salpeter and a Chabrier stellar IMF. As for the projected total mass, the mass of a dPIE profile with a vanishing core radius enclosed within a two-dimensional radius R is (Elíasdóttir et al. 2007)

$$\begin{array}{c}\hfill M(<R)=\frac{\pi {\sigma }_0^2}{G}(R-\sqrt{r_{\mathrm{t}}^2+R^2}+r_{\mathrm{t}}).\end{array}$$(7)

For r_t = 2.3 R_e, Eqs. (5) and (7) give M(R_e) ≈ 0.344 M, where M is the total mass. As Fig. 6 shows, we find that our members follow well the relation found by Grillo (2010). This is significant, as we are considering in both cases massive elliptical galaxies, but the values of their stellar fractions are obtained with different methods. The observed agreement suggests that the inner (i.e. within R_e) mass structure of the massive cluster members of AS1063 and of massive early-type galaxies living in different environments might be very similar. It should be noted that the relation shown in Fig. 6 has been adapted from Eq. (7) for the case of stellar mass values determined using a Salpeter IMF.

Fig. 6. Logarithmic values of the total mass as a function of the logarithmic values of the stellar mass, both projected within the effective radius, of all members with σ0, FP > 80 km s−1. The plot also shows the relation found by Grillo (2010).

Looking instead at the diffuse cluster-scale components of the two models, from the Monte Carlo sampling of the marginalised posterior probability we see a drop of more than 30% of the uncertainty on the value of the core radius of the main DM halo (see Fig. 7). The scatter between the core radius and the velocity dispersion is also considerably reduced. This indicates that a more accurate estimate of the total mass of the member galaxies has resulted in distributions of the parameters with a lower statistical uncertainty.

Fig. 7. Posterior probability distributions of the core radius, rc, and fiducial velocity dispersion, σLT, values of the main DM halo and degeneracy between them for our model (in blue) and that from B19 (in red). The 16th, 50th, and 84th percentiles of the marginalised probability distributions are marked with dashed vertical lines for both models. The values of these percentiles for our model are reported above the histograms.

Monte Carlo Markov chains also show that the core of the main diffuse DM halo has a radius of 86 ± 2 kpc. Our result is in line with those from SL models of similar galaxy clusters (e.g., B19, Grillo et al. 2015; Annunziatella et al. 2017; Lagattuta et al. 2017; Caminha et al. 2019; Diego et al. 2020), which found that cluster haloes are very well described by cored isothermal mass density profiles. Figure 7 also shows a clear degeneracy, with a positive correlation, between the velocity dispersion and the core radius of the main DM halo. This is expected, as different combinations of these two parameters can determine the same value of the projected mass enclosed within a certain radius from the cluster centre. In particular, a higher value of central density scale, here included in the definition of σ_LT, can be compensated by a larger core in the density profile, and vice versa, to obtain the same mass within the multiple images of the same family.

Finally, from Table 2, we notice that the values of our and B19 best-fit parameters of the diffuse DM component are very similar. The only significant discrepancy is the value of the core radius of the main DM halo of the cluster. This can be explained by considering that, far from the cluster centre, the total mass distribution of the cluster is dominated by the cluster-scale DM halo component. The positions of the multiple images observed in these regions provide strong constraints on the values of all the parameters of this halo. The core radius of the cluster DM halo is more sensitive to the total mass profile of the cluster close to its centre, which is more directly related to the mass distribution of the cluster members located in this region. A change in the mass contribution of the cluster members is therefore reflected in a small change of the DM halo core radius.

5. Mass profile decomposition

In this section we obtain and present the cumulative projected mass profiles of all the cluster mass components: the stellar mass and the DM of cluster members, the ICM, and the diffuse, cluster-scale DM. We only consider the mass distribution out to a projected radius R of 350 kpc from the centre of the BCG, to avoid including regions not covered by our photometric data, where the cluster member catalogue is incomplete.

5.1. The baryonic mass components

To measure the stellar mass profile of the cluster, we started from the surface brightness best-fit model of each of the cluster members in the HST F814W band, shown in Fig. 8, which we already used to derive their structural parameters. Figure 8 shows that our procedure does not allow the intra-cluster light distribution, which also contributes to the cluster stellar mass budget, to be completely accounted for. However, the intra-cluster light in AS1063, especially far from its BCG, is less important than in other clusters of comparable mass. Similarly to Annunziatella et al. (2017), we estimated an average value for the stellar mass-to-light ratio (M_*/L) of all cluster members within a projected radius of 350 kpc from their total luminosity and total stellar mass. The median value of the stellar mass-to-light ratio for these members is 1.0 M_⊙/L_⊙, estimated in the HST F814W band. The choice of considering a constant ratio is an approximation, as the value may change in the different regions of each galaxy, as well as among galaxies. However, several galaxies are included in each radial bin above 50 kpc, and this should reduce the impact of this approximation. We then used this ratio to convert the cumulative luminosity profile into a cumulative stellar mass profile. As above, we considered an error of 0.22 on log(M_*) for each galaxy, and we applied the same uncertainty to the total value in each radial bin of the projected profile.

Fig. 8. Reconstructed surface brightness distribution of all cluster members used in the SL model in the HST F814W band. The field of view is 700 kpc wide and centred on the BCG to cover the entire region where the cumulative mass profile is obtained.

We take the cumulative projected ICM mass profile from Bonamigo et al. (2018): following their result, we impose a relative error of 1% on each bin. The total projected mass of the three haloes used to model the hot gas, within a bi-dimensional radius of 350 kpc, is reported in Table 4.

5.2. The dark matter mass distribution

In Sect. 3, we described our modelisation of the DM mass distribution: an ellipsoidal, pseudo-isothermal, non-truncated clump, and a less massive, spherical one represent the cluster-scale DM haloes. On the other hand, 222 spherical, truncated, isothermal haloes were introduced to describe the total mass of cluster members. To disentangle the mass distribution of DM sub-haloes, we therefore need to subtract the stellar mass profile to it.

The parameter space sampling obtained with the Monte Carlo method can be used to estimate the error on the DM cumulative mass profiles. We extracted 100 random sets of parameter values from the chains and derived the cumulative mass profiles of each corresponding model. In Fig. 9 and in Table 5, the error bars show the 16th and 84th percentiles of the mass profiles associated with all the models extracted from the chains. We performed this analysis for our model and repeated it for the model by B19. In the first case, the total mass of the cluster members was fixed, so the uncertainty on the DM sub-halo mass profile is only associated with the error on the stellar mass. In the second case, instead, the mass of the cluster members was also optimised: Monte Carlo Markov chains provided us with an uncertainty on their mass, which was combined with that on M_* to obtain the error on the DM sub-halo cumulative mass profile.

Fig. 9. Cumulative mass and mass fraction profiles of the cluster components out to a projected radius of 350 kpc. Left panel: cumulative two-dimensional mass profiles for all baryonic and dark components of AS1063. The total mass profile and those of the diffuse DM haloes and DM sub-haloes obtained from the model in B19, which differ from ours, are represented with dashed lines. Right panel: cumulative stellar-, gas-, and baryonic-to-total mass fractions from our model. The value of the cosmological baryon fraction from the Planck Collaboration VI (2020) is also shown.

We notice that the total mass of the cluster and the mass of the diffuse DM component are very well constrained by SL: the two models find very similar values and their statistical uncertainties are around 1%. As a matter of fact, the two profiles are almost completely superimposed in Fig. 9, and the error bars are very small. The mass of the DM sub-haloes, instead, slightly differs between the two models, especially in the outer regions, where our model predicts a lower value. A reason for this discrepancy might be in the different scaling laws that determine r_t in the two models. Comparing our results with the total mass profiles obtained, for the same cluster, with a full dynamical reconstruction and with the X-ray hydrostatic method (Sartoris et al. 2020), and with combined weak and SL (Umetsu et al. 2016), we find very good agreement in the radial range considered. These two profiles are shown in Fig. 7 by Sartoris et al. (2020).

The profiles measured from our model also allow us to obtain the two-dimensional, cumulative, stellar-, gas-, and baryonic-to-total mass fractions as a function of R, out to 350 kpc. In Fig. 9, we report the results and compare the baryonic fraction profile with the cosmological baryon fraction value reported in the Planck 2018 results: 0.156 ± 0.002, obtained from the ratio between Ω_bh² and Ω_mh², with information from the cosmic microwave background power spectrum and lensing reconstruction (see Planck Collaboration VI 2020). The error bars of the profiles result from the propagation of the uncertainty on the cumulative values of the mass components. The fraction obtained at R = 350 kpc is 0.147 ± 0.002. This value is similar to the cosmological value and consistent with what Bonamigo et al. (2018) reported for the same cluster at a comparable radius.

6. Comparison with halo occupation distribution studies

In this section we compare the values of the total mass predicted for cluster members by our model with those obtained from the measured stellar mass values using halo occupation distribution (HOD) studies, based on cosmological N-body simulations. In the previous section, we found that both our model and that from B19 lead to median values of the stellar-to-total mass ratio of massive member galaxies of 0.25 at least. As we will see, this result differs significantly from the predictions of statistical (sub)-halo abundance matching studies. We compare our results with the recent stellar-to-halo mass relation (SHMR) presented in Girelli et al. (2020). They considered the halo mass function derived from the cosmological N-body collisionless ΛCDM DUSTGRAIN-pathfinder simulation (Giocoli et al. 2018), which is complete for halo masses M_h > 10^12.5 M_⊙, and extended it to lower masses. The simulation is built from a sample that also includes some cluster members, especially in the high-mass range. Then, taking the observed galaxy stellar mass of the COSMOS field (Scoville et al. 2007), they performed an abundance matching via a deterministic relation between halo and stellar mass. The SHMR is found by Girelli et al. (2020) to be in good agreement both with the observed clustering of massive galaxies and with the SHMR of clusters selected from spectroscopic data from the SDSS. The inclusion of some scatter on the value of the stellar mass leads to consistent results for the relation, so we do not consider it in the following analysis.

The obtained SHMR is valid from stellar mass values between approximately 10⁸ M_⊙ and 10¹² M_⊙, from z = 0 to z = 4. The relation, using the parametrisation proposed by Moster et al. (2010), can be written as

$$\begin{array}{c}\hfill \frac{M_{\ast }}{M_{\mathrm{h}}}(z)=2A(z){[{\left(\frac{M_{\mathrm{h}}}{M_{\mathrm{A}}(z)}\right)}^{-\beta (z)}+{\left(\frac{M_{\mathrm{h}}}{M_{\mathrm{A}}(z)}\right)}^{\gamma (z)}]}^{-1},\end{array}$$(8)

where A(z) is a normalisation factor, M_A(z) a characteristic halo mass, and β(z) and γ(z) are the slopes of the relation in the low- and high-mass ranges, respectively. We adopted the values of the parameters found by Girelli et al. (2020) at the cluster redshift. We then proceeded by using Eq. (8) to compute the total halo mass of the members and, in turn, the values of their truncation radii from Eq. (5). We finally compared these values of the truncation radii with those we obtained in our model from the effective radii of the members. We did not consider any uncertainty on the parameter and stellar mass values, as we will show that the differences are large and cannot be reconciled by the statistical errors. We notice that the SHMR presented by Girelli et al. (2020) was derived adopting a Chabrier stellar IMF. To make use of the relation to estimate the halo mass of the cluster members from their stellar mass, we adopted the same conversion factor of 0.585 between the Salpeter and Chabrier IMFs that we previously used.

Using the SHMR, the projected mass of the DM sub-haloes within a bidimensional radius of 350 kpc grows to $6.1_{-0.1}^{+0.1}\times {10}^{13}{M}_{\odot }$, larger by a factor of almost 14. The ratio $r=\frac{M_{\mathrm{SHMR}}}{M_{\mathrm{model}}}$, where M_SHMR is the total mass obtained for a member from SHMR, and M_model is the corresponding total mass as in the SL model, has a median value of 8.6 for the members with σ_0, FP > 80 km s⁻¹. For these selected high-mass members, as shown in Fig. 10, the median value of the stellar fraction decreases from 0.47 to 0.055. The value of r is also that of the ratio between the value of the truncation radius of a sub-halo obtained from simulations and that adopted in our model.

Fig. 10. Comparison between the values of the stellar over total mass fraction of the cluster members obtained with our lensing model (in blue) and with the SHMR from Girelli et al. (2020) (in orange). Only the 63 cluster members with σ0, FP > 80 km s−1 are shown. Left panel: stellar over total mass fraction of members as a function of their stellar mass. Right panel: SHMR.

We performed a new run of the model, fixing the total mass of each cluster member to the value obtained from the SHMR. Optimising the free parameters of the two diffuse DM haloes, the value of the rms between the observed and predicted positions of the 55 multiple images grows by more than 25% (to 0.76″), despite having the same number of degrees of freedom.

These results seem to indicate a significant discrepancy between SL models and the predictions of an SHMR derived from N-body cosmological simulations. The total halo mass values predicted by the SHMR are significantly higher than what is predicted by both our SL model and that from B19, and their inclusion in an SL model leads to a significant decrease in its accuracy. An important reason for this discrepancy may be the fact that the ΛCDM DUSTGRAIN-pathfinder simulation only includes collisionless particles. Thus, it does not consider the effects of the interplay between DM and baryons, which are very important during galaxy formation. Several complex phenomena taking place on small spatial scales, such as star formation, radiative gas cooling, and active galactic nuclei feedback can have a non-negligible effect on the physical properties of galactic haloes (see Kravtsov & Borgani 2012; Planelles et al. 2014; Barnes et al. 2017). Moreover, the SHMR presented by Girelli et al. (2020) describes both cluster and field galaxies with the same law. During the formation of a cluster, galaxies fall into the cluster potential well. At first, they undergo a phase of star formation triggered by shocks in their gas component. At a later stage, both the baryonic and DM components of the members are stripped by ram pressure and gravitational tidal interactions (see, for instance, McCarthy et al. 2008; van den Bosch et al. 2008). This phase could determine an SHMR different from that of field galaxies: Smith et al. (2016), for instance, found that the DM stripping happens earlier than star stripping, leading to an increase in the ratio between the effective and virial radii of the sub-haloes. This could justify a higher value of the stellar mass fraction of cluster members when compared to field galaxies. While the most massive galaxies considered by Girelli et al. (2020) to calibrate the SHMR are mostly hosted in clusters and groups, the low-mass galaxies considered are mainly field galaxies, which have evolved in an environment different from that of the cluster members we focus on here. To address both the role of baryons and DM interactions and the influence of the dense environment of clusters in the evolution of galaxies, Sirks et al. (2021) used the Hydrangea/C-EAGLE suite (Bahé et al. 2017; Barnes et al. 2017) cosmological hydrodynamical simulations to build two distinct SHMRs, introducing a distinction between cluster and field galaxies. As clear from Fig. 5 in Sirks et al. (2021), the SHMR of cluster members yields stellar-to-total mass fraction values that are almost an order of magnitude higher than those of field galaxies, at a fixed stellar mass. As we see in the next section, such an effect will alleviate the reported discrepancies between the stellar fraction values predicted by SHMRs and reconstructed from our SL model.

7. Comparison with cosmological hydrodynamical simulations

In this section we contrast the features of the sub-haloes in our new SL model with the predictions of the most recent cosmological hydrodynamical simulations. In the first subsection we consider the stellar over total mass fraction of sub-haloes. In the second subsection we study the relation between their maximum circular velocity and total mass, which is an indication of their compactness (Munari et al. 2016).

7.1. Stellar fraction of cluster members

As reported in Sect. 4, our model predicts a median stellar-to-total mass fraction value of 0.47, with a standard deviation of 0.34, for the 63 cluster members with σ_0, FP > 80 km s⁻¹. We have seen, in Sect. 6, that an SHMR obtained with an HOD procedure on the results of N-body cosmological simulations underestimates this value with respect to our model by almost an order of magnitude. We then repeated the comparison with high-resolution cosmological hydrodynamical simulations (Planelles et al. 2014; Rasia et al. 2015) carried out with the GADGET-3 code (Springel 2005). To perform a comparison between the features of sub-haloes as predicted by the simulations and by our model, we considered 18 simulated clusters with M > 10¹⁵ M_⊙ at redshift z = 0.38 and performed three spatial projections of their mass distributions. We then extracted their sub-haloes from the SUBFIND algorithm (Springel et al. 2001; Dolag et al. 2009) catalogue¹. We imposed a minimum value of the velocity dispersion of 80 km s⁻¹, as in the previous analyses carried out in this paper, and a maximum projected distance from the cluster centre of R = 0.15 r_vir, where r_vir is the cluster virial radius. We thus obtained a catalogue of 1756 sub-haloes, and we computed the stellar mass fraction for each of them. Again, as the simulations considered use a Chabrier IMF to obtain the stellar mass of sub-haloes, we used the previously presented conversion factor to compare them with our measured stellar mass values, for which Mercurio et al. (2021) used a Salpeter IMF. We obtain a median stellar mass fraction value of 0.37, with a standard deviation of 0.18, in good agreement with the value suggested by the lensing model. Figure 11 also shows that the simulated sample has a similar relation between the stellar mass fraction and the total sub-halo mass with respect to our model. We note that the sub-haloes extracted from the simulated catalogue have a slightly higher (around 3 × 10¹⁰ M_⊙) lower limit for the total mass compared to those extracted from the SL model adopting the same criterion. This seems to suggest that at a fixed velocity dispersion simulated sub-haloes have, on average, a higher total mass. In the next subsection we address this issue in light of the recent studies that have focused on it.

Fig. 11. Comparison between the values of the stellar over total mass fraction of sub-haloes with velocity dispersion greater than 80 km s−1 as a function of their total mass. The values reconstructed from our model are shown in blue and the values predicted via hydrodynamical simulations by Planelles et al. (2014) in orange.

7.2. Compactness of the cluster sub-haloes

Significant differences between the sub-halo compactness obtained from lensing models and cosmological simulations have recently been reported by M20. They considered 11 SL models of massive galaxy clusters, including the model of AS1063 by B19, and a set of 25 clusters of similar mass and redshift from the cosmological simulations presented in Planelles et al. (2014) and Rasia et al. (2015), which incorporate gas cooling, star formation, and energy feedback from supernovae and accreting supermassive black holes. They found that the probability of producing galaxy-scale SL events is for observed clusters around an order of magnitude larger than it is for simulated clusters. To interpret this discrepancy, they considered the relation between the mass of the sub-haloes in the model and their maximum circular velocity, defined as ${\mathit{v}}_{\max }=\max \sqrt{\frac{\mathit{GM}}{r}}$ (for an isothermal model, ${\mathit{v}}_{\max }=\sqrt{2}{\sigma }_0$). M20 found that, at fixed total sub-halo mass, the value of the maximum circular velocity of cluster sub-haloes is larger than what cosmological simulations predict. This seems to indicate that simulated cluster members are, on average, less compact than observed.

Some other recent results, using the higher-resolution simulations from the Hydrangea/C-EAGLE suite (Bahé et al. 2017; Barnes et al. 2017), do not claim the same large discrepancy in terms of galaxy scale lensing events and of sub-halo compactness (Robertson 2021; Bahé 2021, hereafter B21). Studying simulated clusters, B21 reported a bi-modal v_max(M) function, with a lower sequence of sub-haloes with a baryon fraction smaller than 0.1 dominating for M < 10¹¹ M_⊙, and a second sequence for higher sub-halo mass values. This last branch yields results consistent with those of the lensing models presented in M20.

All SL models included in M20 used two optimised power-law scaling relations to model the truncation radius and the velocity dispersion of the cluster members from their luminosity (B19, Caminha et al. 2019). As all cluster models represent sub-haloes with isothermal profiles, this also leads to a fixed power-law relation between the maximum circular velocity (which is proportional to the velocity dispersion) of each sub-halo and its mass. On the other hand, our new method to assign mass to sub-haloes again avoids a fixed power law between the two quantities, allowing for the inclusion of some scatter. In Fig. 12, we show the v_max(M) relation obtained from our model for all the members with M > 10¹⁰ M_⊙ compared with that from B19: the two models agree in the entire mass range considered.

Fig. 12. Maximum circular velocity of the cluster members as a function of their total mass. The values predicted by the model presented in B19 are shown in red. Our values are colour-coded according to the predicted stellar over total mass fraction. Only members with M > 1010 M⊙ are shown.

Comparing our results with those shown in Fig. 3 from B21, we notice that our relation matches fairly well with the higher-mass branch of the Hydrangea hydrodynamic simulations, both in terms of circular velocity and of the stellar fraction of the cluster members at a fixed total mass. However, we note that the high-mass branch that agrees with our model in terms of maximum circular velocity is not found, for instance, in the Illustris-TNG300s simulation (Marinacci et al. 2018; Naiman et al. 2018; Nelson et al. 2018; Pillepich et al. 2018; Springel et al. 2018), which have a comparable resolution, as shown in Fig. 4 from B21. Furthermore, Ragone-Figueroa et al. (2018) also found that the Hydrangea simulations produce galaxies with a very high total mass that may not agree with observations. On the other hand, as already visible from Fig. 12, we do not find the large number of sub-haloes with M < 10¹¹ M_⊙, v_max < 100 km s⁻¹, and a low stellar fraction predicted by the lower sequence by B21. As a consequence, our model predicts, on average, higher maximum circular velocities with respect to simulations for sub-haloes with mass between 10¹⁰ and 10¹¹ M_⊙. In this total mass range, we find just a few members with maximum circular velocity lower than 100 km s⁻¹, while in B21 most of the galaxies lie below this threshold. This seems to confirm the discrepancy reported by M20 in terms of galaxy-scale SL events, as most of these events are caused by haloes in this mass range, at least for AS1063 and the other clusters with a similar mass. We finally note that in B21 the values of the baryon fractions are computed within a radius of 30 kpc, while in our work they are the ratio between the total stellar mass and the total mass of each sub-halo. To compare the values obtained in the two cases, one must keep in mind that our model predicts that, for the high-mass cluster members, a median value of 94% of the total mass is enclosed within a sphere with a radius of 30 kpc.

The comparative lack of low-mass members in our model is reflected by the cumulative sub-halo mass function within a projected radius of 0.15 r_vir, shown in Fig. 13: our model agrees with the results presented in B19, which in turn found a lower cumulative mass function compared to that of the simulations, particularly considering sub-haloes with M < 10¹¹ M_⊙ and v_max < 100 km s⁻¹, as reported by B21 in Fig. 1. In the same mass range, our cumulative sub-halo mass function has a slope similar to that by B19, significantly shallower than that from the cosmological simulations presented in B21. In the radial range considered, our sample of photometric cluster members is highly complete down to m_F160W = 24.

Fig. 13. Cumulative sub-halo mass function for the members within a projected radius of 0.15 rvir predicted by our model (in blue) by that presented in B19 (in red).

8. Conclusions

In this article we have presented an improved SL model of the HFF galaxy cluster AS1063 thanks to new measurements of the structural parameters, namely magnitude, effective radius, and Sérsic index, of the cluster members. We have also exploited the values of the central stellar velocity dispersion of several luminous members, obtained from MUSE integral-field spectroscopy, to calibrate the FP relation for the cluster.

We used the FP to estimate the value of the central velocity dispersion of each cluster member from its observed magnitude and effective radius, and we chose to consider a proportionality relation between the values of the effective radii, R_e, and truncation radii, r_t. This allowed us to assign a value of total mass to every cluster member that does not depend only on the total luminosity of a galaxy. The main results and conclusions of our analysis are summarised as follows:

1. We find that the lowest value of the rms difference between the observed and model-predicted positions of the 55 considered multiple images (from 20 background sources) is obtained when r_t = 2.3 R_e for all cluster members.

2. Our method allows us to determine a relation between the total mass and velocity dispersion of the cluster members that is more realistic than the simple power law used in B19. Specifically, the new relation has some visible scatter and is less steep than that recently found in B19.

3. A more accurate estimate of the total mass of the cluster members allows for a reduction in the scatter on the determination of some parameters of the diffuse component. In particular, the statistical uncertainty on the value of the core radius of the main DM halo decreases by more than 30% compared to previous studies.

4. We confirm that modelling the main DM haloes of massive galaxy clusters with isothermal mass density profiles requires the presence of significantly large core radii. The value of the core radius of the main DM halo of AS1063 is 86 ± 2 kpc.

5. From our model and from the sampling of the posterior probability distribution of the parameters, we obtain the cumulative projected mass profiles of the various cluster components. Thanks to Chandra X-ray observations and new measurements of the stellar mass values of the cluster members, we can disentangle the hot gas and stellar mass profiles. At a distance from the cluster centre of 350 kpc, the cumulative projected stellar, hot-gas, and baryonic fractions are, respectively, 0.6%, 14.1%, and 14.7%; the last value is close to the cosmological value found by Planck Collaboration VI (2020).

6. For the high-mass members included in our model, the median value of the stellar over total mass fraction is 0.47, with a standard deviation of 0.34. Considering, instead, the ratio between the stellar and total mass projected within the effective radius, we find good agreement with the relation obtained by Grillo (2010) on a wide sample of SDSS early-type galaxies.

7. We used the SHMR from Girelli et al. (2020), based on the ΛCDM DUSTGRAIN-pathfinderN-body simulation, to obtain the halo mass values of the cluster members from their measured stellar mass values. We find that the SHMR predicts a median value of the stellar fraction of the high-mass cluster members that is smaller by almost an order of magnitude with respect to our model.

8. We compared the high-mass cluster members with a catalogue of simulated sub-haloes. We extracted the sub-haloes from the projected cores of 18 clusters from recent cosmological hydrodynamical simulations (Planelles et al. 2014; Rasia et al. 2015), imposing the same selection criterion on their velocity dispersion. We report a good agreement between the values of the stellar over total mass fraction of the two samples.

9. We studied the relation between the maximum circular velocity and total mass of the sub-haloes, which is an indication of their compactness. Comparing our results with those obtained by B21 from the simulated clusters of the Hydrangea/C-EAGLE suite, we obtain consistent results for sub-haloes with total mass greater than 10¹¹ M_⊙. On the contrary, we do not observe the large number of sub-haloes with mass below 10¹¹ M_⊙ and low baryon fraction predicted by B21. In this sub-halo mass range, our relation predicts, on average, higher v_max at a fixed total mass, in agreement with B19 and M20. This corresponds to haloes that are more compact than those in the hydrodynamical simulations. As a consequence, in the same total mass range our cumulative sub-halo mass function has a shallower slope compared to what was reported by B21 and very similar to that from B19 and M20.

Similar analyses can be performed on other massive clusters, for which MUSE data are available, to test the robustness of the results obtained with the new methodology presented here.

Acknowledgments

We thank the anonymous referee for some useful suggestions that helped to improve the paper. We acknowledge financial support by PRIN-MIUR 2017WSCC32 “Zooming into dark matter and proto-galaxies with massive lens-ing clusters” (P.I.: P. Rosati), INAF “main-stream” 1.05.01.86.20 (P.I.: M. Nonino) and INAF 1.05.01.86.31 (P.I.: E. Vanzella). G.G. thanks Ana Acebron Muñoz and Andrea Bolamperti for their help and suggestions. P.B. acknowledges financial support from ASI through the agreement ASI-INAF n. 2018-29-HH.0. M.M. acknowledges support from ASI through the grants ASI-INAF n. 2018-23-HH.0 and ASI-INAF n. 2017-14-H.0. G.B.C. acknowledges the Max Planck Society for financial support through the Max Planck Research Group for S.H. Suyu and the academic support from the German Centre for Cosmological Lensing.

References

All Tables

All Figures

	Fig. 1. Colour composite image of the core of AS1063 created from HFF photometry. Blue, green, and red channels are obtained by combining the filters F435W, F606W + F814W, and F105W + F125W + F140W + F160W, respectively.
In the text

	Fig. 2. Distribution around the FP of the 30 members selected to determine its parameters, with the 1σ uncertainties obtained through the error propagation. For the physical quantities, we use the units listed in the main text.
In the text

	Fig. 3. Values of the rms difference between the observed and model-predicted positions of the 55 multiple images used in our model as a function of the ratio between the truncation and effective radii of the cluster members.
In the text

	Fig. 4. Total mass of the cluster members as a function of their velocity dispersion as predicted by our model (in blue) compared to values predicted by the model presented in B19 (in red).
In the text

	Fig. 5. Stellar fraction of the cluster members as a function of their central velocity dispersion for the 63 galaxies with σ0, FP > 80 km s−1: comparison between the values obtained with our model (in blue) and those from the model in B19 (in red).
In the text

	Fig. 6. Logarithmic values of the total mass as a function of the logarithmic values of the stellar mass, both projected within the effective radius, of all members with σ0, FP > 80 km s−1. The plot also shows the relation found by Grillo (2010).
In the text

Fig. 7. Posterior probability distributions of the core radius, rc, and fiducial velocity dispersion, σLT, values of the main DM halo and degeneracy between them for our model (in blue) and that from B19 (in red). The 16th, 50th, and 84th percentiles of the marginalised probability distributions are marked with dashed vertical lines for both models. The values of these percentiles for our model are reported above the histograms.

In the text

	Fig. 8. Reconstructed surface brightness distribution of all cluster members used in the SL model in the HST F814W band. The field of view is 700 kpc wide and centred on the BCG to cover the entire region where the cumulative mass profile is obtained.
In the text

Fig. 9. Cumulative mass and mass fraction profiles of the cluster components out to a projected radius of 350 kpc. Left panel: cumulative two-dimensional mass profiles for all baryonic and dark components of AS1063. The total mass profile and those of the diffuse DM haloes and DM sub-haloes obtained from the model in B19, which differ from ours, are represented with dashed lines. Right panel: cumulative stellar-, gas-, and baryonic-to-total mass fractions from our model. The value of the cosmological baryon fraction from the Planck Collaboration VI (2020) is also shown.

In the text

	Fig. 10. Comparison between the values of the stellar over total mass fraction of the cluster members obtained with our lensing model (in blue) and with the SHMR from Girelli et al. (2020) (in orange). Only the 63 cluster members with σ0, FP > 80 km s−1 are shown. Left panel: stellar over total mass fraction of members as a function of their stellar mass. Right panel: SHMR.
In the text

	Fig. 11. Comparison between the values of the stellar over total mass fraction of sub-haloes with velocity dispersion greater than 80 km s−1 as a function of their total mass. The values reconstructed from our model are shown in blue and the values predicted via hydrodynamical simulations by Planelles et al. (2014) in orange.
In the text

	Fig. 12. Maximum circular velocity of the cluster members as a function of their total mass. The values predicted by the model presented in B19 are shown in red. Our values are colour-coded according to the predicted stellar over total mass fraction. Only members with M > 1010 M⊙ are shown.
In the text

	Fig. 13. Cumulative sub-halo mass function for the members within a projected radius of 0.15 rvir predicted by our model (in blue) by that presented in B19 (in red).
In the text