© The Authors 2026

1. Introduction

Galaxy clusters serve as a crucial link between astrophysics and cosmology, acting as key tracers of structure formation through hierarchical assembly and as laboratories for studying galaxy evolution in dense environments. As smaller haloes are accreted onto larger structures during cluster formation, they undergo a range of physical processes (Kravtsov & Borgani 2012). Tidal stripping and dynamical heating can remove substantial mass, while baryonic processes, such as radiative cooling and feedback from star formation and active galactic nuclei, reshape their internal structure and regulate their star formation activity. Consequently, the dark matter (DM) distribution in clusters retains signatures of both the large-scale assembly history and the micro-physical interactions between DM and baryons (Giocoli et al. 2008, 2010).

Strong lensing has emerged as a prime technique to map the total mass distribution of massive galaxy clusters. When combined with observational probes of the mass of baryons, it lets us investigate the physical properties and mass structure of DM haloes (e.g. Natarajan & Kneib 1997; Caminha et al. 2019; Granata et al. 2022, 2023). However, several studies have shown that the recovery of the mass distribution of galaxy-scale DM haloes (or sub-haloes) is hindered by parametric degeneracies whose impact can be significantly reduced with independent observational priors obtained from the stellar kinematics of the cluster members (Bergamini et al. 2019; Granata et al. 2022). The stellar line-of-sight velocity distribution of early-type galaxies (ETGs) is a direct and relatively unbiased tracer of the total gravitational potential, less affected by systematics or a priori assumptions than luminosity- or stellar-mass-based estimates. The central stellar velocity dispersion (hereafter indicated with σ), in particular, serves as a robust probe of the gravitational potential well in ETGs. Its distribution, or function, for a galaxy population offers a way to link its properties with theoretical models of galaxy formation and evolution.

The velocity dispersion function (VDF) of ETGs has been extensively characterised thanks to large spectroscopic surveys such as the Sloan Digital Sky Survey (SDSS) and the SDSS Baryon Oscillation Spectroscopic Survey (BOSS; e.g. Sheth et al. 2003; Mitchell et al. 2005; Choi et al. 2007; Bernardi et al. 2010; Montero-Dorta et al. 2017; Sohn et al. 2017b), and in the NMBS COSMOS and UKIDSS UDS deep fields (Bezanson et al. 2012). These efforts have revealed the dependence of the velocity dispersion function on galaxy morphology, environment, and redshift. Following Sheth et al. (2003), the VDF of red galaxies is typically modelled by a modified Schechter function:

$$\begin{array}{c}\hfill \varphi (\sigma )\mathrm{d}\sigma ={\varphi }^{\ast }{\left(\frac{\sigma }{{\sigma }^{\ast }}\right)}^{\alpha }\frac{exp[-{\left(\frac{\sigma }{{\sigma }^{\ast }}\right)}^{\beta }]}{\mathrm{\Gamma }(\alpha /\beta )}\beta \frac{\mathrm{d}\sigma }{\sigma }·\end{array}$$(1)

Here ϕ^* is the normalisation, σ^* a characteristic velocity dispersion value, and α, β control the low- and high-σ slopes. This form originates from the connection between galaxy luminosity and velocity dispersion via the Faber-Jackson relation (Faber & Jackson 1976) and reflects the broader success of the Schechter formalism in describing galaxy demographics (Schechter 1976). Variations in the parameters derived between studies arise primarily from the differing selection criteria for ETGs and sample completeness. Generally, these functions exhibit a downturn at log σ[km s⁻¹]< 2.35, governed by positive α values. Notably, Bezanson et al. (2012) find minimal evolution at the high-σ end, suggesting that massive quiescent galaxies formed early, while the increase at lower σ reflects ongoing quenching in intermediate-mass systems.

Measuring the VDF of cluster members, however, requires spatially dense spectroscopic surveys and very deep observations to achieve the necessary spectral signal-to-noise ratio (S/N) for a large, complete, and unbiased samples of member galaxies. Sohn et al. (2017a) calibrated the VDF for the quiescent galaxies in Coma and Abell 2029 using data from the wide-field spectrograph MMT/Hectospec and found clear deviations from the field VDF. In particular, the cluster VDF is flatter at low dispersions, though this difference may reflect incompleteness of the field sample (Sohn et al. 2017a). An excess at log σ[km s⁻¹]> 2.4 is also seen, likely due to the presence of BCGs in cluster cores. Sohn et al. (2020) extended these results to nine strong lensing clusters at 0.18 < z < 0.29, using data from the same instrument, and confirming both the high-σ excess and compatibility with Coma and Abell 2029, down to log σ[km s⁻¹]∼2.05.

Our aim is to extend the study of the kinematic properties of red cluster members to higher-redshift and lower-mass galaxies by taking advantage of deep observations from the integral-field Multi Unit Spectroscopic Explorer (MUSE, Bacon et al. 2012) at the Very Large Telescope (VLT), which has allowed systematic spectroscopic coverage of the cores of strong lensing galaxy clusters, with exposure times occasionally above 10 h. The depth of the observations extends kinematic measurements to fainter cluster members (see Bergamini et al. 2021), and the integral-field nature of the instrument mitigates the impact of the high galactic density in cluster cores. For this purpose, we have developed a pipeline to systematically measure the central stellar velocity dispersion values of the largest possible set of cluster members observed by MUSE. The relatively small MUSE field of view (FoV) of 1′×1′ limits the region explored to the cluster core (i.e. a few hundred kiloparsecs from the cluster centre). However, we can calibrate scaling laws, such as the Fundamental Plane relation (Dressler et al. 1987; Djorgovski & Davis 1987; Bender et al. 1992, hereafter FP), that link the kinematic properties of cluster members to their structural parameters, and then use these laws to study the kinematic properties of all members within the larger FoV of the photometric observations available. We measure the stellar VDF of red cluster members and discuss the results in terms of the observational picture of galaxy formation and evolution in cluster cores.

The paper is organised as follows. In Sect. 2 we give details on the sample of four galaxy clusters studied and the observational dataset on which this work is based. In Sect. 3 we present our velocity dispersion measurement pipeline, and the catalogue of velocity dispersion values we measured. In Sect. 4 we calibrate the FP for the ETGs in the four clusters. In Sect. 5 we use them to extend the kinematic study to all red galaxies for which we have measured structural parameters, and we calibrate the VDF for the four clusters. Finally, in Sect. 6 we summarise our results. Throughout this work, we use a flat ΛCDM cosmology with Ω_m = 0.3 and H₀ = 70 km s⁻¹ Mpc⁻¹.

2. The galaxy cluster sample

We selected a sample of four massive strong lensing galaxy clusters at 0.31 < z < 0.54 with deep, homogeneous photometric and spectroscopic coverage, and a state-of-the-art strong lensing model available: Abell 2744 (z = 0.307, hereafter A2744, Allen 1998; Ebeling et al. 2010), Abell S1063 (z = 0.346, hereafter AS1063, Abell et al. 1989; Karman et al. 2015; Mercurio et al. 2021), MACS J0416.1−2403 (z = 0.397, hereafter M0416, Ebeling et al. 2001; Balestra et al. 2016), and MACS J1149.6+2223 (z = 0.542, hereafter M1149, Smith et al. 2009; Grillo et al. 2016).

Optical and near-infrared imaging for the clusters was guaranteed by the Hubble Space Telescope (HST). All four clusters were part of the Hubble Frontier Fields (HFF, Lotz et al. 2017) programme, and targeted with deep exposures (140 HST orbits) in seven ACS and WFC3 filters (F435W, F606W, F814W, F105W, F125W, F140W, and F160W). The HFF data products were reduced with the HST science data products pipeline. In Fig. 1, we show the layout of the four clusters with RGB images obtained from the HFF observations. The cores of the four clusters were covered with MUSE, within its 1′×1′ FoV, with a 0.2″ × 0.2″ pixel-scale, a spectral range between 4650 Å and 9300 Å, a spectral resolution between R ≈ 1750 at 4650 Å and R ≈ 3750 at 9300 Å, and a spectral sampling of 1.25 Å/pixel. In our work, MUSE data was reduced with the procedure presented in Caminha et al. (2017a,b, 2019), based on the standard ESO reduction pipeline (Weilbacher et al. 2020), complemented with the Zurich Atmosphere Purge (ZAP, Soto et al. 2016) sky subtraction procedure. The MUSE footprint within the four cluster cores is shown in red in Fig. 1. Further spectroscopic follow-up campaigns were carried out with multi-object spectrographs, such as the CLASH-VLT programme (Rosati et al. 2014) with the VIsible MultiObject Spectrograph (VIMOS) on the VLT, which allowed the measurement of several thousand spectroscopic redshifts within a 20′ field centred on several strong lensing clusters, including AS1063 and M0416. In the following sub-sections, we provide some further details on the spectroscopic dataset on which this work is based, which was crucial for the construction of state-of-the-art strong lensing models for the four clusters, based on extremely large sets of spectroscopically confirmed multiple images (see e.g. Richard et al. 2021; Caminha et al. 2016, 2017a; Grillo et al. 2016; Schuldt et al. 2024), and on a detailed description of the mass structure of the cluster members enhanced by their observed kinematics (see Bergamini et al. 2019, 2023b,a; Granata et al. 2022).

Fig. 1. RGB images of the cores of the four clusters from HFF photometry. The red, green, and blue channels were obtained by combining the filters F105W + F125W + F140W + F160W, F606W + F814W, and F435W, respectively. We show the FoV of the MUSE observations with red contours, and we mark all the cluster members for which we have measured a velocity dispersion value with green circles. Finally, we show the region intersection of all HFF bands for each cluster with cyan contours.

2.1. Abell 2744

The MUSE pointings of A2744 used in this work were obtained within the GTO Programme 094.A-0115 (P.I. Richard). The four main quadrants were covered with exposure times between 3.5 h and 5 h, and a further 2 h pointing was obtained for the cluster central region. The detailed exposure times are presented in Mahler et al. (2018), where the authors measured a point-spread-function (PSF) full width at half maximum (FWHM) of 0.58″ for the combined MUSE image. The spectroscopic redshift catalogue for the field of A2744 was completed with ancillary observations. Specifically, A2744 was targeted for 4.4 h with VIMOS (PI: Böhringer, spectroscopic catalogue in Braglia et al. 2009), with the AAOmega multi-object spectrograph on the AngloAustralian Telescope (AAT; spectroscopic catalogue in Owers et al. 2011), and with the Grism LensAmplified Survey from Space survey (GLASS; Treu et al. 2015; Schmidt et al. 2014).

2.2. Abell S1063

The MUSE observations of Abell S1063 are the combination of two programmes: the first covered the SW half of the cluster core (ESO SV Programme 60.A-9345, P.I. Clement & Caputi), for a total exposure time of 3.1 h and 1.1″ seeing (Karman et al. 2015), while the second covered the NE region of the core (ESO Programme 095.A-0653, P.I. Caputi) with an exposure time of 4.8 h and 0.9″ seeing (Karman et al. 2017). Abell S1063 was also part of the CLASH-VLT sample; the MUSE and VIMOS spectroscopic redshifts were then combined with measurements from Gómez et al. (2012) and from GLASS (Treu & GLASS Team 2016) within the spectroscopic catalogues presented in Mercurio et al. (2021), which includes 1234 cluster members.

2.3. MACS J0416.1−2403

The MUSE observations of MACS J0416 have separately focussed on the NE and SW regions of its core, respectively. The former was targeted by the GTO programme 094.A0115B (P.I. Richard) and the ultra-deep observations 0100.A0763(A) (P.I. Vanzella), for a total integration time of 17.1 h in most of the field and a seeing of approximately 0.6″ (Vanzella et al. 2021), the deepest MUSE pointing for a strong lensing cluster. The latter was targeted in the 094.A0525(A) programme (P.I. Bauer), for 11 h of integration time and 1.0″ seeing. Similarly to AS1063, M0416 was targeted with VIMOS within the CLASH-VLT programme (Balestra et al. 2016).

2.4. MACS J1149.6+2223

The MUSE coverage of the central regions of M1149 was guaranteed by the DDT programme 294.A-5032 (P.I. Grillo), for a total integration time of 6h and seeing of less than 1.1″ in almost all exposures (Grillo et al. 2016), complemented by 5.5 h of observations of the north-west region of the cluster within the programme 105.20P5 (PI: Mercurio), presented in Schuldt et al. (2024). Further spectroscopic data was obtained within the GLASS programme (Treu & GLASS Team 2016).

2.5. The cluster member catalogue

The spectroscopic catalogues presented in the previous section and in the papers therein referenced, built primarily with MUSE and VIMOS data, were the starting point for the selection of the sample of cluster members that we will include in our calibrated VDF. We selected spectroscopic cluster members through analysis of the cluster 1D velocity distribution, following Bergamini et al. (2023a), Caminha et al. (2016), Caminha et al. (2017a), and Grillo et al. (2016) for A2744, AS1063, M0416, and M1149, respectively. These spectroscopic cluster members catalogues can be reliably extended even to regions not covered by spectroscopy when deep multi-band photometry is available. Angora et al. (2020) developed a convolutional neural network to identify cluster members on the basis of pan-chromatic imaging alone, achieving a purity-completeness rate larger than 90% down to a HST F814W band magnitude m₈₁₄ = 25 when all HFF bands are available. We adopted their completed cluster member catalogues, presented in Angora et al. (2020) for AS1063, M0416, and M1149, and in Bergamini et al. (2023a) for A2744, as the basis for our work, limiting ourselves to the intersection between the footprints of all HFF band for each cluster to ensure uniform completeness of the photometric samples. We show this region in Fig. 1 with cyan contours for all four clusters. The resulting cluster member catalogues contain 221 (177), 221 (146), 218 (176), and 253 (191) galaxies overall (spectroscopically confirmed) for A2744, AS1063, MACS0416, and MACS1149, respectively.

3. The velocity dispersion catalogue

The calibration of the VDF of red cluster members requires a large catalogue of robustly measured velocity dispersion values. Thanks to the deep spatially resolved MUSE spectroscopy, the line-of-sight velocity distribution (LOSVD) of the member galaxies can be inferred from the shift and broadening of their spectral features. Since our main point of interest is the determination of the velocity dispersion of the cluster members (often indicated as σ), we focussed on constraining the first two moments of the LOSVD, using the public spectral fitting code pPXF¹ (penalised pixel fitting Cappellari & Emsellem 2004; Cappellari 2017, 2023). We performed a full-spectrum fit of the continuum and of the absorption lines in the rest-frame wavelength range [3700 − 5250] Å. This spectral range contains several absorption features that can be used to constrain the velocity dispersion with high accuracy for passive members within all four clusters, as shown in Fig. 2. The optimisation simultaneously fits both the stellar population and the LOSVD of the observed galaxy by comparing its spectrum with a combination of stellar templates. These templates were selected from the high-resolution UVB stellar spectra in the X-shooter Spectral Library (XSL) DR3 (Verro et al. 2022) and were convolved with a LOSVD. Our template selection was informed by the work of Knabel et al. (2025), who identified template choice as the primary source of systematics in stellar kinematic analyses of ETGs. They reviewed all 830 XSL DR3 spectra, removing those lacking flux calibration or corrections for slit losses and extinction. The remaining templates were then visually inspected, with any suspected issues, such as poor extinction or telluric correction, high noise levels, or the presence of emission lines, flagged along with the affected wavelength regions. The final selection, aimed at ensuring reliable templates on the wavelength range we consider, includes 462 spectra, which we used for our analyses. The fits were performed using a continuum shape correction with 12th-degree additive Legendre polynomials and convolving the spectra with a Gaussian LOSVD, given that we only fitted its first two moments. On the other hand, we did not use any multiplicative polynomial to correct for dust reddening, as we were exclusively focussing on the LOSVD parameters and did not wish to carry out stellar population studies on the spectra. The results of the fit of a high S/N member of M0416 are shown in Fig. 2, as an example. We tested the reliability of the velocity dispersion measurements by building a sample of 16 000 simulated spectra of cluster members, presented in Appendix A. This set-up allows accurate measurements for S/N > 10 and down to velocity dispersion values of 50 km s⁻¹, comparable with the spectral dispersion of MUSE, as shown in Appendix A. We also used the simulated spectra to obtain a robust estimate of the uncertainty on the measured value of σ. On one hand, pPXF provides a formal uncertainty on σ based on the covariance matrix of the fitted parameters, which we re-scaled under the assumption that the best-fit χ² corresponds to the number of degrees of freedom in the fit. While this approach can yield a reasonable order-of-magnitude estimate when the fit quality is high, it tends to significantly underestimate the uncertainty in cases of noisier spectra or low σ values (Cappellari & Emsellem 2004), and in such regimes, simulations offer a more reliable means of assessing the accuracy of the recovered values. In Appendix A, we quantify the relative uncertainty on σ as a function of S/N, based on the scatter of the recovered σ values around the input velocity dispersion for the simulated spectra. Following a similar approach to Bergamini et al. (2019), we fitted a relation between the relative uncertainty δσ_sim emerging from the simulation and S/N, presented in Eq. (A.1). As δσ_sim goes to zero for S/N ≈ 80, we disregarded it for S/N > 80. To be conservative, we combined this simulation-based error in quadrature with the formal uncertainty from pPXF (δσ_pPXF) to derive our final error estimates on σ:

$$\begin{array}{c}\hfill \delta \sigma =\sqrt{{\left(\delta {\sigma }_{\mathtt{pPXF}}\right)}^2+{\left(\delta {\sigma }_{\mathrm{sim}}\right)}^2}.\end{array}$$(2)

Fig. 2. Line-of-sight velocity dispersion fit with pPXF and main spectral absorption features contributing to the measurement. We chose the cluster member ID 82063 of M0416 (RA = 64.0401169, Dec = −24.0659218, S/N = 133). The regions masked as potentially affected by issues with the subtraction of sky features and emission lines are shaded in grey. We note that the mask in the spectral region around 4200 Å is intended to avoid the effects of the adaptive-optics laser emission used for some MUSE exposures of M0416, and as such was not adopted for the other three clusters. We show the observed spectrum in black, the best fit in red, and the fit-measured spectral residuals in green. The main absorption lines are marked with red vertical dashed lines.

As is clear from the results of our simulations, a sufficient spectral S/N is the most critical aspect to guarantee a reliable velocity dispersion measurement. In order to increase the sample of galaxies with S/N > 10, we weighted the MUSE cube for the observed surface brightness of the cluster members (as first done by Horne 1986). We can thus extract the spectra of the member from a larger aperture around their centre, reducing the flux loss that would be given by a smaller aperture, while the brighter pixels, with the highest S/N, contribute the most to the overall spectrum. Specifically, we used the HFF surface brightness in the F814W band, given that it corresponds to a wavelength range where several of the absorption features used for the pPXF fit of the galaxy spectra are found, at the redshift of the four clusters. We measured the HFF and MUSE point-spread-functions (PSFs), and degraded the HFF F814W image to the MUSE resolution and pixel-scale, in order to obtain a weight map for the MUSE cube.

The spectra of all cluster members within the MUSE footprint were extracted summing all pixels within a 1.5″ radius from their centre of light, weighted as described earlier. To guarantee a reliable velocity dispersion measurement, all spectra with S/N < 10 (equivalent to S/N < 8.94 Å⁻¹, given MUSE spectral sampling of 1.25 Å) were removed from the sample. Given the rather large extraction apertures, all members with a source falling within 3″ from their centre were visually checked and excluded if they showed any significant photometric or spectroscopic contamination from nearby objects. All fits were also visually checked to make sure that the spectral absorption features had been correctly recognised and well fitted. These selections leave us with a catalogue of 76 measured velocity dispersion values out of 176 members with a MUSE spectroscopic redshift for A2744, 34 out of 94 for AS1063, 52 out of 158 for M0416, and 51 out of 128 for M1149. In total, we have a catalogue of 213 measured velocity dispersion values, which we report in Appendix C. The cluster members for which we have a reliable velocity dispersion measurements are shown with green circles in Fig. 1.

The weighting technique for spectral light is crucial in collecting such a large sample, as it allows a significant S/N increase. We compare the S/N of our spectra obtained from the weighted cube with that found extracting them directly from the MUSE cube within an aperture of 0.8″ from the galaxy centre of light, which has been found to maximise the S/N (see Bergamini et al. 2019, 2023a). As is clear from Fig. 3, our procedure allows a significant increase in the S/N, typically in the range 10 − 30% for S/N = 10 − 20, allowing a larger number of objects with S/N > 10.

Fig. 3. Relative increase in the value of the S/N of the spectra of cluster member galaxies when weighting with their surface brightness in the HST F814W band. We compare the S/N of the spectrum extracted within an aperture of 1.5″ from the galaxy centre of light from the weighted MUSE cubes with the S/N of the spectrum extracted for the same galaxy within an aperture of 0.8″ directly from the (non-weighted) MUSE cube.

4. The Fundamental Plane relation for the four clusters

Observations show that ETGs follow a tight scaling relation linking the kinematic properties, which descend from their total mass distribution, to the luminosity, which depends on the distribution of the stars. This empiric law is known as the FP (Dressler et al. 1987; Djorgovski & Davis 1987; Bender et al. 1992),

$$\begin{array}{c}\hfill log{R}_e=\alpha log{\sigma }_0+\beta {\mathrm{SB}}_e+\gamma ,\end{array}$$(3)

where R_e is the effective radius (i.e. the projected radius enclosing half of the total luminosity), σ₀ is the central stellar velocity dispersion (i.e. measured within an aperture with a projected radius R_e/8, or σ₀ = σ_ap(R_e/8)), and SB_e is the average surface brightness within R_e. For the calibration of the FP, we complemented the measured velocity dispersion values described in Sect. 3 with measured structural parameters. This allowed us to leverage the accurate information we obtained on the morphology of all the members towards accurate insights into their kinematics, in spite of the lack of a MUSE spectrum.

We measured the structural parameters for all the members in all HFF bands using the Python package MORPHOFIT (Tortorelli & Mercurio 2023; Tortorelli et al. 2018, 2023), based on SExtractor (Bertin & Arnouts 1996) and Galfit (Peng et al. 2010a), and especially developed to automatically fit the surface-brightness profile of a large number of galaxies in a dense, crowded environment. To properly account for the effect of the intra-cluster light in the cluster core, after an initial SExtractor run to inform the starting points for the light profiles parameters, the code proceeds iteratively in cutouts of increasing size, from stamps around the single galaxies to the full HFF images. All measured structural parameters are available at the CDS and presented in Appendix B. The structural parameters of the cluster members spectroscopically confirmed by MUSE in AS1063, M0416, and M1149 were presented by Tortorelli et al. (2023), whilst the remaining have been measured for the first time in this work. In Sect. 5 our structural parameter catalogue is crucial in order to identify red cluster members and to assign a velocity dispersion value for the object whose spectrum does not have a sufficient S/N for a measurement, and for those that fall outside of the MUSE pointing.

As we wish to compare the FP relations for the four clusters, we need to use magnitude and effective radius measurements in the same rest-frame band. We chose the rest-frame SDSS r band since it is mostly covered, especially in the case of A2744 and AS1063, by the HST F814W band, which guarantees a good S/N down to magnitudes of around 28 and has a small PSF FWHM. Firstly, we performed a K correction of the observed F814W magnitudes to the rest-frame SDSS r band, defined as (e.g. Hogg et al. 2002)

$$\begin{array}{c}\hfill K_{814,r}(z)=\frac{1}{1+z}\frac{\int \mathrm{d}{\lambda }_o{\lambda }_{o}L\left(\frac{{\lambda }_o}{1+z}\right)f_{814}({\lambda }_0)}{\int \mathrm{d}{\lambda }_o{\lambda }_{o}g({\lambda }_o)f_{814}({\lambda }_0)}\frac{\int \mathrm{d}{\lambda }_e{\lambda }_{e}g({\lambda }_e)f_r({\lambda }_e)}{\int \mathrm{d}{\lambda }_e{\lambda }_{e}L({\lambda }_e)f_r({\lambda }_e)},\end{array}$$(4)

where z is the cluster redshift, λ_o and λ_e are the observed- and the rest-frame wavelength, L(λ) is the galaxy flux as a function of wavelength, and g(λ) is the standard source for AB magnitudes. We calibrated the correction using L(λ) from a template spectrum built combining ETGs from our MUSE observations. Similarly to Grillo (2010), we obtained the effective radius in the SDSS r band as a linear interpolation of those in the HST F814W and F105W bands,

$$\begin{array}{c}\hfill R_{e,r}=(1+z)(\frac{{\lambda }_r-{\lambda }_{105}/(1+z)}{{\lambda }_{814}-{\lambda }_{105}}{R}_{e,814}+\frac{{\lambda }_r-{\lambda }_{814}/(1+z)}{{\lambda }_{814}-{\lambda }_{105}}{R}_{e,105}),\end{array}$$(5)

where λ_r, λ₈₁₄, and λ₁₀₅ are the effective wavelengths of the SDSS r band, HST F814W and F105W bands, respectively. We calibrated the FP, as defined in Eq. (3), using R_e measured in kiloparsecs, as described in Eq. (5), whilst SB_e is measured in mag arcsec⁻² as

$$\begin{array}{c}\hfill {\mathrm{SB}}_e=m_{814}+K_{814,r}(z)+2.5log(2\pi R_e^2)-10log(1+z),\end{array}$$(6)

where K_814, r(z) is defined as in Eq. (4), the corrected R_e is in arcsec, and the last term takes into account cosmological dimming.

We selected the ETGs used to calibrate the FP from their red colour and early-type morphology. Following Tortorelli et al. (2018, 2023), we calibrated a colour-magnitude relation between the total magnitude in the F814W band and the F606W−F814W colour, for all galaxies included in the photometrically completed cluster member catalogue presented in Sect. 2.5:

$$\begin{array}{c}\hfill m_{606}-m_{814}=a+b{m}_{814}.\end{array}$$(7)

We selected red galaxies as those redder than the best-fit CM relation minus 1σ, with the exception of A2744, where we included galaxies redder than the best-fit CM relation minus 2σ, as the scatter provided by the fitting procedure is significantly smaller than all other clusters and appears underestimated. The best-fit parameters of the colour-magnitude relation for the four clusters are presented in Table 1, and are compatible with those found by Tortorelli et al. (2018) at 1σ level for AS1063 and at 2σ level for M1149. We also followed the criterion adopted by Tortorelli et al. (2018) to identify ETGs from their morphology, requiring a Sérsic index n > 2.5. This left us with a sample of 49, 25, 34, and 25 red early-type members for A2744, AS1063, M0416, and M1149, respectively, that can be used as a basis for the calibration of the FP.

As for the velocity dispersion, we need to correct the measured values from the weighted spectra, which we indicate as σ_w to the corresponding central value σ₀ adopted to calibrate the FP. We adopted a power-law aperture correction (Jorgensen et al. 1996)

$$\begin{array}{c}\hfill {\sigma }_{\mathrm{ap}}(R)={\sigma }_{\mathrm{ap}}(R\prime ){(R/R\prime )}^{{\alpha }_{\mathrm{circ}}},\end{array}$$(8)

where the exponent values α_circ were determined following Zhu et al. (2023), and depend on the r-band absolute magnitude, g − i colour, and Sérsic index of the member. The first step required is understanding how to interpret σ_w in terms of a physical aperture $\widehat{R}$ such that ${\sigma }_w={\sigma }_{\mathrm{ap}}(\widehat{R})$. As is clear from Fig. 4, we find that, on average, σ_w is a good approximation of the velocity dispersion within the effective radius σ_e = σ_ap(R_e), where σ_e has been measured correcting the measured σ_ap(0.8″) through Eq. (8). Specifically, σ_w − σ_e is, on average, for the four clusters, −1.1 km s⁻¹ with a standard deviation of 7.9 km s⁻¹. The weighting technique for spectral light adopted to increase S/N gives more weight to the central pixels, and the resulting weighted combination of their spectra is naturally representative of the region within the effective radius. This represents a second advantage of the procedure, as this region is usually too small to be probed directly from the MUSE cube: R_e typically encircles just a few MUSE pixels and the resulting spectral S/N is very small. In conclusion, we assume σ_e = σ_w for all members and use Eq. (8) to derive

$$\begin{array}{c}\hfill {\sigma }_0={\sigma }_{\mathrm{ap}}(R_e/8)={\sigma }_e\times {(1/8)}^{{\alpha }_{\mathrm{circ}}}={\sigma }_w\times {(1/8)}^{{\alpha }_{\mathrm{circ}}},\end{array}$$(9)

Fig. 4. Velocity dispersion values of the cluster members from weighted spectra (σw) compared with the velocity dispersion within the effective radius (σe). The weighted spectrum was extracted within an aperture of 1.5″ from the galaxy centre of light from the MUSE cubes weighted with the HST F814W surface brightness.

measured in km s⁻¹.

Due to the large number of members of the cluster with a measurement of σ₀, R_e, and SB_e, we can calibrate the FP for the four clusters with an optimisation based on the code ltsfit (Cappellari et al. 2013). In Table 2, we report the best-fit values of α, β, and γ, as defined in Eq. (3), while in Fig. 5 we show the best-fit FP relations and the scatter around them of the members used for their calibration, for the four clusters. We note that, after the corrections performed on the magnitudes and effective radii of the members to obtain comparable values referred to a rest-frame r band, the values of α and β are compatible across the four clusters, with the exception of a lower β value for M1149. We find larger uncertainties on the FP parameters for AS1063 and M1149, for which our sample is less inclusive of fainter galaxies, due to the shallower MUSE data. Our α and β values are also compatible with the FP relation calibrated in the same band (r) for large samples of ETGs observed by the SDSS and the Dark Energy Spectroscopic Instrument (DESI) Peculiar Velocity Survey. Specifically, SDSS and DESI report α values of 1.274 and 1.177, and β values of −0.841 and −0.793, for a sample of 34059 and 3110 ETGs, respectively (see Howlett et al. 2022; Said et al. 2025), and with the FP relation found by D’Addona et al. (2025) in the HST F814W band for a strong lensing cluster at a similar redshift (PLCK G287.0+32.9, z = 0.3833).

Fig. 5. Distribution around the FP relations of the galaxies used for the calibration for the four clusters. The FP relations are shown on the plot diagonal. The values of α, β, and γ for the four clusters are reported in Table 2.

The relation between the zero-point of the FP and redshift can be interpreted as a tracer of the evolution of the mass-to-light ratio (M/L) of galaxies, under the assumption that ellipticals are a passively evolving homologous population, so that the M/L evolution with respect to z = 0 is linked to the variation of γ with respect to a FP at z = 0 by Δlog M/L (z) = Δγ(z) (see Faber et al. 1987; van Dokkum & Stanford 2003; Holden et al. 2005; van Dokkum & van der Marel 2007; van de Sande et al. 2014). Similarly to van de Sande et al. (2014) we fitted the FP relations assuming α = 1.20 and β = −0.83, as shown in Fig. 6. We find Δγ = 0.28 between z = 0.307 and z = 0.542, which is larger than the amount predicted by the relation found by van de Sande et al. (2014), who reported Δlog M/L (z) = (−0.49 ± 0.03)z. Whilst we note that Tortorelli et al. (2018) also found that the zero-point evolution for the Kormendy relation between AS1063 and M1149 is larger than expected when considering only passive luminosity evolution, the shift between planes is of the order of the 1σ scatter of the galaxies used for their calibration, limiting the robustness of any conclusion. The study of the evolution of the FP is beyond the scope of this work, due to the limited redshift interval considered, and should be addressed with a larger sample, spanning a wider z range.

Fig. 6. Evolution of the zero-point of the FP relation. We show the FP relations fit for the four clusters assuming α = 1.20 and β = −0.83, and the 1σ scatter of the galaxies used for its calibration about the plane.

5. The velocity dispersion function of red cluster member galaxies

The completeness of the velocity dispersion catalogue presented in Sect. 3 is limited by the smaller MUSE FoV, which only covers the central regions of the cluster core, as well as by the exclusion of members in the region covered by MUSE, but whose spectrum is contaminated by a nearby object. In this section, we take advantage of the FP relations calibrated in Sect. 4, and of the larger sets of cluster members described in Sect. 2.5, obtained with a combination of spectroscopic and photometric data, to extend the velocity dispersion catalogue to the red members falling within the HFF footprint, including the area not covered by MUSE. As detailed in Sect. 2.5, our sample of cluster members is complete across the intersection of all HFF bands, marked in cyan in Fig. 1. Whilst the FP strictly only applies to morphologically selected ETGs, we found that including all red galaxies in the FP calibration leads to compatible parameter values for all four clusters, albeit with a larger scatter than with a stricter morphological selection. This is analogous to what found by Tortorelli et al. (2018) for the Kormendy relation. This suggests that the FP describes in good approximation all red cluster galaxies, and as such we used the calibrated FP relation for each cluster to infer the central velocity dispersion σ₀ for all the red cluster members falling within the HFF images and use the scatter about the plane to estimate the uncertainty on the recovered values. Similarly to what found by Tortorelli et al. (2018), our cluster member sample selected by colour is made up almost exclusively by passively evolving galaxies: fewer than 3% of the red cluster members with a high S/N MUSE spectra (S/N > 10) shows clear emission lines.

5.1. Calibration of the velocity dispersion function

We selected red members with the same criteria described in Sect. 4, using the colour-magnitude relations for the four clusters presented in Eq. (7) and Fig. 1, for a total of 194, 169, 172, and 188 red members for A2744, AS1063, M0416, and M1149, respectively. For these objects we derived log σ₀ from log R_e and SB_e through Eq. (3). The measured values of R_e and SB_e are converted into effective values in the rest-frame r band using Eqs. (4) and (5). As detailed, our sample of cluster members is highly complete up to m₈₁₄ = 25, so our VDF may suffer from biases due to incompleteness in the bins where members with m₈₁₄ > 25 fall. We find that setting a threshold log σ₀ [km s⁻¹] > 1.5 ensures a sample completely below m₈₁₄ = 25 – even for M1149, the highest-redshift cluster in our sample. This limit is also very close to the minimum measured velocity dispersion value of our catalogue, thus avoiding extrapolating the FP relations beyond the velocity dispersion range used for their calibration.

Having estimated σ₀ for the red members, we focus on the VDF for each cluster. As detailed, we find an increased scatter about the FP for the full population of red cluster galaxies compared to the sample selected to calibrate it. We used this scatter to evaluate the uncertainty on the σ₀ values inferred through the FP and how it propagates on the VDF. For each cluster, we took the 1σ scatter of all red galaxies about the FP in terms of log σ₀ and used it to extract 100 values of σ₀ from a Gaussian distribution for all red member galaxies. After choosing some bins in log σ₀, this provided us, for each cluster, with 100 realisations of the number counts of red members in a given velocity dispersion bin. Finally, with the median, and with the 16th and 84th percentiles of the distribution of number counts, we estimated the best value for the VDF in each bin and its 1σ uncertainty.

Following Sohn et al. (2020), we describe the VDF of red cluster galaxies with a Schechter function (Schechter 1976)

$$\begin{array}{c}\hfill \varphi (V)={\varphi }^{\ast }{10}^{\alpha (V-log{\sigma }^{\ast })}exp(-{10}^{(V-log{\sigma }^{\ast })}),\end{array}$$(10)

where ϕ(V) is the number of galaxies falling in a bin centred in V = log σ₀, α is the slope of the power-law regime describing the faint end of the function (for V ≪ log σ^*) and log σ^* marks the transition to an exponentially decreasing regime for the bright end. Similarly to Sohn et al. (2020), we find that we cannot simultaneously constrain all parameters of a modified Schechter function parametrisation (Eq. 1) without parametric degeneracy. Compared with Eq. (1), we are thus setting β = 1, and evaluating the function on finite V bins. The normalisation parameter ϕ^* depends on the richness of the cluster, as well as on the ratio between the cluster region covered by all HFF bands, where our cluster member catalogue is complete, and the area within the cluster virial radius R₂₀₀. Focussing on AS1063, M0416, and M1149, Angora et al. (2020) found an approximate self-similarity between the radial distribution of the cluster members when the number counts for a radial bin at radius R are normalised by the total number of members within 0.15 R₂₀₀, N₀, and by the ratio between the areas enclosed within the radius R and within R₂₀₀. In other words, the number of members in a given area within the cluster core is approximately the same when normalised by the fraction of the virial area covered and the richness of the cluster core. We adopted this approach to allow a comparison between the velocity functions of the four clusters, normalising the VDF in each cluster by

$$\begin{array}{c}\hfill \stackrel{\sim }{N}=N_0\frac{A_{\mathrm{HFF}}}{\pi R_{200}^2},\end{array}$$(11)

where A_HFF is the area covered by the intersection of all HFF bands. For clarity, we use $\stackrel{\sim }{\varphi }$ to refer to the value of ϕ^* for the normalised Schechter function. We used the values of R₂₀₀ by Medezinski et al. (2016) for A2744 and those by Umetsu et al. (2018) for the remaining three clusters, following Angora et al. (2020). Given that the region covered by HFF does not extend fully to 0.15 R₂₀₀ for all clusters, we chose the number of red members within 0.1 R₂₀₀ as a proxy of richness (N₀). We report the values of $\stackrel{\sim }{N}$ in Table 3. We remark that this self-similarity only holds approximately and is potentially influenced by the morphology and projected structure of the clusters. We fitted Eq. (10) with an affine-invariant Markov chain Monte Carlo (MCMC; Goodman & Weare 2010; Foreman-Mackey et al. 2013) algorithm. We explored the posterior probability distribution for the values of $\stackrel{\sim }{\varphi }$, α, and log σ^* as defined in Eq. (10) and in the text, expressed as a product between a likelihood

$$\begin{array}{c}\hfill ln(\mathcal{L}(\varphi (V)|V,\delta \varphi (V),\stackrel{\sim }{\varphi },\alpha ,log{\sigma }^{\ast }))=-\frac{1}{2}{\displaystyle \sum _{i=1}^{N_{\mathrm{bins}}}}\frac{{(\varphi (V_i)-\widehat{\varphi }(V_i))}^2}{{(\delta \varphi (V_i))}^2},\end{array}$$(12)

and flat priors for all the free parameters, chosen to limit the parameter space. In Eq. (12), ϕ(V_i) and δϕ(V_i) are the median and 1σ uncertainty of the distribution of number counts in a given bin V_i in log σ₀ space, normalised by the factor $\stackrel{\sim }{N}$, introduced in Eq. (11), whilst $\widehat{\varphi }(V_i)$ are the model-predicted values in the same bin for a Schechter function model with parameters $\stackrel{\sim }{\varphi }$, α, and log σ^*. We notice that some bins in the high-V regime only contains a few objects, and that this causes the recovered Schechter function parameters to be sensitive to the binning choice. To take this source of uncertainty into account, we took nine realisation of the bins {V_i}_i = 1^N_bins, shifted by {−0.5;−0.375;−0.25;−0.125;0;0.125;0.25;0.375;0.5} × ΔV, where ΔV is the bin width. Then we repeated the MCMC optimisation for each bin shift. The final posterior probability distributions for the parameters are the sum of the nine distributions obtained for each of the realisation. We took the median, 16th, and 84th percentiles of the sampled posterior probability distributions to estimate the best-fit values of the optimised parameters for all clusters and their uncertainties, reported in Table 3. The posterior probability distribution of the three parameters for the four clusters, as sampled by MCMC chains, is instead shown in Fig. 7. The values of α and log σ^* are always compatible at 1σ or 2σ level, and Fig. 7 shows a clear overlap and alignment of the corner-plots between them, suggesting a similar correlation for the two parameters across the four clusters. In general, the results favour a positive α slope in the range 0.55 − 1.60 and log σ^* [km s⁻¹] between 2.18 and 2.47, corresponding to σ^* values of 151 − 295 km s⁻¹. On the other hand, the distributions of $\stackrel{\sim }{\varphi }$ are brought closer by the normalisation adopted (see Eq. 11), but the scatter between the clusters remains more significant, which points towards intrinsic differences in the cluster richness. We extracted correlated sets of parameters from the posterior probability distributions shown in Fig. 7 to plot the VDFs with their propagated uncertainty in Fig. 8. The bin shifts in the optimisation procedure are reflected in the bin width in Fig. 8. We note that the posterior probability distributions for M0416 appear broader and occasionally double-peaked. This behaviour arises because, for some bin configurations, the highest-V bin becomes empty and is consequently excluded from the fit – giving rise to a second peak in the parameter space for the posterior probability distribution. As such, uncertainty on the parameters for M0416 may be overestimated. We also note that M0416 is a distinctly bi-modal cluster (see Fig. 1), with two sub-clusters clearly emerging, which might determine a different evolution history for its member galaxies, affecting the VDF.

Fig. 7. Posterior probability distribution for the parameters of the Schechter function fit of the VDF for the four clusters, as sampled by MCMC chain. The parameters are defined in Eq. (10); the colour scheme identifies the four clusters following Fig. 8.

Fig. 8. Velocity dispersion function for the four clusters. We show the data bins and the Schechter function fit with their uncertainties, as defined in Eq. (10).

5.2. Comparison between the calibrated functions and with other works

As pointed out in the previous sub-section, the best-fit values of α and log σ^* for a Schechter function-fit of the normalised VDFs are compatible across the four clusters. However, some differences between the four VDFs emerge when looking at the higher-velocity-dispersion regime. Whilst for log σ₀ [km s⁻¹] < 2.2 data bins for the four functions are compatible and very close to one another, M1149 shows some excess for log σ₀ [km s⁻¹] > 2.2. Whilst any conclusion is limited by the small cluster sample, it is interesting to note how M1149 has a significantly higher virial mass (see Table 3) and redshift compared to the other three. This may be a trace of major dry dissipationless mergers (see e.g. Khochfar & Burkert 2003; Boylan-Kolchin et al. 2006; De Lucia & Blaizot 2007; Naab et al. 2009) reducing the number of high-velocity-dispersion galaxies in lower-z clusters. This hypothesis is also suggested by the higher BCG velocity dispersion values we find for the lower-z clusters A2744 and AS1063 (see Tables C.1, C.2, C.3, C.4), which can also be a hint of major dry mergers affecting the BCG structure (Cox et al. 2006).

The integral-field nature of MUSE allows us to simultaneously measure the spectra of all cluster members in the core with very few pointings, efficiently allowing deep exposures. As shown in Fig. 9, we can reliably extend our catalogue to values well below 100 km s⁻¹, corresponding to a total rest-frame SDSS r-band luminosity ≳10⁹ L_⊙, probing galaxy mass and luminosity regimes unexplored at these redshifts. Figure 9 also allows us to probe the slope of power-law scaling relations between σ and the total luminosity, akin to those adopted in SL, albeit in a different band: we find that a fit of σ ∝ L_r^α gives α values in the range 0.30 − 0.38 for the four clusters, compatible with those found by cluster-scale SL models (see Granata et al. 2025). Such low-luminosity targets were inaccessible to the velocity dispersion survey from Hectospec spectra used to calibrate the VDF in Sohn et al. (2020) for a set of nine strong lensing clusters in the redshift range 0.18 − 0.29 (see Fig. 10). Looking at the VDF parameters, Sohn et al. (2020) found higher values of α, between 1.75 and 4.40, and lower values of log σ^*, in the range 1.46 − 1.72, compared to our work. Due to the stricter lower limit log σ [km s⁻¹] > 2.05 on the velocity dispersion values they consider, the recovered parameter values may be influenced by the fact that log σ^* falls outside the velocity dispersion range considered. In our work, on the other hand, the depth of the MUSE observations has allowed us to extend the VDF to a wider velocity dispersion range, accurately sampling its power-law regime. Sohn et al. (2020) also found that higher-mass galaxy clusters have higher σ values for the BCG, but fewer high-σ members, which may be a trace of major dry mergers, similarly to what we found for lower z. We do not find any such trend with the virial mass values of the clusters (2.06 × 10¹⁵ M_⊙, 1.98 × 10¹⁵ M_⊙, 1.14 × 10¹⁵ M_⊙, and 2.89 × 10¹⁵ M_⊙ for A2744, AS1063, M0416, and M1149, respectively) but we should remark that any trend emerging from a comparison between the four clusters is significantly limited by the small sample.

Fig. 9. Measured velocity dispersion as a function of rest-frame SDSS r-band luminosity for our work.

Fig. 10. Relative uncertainty on the measurement of the velocity dispersion. Comparison between the σw values presented in this work and the Hectospec observations of nine strong lensing clusters (Sohn et al. 2020, in cyan).

6. Summary and results

We have presented a detailed and comprehensive study of the kinematic properties of the red member galaxies in the cores of four strong lensing clusters at z = 0.307 − 0.542 chosen from the HFF sample, for which exquisite photometry and deep integral-field spectroscopy from MUSE are available. Taking advantage of the MUSE data, we built a kinematic pipeline based on the spectral fitting code pPXF aimed at reliably measuring the velocity dispersion of cluster members falling within the instrument FoV. We tested our pipeline with a library of 16 000 simulated spectra (see Appendix A), with the objective of a more accurate assessment of the uncertainty on our measurements. Using a library of 462 XSL DR3 spectra, selected in Knabel et al. (2025), we found that pPXF allows us to measure σ with a statistical uncertainty consistently below 5% when the spectral S/N is higher than 10. For higher S/N spectra, even lower uncertainties can be achieved for the brightest galaxies, potentially of the order of 1%.

We introduced a spectral light weighting technique that exploits the HFF images to increase the S/N and introduce a proxy for the galaxy light distribution, which acts as an approximate aperture correction to the effective radius of the member galaxy. With the pipeline in place, we built an extensive catalogue of 213 measured velocity dispersion values across the four clusters, which we report in Appendix C. This sample is unprecedented in terms of combined depth and accuracy for cluster galaxies at this redshift. The pipeline and the procedure presented in this work can be applied to gain kinematic priors to calibrate the scaling relations that describe the total mass distribution of cluster members in strong lensing models of galaxy clusters and has been used in several studies for this purpose (Bolamperti et al. 2024; Granata et al. 2025; Ertl et al. 2025; Acebron et al. 2025; Abriola et al. 2025; Bianchetti et al., in prep.).

We took advantage once again of the deep multi-band HFF photometry to extend the catalogue of measured velocity dispersion to all red members, including those whose spectrum is contaminated by a nearby source or has an insufficient S/N, and those falling outside of the footprint of the MUSE data cubes but within that of the HFF observations. We measured the structural parameters using MORPHOFIT for an extensive sample of 723 red cluster members within the four clusters, selected either spectroscopically or photometrically using a convolutional neural network. The sample is highly pure and complete down to m_F814W = 25. In order to do so, we calibrated the FP relation for all four clusters in the rest-frame r band, using K-corrected magnitudes and effective radii obtained from the multi-band structural parameters. The best-fit values of the parameters α and β of the FP relations for the four clusters are compatible with each other (see Table 2), and with FP relations calibrated for samples of ETGs in the same band. The third parameter γ shows the expected signs of evolution with redshift. The offset we find across the redshift range considered is larger than expected under the effect of passive luminosity evolution only, similarly to what was found for the Kormendy relation zero-point (Tortorelli et al. 2018), and could point to hints of size evolution. However, the small sample and the z range considered prevent us from drawing any strong conclusions.

With the calibrated FP relations and their scatter, we assigned a velocity dispersion to all 723 red cluster members and built a VDF for each cluster, extended down to log σ₀ [km s⁻¹] = 1.5. We used the number of members within R = 0.1R₂₀₀ and the ratio of the area covered by HFF observations to the virial area to normalise the functions, and we fitted a Schechter function to the relations; we found that they are well described by a positive α slope with values in the range 0.55 − 1.60 and log σ^* [km s⁻¹] between 2.18 and 2.47. The two parameters α and log σ^* are compatible across the four clusters; however, the highest-z cluster of the sample, M1149, shows an excess for log σ₀ [km s⁻¹] > 2.2, which might be explained by major dry mergers reducing the number of high-σ₀ member galaxies at lower redshifts, whilst increasing the velocity dispersion of the BCG. On the other hand, we find significant differences with α and log σ^*, as found for nine strong lensing clusters at z = 0.18 − 0.29 by Sohn et al. (2020), which may depend on limitations descending from their lower bound log σ₀ [km s⁻¹] > 2.05. Comparing our results with the local massive clusters Coma and Abell 2029 (Sohn et al. 2017a) to probe the redshift evolution of the cluster-core VDF is not straightforward, given the different parametrisation chosen. However, a side-by-side comparison reveals broadly similar slopes for the VDFs in the low-σ power-law regime. This suggests that the abundance of intermediate-mass early-type galaxies in cluster cores does not evolve strongly between z ∼ 0 and the redshift range probed by the HFF clusters. Our conclusions on the VDF are, in any case, limited by the small sample, as highlighted by the cluster-to-cluster scatter within our sample, and by the limited coverage of the HFF observations, extended to R ≲ 0.15R₂₀₀, but explore velocity dispersion regimes mostly uncharted in previous studies.

The deep pan-chromatic observations and MUSE integral-field spectroscopy now available for the cores of massive strong lensing clusters have allowed us to carry out an extremely detailed and comprehensive study of the kinematics of red cluster galaxies, accurately probing the depth of the potential wells of the sub-haloes that host them. Luminosity-function-based studies depend on stellar population properties (age, metallicity, star formation history) and dust attenuation, and calibrating the stellar mass function of cluster galaxies requires assumptions on their mass-to-light ratios and initial mass functions (e.g. Bernardi et al. 2010; Montero-Dorta et al. 2017). Furthermore, the luminosities and mass-to-light ratios of cluster galaxies are affected by quenching effects arising from the dense cluster environment (e.g. ram-pressure stripping and starvation, see Gunn & Gott 1972; Larson et al. 1980; Peng et al. 2010b), which do not substantially modify their total gravitational potential, thus making the luminosity and stellar-mass functions environment dependent tracers of total mass. In contrast, the central stellar velocity dispersion of galaxies is inferred directly from the width of their spectral absorption lines, and closely traces the depth of the gravitational potential in their inner regions, providing a much more straightforward tracer of the total mass distribution of cluster galaxies, and thus of the properties of the DM haloes that host them. This allows a direct comparison with the mass properties of cluster galaxies as inferred by strong lensing models, which are likewise sensitive to the structure of their inner potential wells. This will crucially inform attempts at resolving the persistent excess of galaxy-galaxy strong lensing events in galaxy clusters with respect to cosmological hydrodynamical simulations (Meneghetti et al. 2020, 2022), a discrepancy that also depends on the detailed mass structure of cluster sub-haloes. Whilst this work focussed on a small sample with homogeneous photometry, extending it to a larger sample, spanning a wider redshift range, will allow us to test the robustness of our conclusions and meaningfully probe emerging evolutionary trends, thus constraining galaxy formation and evolution scenarios in clusters.

Data availability

The catalogues of structural parameters presented in Appendix B are available at the CDS via https://cdsarc.cds.unistra.fr/viz-bin/cat/J/A+A/709/A254

Acknowledgments

We thank Bodo Ziegler for the useful feedback and comments that helped to improve the quality of the manuscript. We acknowledge support from the Italian Ministry of University and Research through grant PRIN-MIUR 2020SKSTHZ. PB acknowledges support from the Italian Space Agency (ASI) through contract “Euclid – Phase E”, INAF Grants “The Big-Data era of cluster lensing” and “Probing Dark Matter and Galaxy Formation in Galaxy Clusters through Strong Gravitational Lensing”. This work uses the following software packages: Astropy (Astropy Collaboration 2013, 2018), matplotlib (Hunter 2007), NumPy (van der Walt et al. 2011; Harris et al. 2020), pPXF (Cappellari & Emsellem 2004; Cappellari 2023), MORPHOFIT (Tortorelli & Mercurio 2023; Tortorelli et al. 2023), Python (Van Rossum & Drake 2009), Scipy (Virtanen et al. 2020).

References

Appendix A: Testing the accuracy of the velocity dispersion measurements with simulated spectra

In this appendix we test the reliability of the kinematic measurements performed with pPXF on a set of simulated MUSE spectra of member galaxies. Here and in the main body of the paper, we refer to the S/N measured as an average over the studied wavelength range of the ratio between the observed spectrum and its residuals with respect to the best fit found by pPXF, computed within 150 Å bins.

We generated simulated spectra of cluster galaxies using our MUSE observations as empirical templates. Specifically, we adopted the best-fit combinations of X-Shooter stellar templates derived with pPXF-fitting for a sample of 40 high S/N members of M0416 for which robust velocity dispersion measurements are available. Figure A.1 shows the quality of our spectra and the available absorption features. These spectra serve as reference bare galaxy spectra at σ = 0, with a resolution matching that of the templates (R ≈ 10000).

Fig. A.1. VLT/MUSE mean stacked normalised spectra of 14, 6, 17, and 27 member galaxies of A2744, AS1063, M0416, and M1149, respectively, with σ > 180 km s−1 and S/N > 15. The shaded regions are the standard deviation of each spectral pixel. The spectra are smoothed with a Gaussian kernel with a standard deviation of 3.75 Å. The red dashed vertical lines locate the absorption features detected in the spectra. For clearer visualisation, the normalised flux values from AS1063, M0416, and M1149 have been increased by 1, 2, and 3, respectively.

To construct the simulated sample, we randomly selected one of 40 bare spectra and re-sampled it onto a wavelength grid ten times finer to avoid interpolation artifacts. We then convolved each oversampled spectrum with a Gaussian LOSVD, parametrised by a velocity v and dispersion σ, drawn from distributions shown in Table A.1. The spectrum was redshifted to z = 0.4, degraded to the MUSE resolution (FWHM = 2.6 Å), and re-sampled to the MUSE pixel scale (Δλ = 1.25 Å). Random Gaussian noise was added to match a specified signal-to-noise ratio. Velocity values were sampled uniformly in the range v ∈ [ − 50, 50] km s⁻¹, while σ and S/N were drawn from the intervals listed in Table A.1.

Given that the sample of galaxies selected for our velocity dispersion catalogues does not push below 30 km s⁻¹, we adopted it as lower limit for our simulation sample and for the range on which we tested the reliability of our pipeline. The mock galaxy spectra were fitted in the same rest-frame wavelength range of [3700 − 5250] Å and the same pPXF set-up described in Sect. 3. As anticipated in the text, we tested different degrees for the additive Legendre polynomials used to correct for the continuum shape, and find the fits to be largely insensitive to the choice (similarly to Cappellari 2017). We fixed the degree to 12 to avoid possible oscillations of the recovered values for low-degree polynomials (as found by D’Ago et al. 2023). Both these fits and those of observed cluster members show that for S/N < 10 the fit often completely fails, or struggles to recover the observed shape of the absorption features. As such, we disregarded cluster members with a spectral S/N < 10 for our velocity dispersion measurements.

As a first test, we study the difference between the measured velocity dispersion σ_m and the input value of σ. As is clear from the left panel of Fig. A.2, the average offset is very close to zero, and on average ≈0.03 km s⁻¹ for the members with a reliable fit (S/N > 10). Unlike Bergamini et al. (2019), we find that we do not need to correct for any bias to the value of σ, even in the range σ < 100 km s⁻¹.

Fig. A.2. Velocity dispersion measurements for simulated spectra. Left panel: Relative offset of the recovered velocity dispersion with respect to the true value as a function of the S/N. We show the median relative uncertainty from Eq. (A.1) as a black line. Right panel: Absolute offset between the recovered and the true velocity dispersion values as a function of the input velocity dispersion. The points are colour-coded by the corresponding spectral S/N.

As a following step, we tested the impact of the spectral noise on the uncertainty on the velocity dispersion value. Even for S/N > 10, the noise determines random offsets between the recovered and true values of the velocity dispersion. During the LOSVD fit, pPXF determines the statistical uncertainty on the best-fit value of σ required to fit the stellar templates to the observed spectrum. As found by Cappellari & Emsellem (2004), this method can yield a significantly underestimated uncertainty value for noisy spectra and/or lower σ values. As detailed in Sect. 3, we increased this value by adding in quadrature another contribution from the random offsets of σ determined by the spectral noise. To estimate this contribution, we used our simulated spectra to study the relative offsets between the measured and observed values of σ as a function of S/N. We choose a set of S/N bins and take the standard deviation of the difference between the measured and input velocity dispersion in each bin, which we indicate as |δσ_sim|. As shown in Fig. A.2, the relative offset |δσ|/σ is typically smaller than 5% for S/N > 10 and goes to zero for S/N ≈ 10. A polynomial fit of the relative uncertainty δσ/σ as a function of the value of S/N gives the following relation

$$\begin{array}{c}\hfill |\delta {\sigma }_{\mathrm{sim}}|/\sigma =0.0705-5.1\times {10}^{-3}S/N+1.38\times {10}^{-4}{(S/N)}^2\\ \hfill -1.57\times {10}^{-6}{(S/N)}^3+6.29\times {10}^{-9}{(S/N)}^4,\end{array}$$(A.1)

where σ is the measured value of the LOS velocity dispersion provided by pPXF. In Sect. 3 we detail how we take advantage of this relation to estimate the total uncertainty on the values of the stellar velocity dispersion we report in Appendix C.

Finally, we test the impact of the optimisation set-up chosen for this work and presented in Sect. 3. We confirm that the accuracy in recovering the value of the LOSVD of the simulated spectra presented in Table A.1 is not affected by the choice of a different degree for the additive polynomials adopted for the spectral fit nor improved by the inclusion of multiplicative polynomials. We verify that masking some pixels of the wavelength range, which is necessary to avoid including contaminated regions of the observed spectra, does not bias the measured values of the velocity dispersion. This is also the case when some absorption lines are masked, as is often the case (see Fig. 2), although we notice an increase in the scatter on the recovered values when either the CaII H or K lines are not included in the spectral fit.

Appendix B: Catalogues of structural parameters

The catalogues of structural parameters for the four clusters A2744, AS1063, M0416, and M1149 are available at the CDS via anonymous ftp to cdsarc.cds.unistra.fr or via https://cdsarc.cds.unistra.fr/viz-bin/cat/J/A+A/. The table column names (in typewritten font) and their descriptions are as follows:

• ID: ID in serial order.

• RA: right ascension in degrees.

• DEC: declination in degrees.

• MAG_x: total magnitude in the x-band.

• MAG_ERR_x: total magnitude error in the x-band.

• RE_x: circularised effective radius in kiloparsecs in the x-band.

• RE_ERR_x: circularised effective radius error in kiloparsecs in the x-band.

• MU_x: average surface brightness within the effective radius in mag arcsec⁻² in the x-band.

• MU_ERR_x: average surface brightness within the effective radius error in mag arcsec⁻² in the x-band.

• N_x: Sérsic index in the x-band.

• N_ERR_x: Sérsic index error in the x-band.

• AR_x: axis ratio in the x-band.

• AR_ERR_x: axis ratio error in the x-band.

• PA_x: position angle in degrees in the x-band.

• PA_ERR_x: position angle error in degrees in the x-band.

MAG_x, N_x, AR_x, PA_x and their respective errors are the final structural parameters obtained from the fit on the full images. The weighted means of the semi-major axes in pixels and of the axis ratios in the x-band (r_e_x and AR_x, and their errors r_e_ERR_x and AR_ERR_x) are used to compute the circularised effective radii in kiloparsecs,

$$\begin{array}{c}\hfill \mathtt{RE}\_\mathtt{x}={\mathtt{r}}_{\mathtt{e}}\_\mathtt{x}\times \sqrt{\mathtt{AR}\_\mathtt{x}}\times \mathrm{pixel}\_\mathrm{scale}\times \mathrm{kpc}\_\mathrm{per}\_\mathrm{arcsec},\end{array}$$(B.1)

and their errors RE_ERR_x. These are combined with the magnitudes in the x-band MAG_x and their errors MAG_ERR_x to compute the average surface brightnesses within that radius,

$$\begin{array}{cc}\hfill \mathtt{MU}\_\mathtt{x}& =\mathtt{MAG}\_\mathtt{x}+2.5log(2\pi )+5log(\mathtt{RE}\_\mathtt{x}/\mathrm{kpc}\_\mathrm{per}\_\mathrm{arcsec})\hfill \\ \hfill & -10log(1+z_{\mathrm{cluster}}),\hfill \end{array}$$(B.2)

and their errors MU_ERR_x.

Appendix C: Catalogues of measured velocity dispersion values

In this appendix we present the final kinematic catalogues for the four clusters. We include 76 members of A2744, 34 of AS1063, 52 of M0416, and 51 of M1149. The procedure for the kinematic measurements is described in Sect. 3. The IDs and coordinates of the members follow those reported in Appendix B.

All Tables

All Figures

Fig. 1. RGB images of the cores of the four clusters from HFF photometry. The red, green, and blue channels were obtained by combining the filters F105W + F125W + F140W + F160W, F606W + F814W, and F435W, respectively. We show the FoV of the MUSE observations with red contours, and we mark all the cluster members for which we have measured a velocity dispersion value with green circles. Finally, we show the region intersection of all HFF bands for each cluster with cyan contours.

In the text

Fig. 2. Line-of-sight velocity dispersion fit with pPXF and main spectral absorption features contributing to the measurement. We chose the cluster member ID 82063 of M0416 (RA = 64.0401169, Dec = −24.0659218, S/N = 133). The regions masked as potentially affected by issues with the subtraction of sky features and emission lines are shaded in grey. We note that the mask in the spectral region around 4200 Å is intended to avoid the effects of the adaptive-optics laser emission used for some MUSE exposures of M0416, and as such was not adopted for the other three clusters. We show the observed spectrum in black, the best fit in red, and the fit-measured spectral residuals in green. The main absorption lines are marked with red vertical dashed lines.

In the text

Fig. 3. Relative increase in the value of the S/N of the spectra of cluster member galaxies when weighting with their surface brightness in the HST F814W band. We compare the S/N of the spectrum extracted within an aperture of 1.5″ from the galaxy centre of light from the weighted MUSE cubes with the S/N of the spectrum extracted for the same galaxy within an aperture of 0.8″ directly from the (non-weighted) MUSE cube.

In the text

	Fig. 4. Velocity dispersion values of the cluster members from weighted spectra (σw) compared with the velocity dispersion within the effective radius (σe). The weighted spectrum was extracted within an aperture of 1.5″ from the galaxy centre of light from the MUSE cubes weighted with the HST F814W surface brightness.
In the text

	Fig. 5. Distribution around the FP relations of the galaxies used for the calibration for the four clusters. The FP relations are shown on the plot diagonal. The values of α, β, and γ for the four clusters are reported in Table 2.
In the text

	Fig. 6. Evolution of the zero-point of the FP relation. We show the FP relations fit for the four clusters assuming α = 1.20 and β = −0.83, and the 1σ scatter of the galaxies used for its calibration about the plane.
In the text

	Fig. 7. Posterior probability distribution for the parameters of the Schechter function fit of the VDF for the four clusters, as sampled by MCMC chain. The parameters are defined in Eq. (10); the colour scheme identifies the four clusters following Fig. 8.
In the text

	Fig. 8. Velocity dispersion function for the four clusters. We show the data bins and the Schechter function fit with their uncertainties, as defined in Eq. (10).
In the text

	Fig. 9. Measured velocity dispersion as a function of rest-frame SDSS r-band luminosity for our work.
In the text

	Fig. 10. Relative uncertainty on the measurement of the velocity dispersion. Comparison between the σw values presented in this work and the Hectospec observations of nine strong lensing clusters (Sohn et al. 2020, in cyan).
In the text

Fig. A.1. VLT/MUSE mean stacked normalised spectra of 14, 6, 17, and 27 member galaxies of A2744, AS1063, M0416, and M1149, respectively, with σ > 180 km s−1 and S/N > 15. The shaded regions are the standard deviation of each spectral pixel. The spectra are smoothed with a Gaussian kernel with a standard deviation of 3.75 Å. The red dashed vertical lines locate the absorption features detected in the spectra. For clearer visualisation, the normalised flux values from AS1063, M0416, and M1149 have been increased by 1, 2, and 3, respectively.

In the text

Fig. A.2. Velocity dispersion measurements for simulated spectra. Left panel: Relative offset of the recovered velocity dispersion with respect to the true value as a function of the S/N. We show the median relative uncertainty from Eq. (A.1) as a black line. Right panel: Absolute offset between the recovered and the true velocity dispersion values as a function of the input velocity dispersion. The points are colour-coded by the corresponding spectral S/N.

In the text