# A new step forward in realistic cluster lens mass modelling: Analysis of Hubble Frontier Field Cluster Abell S1063 from joint lensing, X-ray and galaxy kinematics data

Beauchesne Benjamin, Benjamin Clément, Pascale Hibon, Marceau Limousin, Dominique Eckert, Jean-Paul Kneib, Johan Richard, Priyamvada Natarajan, Mathilde Jauzac, Mireia Montes, et al.

► To cite this version:

Beauchesne Benjamin, Benjamin Clément, Pascale Hibon, Marceau Limousin, Dominique Eckert, et al. A new step forward in realistic cluster lens mass modelling: Analysis of Hubble Frontier Field Cluster Abell S1063 from joint lensing, X-ray and galaxy kinematics data. Monthly Notices of the Royal Astronomical Society, 2024, 527 (2), pp.3246-3275. 10.1093/mnras/stad3308. hal-03979897

HAL Id: hal-03979897

https://hal.science/hal-03979897

Submitted on 24 May 2024

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# A new step forward in realistic cluster lens mass modelling: Analysis of Hubble Frontier Field Cluster Abell S1063 from joint lensing, X-ray and galaxy kinematics data

Beauchesne Benjamin, Benjamin Clément, Pascale Hibon, Marceau Limousin, Dominique Eckert, Jean-Paul Kneib, Johan Richard, Priyamvada Natarajan, Mathilde Jauzac, Mireia Montes, et al.

► To cite this version:

Beauchesne Benjamin, Benjamin Clément, Pascale Hibon, Marceau Limousin, Dominique Eckert, et al.. A new step forward in realistic cluster lens mass modelling: Analysis of Hubble Frontier Field Cluster Abell S1063 from joint lensing, X-ray and galaxy kinematics data. Monthly Notices of the Royal Astronomical Society, 2023, 527 (2), pp.3246-3275. 10.1093/mnras/stad3308. hal-03979897

HAL Id: hal-03979897

https://hal.science/hal-03979897v1

Submitted on 24 May 2024

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
<figure><img src="image_2.png" /></figure>

# A new step forward in realistic cluster lens mass modelling: analysis of Hubble Frontier Field Cluster Abell S1063 from joint lensing, X-ray, and galaxy kinematics data

Benjamin Beauchesne<sup>1,2*</sup> Benjamin Clément,<sup>1</sup> Pascale Hibon<sup>1,2</sup> Marceau Limousin,<sup>3</sup>

Dominique Eckert,<sup>4</sup> Jean-Paul Kneib,<sup>1,3</sup> Johan Richard<sup>1,4</sup> Priyamvada Natarajan<sup>1,5,6</sup>

Mathilde Jauzac<sup>1,7,8,9,10</sup> Mireia Montes<sup>1,11,12</sup> Guillaume Mahler<sup>1,7,8</sup> Adélaïde Claeysens<sup>1,13</sup>

Alexandre Jeanneau,<sup>14</sup> Anton M. Koekemoer<sup>1,15</sup> David Lagattuta,<sup>7,8</sup> Amanda Pagul<sup>1,15,16</sup>
and Javier Sánchez<sup>1,15,17</sup>

*Affiliations are listed at the end of the paper*

Accepted 2023 October 23. Received 2023 June 19; in original form 2023 January 26

## ABSTRACT

We present a new method to simultaneously and self-consistently model the mass distribution of galaxy clusters that combines constraints from strong lensing features, X-ray emission, and galaxy kinematics measurements. We are able to successfully decompose clusters into their collisionless and collisional mass components thanks to the X-ray surface brightness, as well as use the dynamics of cluster members, to obtain more accurate masses exploiting the fundamental plane of elliptical galaxies. Knowledge from all observables is included through a consistent Bayesian approach in the likelihood or in physically motivated priors. We apply this method to the galaxy cluster Abell S1063 and produce a mass model that we publicly release with this paper. The resulting mass distribution presents different ellipticities for the intra-cluster gas and the other large-scale mass components as well as deviation from elliptical symmetry in the main halo. We assess the ability of our method to recover the masses of the different elements of the cluster using a mock cluster based on a simplified version of our Abell S1063 model. Thanks to the wealth of mutliwavelength information provided by the mass model and the detected X-ray emission, we also found evidence for an ongoing merger event with gas sloshing from a smaller infalling structure into the main cluster. In agreement with previous findings, the total mass, gas profile, and gas mass fraction are all consistent with small deviations from the hydrostatic equilibrium. This new mass model for Abell S1063 is publicly available, as the LENSTOOL extension used to construct it.

**Key words:** gravitational lensing: strong–galaxies: clusters: general–galaxies: clusters: individual: Abell S1063–X-rays: galaxies: clusters.

## 1 INTRODUCTION

Galaxy clusters are amongst the largest structures bound by gravity and, as such, represent a formidable probe at the cross-roads between cosmology and astrophysics. From their formation and growth to the evolution of the galaxies within them as well as the physics of plasma with the intra-cluster gas heated by the strong gravitational interaction ruling them, galaxy clusters provide a laboratory for astrophysics on all scales. It is possible to use them to constrain models of structure formation and evolution (see Allen, Evrard & Mantz 2011; Kravtsov & Borgani 2012, for review) or probe possible deviations from General relativity (Clowe et al. 2006; Lam et al. 2012; Cataneo & Rapetti 2018) as well as estimators of cosmological parameters (Jullo et al. 2010; Acebron et al. 2017). In the context of the Λ Cold Dark Matter (ΛCDM) paradigm, they represent one of the best tools to study dark matter (DM) properties (Clowe, Gonzalez & Markevitch 2004; Bradač et al. 2008; Natarajan et al.

2017) and divergences to the collisionless model (Harvey et al. 2015; Meneghetti et al. 2020; Limousin, Beauchesne & Jullo 2022). Such analyses rely on accurate measurements of the cluster mass distribution through different physical phenomena, amongst which gravitational lensing has an unquestionable place.

In the framework of General Relativity, gravitational lensing refers to the bending of light rays passing near massive objects such as galaxy clusters called lenses. The most prominent advantage of using this phenomenon is its ability to account for the total mass, making the DM indirectly visible (see review by Kneib & Natarajan 2011). It also requires fewer assumptions on the studied cluster, besides the thin lens approximation, compared to other methods based on hydrostatic equilibrium (Ettori et al. 2010) or galaxy kinematics (Mamon, Biviano & Boué 2013). The lensing effect can be divided into two regimes of different intensities. The weak regime refers to a statistical effect on the apparent shapes of background galaxies, which is mainly used to study the outskirts of galaxy clusters (Jauzac et al. 2012; Jullo et al. 2014). Closer to the cluster core, the lensing effect intensifies, and features in the strong regime appear, such as lensed galaxies distorted into giant arcs or an apparent presence of the same object

<footnote>*E-mail: benjamin.beauchesne@epfl.ch</footnote>
<figure><img src="image_2.png" /></figure>

# A new step forward in realistic cluster lens mass modelling: analysis of Hubble Frontier Field Cluster Abell S1063 from joint lensing, X-ray, and galaxy kinematics data

Benjamin Beauchesne<sup>1,2*</sup> Benjamin Clément,<sup>1</sup> Pascale Hibon<sup>1,2</sup> Marceau Limousin,<sup>3</sup>

Dominique Eckert,<sup>4</sup> Jean-Paul Kneib,<sup>1,3</sup> Johan Richard<sup>1,4</sup> Priyamvada Natarajan<sup>1,5,6</sup>

Mathilde Jauzac<sup>1,7,8,9,10</sup> Mireia Montes<sup>1,11,12</sup> Guillaume Mahler<sup>1,7,8</sup> Adélaïde Claeysens<sup>1,13</sup>

Alexandre Jeanneau,<sup>14</sup> Anton M. Koekemoer<sup>1,15</sup> David Lagattuta,<sup>7,8</sup> Amanda Pagul<sup>1,15,16</sup>
and Javier Sánchez<sup>1,15,17</sup>

*Affiliations are listed at the end of the paper*

Accepted 2023 October 23. Received 2023 June 19; in original form 2023 January 26

## ABSTRACT

We present a new method to simultaneously and self-consistently model the mass distribution of galaxy clusters that combines constraints from strong lensing features, X-ray emission, and galaxy kinematics measurements. We are able to successfully decompose clusters into their collisionless and collisional mass components thanks to the X-ray surface brightness, as well as use the dynamics of cluster members, to obtain more accurate masses exploiting the fundamental plane of elliptical galaxies. Knowledge from all observables is included through a consistent Bayesian approach in the likelihood or in physically motivated priors. We apply this method to the galaxy cluster Abell S1063 and produce a mass model that we publicly release with this paper. The resulting mass distribution presents different ellipticities for the intra-cluster gas and the other large-scale mass components as well as deviation from elliptical symmetry in the main halo. We assess the ability of our method to recover the masses of the different elements of the cluster using a mock cluster based on a simplified version of our Abell S1063 model. Thanks to the wealth of mutliwavelength information provided by the mass model and the detected X-ray emission, we also found evidence for an ongoing merger event with gas sloshing from a smaller infalling structure into the main cluster. In agreement with previous findings, the total mass, gas profile, and gas mass fraction are all consistent with small deviations from the hydrostatic equilibrium. This new mass model for Abell S1063 is publicly available, as the LENSTOOL extension used to construct it.

**Key words:** gravitational lensing: strong – galaxies: clusters: general – galaxies: clusters: individual: Abell S1063 – X-rays: galaxies: clusters.

## 1 INTRODUCTION

Galaxy clusters are amongst the largest structures bound by gravity and, as such, represent a formidable probe at the cross-roads between cosmology and astrophysics. From their formation and growth to the evolution of the galaxies within them as well as the physics of plasma with the intra-cluster gas heated by the strong gravitational interaction ruling them, galaxy clusters provide a laboratory for astrophysics on all scales. It is possible to use them to constrain models of structure formation and evolution (see Allen, Evrard & Mantz 2011; Kravtsov & Borgani 2012, for review) or probe possible deviations from General relativity (Clowe et al. 2006; Lam et al. 2012; Cataneo & Rapetti 2018) as well as estimators of cosmological parameters (Jullo et al. 2010; Acebron et al. 2017). In the context of the Λ Cold Dark Matter (ΛCDM) paradigm, they represent one of the best tools to study dark matter (DM) properties (Clowe, Gonzalez & Markevitch 2004; Bradač et al. 2008; Natarajan et al.

2017) and divergences to the collisionless model (Harvey et al. 2015; Meneghetti et al. 2020; Limousin, Beauchesne & Jullo 2022). Such analyses rely on accurate measurements of the cluster mass distribution through different physical phenomena, amongst which gravitational lensing has an unquestionable place.

In the framework of General Relativity, gravitational lensing refers to the bending of light rays passing near massive objects such as galaxy clusters called lenses. The most prominent advantage of using this phenomenon is its ability to account for the total mass, making the DM indirectly visible (see review by Kneib & Natarajan 2011). It also requires fewer assumptions on the studied cluster, besides the thin lens approximation, compared to other methods based on hydrostatic equilibrium (Ettori et al. 2010) or galaxy kinematics (Mamon, Biviano & Boué 2013). The lensing effect can be divided into two regimes of different intensities. The weak regime refers to a statistical effect on the apparent shapes of background galaxies, which is mainly used to study the outskirts of galaxy clusters (Jauzac et al. 2012; Jullo et al. 2014). Closer to the cluster core, the lensing effect intensifies, and features in the strong regime appear, such as lensed galaxies distorted into giant arcs or an apparent presence of the same object

<footnote>*E-mail: benjamin.beauchesne@epfl.ch</footnote>
multiple times in the case of a multiply imaged system. Different
modelling techniques have been developed to use the strong lensing
regime to recover the mass distribution, particularly with multiple
images positions and/or shapes (LENSTOOL, GLAFIC, LTM, and WSLAP
; Diego et al. 2007; Jullo et al. 2007; Zitrin et al. 2009; Oguri 2010).
All of them have different biases given their specific assumptions.
Still, these independent modelling methods obtain total mass profiles
of clusters with only a few per cent error (Jauzac et al. 2014;
Meneghetti et al. 2017) and yet offer one of the most robust probes of
the ΛCDM paradigm (Robertson et al. 2019; Meneghetti et al. 2020).

Though lensing analyses are now almost unavoidable to study galaxy clusters in detail, they are not yet providing exhaustive pictures of clusters as additional multiwavelength information that is available needs to be and can be incorporated. Indeed, even if the DM largely dominates the mass content, that profile can only be extracted accurately with further knowledge of the baryonic matter, i.e. the intra-cluster gas and the stars belonging to the cluster galaxies. The intra-cluster gas, typically expected to be in hydrostatic equilibrium in the cluster potential, is heated to high temperature and emits X-ray photons through a Bremsstrahlung emission that allows us to obtain its distribution. Thanks to the optical and spectroscopic observations needed to acquire lensing information, it is possible to recover the mass of cluster members. Therefore, together with the X-ray observations, we can almost fully disentangle the masses of the various cluster components – DM, hot gas and star content. Such decompositions are able to probe in detail the behaviour of the different components, allowing us to see any displacement between them and constraining potential self-interactions of DM particles, like in the case of merging clusters (Markevitch et al. 2004; Harvey et al. 2015) and test the assumption of hydrostatic equilibrium (Cerini, Cappelluti & Natarajan 2023). Combining different observables to estimate the mass also reduces systematic effects arising inherently from each method and allows us to increase the complexity of the mass reconstructions as more information can be included (Bonamigo et al. 2018; Granata et al. 2022).

Recent work by Bonamigo et al. (2017, 2018) incorporated an independent gas component from X-ray observables to the lensing modelling software LENSTOOL (Jullo et al. 2007). This component is added with fixed parameters, and the rest of the mass components are optimized with the systems of multiple images and this extra halo. In particular, this adds asymmetry to the reconstruction that would be difficult to obtain with strong lensing alone due to the sparsity of the data, and the low contribution of the gas to the overall mass budget (Bonamigo et al. 2018). However, this two-step approach does not allow us to include the full information provided by the X-ray emission as only the best-fitting model for the gas is explicitly included. This precludes the assessment of possible degeneracies between the different components and obtaining overall errors for the entire model through a single posterior distribution sampling. This method does not therefore provide a balance between the different pieces of information; the best-fitting model with gas-only is a priori not the best-fitting model according to the combination of the X-ray and lensing observables. This prior work has paved the way by developing the necessary tools and set the stage for merging the two steps with a single, simultaneous joint analysis that would address the previous weaknesses and allow comprehensive studies of galaxy clusters truly combining X-ray and lensing.

Continuing these previous improvements to LENSTOOL, the mod-
elling of cluster members has also been enhanced by explicitly
including spectroscopic information to calibrate the scaling relations
based on the Faber & Jackson law (Bergamini et al. 2019, hereafter
B19) and combining it with more complex scaling relations such

as the so-called fundamental plane of elliptical galaxies (Granata
et al. 2022, hereafter G22). The analysis of G22 enables the addition
of a scatter in the possible mass of each galaxy based on their
optical morphology and star kinematics. Notably, both of these
works include the gas component from Bonamigo et al. (2018)
and thus provide a more accurate decomposition of the cluster
components. In comparison to G22, B19 were able to include the
full information from the Faber & Jackson relation calibration in the
form of Bayesian priors because this scaling relation was the only one
currently implemented in LENSTOOL. Indeed, as for Bonamigo et al.
(2017), the fundamental plane measurement was only represented
by the best-fitting values, and the scaling between the light and mass
profiles was adjusted outside the parameter sampling. It was fixed
for a specific run of the sampler, and the final value was chosen by
selecting the one producing the best-fitting model amongst all runs. It
is then difficult to assess the proper influence of these assumptions on
other parts of the model as this amounts to slicing the likelihood on
some parameters instead of marginalizing on it as a fully consistent
Bayesian method would do.

Taking inspiration from these previous works, we present in this
paper a new method aimed at tackling the previously mentioned
drawbacks on the modelling of the gas and cluster galaxy members by
integrating both components in a single self-consistent homogeneous
process. Our method can be applied to all other strong lensing mass
reconstruction methods, but we have implemented it in the publicly
available software package LENSTOOL to make it broadly available to
the rest of the community. Thus, we start the paper with details of the
observational data of the cluster Abell S1063 in Section 2. We explain
in Section 3 the analysis of the X-ray and galaxy spectroscopic data
that are required as part of this new method. We, then, detail the
building blocks of the mass model in Section 4, followed by the
definition of our joint X-ray and lensing likelihood in Section 5. In
Sections 6 and 7, we analyse the mass distribution of Abell S1063
and the reconstruction of a mock cluster for comparison. We finally
discuss the deviation from the hydrostatic equilibrium; a possible
merging event and interrogate our modelling hypotheses in Section 8
and concluding with prospects for the application of this new method
in Section 9.

We adopt a flat ΛCDM cosmology with Ω_Λ = 0.7, Ω_m = 0.3, and H_0 = 70 km s⁻¹Mpc⁻¹, and we use magnitudes quoted in the AB system throughout this paper. Regarding the statistical treatment of all the analyses, we compute the uncertainties with the median-centred credible intervals (CI) based on the posterior distribution of the considered random variable. We choose the size of these intervals in analogy with the σ levels of a normal distribution; thus, uncertainties are expressed in intervals that contain 100 × erf (n/√2) per cent of the posterior with n an integer. To avoid heavy notations and to make clear that we are considering CI, we denote them "nσ"CI. We only use σ to refer to the actual standard deviation.

## 2 ABELL S1063: OBSERVATIONAL DATA

Abell S1063 (AS1063; also known as RXC J2248.7-4431 and SPT-CL J2248-4431) is a massive galaxy cluster at z = 0.3475 first identified by Abell, Corwin & Olowin (1989). It is a bright X-ray source with one of the highest X-ray temperatures measured, L<sub>X</sub> ≈ 2.5 × 10<sup>45</sup> erg s<sup>-1</sup> in the [0.1,2.4] keV band and T ≈ 12 keV, respectively (Rahaman et al. 2021). It also presents a strong Sunyaev-Zel'dovich (SZ) effect detection in the South Pole Telescope survey (Williamson et al. 2011), with an estimated SZ mass at the virial
multiple times in the case of a multiply imaged system. Different
modelling techniques have been developed to use the strong lensing
regime to recover the mass distribution, particularly with multiple
images positions and/or shapes (LENSTOOL, GLAFIC, LTM, and WSLAP
; Diego et al. 2007; Jullo et al. 2007; Zitrin et al. 2009; Oguri 2010).
All of them have different biases given their specific assumptions.
Still, these independent modelling methods obtain total mass profiles
of clusters with only a few per cent error (Jauzac et al. 2014;
Meneghetti et al. 2017) and yet offer one of the most robust probes of
the ΛCDM paradigm (Robertson et al. 2019; Meneghetti et al. 2020).

Though lensing analyses are now almost unavoidable to study galaxy clusters in detail, they are not yet providing exhaustive pictures of clusters as additional multiwavelength information that is available needs to be and can be incorporated. Indeed, even if the DM largely dominates the mass content, that profile can only be extracted accurately with further knowledge of the baryonic matter, i.e. the intra-cluster gas and the stars belonging to the cluster galaxies. The intra-cluster gas, typically expected to be in hydrostatic equilibrium in the cluster potential, is heated to high temperature and emits X-ray photons through a Bremsstrahlung emission that allows us to obtain its distribution. Thanks to the optical and spectroscopic observations needed to acquire lensing information, it is possible to recover the mass of cluster members. Therefore, together with the X-ray observations, we can almost fully disentangle the masses of the various cluster components – DM, hot gas and star content. Such decompositions are able to probe in detail the behaviour of the different components, allowing us to see any displacement between them and constraining potential self-interactions of DM particles, like in the case of merging clusters (Markevitch et al. 2004; Harvey et al. 2015) and test the assumption of hydrostatic equilibrium (Cerini, Cappelluti & Natarajan 2023). Combining different observables to estimate the mass also reduces systematic effects arising inherently from each method and allows us to increase the complexity of the mass reconstructions as more information can be included (Bonamigo et al. 2018; Granata et al. 2022).

Recent work by Bonamigo et al. (2017, 2018) incorporated an independent gas component from X-ray observables to the lensing modelling software LENSTOOL (Jullo et al. 2007). This component is added with fixed parameters, and the rest of the mass components are optimized with the systems of multiple images and this extra halo. In particular, this adds asymmetry to the reconstruction that would be difficult to obtain with strong lensing alone due to the sparsity of the data, and the low contribution of the gas to the overall mass budget (Bonamigo et al. 2018). However, this two-step approach does not allow us to include the full information provided by the X-ray emission as only the best-fitting model for the gas is explicitly included. This precludes the assessment of possible degeneracies between the different components and obtaining overall errors for the entire model through a single posterior distribution sampling. This method does not therefore provide a balance between the different pieces of information; the best-fitting model with gas-only is a priori not the best-fitting model according to the combination of the X-ray and lensing observables. This prior work has paved the way by developing the necessary tools and set the stage for merging the two steps with a single, simultaneous joint analysis that would address the previous weaknesses and allow comprehensive studies of galaxy clusters truly combining X-ray and lensing.

Continuing these previous improvements to LENSTOOL, the mod-
elling of cluster members has also been enhanced by explicitly
including spectroscopic information to calibrate the scaling relations
based on the Faber & Jackson law (Bergamini et al. 2019, hereafter
B19) and combining it with more complex scaling relations such

as the so-called fundamental plane of elliptical galaxies (Granata
et al. 2022, hereafter G22). The analysis of G22 enables the addition
of a scatter in the possible mass of each galaxy based on their
optical morphology and star kinematics. Notably, both of these
works include the gas component from Bonamigo et al. (2018)
and thus provide a more accurate decomposition of the cluster
components. In comparison to G22, B19 were able to include the
full information from the Faber & Jackson relation calibration in the
form of Bayesian priors because this scaling relation was the only one
currently implemented in LENSTOOL. Indeed, as for Bonamigo et al.
(2017), the fundamental plane measurement was only represented
by the best-fitting values, and the scaling between the light and mass
profiles was adjusted outside the parameter sampling. It was fixed
for a specific run of the sampler, and the final value was chosen by
selecting the one producing the best-fitting model amongst all runs. It
is then difficult to assess the proper influence of these assumptions on
other parts of the model as this amounts to slicing the likelihood on
some parameters instead of marginalizing on it as a fully consistent
Bayesian method would do.

Taking inspiration from these previous works, we present in this
paper a new method aimed at tackling the previously mentioned
drawbacks on the modelling of the gas and cluster galaxy members by
integrating both components in a single self-consistent homogeneous
process. Our method can be applied to all other strong lensing mass
reconstruction methods, but we have implemented it in the publicly
available software package LENSTOOL to make it broadly available to
the rest of the community. Thus, we start the paper with details of the
observational data of the cluster Abell S1063 in Section 2. We explain
in Section 3 the analysis of the X-ray and galaxy spectroscopic data
that are required as part of this new method. We, then, detail the
building blocks of the mass model in Section 4, followed by the
definition of our joint X-ray and lensing likelihood in Section 5. In
Sections 6 and 7, we analyse the mass distribution of Abell S1063
and the reconstruction of a mock cluster for comparison. We finally
discuss the deviation from the hydrostatic equilibrium; a possible
merging event and interrogate our modelling hypotheses in Section 8
and concluding with prospects for the application of this new method
in Section 9.

We adopt a flat ΛCDM cosmology with Ω_Λ = 0.7, Ω_m = 0.3, and H_0 = 70 km s⁻¹Mpc⁻¹, and we use magnitudes quoted in the AB system throughout this paper. Regarding the statistical treatment of all the analyses, we compute the uncertainties with the median-centred credible intervals (CI) based on the posterior distribution of the considered random variable. We choose the size of these intervals in analogy with the σ levels of a normal distribution; thus, uncertainties are expressed in intervals that contain 100 × erf (n/√2) per cent of the posterior with n an integer. To avoid heavy notations and to make clear that we are considering CI, we denote them "nσ"CI. We only use σ to refer to the actual standard deviation.

## 2 ABELL S1063: OBSERVATIONAL DATA

Abell S1063 (AS1063; also known as RXC J2248.7-4431 and SPT-CL J2248-4431) is a massive galaxy cluster at z = 0.3475 first identified by Abell, Corwin & Olowin (1989). It is a bright X-ray source with one of the highest X-ray temperatures measured, L<sub>X</sub> ≈ 2.5 × 10<sup>45</sup> erg s<sup>-1</sup> in the [0.1,2.4] keV band and T ≈ 12 keV, respectively (Rahaman et al. 2021). It also presents a strong Sunyaev-Zel'dovich (SZ) effect detection in the South Pole Telescope survey (Williamson et al. 2011), with an estimated SZ mass at the virial
<figure><img src="image_3.png" /><figcaption>**Figure 1.** Composite colour image of AS1063 made with HST filters F435W (blue), F606W (green), and F814W (red). Magenta and cyan dots highlight the positions of the modelled cluster members and multiple images, respectively. The two dashed white boxes show the position of the two main mass haloes as detailed in Section 4.1. Finally, dashed red contours and dashed green box present the X-ray surface brightness seen by the *Chandra* observatory and the footprint of the available VLT/MUSE observations, respectively.</figcaption></figure>

radius of $M_{200} = 2.9 \times 10^{15} M_{\odot}$. Fig. 1 shows a colour image of the cluster core made with the *Hubble Space Telescope* (HST) broadband filters F435W (blue), F606W (green), and F814W (red).

AS1063 is well suited for our analysis as it has been observed in multiple wavelengths – X-ray and optical – from space, as well as with an integral field spectrograph from the ground. Thus, we have all the necessary observational data to constrain the gas distribution and the galaxy scaling laws in addition to the lensing observables. Its morphology is also simpler in comparison to other cluster lenses with the same wealth of data. For starters, the cluster mass distribution comprises a single peak, i.e. it is unimodal. In addition, although it appears to be relaxed at first sight, previous analyses have favoured

a slightly perturbed cluster. As proposed by Gómez et al. (2012) and supported by Rahaman et al. (2021), de Oliveira, Jiménez-Teja & Dupke (2022), Xie et al. (2020), and Mercurio et al. (2021), AS1063 X-ray emission; intra-cluster light (ICL Montes 2022, for a review); the presence of a giant radio halo, and galaxy kinematics seem to show that the cluster had undergone a recent off-axis merger.

## 2.1 X-ray observations

As we are focused on a strong lensing analysis, we will model the core of AS1063. Thus, we choose observations from the *Chandra X-ray observatory* as it provides the best spatial resolution amongst
<figure><img src="image_3.png" /><figcaption>**Figure 1.** Composite colour image of AS1063 made with HST filters F435W (blue), F606W (green), and F814W (red). Magenta and cyan dots highlight the positions of the modelled cluster members and multiple images, respectively. The two dashed white boxes show the position of the two main mass haloes as detailed in Section 4.1. Finally, dashed red contours and dashed green box present the X-ray surface brightness seen by the *Chandra* observatory and the footprint of the available VLT/MUSE observations, respectively.</figcaption></figure>

radius of $M_{200} = 2.9 \times 10^{15}$ M☉. Fig. 1 shows a colour image of the cluster core made with the *Hubble Space Telescope* (HST) broadband filters F435W (blue), F606W (green), and F814W (red).

AS1063 is well suited for our analysis as it has been observed in multiple wavelengths – X-ray and optical – from space, as well as with an integral field spectrograph from the ground. Thus, we have all the necessary observational data to constrain the gas distribution and the galaxy scaling laws in addition to the lensing observables. Its morphology is also simpler in comparison to other cluster lenses with the same wealth of data. For starters, the cluster mass distribution comprises a single peak, i.e. it is unimodal. In addition, although it appears to be relaxed at first sight, previous analyses have favoured

a slightly perturbed cluster. As proposed by Gómez et al. (2012) and supported by Rahaman et al. (2021), de Oliveira, Jiménez-Teja & Dupke (2022), Xie et al. (2020), and Mercurio et al. (2021), AS1063 X-ray emission; intra-cluster light (ICL Montes 2022, for a review); the presence of a giant radio halo, and galaxy kinematics seem to show that the cluster had undergone a recent off-axis merger.

## 2.1 X-ray observations

As we are focused on a strong lensing analysis, we will model the core of AS1063. Thus, we choose observations from the *Chandra X-ray observatory* as it provides the best spatial resolution amongst
X-ray telescopes with a field of view (FoV) covering the entire area of interest. In particular, we use three combined archival pointings of this cluster with a total exposure time of 123 ks in very faint<fnref n="1" /> (VFAINT) mode only, with data from proposals IDs: 4966 (PI: Romer, 2004), 18 818 and 18 611 (PI: Kraft, 2016).

As the cluster covers a large part of the FoV of each observation, we also retrieved the *Blank-sky* backgrounds (Hickox & Markevitch 2006) associated with them to estimate the emission from the background sky. All of these data have been obtained with the Advanced CCD Imaging Spectrometer (Garmire et al. (2003) onboard *Chandra*.

We used a PYTHON wrapper of CIAO<fnref n="2" /> 4.13 (Fruscione et al. 2006) and CALDB 4.9.6 provided in the LENSTOOL public repository<fnref n="3" /> to reduce the data. We proceed first by removing point sources through a combination of visual inspection and the WAVEDETECT routine. Then, we correct for background flares by successively cleaning them in the [9.5,12] keV bands and in the whole energy range with the DEFLARE tool. We finally remove the area covered by the point sources from the background observations to avoid having negative number counts when subtracting.

These data are then used to create four images binned to four times the initial spatial resolution in the [4.0, 7.0] keV band (27531 observed and 12 512 background photon counts in the considered area) but also in the soft ([0.5, 1.2] keV; 29 256 photon counts), medium ([1.2, 2.0] keV; 50 977 photon counts), and hard ([2.0, 7.0] keV; 69 660 photon counts) science energy bands as defined in the *Chandra* Source Catalog (Evans et al. 2010). The former is used in the analysis of the X-ray observations to produce a tessellation of the FoV to map the thermodynamic properties of the intra-cluster gas as detailed in Section 3.1. As for the images in the soft, medium, and hard bands, they will be fed to LENSTOOL to perform the fit of the X-ray emitting gas, which is outlined in Section 4. The background level amongst these three bands (or the broad-band) is of 32 838 counts after re-scaling it to the observations.

## 2.2 Photometric data

AS1063 has been widely observed with HST through the Cluster Lensing and Supernova Survey [CLASH; Postman et al. (2012)]; the Hubble Frontier Fields [HFF; Lotz et al. (2017)] and the Beyond Ultra-deep Frontier Fields And Legacy Observations (BUFFALO; Steinhardt et al. (2020)). We use fully reprocessed mosaic images produced by the BUFFALO collaboration combining all available observations and including the HFF data set. For more details about the different available filters, FoV and depth of the HST observations, we refer the reader to Steinhardt et al. (2020).

We use the publicly available measurements on galaxies in the FoV made on the HFF data for this cluster from three different catalogues. In particular, we use the multiwavelength photometric catalogue produced by Pagul et al. (2023), combined with two catalogues of galaxy structural properties from Tortorelli et al. (2018) and Nedkova et al. (2021). We primarily use the morphology measurements from Tortorelli et al. (2018) that are made by fitting a single Sérsic light profile with GALAPAGOS and GALFIT [i.e. the single-band version; Peng et al. (2010)] in a two-step process. As per this process, profiles are first optimized on postage-stamp size images of each object. In the second step, all galaxies are fitted at the same time to the whole image, which improves the estimation of the local background, in particular

in the central region where the postage-stamp images are dominated by the light from the bright cluster members. The second structural catalogue presented in Nedkova et al. (2021) has been produced as part of the DeepSpace project and relies on a similar method to the one described above, but with only postage-stamp images included in the optimization with the multiwavelength fitting tools GALFITM and GALAPAGOS-2, developed as part of the MegaMorph project (Häußler et al. 2013; Vika et al. 2013). Hence, as it has a lower robustness to the intra-cluster light and the bright cluster members, we only use these data if they are not available from the first morphology catalogue to maximize our available knowledge at galaxy-scales.

## 2.3 Spectroscopic data

AS1063 was observed by the Multi-Unit Spectroscopic Explorer (MUSE) mounted on the Very Large Telescope (VLT; Bacon et al. 2010) with two pointings covering the south-west (SW) and north-east (NE) regions of the cluster (see the green dashed box in Fig. 1 for the combined FoV of the MUSE observations) with the following proposals IDs:

60.A-9345(A) (PI: Caputi & Grillo) and 095.A-0653(A) (PI: Caputi)

These data have already been reduced and analysed in Karman et al. (2015, 2017) and Caminha et al. (2016), wherein a redshift catalogue of the cluster core has been presented. However, we reduce and re-analyse the data set with the improved pipeline detailed in Richard et al. (2021), which has been specifically developed with a focus on cluster fields, to take into account more accurately the ICL and the edge of each integral field unit. Going forward, we use this new redshift catalogue and new MUSE data cube for the rest of this work. We note 10 per cent of the objects have their first redshift measurement from this re-analysis.

In addition, AS1063 is a target of the CLASH-VLT program (ESO ID: 186.A-0798; PI: P. Rosati) that obtained 200 h of observations on clusters from the CLASH sample with the VISIBLE Multi-Object Spectrograph previously installed on the VLT before its decommissioning. A catalogue containing almost 4000 redshifts has been produced by Mercurio et al. (2021) and released publicly. Thus, we use it to complete the MUSE catalogue mentioned previously on the HFF footprint.

# 3 PRE-MODELLING ANALYSIS

The approach used for this work is based on the one developed for the LENSTOOL software (Jullo et al. 2007) enhanced with the following refinements:

(i) Modelling of the X-ray emitting intra-cluster gas.

(ii) Inclusion of the kinematics of Galaxy-scale perturbers.

(iii) Addition of a perturbative surface of B-spline on the lensing potential (see Beauchesne et al. 2021, and Section 4.4 for details).

Hence, in addition to the usual modelling analysis presented in Section 4, we need measurements of the temperature and metallicity of the intra-cluster gas to fully define our plasma emission model. This process is detailed in Section 3.1. We also recover the kinematics of cluster members described in Section 3.2.

## 3.1 X-ray data analysis

We are interested in modelling the intra-cluster gas mass distribution using the X-ray surface brightness ($S_X$) maps that are presented in

<footnote n="1">¹https://cxc.cfa.harvard.edu/ciao/why/aciscleanvf.html</footnote>

<footnote n="2">²https://cxc.cfa.harvard.edu/ciao/</footnote>

<footnote n="3">³https://git-cral.univ-lyon1.fr/lenstool/lenstool</footnote>
X-ray telescopes with a field of view (FoV) covering the entire area of interest. In particular, we use three combined archival pointings of this cluster with a total exposure time of 123 ks in very faint<fnref n="1" /> (VFAINT) mode only, with data from proposals IDs: 4966 (PI: Romer, 2004), 18 818 and 18 611 (PI: Kraft, 2016).

As the cluster covers a large part of the FoV of each observation, we also retrieved the Blank-sky backgrounds (Hickox & Markevitch 2006) associated with them to estimate the emission from the background sky. All of these data have been obtained with the Advanced CCD Imaging Spectrometer (Garmire et al. 2003) onboard Chandra.

We used a PYTHON wrapper of CIAO<fnref n="2" /> 4.13 (Fruscione et al. 2006) and CALDB 4.9.6 provided in the LENSTOOL public repository<fnref n="3" /> to reduce the data. We proceed first by removing point sources through a combination of visual inspection and the WAVEDETECT routine. Then, we correct for background flares by successively cleaning them in the [9.5,12] keV bands and in the whole energy range with the DEFLARE tool. We finally remove the area covered by the point sources from the background observations to avoid having negative number counts when subtracting.

These data are then used to create four images binned to four times the initial spatial resolution in the [4.0, 7.0] keV band (27531 observed and 12 512 background photon counts in the considered area) but also in the soft ([0.5, 1.2] keV; 29 256 photon counts), medium ([1.2, 2.0] keV; 50 977 photon counts), and hard ([2.0, 7.0] keV; 69 660 photon counts) science energy bands as defined in the Chandra Source Catalog (Evans et al. 2010). The former is used in the analysis of the X-ray observations to produce a tessellation of the FoV to map the thermodynamic properties of the intra-cluster gas as detailed in Section 3.1. As for the images in the soft, medium, and hard bands, they will be fed to LENSTOOL to perform the fit of the X-ray emitting gas, which is outlined in Section 4. The background level amongst these three bands (or the broad-band) is of 32 838 counts after re-scaling it to the observations.

## 2.2 Photometric data

AS1063 has been widely observed with HST through the Cluster Lensing and Supernova Survey [CLASH; Postman et al. (2012)]; the Hubble Frontier Fields [HFF; Lotz et al. (2017)] and the Beyond Ultra-deep Frontier Fields And Legacy Observations (BUFFALO; Steinhardt et al. (2020)). We use fully reprocessed mosaic images produced by the BUFFALO collaboration combining all available observations and including the HFF data set. For more details about the different available filters, FoV and depth of the HST observations, we refer the reader to Steinhardt et al. (2020).

We use the publicly available measurements on galaxies in the FoV made on the HFF data for this cluster from three different catalogues. In particular, we use the multiwavelength photometric catalogue produced by Pagul et al. (2023), combined with two catalogues of galaxy structural properties from Tortorelli et al. (2018) and Nedkova et al. (2021). We primarily use the morphology measurements from Tortorelli et al. (2018) that are made by fitting a single Sérsic light profile with GALAPAGOS and GALFIT [i.e. the single-band version; Peng et al. (2010)] in a two-step process. As per this process, profiles are first optimized on postage-stamp size images of each object. In the second step, all galaxies are fitted at the same time to the whole image, which improves the estimation of the local background, in particular

in the central region where the postage-stamp images are dominated by the light from the bright cluster members. The second structural catalogue presented in Nedkova et al. (2021) has been produced as part of the DeepSpace project and relies on a similar method to the one described above, but with only postage-stamp images included in the optimization with the multiwavelength fitting tools GALFITM and GALAPAGOS-2, developed as part of the MegaMorph project (Häußler et al. 2013; Vika et al. 2013). Hence, as it has a lower robustness to the intra-cluster light and the bright cluster members, we only use these data if they are not available from the first morphology catalogue to maximize our available knowledge at galaxy-scales.

## 2.3 Spectroscopic data

AS1063 was observed by the Multi-Unit Spectroscopic Explorer (MUSE) mounted on the Very Large Telescope (VLT; Bacon et al. 2010) with two pointings covering the south-west (SW) and north-east (NE) regions of the cluster (see the green dashed box in Fig. 1 for the combined FoV of the MUSE observations) with the following proposals IDs:

60.A-9345(A) (PI: Caputi & Grillo) and 095.A-0653(A) (PI: Caputi)

These data have already been reduced and analysed in Karman et al. (2015, 2017) and Caminha et al. (2016), wherein a redshift catalogue of the cluster core has been presented. However, we reduce and re-analyse the data set with the improved pipeline detailed in Richard et al. (2021), which has been specifically developed with a focus on cluster fields, to take into account more accurately the ICL and the edge of each integral field unit. Going forward, we use this new redshift catalogue and new MUSE data cube for the rest of this work. We note 10 per cent of the objects have their first redshift measurement from this re-analysis.

In addition, AS1063 is a target of the CLASH-VLT program (ESO ID: 186.A-0798; PI: P. Rosati) that obtained 200 h of observations on clusters from the CLASH sample with the VISIBLE Multi-Object Spectrograph previously installed on the VLT before its decommissioning. A catalogue containing almost 4000 redshifts has been produced by Mercurio et al. (2021) and released publicly. Thus, we use it to complete the MUSE catalogue mentioned previously on the HFF footprint.

# 3 PRE-MODELLING ANALYSIS

The approach used for this work is based on the one developed for the LENSTOOL software (Jullo et al. 2007) enhanced with the following refinements:

(i) Modelling of the X-ray emitting intra-cluster gas.

(ii) Inclusion of the kinematics of Galaxy-scale perturbers.

(iii) Addition of a perturbative surface of B-spline on the lensing potential (see Beauchesne et al. 2021, and Section 4.4 for details).

Hence, in addition to the usual modelling analysis presented in Section 4, we need measurements of the temperature and metallicity of the intra-cluster gas to fully define our plasma emission model. This process is detailed in Section 3.1. We also recover the kinematics of cluster members described in Section 3.2.

## 3.1 X-ray data analysis

We are interested in modelling the intra-cluster gas mass distribution using the X-ray surface brightness ($S_X$) maps that are presented in

<footnote n="1">¹https://cxc.cfa.harvard.edu/ciao/why/aciscleanvf.html</footnote>

<footnote n="2">²https://cxc.cfa.harvard.edu/ciao/</footnote>

<footnote n="3">³https://git-cral.univ-lyon1.fr/lenstool/lenstool</footnote>
<figure><img src="image_4.png" /><figcaption>**Figure 2.** Flowchart presenting the process applied on X-ray data to recover the intra-cluster gas properties and the mass distribution. Each step is detailed in Section 3.1.</figcaption></figure>

Section 2.1, defined as:

$$S_X(x, y) = \frac{\Lambda(T, Z)}{4\pi(1+z)^4} \int_{\mathbb{R}} n_e(x, y, z) n_p(x, y, z) dz \quad (1)$$

Where $(x, y, z)$ are the observer coordinates with the z-axis aligned with the line of sight (LOS). Parameters $z$, $T$, $Z$, $n_e$, and $n_p$, are the cluster redshift, temperature, metallicity, electron density, and proton density, respectively. $\Lambda(T, Z)$ is the associated cooling function (Boehringer & Hensler 1989; Sutherland & Dopita 1993). To derive $S_X$ from the mass modelling (i.e. $n_e$ and $n_p$), we need to extract the thermodynamic properties of the gas. The method associated with this extraction is outlined in the following two sections, and the multiple steps are shown schematically in the flowchart in Fig. 2.

### 3.1.1 X-ray spectra fitting

To obtain the thermodynamic properties of the gas from its X-ray spectrum, we follow a similar approach to the one detailed in Rossetti et al. (2007). We make a Voronoi tessellation of the X-rays images of the cluster field with the VORBIN python packages (Cappellari & Copin 2003) from the photon counts associated with the cluster in the [4, 7] keV band. We need to rely on a tessellation to obtain enough photons on the X-ray spectra. We remove the signal associated from the background of the total signal, and we use the Gehrels (1986) approximation to obtain the uncertainty as the sum of the two associated variances. A signal-to-noise (S/N) threshold parametrizes this tessellation, for which we use an S/N of 10 as it seems to be the best compromise between the uncertainty and the mapping resolution in the different S/N that we tried. We obtained $9962_{-2747}^{+2489}$ observed

counts per bin on the [0.5, 12.5] keV band amongst all observations with this S/N as well as $16913_{-15387}^{+72685}$ for the blank-sky background before rescaling. We also reproduce the spectral fitting for different S/N to obtain more details on gas properties and insight into the induced bias on the gas distribution in Appendix E.

We assumed the metallicity to be constant across the whole cluster field, and as the temperature is high in AS1063 (<8 keV), this assumption should not bias the result on the mass fitting. In that range of temperature, the emission lines have only a limited contribution to the overall spectra. Other parameters of the gas emission model are mapped in the field based on the tessellation described above. All these quantities are fitted with the Astrophysical Plasma Emission Code (APEC) model,<fnref n="4" /> combined with a photoelectric-absorption (PHABS) model,<fnref n="5" /> to account for the foreground galactic gas. For the APEC model, we assume the abundance ratio provided by Asplund et al. (2009) and keep the other parameters free. The hydrogen column density ($n_H$) needed by the PHABS model is fixed to the value measured by HI4PI Collaboration (2016). Before proceeding to the two successive fitting procedures shown in the flowchart in Fig. 2, we empirically model the instrumental and sky background provided by the blank-sky observations with the B-spline functions. This background modelling approach is detailed in Appendix A.

Our fitting process is based on the fitting environment SHERPA (Burke et al. 2022) combined with the nested sampling method

<footnote n="4">⁴http://atomdb.org/</footnote>

<footnote n="5">⁵https://heasarc.gsfc.nasa.gov/docs/xanadu/xspec/manual/XSmodelPhabs.html</footnote>
<figure><img src="image_4.png" /><figcaption>**Figure 2.** Flowchart presenting the process applied on X-ray data to recover the intra-cluster gas properties and the mass distribution. Each step is detailed in Section 3.1.</figcaption></figure>

Section 2.1, defined as:

$$S_X(x, y) = \frac{\Lambda(T, Z)}{4\pi(1+z)^4} \int_{\mathbb{R}} n_e(x, y, z) n_p(x, y, z) dz \quad (1)$$

Where $(x, y, z)$ are the observer coordinates with the z-axis aligned with the line of sight (LOS). Parameters $z$, $T$, $Z$, $n_e$, and $n_p$, are the cluster redshift, temperature, metallicity, electron density, and proton density, respectively. $\Lambda(T, Z)$ is the associated cooling function (Boehringer & Hensler 1989; Sutherland & Dopita 1993). To derive $S_X$ from the mass modelling (i.e. $n_e$ and $n_p$), we need to extract the thermodynamic properties of the gas. The method associated with this extraction is outlined in the following two sections, and the multiple steps are shown schematically in the flowchart in Fig. 2.

### 3.1.1 X-ray spectra fitting

To obtain the thermodynamic properties of the gas from its X-ray spectrum, we follow a similar approach to the one detailed in Rossetti et al. (2007). We make a Voronoi tessellation of the X-rays images of the cluster field with the VORBIN python packages (Cappellari & Copin 2003) from the photon counts associated with the cluster in the [4, 7] keV band. We need to rely on a tessellation to obtain enough photons on the X-ray spectra. We remove the signal associated from the background of the total signal, and we use the Gehrels (1986) approximation to obtain the uncertainty as the sum of the two associated variances. A signal-to-noise (S/N) threshold parametrizes this tessellation, for which we use an S/N of 10 as it seems to be the best compromise between the uncertainty and the mapping resolution in the different S/N that we tried. We obtained $9962_{-2747}^{+2489}$ observed

counts per bin on the [0.5, 12.5] keV band amongst all observations with this S/N as well as $16913_{-15387}^{+72685}$ for the blank-sky background before rescaling. We also reproduce the spectral fitting for different S/N to obtain more details on gas properties and insight into the induced bias on the gas distribution in Appendix E.

We assumed the metallicity to be constant across the whole cluster field, and as the temperature is high in AS1063 (<8 keV), this assumption should not bias the result on the mass fitting. In that range of temperature, the emission lines have only a limited contribution to the overall spectra. Other parameters of the gas emission model are mapped in the field based on the tessellation described above. All these quantities are fitted with the Astrophysical Plasma Emission Code (APEC) model,<fnref n="4" /> combined with a photoelectric-absorption (PHABS) model,<fnref n="5" /> to account for the foreground galactic gas. For the APEC model, we assume the abundance ratio provided by Asplund et al. (2009) and keep the other parameters free. The hydrogen column density ($n_H$) needed by the PHABS model is fixed to the value measured by HI4PI Collaboration (2016). Before proceeding to the two successive fitting procedures shown in the flowchart in Fig. 2, we empirically model the instrumental and sky background provided by the blank-sky observations with the B-spline functions. This background modelling approach is detailed in Appendix A.

Our fitting process is based on the fitting environment SHERPA (Burke et al. 2022) combined with the nested sampling method

<footnote n="4">⁴http://atomdb.org/</footnote>

<footnote n="5">⁵https://heasarc.gsfc.nasa.gov/docs/xanadu/xspec/manual/XSmodelPhabs.html</footnote>
PYMULTINEST (Buchner et al. 2014; Feroz et al. 2019). We adopt the following Poisson log-likelihood, $\mathcal{L}_{X-ray}$, which is similar to XSPEC CSTAT<fnref n="6" /> without the Sterling approximation:

$$ \mathcal{L}_{X-ray} = \sum_{i} D_{i} \log(M_{i}) - M_{i} - \log(D_{i}!) \quad (2) $$

where $D_i$ and $M_i$ represent the observed and model number of counts in each spectral bin, respectively. The term, $\log(D_i!)$, is computed before and kept in memory. We benefit from the following relation to reduce the computational cost:

$$ \log(D_i!) = \log(D_{i-1}!) + \log(D_i) \quad (3) $$

We first fit the whole cluster field X-ray spectrum to extract the metallicity of the gas. We then fix it to the best-fitting value and perform the same spectral fit on each individual bin. All spectra are taken in the [0.5, 12.0] keV range, the high energy part being dedicated to the normalization of the background model. Hence, the sky background is assumed to follow the instrumental one, as this normalization is mainly constrained on energy levels where the effective area of the telescope goes to zero. We also check the suitability of our assumed $n_H$ values by allowing it free to optimize, which leads to $n_H = (2.285_{-0.681}^{+0.496}) \times 10^{20} \text{ cm}^{-2}$ for the "1$\sigma$"CI. The measured value of $n_H = 1.293 \times 10^{20} \text{ cm}^{-2}$ from HI4PI Collaboration (2016) is included in the "2$\sigma$"CI. Hence, it is a fair assumption to use the measured values for the rest of our analysis.

### 3.1.2 Gas properties

The fitting approach above provides us with the temperature, $T$, the metallicity, $Z$, and the APEC norm, $N$, that we use to obtain the pseudo-pressure, $P$, and the pseudo entropy, $E$, with the following formula [e.g. similar to Rossetti et al. (2007) with different units]:

$$ P = T \sqrt{\frac{10^{-14} \pi \text{nhc}}{4.6656 \times 10^8 (1+z)^2} \frac{N}{\text{Area}_{\text{bin}}} \quad (4) $$

$$ E = T \left( \frac{10^{-14} \pi \text{nhc}}{4.6656 \times 10^8 (1+z)^2} \frac{N}{\text{Area}_{\text{bin}}} \right)^{-1/3} \quad (5) $$

Where $\text{Area}_{\text{bin}}$ is the area covered by each bin in arcsec and nhc is the conversion factor from $n_e$ to $n_p$. Regarding $\Lambda(T, Z)$, it is obtained as a single emission measure for the associated plasma model with a *Sherpa* routine. We compute it for the soft, medium, and hard energy bands that we multiply with the exposure maps associated before summing them. The obtained maps are then fed in to *LENSTOOL* with the total counts amongst all bands and the *blank-sky* backgrounds as detailed in Section 5.1.

## 3.2 Cluster member kinematics

In this section, we detail how we measure the stellar LOS velocity dispersion (hereafter LOSVD) of cluster members from the MUSE data cube. Our approach is similar to the one outlined in B19, with some modifications needed for our modelling choices as presented in Section 4.

### 3.2.1 LOSVD measurements

We extract galaxy spectra from the MUSE data cube with an elliptical aperture, $R_{ap} = R_e$, where $R_e$ is the half-light radius of

the galaxy. To create this elliptical mask, we measure the semimajor axis, A_WORLD, semiminor axis, B_WORLD, and position angles, THETA_J2000, using SExTRACTOR<fnref n="7" /> (Bertin & Arnouts 1996) in the HFF images which have the same astrometry as our data cube. We use $R_e$ from the structural catalogues presented in Section 2.3, using FLUX_RADIUS (i.e. $R_e$ proxy estimation) when the previous ones are not available. Thus, we obtained spectra for all the 107 spectroscopically confirmed cluster members in the MUSE FoV. As the BUFFALO program mainly expands the HFF observations spatially, these measurements should not be affected by the use of the later program images. The extraction of the spectra is performed with the python package MPDAF (Bacon et al. 2016) developed for the analysis of MUSE data. We use the optimal extraction algorithm for CCD spectroscopy implemented (Horne 1986). Then, we limit our wavelengths of analysis to 4850–7160 nm in the observer frame, which corresponds to the range used by B19.

Measurements of the LOSVD are performed with the python package PPXF (Cappellari 2017) modified to use the nested sampling engine pyMultiNest (Buchner et al. 2014; Feroz et al. 2019) as the non-linear optimizer method. LOSVD parameters are obtained from a cross-correlation between the galaxy and a set of star spectra templates. In the original package, parameters are optimized through the minimization of a $\chi^2$ statistic that has been modified here to maximize the associated Gaussian log-likelihood. In particular, we sample the LOS velocity with respect to the cluster redshift, $V$, the velocity dispersion, $\sigma_e$, and the two first Hermite moments, ($h_3, h_4$), from equation 13 of Cappellari (2017). We also use a Legendre polynomial of degree 3 to modify in a multiplicative way the continuum, as showed in equation 11 from the same reference. Thus, we fit seven parameters for each spectrum, where we use Gaussian priors for $V$, $h_3$, and $h_4$ with the following law $\mathcal{N}(0, 250)$, $\mathcal{N}(0, 0.05)$, and $\mathcal{N}(0, 0.05)$, respectively. The Legendre polynomial coefficient and $\sigma_e$, have uniform priors with the [-0.3, 0.3] interval (i.e. default bounds) for the first, and in the range [0,1000] for the second. For the star templates, we use the Indo-US Library of Coudé Feed Stellar Spectra (Valdes et al. 2004), which have a full-width half maximum (FWHM) of 1.35 Å, and a pixel-scale of 0.44 Å pixel$^{-1}$. To automatically correct for the emission-line contamination, we use the CLEAN parameter of PPXF in a two-step approach. We first run PPXF with a geometrical optimizer that is, in our case, the Trust Region Reflective algorithm implemented in SCIPY (Virtanen et al. 2020) with PPXF cleaning mode activated. This mode rejects outliers recursively until it converges to a stable number of masked data points. Then, we use the previous mask to perform the nested sampling. This is a compelling computing optimization, as running the nested sampling iteratively to improve the masking is time expensive.

We restrict our sample to galaxies with $S/N > 3$ for their spectra, and a redshift estimation with more than 80 per cent reliability [i.e. $z_{conf} \ge 2$ for MUSE redshift; see Bacon et al. (2015) for details]. We analyse the bias of these measurements in Appendix B, and we also apply the correction defined in equation (B1).

### 3.2.2 The fundamental plane of elliptical galaxies and the Faber & Jackson relation

Thanks to high-resolution images provided by HST within the HFF program, the structural properties of cluster members have been measured. With the addition of the velocity dispersion, $\sigma_e$,

<footnote n="6">⁶https://cxc.harvard.edu/sherpa4.13/ahelp/cstat.html</footnote>

<footnote n="7">⁷https://sextractor.readthedocs.io/en/latest/Position.html</footnote>
PYMULTINEST (Buchner et al. 2014; Feroz et al. 2019). We adopt the following Poisson log-likelihood, $\mathcal{L}_{X-ray}$, which is similar to XSPEC CSTAT<fnref n="6" /> without the Sterling approximation:

$$ \mathcal{L}_{X-ray} = \sum_{i} D_{i} \log(M_{i}) - M_{i} - \log(D_{i}!) \quad (2) $$

where $D_i$ and $M_i$ represent the observed and model number of counts in each spectral bin, respectively. The term, $\log(D_i!)$, is computed before and kept in memory. We benefit from the following relation to reduce the computational cost:

$$ \log(D_i!) = \log(D_{i-1}!) + \log(D_i) \quad (3) $$

We first fit the whole cluster field X-ray spectrum to extract the metallicity of the gas. We then fix it to the best-fitting value and perform the same spectral fit on each individual bin. All spectra are taken in the [0.5, 12.0] keV range, the high energy part being dedicated to the normalization of the background model. Hence, the sky background is assumed to follow the instrumental one, as this normalization is mainly constrained on energy levels where the effective area of the telescope goes to zero. We also check the suitability of our assumed $n_H$ values by allowing it free to optimize, which leads to $n_H = (2.285_{-0.681}^{+0.496}) \times 10^{20} \text{ cm}^{-2}$ for the "1$\sigma$"CI. The measured value of $n_H = 1.293 \times 10^{20} \text{ cm}^{-2}$ from HI4PI Collaboration (2016) is included in the "2$\sigma$"CI. Hence, it is a fair assumption to use the measured values for the rest of our analysis.

### 3.1.2 Gas properties

The fitting approach above provides us with the temperature, T, the metallicity, Z, and the APEC norm, N, that we use to obtain the pseudo-pressure, P, and the pseudo entropy, E, with the following formula [e.g. similar to Rossetti et al. (2007) with different units]:

$$ P = T \sqrt{\frac{10^{-14} \pi \text{nhc}}{4.6656 \times 10^8 (1+z)^2} \frac{N}{\text{Area}_{\text{bin}}} \quad (4) $$

$$ E = T \left( \frac{10^{-14} \pi \text{nhc}}{4.6656 \times 10^8 (1+z)^2} \frac{N}{\text{Area}_{\text{bin}}} \right)^{-1/3} \quad (5) $$

Where $\text{Area}_{\text{bin}}$ is the area covered by each bin in arcsec and nhc is the conversion factor from $n_e$ to $n_p$. Regarding $\Lambda(T, Z)$, it is obtained as a single emission measure for the associated plasma model with a *Sherpa* routine. We compute it for the soft, medium, and hard energy bands that we multiply with the exposure maps associated before summing them. The obtained maps are then fed in to *LENSTOOL* with the total counts amongst all bands and the *blank-sky* backgrounds as detailed in Section 5.1.

## 3.2 Cluster member kinematics

In this section, we detail how we measure the stellar LOS velocity dispersion (hereafter LOSVD) of cluster members from the MUSE data cube. Our approach is similar to the one outlined in B19, with some modifications needed for our modelling choices as presented in Section 4.

### 3.2.1 LOSVD measurements

We extract galaxy spectra from the MUSE data cube with an elliptical aperture, $R_{ap} = R_e$, where $R_e$ is the half-light radius of

the galaxy. To create this elliptical mask, we measure the semimajor axis, A_WORLD, semiminor axis, B_WORLD, and position angles, THETA_J2000, using SExTRACTOR<fnref n="7" /> (Bertin & Arnouts 1996) in the HFF images which have the same astrometry as our data cube. We use $R_e$ from the structural catalogues presented in Section 2.3, using FLUX_RADIUS (i.e. $R_e$ proxy estimation) when the previous ones are not available. Thus, we obtained spectra for all the 107 spectroscopically confirmed cluster members in the MUSE FoV. As the BUFFALO program mainly expands the HFF observations spatially, these measurements should not be affected by the use of the later program images. The extraction of the spectra is performed with the python package MPDAF (Bacon et al. 2016) developed for the analysis of MUSE data. We use the optimal extraction algorithm for CCD spectroscopy implemented (Horne 1986). Then, we limit our wavelengths of analysis to 4850–7160 nm in the observer frame, which corresponds to the range used by B19.

Measurements of the LOSVD are performed with the python package PPXF (Cappellari 2017) modified to use the nested sampling engine pyMultiNest (Buchner et al. 2014; Feroz et al. 2019) as the non-linear optimizer method. LOSVD parameters are obtained from a cross-correlation between the galaxy and a set of star spectra templates. In the original package, parameters are optimized through the minimization of a $\chi^2$ statistic that has been modified here to maximize the associated Gaussian log-likelihood. In particular, we sample the LOS velocity with respect to the cluster redshift, V, the velocity dispersion, $\sigma_e$, and the two first Hermite moments, ($h_3, h_4$), from equation 13 of Cappellari (2017). We also use a Legendre polynomial of degree 3 to modify in a multiplicative way the continuum, as showed in equation 11 from the same reference. Thus, we fit seven parameters for each spectrum, where we use Gaussian priors for V, $h_3$, and $h_4$ with the following law $N(0, 250)$, $N'(0, 0.05)$, and $N''(0, 0.05)$, respectively. The Legendre polynomial coefficient and $\sigma_e$, have uniform priors with the [-0.3, 0.3] interval (i.e. default bounds) for the first, and in the range [0,1000] for the second. For the star templates, we use the Indo-US Library of Coudé Feed Stellar Spectra (Valdes et al. 2004), which have a full-width half maximum (FWHM) of 1.35 Å, and a pixel-scale of 0.44 Å pixel$^{-1}$. To automatically correct for the emission-line contamination, we use the CLEAN parameter of PPXF in a two-step approach. We first run PPXF with a geometrical optimizer that is, in our case, the Trust Region Reflective algorithm implemented in SCIPY (Virtanen et al. 2020) with PPXF cleaning mode activated. This mode rejects outliers recursively until it converges to a stable number of masked data points. Then, we use the previous mask to perform the nested sampling. This is a compelling computing optimization, as running the nested sampling iteratively to improve the masking is time expensive.

We restrict our sample to galaxies with $S/N > 3$ for their spectra, and a redshift estimation with more than 80 per cent reliability [i.e. $z_{conf} \ge 2$ for MUSE redshift; see Bacon et al. (2015) for details]. We analyse the bias of these measurements in Appendix B, and we also apply the correction defined in equation (B1).

### 3.2.2 The fundamental plane of elliptical galaxies and the Faber & Jackson relation

Thanks to high-resolution images provided by HST within the HFF program, the structural properties of cluster members have been measured. With the addition of the velocity dispersion, $\sigma_e$,

<footnote n="6">⁶https://cxc.harvard.edu/sherpa4.13/ahelp/cstat.html</footnote>

<footnote n="7">⁷https://sextractor.readthedocs.io/en/latest/Position.html</footnote>
<figure><img src="image_5.png" /></figure>

<figure><img src="image_6.png" /><figcaption>**Figure 3.** Left-hand panel: Graph of the LOSVD measurements as a function of the magnitude in the HST F814W filter of the associated cluster member. LOSVD errors are the standard deviation of the distribution obtained by the nested sampling run, and the magnitude ones are provided by the Pagul et al. (2023) catalogue. The plain black line and the light orange area represent the best-fitting relations as well as the "1σ"CI of the whole sampled relations. Right-panel: Plot of the fundamental plane against the logarithm of the half-light radius. Errors represent the standard deviation and have been taken from the morphological catalogues detailed in Section 2.2 or computed from a sampling run. The black plain line shows the best fit of the fundamental plane, and the light orange area represents the position of the green dots amongst the "1σ"CI of the sample of fundamental plane parameters. For both panels, green dots highlight objects considered in the final fit when the blue ones have been rejected from the calibration process.</figcaption></figure>

measurement, we can calibrate the fundamental plane of elliptical galaxies, as defined in Hyde & Bernardi (2009) as follows:

$$ \log_{10}(R_e) = a + b\mu_e + c \log_{10}(\sigma_e) \quad (6) $$

Where $\mu_e$ is the averaged surface brightness inside $R_e$, in mag arcsec$^{-2}$, and a, b, and c are the parametrization of the plane. Regarding the units, we adopt arcsec and km s$^{-1}$ for $R_e$ and $\sigma_e$, respectively. We use $R_e$ and $\mu_e$ provided by the structural catalogues that we previously obtained, as outlined in Section 2.2.

From the expression of the fundamental plane, we can derive the Faber & Jackson relation (Faber & Jackson 1976), which has the following form:

$$ L \propto \sigma_e^\alpha \quad (7) $$

Where L is the galaxy luminosity and $\alpha$ is associated with the coefficient c. Indeed, we have:

$$ L \propto R_e^2 10^{-\mu_e/2.5} \quad (8) $$

$$ \Rightarrow L \propto 10^{-\mu_e/2.5} 10^{2(a+b\mu_e)} \sigma_e^{2c} \quad (9) $$

$$ \Rightarrow L \propto \sigma_e^{2c}, \text{ assuming } \mu_e = \text{Cst.} \quad (10) $$

Following B19, we calibrate the Faber & Jackson relation through the associated scaling relation usually defined in LENSTOOL (Natarajan & Kneib 1997; Richard et al. 2010), which is:

$$ \sigma_e = \sigma_e^* \left( \frac{L}{L^*} \right)^{1/2c} \quad (11) $$

$$ \Rightarrow \log(\sigma_e) = \log(\sigma_e^*) + \frac{1}{2c} \log \left( \frac{L}{L^*} \right) \quad (12) $$

Where $L^*$ and $\sigma_e^*$ are the luminosity and the velocity dispersion of a galaxy that is typical of the galaxy population at that redshift. $L^*$ is usually chosen where the elliptical galaxy luminosity function cuts off (Schechter 1976). Thus, the calibration of this relation is identical to fitting a line in the log space. Joint optimization of both relations can be done with four parameters, three for the fundamental plane and one for the Faber & Jackson relation.

<table><caption>**Table 1.** Parameters of the Faber & Jackson scaling relations and the fundamental plane of elliptical galaxies used or sampled in this work.</caption><thead><tr><th>Parameters</th><th>Mean</th><th>σ</th></tr></thead><tbody><tr><td>a</td><td>-0.543</td><td>0.029</td></tr><tr><td>b</td><td>0.335</td><td>0.028</td></tr><tr><td>c</td><td>1.365</td><td>0.122</td></tr><tr><td>σ<sub>e</sub><sup>*</sup></td><td>195.4</td><td>6.8</td></tr><tr><td>b<sub>FB</sub></td><td>[0.142634]</td><td>—</td></tr><tr><td>b<sub>FP</sub></td><td>[0.159140]</td><td>—</td></tr><tr><td>log(σ<sub>e</sub>)</td><td>[2.08949]</td><td>—</td></tr><tr><td>M<sub>AB</sub></td><td>[19.62834]</td><td>—</td></tr></tbody></table>

To perform the calibration of the plane and the Faber & Jackson relation, we use a combined approach inspired by the method outlined in Cappellari et al. (2013). We use the same $\chi^2$ statistics that are defined to fit a plane or a line, taking into account the scatter on all data sets, but we sum the two associated likelihoods to perform the joint optimization with PYMULTINEST. We assume the statistical independence of both relations in the combined likelihood, which is a priori not true, but allows us to derive a consistent slope between both relations. As shown in G22 and B19, the separate fit can lead to different slopes in the mass versus $\sigma$ plane, and thus inconsistent modelling of cluster members depending on which relation is used. Regarding the outliers, we use an iterative approach applied on the fundamental plane only, where all data points that are at more than 2.5$\sigma$ than the associated best-fitting relations are rejected. This rejection is performed until there are no more data points to exclude. However, this method results in an increased intrinsic scatter on this second relation. Regarding this scatter, it is estimated at each step of the rejection by running the nested sampler a first time with scatter fixed to zero and choosing it in a way that the associated reduced $\chi^2_{best fit} = 1$. The actual sampling is then done with this value.

Fig. 3 presents the fundamental plane (left-hand panel) as well as the Faber & Jackson law (right-panel panel) sampling results, and the parameter posterior distribution statistics are shown in Table 1. The left plot shows that only two sources have been rejected from the fundamental plane (i.e. red dots), signalling that our selection
<figure><img src="image_5.png" /></figure>

<figure><img src="image_6.png" /><figcaption>**Figure 3.** Left-hand panel: Graph of the LOSVD measurements as a function of the magnitude in the HST F814W filter of the associated cluster member. LOSVD errors are the standard deviation of the distribution obtained by the nested sampling run, and the magnitude ones are provided by the Pagul et al. (2023) catalogue. The plain black line and the light orange area represent the best-fitting relations as well as the "1σ"CI of the whole sampled relations. Right-panel: Plot of the fundamental plane against the logarithm of the half-light radius. Errors represent the standard deviation and have been taken from the morphological catalogues detailed in Section 2.2 or computed from a sampling run. The black plain line shows the best fit of the fundamental plane, and the light orange area represents the position of the green dots amongst the "1σ"CI of the sample of fundamental plane parameters. For both panels, green dots highlight objects considered in the final fit when the blue ones have been rejected from the calibration process.</figcaption></figure>

measurement, we can calibrate the fundamental plane of elliptical galaxies, as defined in Hyde & Bernardi (2009) as follows:

$$ \log_{10}(R_e) = a + b\mu_e + c \log_{10}(\sigma_e) \quad (6) $$

Where $\mu_e$ is the averaged surface brightness inside $R_e$, in mag arcsec$^{-2}$, and a, b, and c are the parametrization of the plane. Regarding the units, we adopt arcsec and km s$^{-1}$ for $R_e$ and $\sigma_e$, respectively. We use $R_e$ and $\mu_e$ provided by the structural catalogues that we previously obtained, as outlined in Section 2.2.

From the expression of the fundamental plane, we can derive the Faber & Jackson relation (Faber & Jackson 1976), which has the following form:

$$ L \propto \sigma_e^\alpha \quad (7) $$

Where L is the galaxy luminosity and $\alpha$ is associated with the coefficient c. Indeed, we have:

$$ L \propto R_e^2 10^{-\mu_e/2.5} \quad (8) $$

$$ \Rightarrow L \propto 10^{-\mu_e/2.5} 10^{2(a+b\mu_e)} \sigma_e^{2c} \quad (9) $$

$$ \Rightarrow L \propto \sigma_e^{2c}, \text{ assuming } \mu_e = \text{Cst.} \quad (10) $$

Following B19, we calibrate the Faber & Jackson relation through the associated scaling relation usually defined in LENSTOOL (Natarajan & Kneib 1997; Richard et al. 2010), which is:

$$ \sigma_e = \sigma_e^* \left( \frac{L}{L^*} \right)^{1/2c} \quad (11) $$

$$ \Rightarrow \log(\sigma_e) = \log(\sigma_e^*) + \frac{1}{2c} \log \left( \frac{L}{L^*} \right) \quad (12) $$

Where $L^*$ and $\sigma_e^*$ are the luminosity and the velocity dispersion of a galaxy that is typical of the galaxy population at that redshift. $L^*$ is usually chosen where the elliptical galaxy luminosity function cuts off (Schechter 1976). Thus, the calibration of this relation is identical to fitting a line in the log space. Joint optimization of both relations can be done with four parameters, three for the fundamental plane and one for the Faber & Jackson relation.

<table><caption>**Table 1.** Parameters of the Faber & Jackson scaling relations and the fundamental plane of elliptical galaxies used or sampled in this work.</caption><thead><tr><th>Parameters</th><th>Mean</th><th>σ</th></tr></thead><tbody><tr><td>a</td><td>-0.543</td><td>0.029</td></tr><tr><td>b</td><td>0.335</td><td>0.028</td></tr><tr><td>c</td><td>1.365</td><td>0.122</td></tr><tr><td>σ<sup>*</sup><sub>e</sub></td><td>195.4</td><td>6.8</td></tr><tr><td>b<sub>FB</sub></td><td>[0.142634]</td><td>—</td></tr><tr><td>b<sub>FP</sub></td><td>[0.159140]</td><td>—</td></tr><tr><td>log(σ<sub>e</sub>)</td><td>[2.08949]</td><td>—</td></tr><tr><td>M<sub>AB</sub></td><td>[19.62834]</td><td>—</td></tr></tbody></table>

To perform the calibration of the plane and the Faber & Jackson relation, we use a combined approach inspired by the method outlined in Cappellari et al. (2013). We use the same $\chi^2$ statistics that are defined to fit a plane or a line, taking into account the scatter on all data sets, but we sum the two associated likelihoods to perform the joint optimization with PYMULTINEST. We assume the statistical independence of both relations in the combined likelihood, which is a priori not true, but allows us to derive a consistent slope between both relations. As shown in G22 and B19, the separate fit can lead to different slopes in the mass versus $\sigma$ plane, and thus inconsistent modelling of cluster members depending on which relation is used. Regarding the outliers, we use an iterative approach applied on the fundamental plane only, where all data points that are at more than 2.5$\sigma$ than the associated best-fitting relations are rejected. This rejection is performed until there are no more data points to exclude. However, this method results in an increased intrinsic scatter on this second relation. Regarding this scatter, it is estimated at each step of the rejection by running the nested sampler a first time with scatter fixed to zero and choosing it in a way that the associated reduced $\chi^2_{bestfit} = 1$. The actual sampling is then done with this value.

Fig. 3 presents the fundamental plane (left-hand panel) as well as the Faber & Jackson law (right-panel panel) sampling results, and the parameter posterior distribution statistics are shown in Table 1. The left plot shows that only two sources have been rejected from the fundamental plane (i.e. red dots), signalling that our selection
of cluster members lies well into that plane, unlike the Faber & Jackson relations where some galaxies show a larger scatter from the relation. In comparison to G22 and B19, we have a higher estimation of c at 3.1σ, and a similar one at less than 1σ, respectively. However, we used different apertures than B19 and G22 to extract cluster member spectra that can partially explain the difference with G22 as well as in the joint optimization. Both plane estimations are in agreement with blank field results at less than 2σ in the R-band (Hyde & Bernardi 2009). Regarding our b estimation, our results agree with the two other studies at less than 2σ, depending on the case. Therefore, we conclude that our calibration process is able to recover estimates consistent with the previously derived values available in the literature.

### 3.2.3 Parameters estimation through the fundamental plane

The quality of our fundamental plane calibration can be assessed through its prediction, $R_{e,FP}$, against the measured one, $R_e$, and the same can be done for $\sigma_{e,FP}$ and $\sigma_e$. Except for a few outliers, most predictions are in good agreement with the measured values. The $R_e$ and $\sigma_e$ predictions show a mean relative error on the considered sample of 4.0 and -3.6 per cent, respectively. However, the uncertainties based on the calibrations show larger uncertainties. In fact, it is a factor of four of the observed errors for $\sigma_e$ and twice for $R_e$. As we are using the fundamental plane to recover the mass of cluster members, we are more interested in the accuracy and precision of the relations for a combination of $\sigma_e$ and $R_e$ rather than them independently as the combination appears in the mass estimate. In particular, the total mass of a model following a mass profile of a dual Pseudo-Isothermal Elliptical (dPIE; Limousin, Kneib & Natarajan 2005) follows $M_{tot,dPIE} \propto R_e\sigma_e^2$, if we assume that the central velocity, $\sigma_0$, and the cut radius, $r_{cut}$, are proportional to $\sigma_e$ and $R_e$, respectively.

We need two quantities from the plane to predict the third one, though only $R_e$ and $\sigma_e$ are relevant for the mass modelling. Thus, we use the combination of $\sigma_e$ and $\mu_e$ to predict $R_e$ (i.e. $R_{e,FP}$) and combine $R_e$ and $\mu_e$ to obtain $\sigma_{e,FP}$. From these two uses of the fundamental plane, we create two mass proxies. We estimate $R_e$ and $\sigma_e$ from the fundamental plane, which leads to the two mass proxies, $R_{e,FP}\sigma_e^2$ and $R_e\sigma_{e,FP}^2$. Both proxies agree well with the equivalent based on measured values only. Notably, the averaged relative error for $R_e\sigma_e^2$, and $R_{e,FP}\sigma_e^2$, are of -4.4 per cent and 8.0 per cent. Regarding the uncertainties, we assume that for the measurements-based proxies, $R_e$ and $\sigma_e$ are statistically independent, and we sample them according to their measurement uncertainties. The fundamental plane-based proxies uncertainties are assumed to be only from $\sigma_{e,FP}$ or $R_{e,FP}$, as $R_e$ or $\sigma_e$ will be fixed. Hence, for both proxies, the propagated uncertainties are ~4-5 times smaller than the observed ones. Though, these are large enough to be in agreement at less than the 0.6σ level with their measured counterparts. We note here that a proper comparison with the observational uncertainty will require the correlation between $R_e$ and $\sigma_e$ but it is beyond the scope of this paper.

## 4 MASS MODELLING

As previously mentioned, our modelling incorporates the usual cluster-scale model represented mainly by the DM, the intra-cluster gas, cluster members through the fundamental plane of elliptical galaxies, and the Faber & Jackson law and an additional perturbative component. We detail the modelling of each in order as Section 4.1 for the cluster-scale DM and Section 4.2 for the gas. We finish with

the galaxy-scale elements and the perturbative piece in Sections 4.3 and 4.4, respectively.

Except for the modelling of the additional perturbation, all the building blocks of our models are dPIE potentials (Limousin, Kneib & Natarajan 2005). Each of the different components is composed of a sum of this analytical profile with specific assumptions on their parameters. Such potentials are defined by seven parameters, the central coordinates, $(x_c, y_c)$, the position angle, $\theta$, the ellipticity, $\epsilon$, a central velocity dispersion, $\sigma_0$, a core radius, $r_{core}$, and a cut radius, $r_{cut}$. The 3D mass density, $\rho$, follows the relation:

$$ \rho(r) = \frac{r_{core} + r_{cut}}{2\pi G r_{core}^2 r_{cut}} \times \frac{\sigma_0^2}{\left(1 + \left(\frac{r}{r_{core}}\right)^2\right)\left(1 + \left(\frac{r}{r_{cut}}\right)^2\right)} \quad (13) $$

With $r$, the elliptical radius in the coordinates of the potential, and $G$, the universal constant of gravitation. Now, we outline how we parametrize these dPIE potentials for each of the cluster components, followed by the definition of the small B-spline perturbation that we add on top of the overall parametric modelling of the various components.

### 4.1 Large-scale dark matter distribution

Similarly to previous studies on AS1063 (Caminha et al. 2016; Limousin, Beauchesne & Jullo 2022), we use two haloes to represent the smooth DM component of the cluster. The main potential is associated with the brightest cluster galaxy (BCG) to account for most of the cluster DM, while the second one is placed near a higher concentration of galaxies in the NE of the cluster. This placement is shown as the dashed white diamonds in Fig. 1. The main cluster halo has all of its parameters free to vary except its position, and $r_{cut}$, which are fixed to the BCG centre, and to 3 Mpc (beyond the virial radius of the cluster), respectively. The parameter $r_{cut}$, is, in fact, ill-constrained due to the lack of strong lensing constraints in the cluster outskirts. Regarding the potential in the NE, we use a Bayesian regularization on its position. We assume it to be centred around the light distribution with Gaussian priors on its centre coordinate. It is similar to the assumptions made in Limousin, Beauchesne & Jullo (2022) and avoids dPIE halo positions to be inconsistent with the luminous distribution. We assume a standard deviation of 2 arcsec around the centre of the group of three galaxies. As there are no spectroscopically confirmed multiple images near the NE halo, its $r_{core}$ is ill-constrained and fixed to 0.5 arcsec. The remaining parameters of this halo are assigned uniform priors.

### 4.2 Intra-cluster gas

To model the X-ray emitting gas in the cluster, we use a similar approach to Bonamigo et al. (2017, 2018) for the definition of the gas mass component, which consists of a free-form modelling of the distribution with elliptical dPIE potentials. We deviate significantly from these previous works in the other aspects of the modelling, such as during parameter optimization or in model discrimination described in Sections 5.3 and 6.1, respectively. The shapes of these potentials are constrained by the X-ray surface brightness in the form of the number of photon counts and by strong lensing through their contributions to the total mass. We will detail how we discriminate the number of potentials to be used in Section 6.1.

All parameters of the dPIE potentials are free to vary except for the cut radius, $r_{cut}$. When using more than one potential, this parameter is degenerated, and our optimization process does not converge. Thus, we fix it to the large value of 1250 kpc for all of the potentials. This
of cluster members lies well into that plane, unlike the Faber & Jackson relations where some galaxies show a larger scatter from the relation. In comparison to G22 and B19, we have a higher estimation of c at 3.1σ, and a similar one at less than 1σ, respectively. However, we used different apertures than B19 and G22 to extract cluster member spectra that can partially explain the difference with G22 as well as in the joint optimization. Both plane estimations are in agreement with blank field results at less than 2σ in the R-band (Hyde & Bernardi 2009). Regarding our b estimation, our results agree with the two other studies at less than 2σ, depending on the case. Therefore, we conclude that our calibration process is able to recover estimates consistent with the previously derived values available in the literature.

### 3.2.3 Parameters estimation through the fundamental plane

The quality of our fundamental plane calibration can be assessed through its prediction, $R_{e,FP}$, against the measured one, $R_e$, and the same can be done for $\sigma_{e,FP}$ and $\sigma_e$. Except for a few outliers, most predictions are in good agreement with the measured values. The $R_e$ and $\sigma_e$ predictions show a mean relative error on the considered sample of 4.0 and -3.6 per cent, respectively. However, the uncertainties based on the calibrations show larger uncertainties. In fact, it is a factor of four of the observed errors for $\sigma_e$ and twice for $R_e$. As we are using the fundamental plane to recover the mass of cluster members, we are more interested in the accuracy and precision of the relations for a combination of $\sigma_e$ and $R_e$ rather than them independently as the combination appears in the mass estimate. In particular, the total mass of a model following a mass profile of a dual Pseudo-Isothermal Elliptical (dPIE; Limousin, Kneib & Natarajan 2005) follows $M_{tot,dPIE} \propto R_e\sigma_e^2$, if we assume that the central velocity, $\sigma_0$, and the cut radius, $r_{cut}$, are proportional to $\sigma_e$ and $R_e$, respectively.

We need two quantities from the plane to predict the third one, though only $R_e$ and $\sigma_e$ are relevant for the mass modelling. Thus, we use the combination of $\sigma_e$ and $\mu_e$ to predict $R_e$ (i.e. $R_{e,FP}$) and combine $R_e$ and $\mu_e$ to obtain $\sigma_{e,FP}$. From these two uses of the fundamental plane, we create two mass proxies. We estimate $R_e$ and $\sigma_e$ from the fundamental plane, which leads to the two mass proxies, $R_{e,FP}\sigma_e^2$ and $R_e\sigma_{e,FP}^2$. Both proxies agree well with the equivalent based on measured values only. Notably, the averaged relative error for $R_e\sigma_e^2$, and $R_{e,FP}\sigma_e^2$, are of -4.4 per cent and 8.0 per cent. Regarding the uncertainties, we assume that for the measurements-based proxies, $R_e$ and $\sigma_e$ are statistically independent, and we sample them according to their measurement uncertainties. The fundamental plane-based proxies uncertainties are assumed to be only from $\sigma_{e,FP}$ or $R_{e,FP}$, as $R_e$ or $\sigma_e$ will be fixed. Hence, for both proxies, the propagated uncertainties are ~4-5 times smaller than the observed ones. Though, these are large enough to be in agreement at less than the 0.6σ level with their measured counterparts. We note here that a proper comparison with the observational uncertainty will require the correlation between $R_e$ and $\sigma_e$ but it is beyond the scope of this paper.

## 4 MASS MODELLING

As previously mentioned, our modelling incorporates the usual cluster-scale model represented mainly by the DM, the intra-cluster gas, cluster members through the fundamental plane of elliptical galaxies, and the Faber & Jackson law and an additional perturbative component. We detail the modelling of each in order as Section 4.1 for the cluster-scale DM and Section 4.2 for the gas. We finish with

the galaxy-scale elements and the perturbative piece in Sections 4.3 and 4.4, respectively.

Except for the modelling of the additional perturbation, all the building blocks of our models are dPIE potentials (Limousin, Kneib & Natarajan 2005). Each of the different components is composed of a sum of this analytical profile with specific assumptions on their parameters. Such potentials are defined by seven parameters, the central coordinates, $(x_c, y_c)$, the position angle, $\theta$, the ellipticity, $\epsilon$, a central velocity dispersion, $\sigma_0$, a core radius, $r_{core}$, and a cut radius, $r_{cut}$. The 3D mass density, $\rho$, follows the relation:

$$ \rho(r) = \frac{r_{core} + r_{cut}}{2\pi G r_{core}^2 r_{cut}} \times \frac{\sigma_0^2}{\left(1 + \left(\frac{r}{r_{core}}\right)^2\right)\left(1 + \left(\frac{r}{r_{cut}}\right)^2\right)} \quad (13) $$

With $r$, the elliptical radius in the coordinates of the potential, and $G$, the universal constant of gravitation. Now, we outline how we parametrize these dPIE potentials for each of the cluster components, followed by the definition of the small B-spline perturbation that we add on top of the overall parametric modelling of the various components.

### 4.1 Large-scale dark matter distribution

Similarly to previous studies on AS1063 (Caminha et al. 2016; Limousin, Beauchesne & Jullo 2022), we use two haloes to represent the smooth DM component of the cluster. The main potential is associated with the brightest cluster galaxy (BCG) to account for most of the cluster DM, while the second one is placed near a higher concentration of galaxies in the NE of the cluster. This placement is shown as the dashed white diamonds in Fig. 1. The main cluster halo has all of its parameters free to vary except its position, and $r_{cut}$, which are fixed to the BCG centre, and to 3 Mpc (beyond the virial radius of the cluster), respectively. The parameter $r_{cut}$, is, in fact, ill-constrained due to the lack of strong lensing constraints in the cluster outskirts. Regarding the potential in the NE, we use a Bayesian regularization on its position. We assume it to be centred around the light distribution with Gaussian priors on its centre coordinate. It is similar to the assumptions made in Limousin, Beauchesne & Jullo (2022) and avoids dPIE halo positions to be inconsistent with the luminous distribution. We assume a standard deviation of 2 arcsec around the centre of the group of three galaxies. As there are no spectroscopically confirmed multiple images near the NE halo, its $r_{core}$ is ill-constrained and fixed to 0.5 arcsec. The remaining parameters of this halo are assigned uniform priors.

### 4.2 Intra-cluster gas

To model the X-ray emitting gas in the cluster, we use a similar approach to Bonamigo et al. (2017, 2018) for the definition of the gas mass component, which consists of a free-form modelling of the distribution with elliptical dPIE potentials. We deviate significantly from these previous works in the other aspects of the modelling, such as during parameter optimization or in model discrimination described in Sections 5.3 and 6.1, respectively. The shapes of these potentials are constrained by the X-ray surface brightness in the form of the number of photon counts and by strong lensing through their contributions to the total mass. We will detail how we discriminate the number of potentials to be used in Section 6.1.

All parameters of the dPIE potentials are free to vary except for the cut radius, $r_{cut}$. When using more than one potential, this parameter is degenerated, and our optimization process does not converge. Thus, we fix it to the large value of 1250 kpc for all of the potentials. This
<figure><img src="image_7.png" /><figcaption>**Figure 4.** Left-hand panel: colour–magnitude diagram of m<sub>F435W</sub> – m<sub>F606W</sub> versus m<sub>F814W</sub>. The blue and orange dots highlight the spectroscopically confirmed cluster members, and the orange ones represent the objects that were not rejected by the iterative procedure to fit the red sequence. The green and purple dots show the full photometric catalogue. The former represents all objects, while the second represents the galaxies that have been selected by both colour–magnitude diagrams. The red plain and dashed lines represent the best-fitting red sequence and the 3σ bounds used to select galaxies, respectively. The dashed black lines show the magnitude threshold of 0.1L*. Right-panel panel: Same as the left-hand panel but for the colour–magnitude diagram of m<sub>F606W</sub> – m<sub>F814W</sub> versus m<sub>F814W</sub>. There are more purple dots in this panel than in the left one, as some objects do not have a m<sub>F435W</sub> measurement.</figcaption></figure>

choice corresponds to the radius at which the X-ray signal starts to be equivalent to the noise. As we only use a few dPIE profiles, we only model the cluster-scale distribution of the gas. Hence, small-scale variations like gas sloshing, shock front, or micro-physics processes are not taken into account in our current approach. If such processes have enough amplitude, they will leave traces in the count residual between the observations and our models.

## 4.3 Cluster member masses

We select the cluster members using a spectroscopic confirmation combined with a red sequence selection, if spectroscopy is not available. We include all galaxies that have a redshift, z, in the range of [0.327, 0.360], and are in the HFF footprint, independently of their type, similar to Lagattuta et al. (2017). We use this galaxy sample to calibrate the red sequence and select the cluster member without spectroscopic confirmation. We perform a red sequence fitting on the (m<sub>F606W</sub> – m<sub>F814W</sub>, m<sub>F814W</sub>) and the (m<sub>F435W</sub> – m<sub>F606W</sub>, m<sub>F814W</sub>) planes. We use a similar iterative approach as in our fit for the fundamental plane and Faber & Jackson scaling relations. We fit a linear model in the considered plane, rejecting all elements that are at more than 2.5σ of the fitted lines. We reiterate this procedure until no galaxies are rejected. We found that this value of 2.5σ leads to fast convergence of the rejection process. To constitute the red sequence selected sample, we consider all galaxies that are brighter than 0.1L* and are within 3σ of the relation of fitted on both planes. In case some galaxies have no measurement for m<sub>F435W</sub>, we select galaxies with the other relation only. With the addition of these photometrically selected galaxies, we reach a total of 250 elements. Fig. 4 presents the two colour–magnitude diagrams previously mentioned with the spectroscopic and photometric samples. The red plain and dashed lines show the two red sequence relations, representing the best-fitting relation and the 3σ range. As we can see by the number of purple dots (i.e. red sequence selected cluster member) in comparison to the spectroscopic sample (e.g. blue and orange dots), most of the cluster members used in the modelling are spectroscopically confirmed.

We assume, to some extent, a light-traces-mass hypothesis for profiles associated with selected cluster members. Indeed, we follow the same geometry by fixing (x<sub>c</sub>, y<sub>c</sub>), θ and ε, to the values measured

in the light distribution with SExTRACTOR, and use a dPIE profile to represent both. For the remaining parameters, we suppose that light and mass profiles have the same shape with only a scaling as a difference. This allows us to link σ<sub>0</sub> from the measurement of σ<sub>e</sub> with the projection factor c<sub>p</sub> (its calculation is detailed in Appendix C). The scaling is defined as follows on the three relevant parameters:

$$ \sigma_{0,\text{light}} = \sigma_0 \qquad (14) $$

$$ v r_{\text{cut,light}} = r_{\text{cut}} \qquad (15) $$

$$ v r_{\text{core,light}} = r_{\text{core}} \qquad (16) $$

Where v is the scaling factor between mass and light profiles. σ<sub>0,light</sub>, r<sub>cut, light</sub> and r<sub>core, light</sub> are dPIE parameters associated with the light profile. These hypotheses imply a constant mass-to-light ratio of v as well as a proportionality relation between the R<sub>e</sub> and half-mass radius, R<sub>M/2</sub>, in the form of R<sub>M/2</sub> = vR<sub>e</sub>. Regarding the proper scaling of r<sub>cut</sub>, r<sub>core</sub> and σ<sub>0</sub>, we apply different relations depending on the available observational data. We split the cluster member selection into three groups depending on the photometric measurements that are publicly available. If no light profile has been fitted, meaning that we only have the integrated luminosity, L, we use a Faber & Jackson scaling, and if we have it, we take benefit of the fundamental plane with one relation, of our two different scalings.

### 4.3.1 Case 1: if only L is available

In this case, we rely on the Faber & Jackson law with small differences in the power law in comparison to previous works such as Richard et al. (2010) and Bergamini et al. (2019). We use the following formulae:

$$ r_{\text{core}} = r_{\text{core}}^{*} \left( \frac{L}{L^{*}} \right)^{1-1/c} \qquad (17) $$

$$ r_{\text{cut}} = r_{\text{cut}}^{*} \left( \frac{L}{L^{*}} \right)^{1-1/c} \qquad (18) $$

$$ \sigma_0 = \frac{\sigma_e^*}{c_p (r_{\text{cut}}^*/r_{\text{core}}^*, v)} \left(\frac{L}{L^*}\right)^{1/2c} \qquad (19) $$
<figure><img src="image_7.png" /><figcaption>**Figure 4.** Left-hand panel: colour–magnitude diagram of m<sub>F435W</sub> – m<sub>F606W</sub> versus m<sub>F814W</sub>. The blue and orange dots highlight the spectroscopically confirmed cluster members, and the orange ones represent the objects that were not rejected by the iterative procedure to fit the red sequence. The green and purple dots show the full photometric catalogue. The former represents all objects, while the second represents the galaxies that have been selected by both colour–magnitude diagrams. The red plain and dashed lines represent the best-fitting red sequence and the 3σ bounds used to select galaxies, respectively. The dashed black lines show the magnitude threshold of 0.1L*. Right-panel panel: Same as the left-hand panel but for the colour–magnitude diagram of m<sub>F606W</sub> – m<sub>F814W</sub> versus m<sub>F814W</sub>. There are more purple dots in this panel than in the left one, as some objects do not have a m<sub>F435W</sub> measurement.</figcaption></figure>

choice corresponds to the radius at which the X-ray signal starts to be equivalent to the noise. As we only use a few dPIE profiles, we only model the cluster-scale distribution of the gas. Hence, small-scale variations like gas sloshing, shock front, or micro-physics processes are not taken into account in our current approach. If such processes have enough amplitude, they will leave traces in the count residual between the observations and our models.

## 4.3 Cluster member masses

We select the cluster members using a spectroscopic confirmation combined with a red sequence selection, if spectroscopy is not available. We include all galaxies that have a redshift, z, in the range of [0.327, 0.360], and are in the HFF footprint, independently of their type, similar to Lagattuta et al. (2017). We use this galaxy sample to calibrate the red sequence and select the cluster member without spectroscopic confirmation. We perform a red sequence fitting on the (m<sub>F606W</sub> – m<sub>F814W</sub>, m<sub>F814W</sub>) and the (m<sub>F435W</sub> – m<sub>F606W</sub>, m<sub>F814W</sub>) planes. We use a similar iterative approach as in our fit for the fundamental plane and Faber & Jackson scaling relations. We fit a linear model in the considered plane, rejecting all elements that are at more than 2.5σ of the fitted lines. We reiterate this procedure until no galaxies are rejected. We found that this value of 2.5σ leads to fast convergence of the rejection process. To constitute the red sequence selected sample, we consider all galaxies that are brighter than 0.1L* and are within 3σ of the relation of fitted on both planes. In case some galaxies have no measurement for m<sub>F435W</sub>, we select galaxies with the other relation only. With the addition of these photometrically selected galaxies, we reach a total of 250 elements. Fig. 4 presents the two colour–magnitude diagrams previously mentioned with the spectroscopic and photometric samples. The red plain and dashed lines show the two red sequence relations, representing the best-fitting relation and the 3σ range. As we can see by the number of purple dots (i.e. red sequence selected cluster member) in comparison to the spectroscopic sample (e.g. blue and orange dots), most of the cluster members used in the modelling are spectroscopically confirmed.

We assume, to some extent, a light-traces-mass hypothesis for profiles associated with selected cluster members. Indeed, we follow the same geometry by fixing (x<sub>c</sub>, y<sub>c</sub>), θ and ε, to the values measured

in the light distribution with SExTRACTOR, and use a dPIE profile to represent both. For the remaining parameters, we suppose that light and mass profiles have the same shape with only a scaling as a difference. This allows us to link σ<sub>0</sub> from the measurement of σ<sub>e</sub> with the projection factor c<sub>p</sub> (its calculation is detailed in Appendix C). The scaling is defined as follows on the three relevant parameters:

$$ \sigma_{0,\text{light}} = \sigma_0 \qquad (14) $$

$$ \nu r_{\text{cut,light}} = r_{\text{cut}} \qquad (15) $$

$$ \nu r_{\text{core,light}} = r_{\text{core}} \qquad (16) $$

Where ν is the scaling factor between mass and light profiles. σ<sub>0,light</sub>, r<sub>cut,light</sub> and r<sub>core,light</sub> are dPIE parameters associated with the light profile. These hypotheses imply a constant mass-to-light ratio of ν as well as a proportionality relation between the R<sub>e</sub> and half-mass radius, R<sub>M/2</sub>, in the form of R<sub>M/2</sub> = νR<sub>e</sub>. Regarding the proper scaling of r<sub>cut</sub>, r<sub>core</sub> and σ<sub>0</sub>, we apply different relations depending on the available observational data. We split the cluster member selection into three groups depending on the photometric measurements that are publicly available. If no light profile has been fitted, meaning that we only have the integrated luminosity, L, we use a Faber & Jackson scaling, and if we have it, we take benefit of the fundamental plane with one relation, of our two different scalings.

### 4.3.1 Case 1: if only L is available

In this case, we rely on the Faber & Jackson law with small differences in the power law in comparison to previous works such as Richard et al. (2010) and Bergamini et al. (2019). We use the following formulae:

$$ r_{\text{core}} = r_{\text{core}}^* \left( \frac{L}{L^*} \right)^{1-1/c} \qquad (17) $$

$$ r_{\text{cut}} = r_{\text{cut}}^* \left( \frac{L}{L^*} \right)^{1-1/c} \qquad (18) $$

$$ \sigma_0 = \frac{\sigma_e^*}{c_p (r_{\text{cut}}^*/r_{\text{core}}^*, \nu)} \left(\frac{L}{L^*}\right)^{1/2c} \qquad (19) $$
Where $r_{\text{core}}^*$, $r_{\text{cut}}^*$, and $\sigma_e^*$ are the parameters associated with the galaxy representing the population at that redshift defined in Section 3.2.2. In our case, we follow previous studies (Richard et al. 2010; Mahler et al. 2018; Lagattuta et al. 2019), and we fix $r_{\text{core}}^*$ to 0.15 kpc as its effect is negligible on strong lensing constraints used here.

### 4.3.2 Case 2: if $\sigma_e$ and the light-profile are available

If available, we use the $\sigma_e$ measurement, and we rely on the fundamental plane prediction for $R_{e,\text{FP}}$ combined with our assumption on the relationship between light and mass profiles. Thus, we obtain the following relations:

$$r_{\text{core}} = \frac{r_{\text{core}}^*}{r_{\text{cut}}^*} r_{\text{cut}} \qquad (20)$$

$$r_{\text{cut}} = \frac{\nu 10^{a+b\mu_e+c \log \sigma_e}}{\sqrt[3]{1 + \frac{10}{3} \frac{r_{\text{core}}^*}{r_{\text{cut}}^*} + \left(\frac{r_{\text{core}}^*}{r_{\text{cut}}^*}\right)^2}} \qquad (21)$$

$$\sigma_0 = \frac{\sigma_e}{c_p \left( \frac{r_{\text{cut}}^*}{r_{\text{core}}^*}, \nu \right)} \qquad (22)$$

The denominator of the $r_{\text{cut}}$ relation comes from the transformation of $R_{M/2}$ to $r_{\text{cut}}$ with the constant ratio assumed between $r_{\text{core}}$ and $r_{\text{cut}}$. Hence, we do not make direct use of $R_e$ measurement, and we instead use the fundamental plane estimation, which allows us to propagate partially the observational error on these measurements in our MCMC sampling. The proper solution would be to consider $\sigma_e$ and $R_e$ as random variables, but this would overwhelmingly increase the computational cost of the analysis.

### 4.3.3 Case 3: if the light-profile is available

We use a similar scaling as the previous case, and the only difference is that we are using $\sigma_{e,\text{FP}}$ instead as we do not have a measurement in this case. Thus, now the relations are the following:

$$r_{\text{core}} = \frac{r_{\text{core}}^*}{r_{\text{cut}}^*} r_{\text{cut}} \qquad (23)$$

$$r_{\text{cut}} = \frac{\nu R_e}{\sqrt[3]{1 + \frac{10}{3} \frac{r_{\text{core}}^*}{r_{\text{cut}}^*} + \left(\frac{r_{\text{core}}^*}{r_{\text{cut}}^*}\right)^2}} \qquad (24)$$

$$\sigma_0 = \frac{10^{(-a-b\mu_e+\log R_e)/c}}{c_p \left( \frac{r_{\text{cut}}^*}{r_{\text{core}}^*}, \nu \right)} \qquad (25)$$

Between the two relations in cases 2 and 3, we prioritize the case 2, which is motivated by the accuracy of the associated mass proxies in recovering the true parameter value from the fundamental plane, as discussed in Section 3.2.3. We use the same treatment for both of them as we use the same prior based on our calibration for the fundamental plane parameters (i.e. Gaussian priors based on the mean and standard deviation of the calibration). Notably, both proxies present similar trends but with a different amplitude in the discrepancy between the propagated uncertainty and the observed one. We choose not to treat them differently as the insights from our mass proxy analyses are not robust enough at the present time because we lack the covariance between the $R_e$ and $\sigma_e$.

From all the parameters mentioned, only $r_{\text{cut}}^*$ and $\nu$ are optimized with uniform priors as they are not involved in the calibration process. This ensures consistency with the galaxies modelled with the fundamental planes and Faber & Jackson relation, as in both cases, the parameter that we optimize with minimum prior knowledge is $r_{\text{cut}}$.

Regarding the BCG, we apply a special treatment as it is not modelled according to the previous relations. Its associated dPIE has

its shape and position parameters fixed by its luminous distribution as for other cluster members, but its $r_{\text{core}}$, $r_{\text{cut}}$, and $\sigma_0$ are optimized. To define physically motivated priors, we ran a nested sampling algorithm to fit its LOSVD with data points measured at different radii for the BCG. We used the $\sigma_0$ posterior distribution to obtain a Gaussian prior for the sampling using the lensing and X-ray constraint. $r_{\text{core}}$ and $r_{\text{cut}}$ did not converge in the previous run; thus, we assign uniform priors to them.

## 4.4 Additional perturbation

To complete our modelling, we consider a perturbation under the form of a B-spline surface added to the lensing potential (Beauchesne et al. 2021). It allows us to incorporate effects like complex asymmetries that can not be captured with the number of large-scale haloes considered to represent the smooth DM distribution. The surface expression, $\Delta\psi$, is the following:

$$\Delta\psi(x, y) = \frac{D_{ls}}{D_s} \sum_{j,l=1}^{m} C_{j,l} B_{j,p,t_x}(x) B_{l,p,t_y}(y) \qquad (26)$$

where $B_{j,p,t_x}$ are the $j^{th}$ 1D B-spline basis functions of polynomial degree, $p$, associated with the knot vector, $t_x$. $C_{j,l}$ are the coefficients of each B-spline basis, and $m$ is the number of B-spline functions per axis. $D_{ls}$ and $D_s$ are the angular distance between the lens and the source, and the observer and the source, respectively.

This component is added to the model at a different step than the previous ones as we follow a two-step optimization (Beauchesne et al. 2021), where we first run a nested sampling with the priors previously defined on the dPIEs parameters without the perturbation. We use the obtained posterior distributions to define Gaussian priors on each parameter for the second sampling run. We use the best-fitting values as the mean, and three times the empirical standard deviation. For the physically motivated priors related to cluster members, we use the mean of the posterior and not the best-fitting values as it does not take into account the bias from the priors and only one time the empirical standard deviation, as the widths of the posterior distributions are similar to the ones of their priors.

We tested the addition of an external shear (Mahler et al. 2018; Lagattuta et al. 2019) at the same step to take into account the effect of haloes on the cluster outskirts. However, it did not provide improvements on the reconstruction of the lensing constraints, so we choose not to include it in the final method. In the case of cluster lenses such as Abell 2744 (Mahler et al. 2018) where structures surrounding the cluster core are clearly seen, it is important to include this kind of physically supported perturbation.

# 5 JOINT LENSING AND X-RAY LIKELIHOOD

We now outline how we combine the two probes of the mass distribution of the galaxy cluster by detailing what quantities we use as constrain, and how the likelihoods are defined and combined. We start with the X-ray and lensing constraints in Sections 5.1 and 5.2, respectively, and then go on to describe the likelihoods in Section 5.3.

## 5.1 X-ray emission: constraints and modelling

As indicated in Section 3.1, we fit the X-ray counts in the 0.5–7 keV band, which corresponds to the combined range of the soft, medium, and hard science energy bands of Chandra. To define a mask where the X-ray emission is taken into account, we use the radius at which the X-ray signal is equivalent to the associated noise. We obtain this radius by making circularly averaged profiles of the counts from both
Where $r_{\text{core}}^*$, $r_{\text{cut}}^*$, and $\sigma_e^*$ are the parameters associated with the galaxy representing the population at that redshift defined in Section 3.2.2. In our case, we follow previous studies (Richard et al. 2010; Mahler et al. 2018; Lagattuta et al. 2019), and we fix $r_{\text{core}}^*$ to 0.15 kpc as its effect is negligible on strong lensing constraints used here.

### 4.3.2 Case 2: if $\sigma_e$ and the light-profile are available

If available, we use the $\sigma_e$ measurement, and we rely on the fundamental plane prediction for $R_{e\text{FP}}$ combined with our assumption on the relationship between light and mass profiles. Thus, we obtain the following relations:

$$r_{\text{core}} = \frac{r_{\text{core}}^*}{r_{\text{cut}}^*} r_{\text{cut}} \qquad (20)$$

$$r_{\text{cut}} = \frac{\nu 10^{a+b\mu_e+c \log \sigma_e}}{\frac{3}{4}\sqrt{1+\frac{10}{3}\frac{r_{\text{core}}^*}{r_{\text{cut}}^*}+\left(\frac{r_{\text{core}}^*}{r_{\text{cut}}^*}\right)^2}} \qquad (21)$$

$$\sigma_0 = \frac{\sigma_e}{c_p \left( \frac{r_{\text{cut}}^*}{r_{\text{core}}^*}, \nu \right)} \qquad (22)$$

The denominator of the $r_{\text{cut}}$ relation comes from the transformation of $R_{M/2}$ to $r_{\text{cut}}$ with the constant ratio assumed between $r_{\text{core}}$ and $r_{\text{cut}}$. Hence, we do not make direct use of $R_e$ measurement, and we instead use the fundamental plane estimation, which allows us to propagate partially the observational error on these measurements in our MCMC sampling. The proper solution would be to consider $\sigma_e$ and $R_e$ as random variables, but this would overwhelmingly increase the computational cost of the analysis.

### 4.3.3 Case 3: if the light-profile is available

We use a similar scaling as the previous case, and the only difference is that we are using $\sigma_{e,\text{FP}}$ instead as we do not have a measurement in this case. Thus, now the relations are the following:

$$r_{\text{core}} = \frac{r_{\text{core}}^*}{r_{\text{cut}}^*} r_{\text{cut}} \qquad (23)$$

$$r_{\text{cut}} = \frac{\nu R_e}{\frac{3}{4}\sqrt{1 + \frac{10}{3}\frac{r_{\text{core}}^*}{r_{\text{cut}}^*} + \left(\frac{r_{\text{core}}^*}{r_{\text{cut}}^*}\right)^2}} \qquad (24)$$

$$\sigma_0 = \frac{10^{(-a-b\mu_e+\log R_e)/c}}{c_p \left( \frac{r_{\text{cut}}^*}{r_{\text{core}}^*}, \nu \right)} \qquad (25)$$

Between the two relations in cases 2 and 3, we prioritize the case 2, which is motivated by the accuracy of the associated mass proxies in recovering the true parameter value from the fundamental plane, as discussed in Section 3.2.3. We use the same treatment for both of them as we use the same prior based on our calibration for the fundamental plane parameters (i.e. Gaussian priors based on the mean and standard deviation of the calibration). Notably, both proxies present similar trends but with a different amplitude in the discrepancy between the propagated uncertainty and the observed one. We choose not to treat them differently as the insights from our mass proxy analyses are not robust enough at the present time because we lack the covariance between the $R_e$ and $\sigma_e$.

From all the parameters mentioned, only $r_{\text{cut}}^*$ and $\nu$ are optimized with uniform priors as they are not involved in the calibration process. This ensures consistency with the galaxies modelled with the fundamental planes and Faber & Jackson relation, as in both cases, the parameter that we optimize with minimum prior knowledge is $r_{\text{cut}}$.

Regarding the BCG, we apply a special treatment as it is not modelled according to the previous relations. Its associated dPIE has

its shape and position parameters fixed by its luminous distribution as for other cluster members, but its $r_{\text{core}}$, $r_{\text{cut}}$, and $\sigma_0$ are optimized. To define physically motivated priors, we ran a nested sampling algorithm to fit its LOSVD with data points measured at different radii for the BCG. We used the $\sigma_0$ posterior distribution to obtain a Gaussian prior for the sampling using the lensing and X-ray constraint. $r_{\text{core}}$ and $r_{\text{cut}}$ did not converge in the previous run; thus, we assign uniform priors to them.

## 4.4 Additional perturbation

To complete our modelling, we consider a perturbation under the form of a B-spline surface added to the lensing potential (Beauchesne et al. 2021). It allows us to incorporate effects like complex asymmetries that can not be captured with the number of large-scale haloes considered to represent the smooth DM distribution. The surface expression, $\Delta\psi$, is the following:

$$\Delta\psi(x, y) = \frac{D_{ls}}{D_s} \sum_{j,l=1}^{m} C_{j,l} B_{j,p,t_x}(x) B_{l,p,t_y}(y) \qquad (26)$$

where $B_{j,p,t_x}$ are the $j^{th}$ 1D B-spline basis functions of polynomial degree, $p$, associated with the knot vector, $t_x$. $C_{j,l}$ are the coefficients of each B-spline basis, and $m$ is the number of B-spline functions per axis. $D_{ls}$ and $D_s$ are the angular distance between the lens and the source, and the observer and the source, respectively.

This component is added to the model at a different step than the previous ones as we follow a two-step optimization (Beauchesne et al. 2021), where we first run a nested sampling with the priors previously defined on the dPIEs parameters without the perturbation. We use the obtained posterior distributions to define Gaussian priors on each parameter for the second sampling run. We use the best-fitting values as the mean, and three times the empirical standard deviation. For the physically motivated priors related to cluster members, we use the mean of the posterior and not the best-fitting values as it does not take into account the bias from the priors and only one time the empirical standard deviation, as the widths of the posterior distributions are similar to the ones of their priors.

We tested the addition of an external shear (Mahler et al. 2018; Lagattuta et al. 2019) at the same step to take into account the effect of haloes on the cluster outskirts. However, it did not provide improvements on the reconstruction of the lensing constraints, so we choose not to include it in the final method. In the case of cluster lenses such as Abell 2744 (Mahler et al. 2018) where structures surrounding the cluster core are clearly seen, it is important to include this kind of physically supported perturbation.

# 5 JOINT LENSING AND X-RAY LIKELIHOOD

We now outline how we combine the two probes of the mass distribution of the galaxy cluster by detailing what quantities we use as constrain, and how the likelihoods are defined and combined. We start with the X-ray and lensing constraints in Sections 5.1 and 5.2, respectively, and then go on to describe the likelihoods in Section 5.3.

## 5.1 X-ray emission: constraints and modelling

As indicated in Section 3.1, we fit the X-ray counts in the 0.5–7 keV band, which corresponds to the combined range of the soft, medium, and hard science energy bands of Chandra. To define a mask where the X-ray emission is taken into account, we use the radius at which the X-ray signal is equivalent to the associated noise. We obtain this radius by making circularly averaged profiles of the counts from both
the source and the blank-sky background and measuring where the signal is going below two times the noise (2σ). Thus, we use a square mask of 1250 kpc centred on the position of the centroid of the BCG light distribution.

To constrain the gas distribution, we associate the dPIE potentials with the X-ray surface brightness. We start by obtaining, $\int \rho_{gas}^2$, along the LOS with the analytical formulae computed by Bonamigo et al. (2017). To switch from the previous quantity to the definition of the surface brightness outlined in equation 1, we need to transform $\rho_{gas}^2$, with $\rho_{gas}$ the gas mass density, into $n_e n_p$ for a fully ionized gas, which is done in the following way:

$$ \rho_{gas} = \mu(n_e + n_p) \quad (27) $$

$$ \Rightarrow n_e n_p = \frac{\rho^2}{\mu^2 (1 + \text{nhc}) (1 + \frac{1}{\text{nhc}})} \quad (28) $$

With $\mu$, the mean molecular weight per particle in a fully ionized gas, and nhc, the conversion factor from $n_p$ to $n_e$. For both parameters, we assume the value tabulated in Asplund et al. (2009). The total count number is obtained by multiplying $\int n_e n_p$ by the exposure time, the effective area of the telescope, and the cooling function $\Lambda(Z, T, E)$. $\Lambda(Z, T, E)$ used in this analysis has been computed with a binning of the FoV made with an SN threshold of 10 as detailed in Section 3.1, but we analyse its possible bias in Appendix E with SN thresholds of 6 and 14. As we fit the count number on the broad-band of Chandra science band, we estimate the effective area of the CCD pixels at different energies. In particular, we use the following approximation:

$$ \int t_{\text{exp}} A_{\text{eff}}(E) \Lambda(Z, T, E) dE = \sum_{\text{Band}} t_{\text{exp}} A_{\text{eff}}(E_{\text{Band}}) \quad (29) $$

$$ \times \int_{E_{\text{min,Band}}}^{E_{\text{max,Band}}} \Lambda(Z, T, E) dE \quad (30) $$

Where Band refers to the soft, medium, and hard energy bands of Chandra, $t_{exp}$, the exposure time, $A_{eff}$, effective area of the CCD pixels. $E_{Band}$ is the energy value associated with each of these bands for the computation of the exposure map, and $E_{min/max, Band}$ are the associated bounds. We finally add the blank-sky background to obtain the count model.

## 5.2 Multiply imaged systems

Regarding strong lensing constraints, we build our set of multiply imaged systems from the one put in place by the HFF modelling teams. In particular, we use the Gold sample (see section 4.2 from Lagattuta et al. 2019, for the explanation of all HFF labels), which is a set of 50 images over 19 systems all spectroscopically confirmed. We complete it with other additional spectroscopically confirmed systems presented in Caminha et al. (2016), which brings the total to 61 multiple images from 22 systems. Notably, we have three more images and one more system than previous works, which is due to different selection of multiple image since the objects are present in both sets.

Thanks to our new reduction of the MUSE data cube, we are able to measure the spectroscopic redshift of some objects for the first time, which allows us to also identify a new pair of multiple images. One object (system 303.1) is already present in previous spectroscopic catalogues, while the other comes from these new measurements. In addition, we identify a close pair as a multiply imaged system that was previously considered as a single image (SW - 54 in Karman et al. 2015, 2017, and system 304 in this work) and a second pair

from one confirmed image and the MUSE narrow-band data cube (system 305). These new constraints bring our final data set to a total of 67 images for 25 systems that are indicated with red circles in Fig. 1. The spectra and optical images of these new systems are presented in Appendix D as well as their coordinates and redshifts. The method used to find these new systems is detailed in Richard et al. (2021).

## 5.3 Likelihoods

In our optimization process, we simultaneously constrain the mass distribution with the strong lensing and the X-ray surface brightness. We obtain the combined likelihood, $\mathcal{L}_{tot}$, by assuming that both observables are independent, which leads to the following expression:

$$ \mathcal{L}_{\text{tot}} = \mathcal{L}_{\text{SL}} \times \mathcal{L}_{\text{X-ray}} \quad (31) $$

where $\mathcal{L}_{SL}$ and $\mathcal{L}_{X-ray}$ are the likelihoods for the strong lensing only and the X-ray only, respectively. For the latter one, we use a Poisson-Gamma mixture likelihood. For the former, we use a Gaussian likelihood slightly modified in comparison to previous works based on LENSTOOL (Lagattuta et al. 2017; Mahler et al. 2018).

### 5.3.1 Constructing the X-ray likelihood function

As our modelling of the X-ray gas mass solely uses a couple of dPIE potentials, we can only reproduce the large-scale distribution, excluding all microphysical gas phenomena. These small-scale events in the plasma are linked to surface brightness fluctuations (Eckert et al. 2017), which we take into account by adding an uncertainty on the mean, $\lambda_i$, of the Poisson distribution of the observed photon counts in the $i^{th}$ bin. We assume that $\lambda_i$ follows a Gamma distribution, $\Gamma(\mu_i, \sigma_X^2)$, parametrized by its expected value $\mu_i$ and its variance $\sigma_X^2$ instead of the usual shape and scale parameters. $\mu_i$ is our predicted count model excluding the microphysics, whereas $\sigma_X^2$ represents the induced uncertainty assumed to be the same for all bins. As we only intend to model the count model uncertainty statistically, we marginalize on the realization of the Gamma distribution, which leads to the following likelihood in the $i^{th}$ bin to observe $k_i$ count:

$$ \mathcal{L}_{\text{X-ray},i} = \int_{0}^{\infty} \mathbb{P}_{i}(X_{\lambda} = k) f_{\Gamma}(\lambda, \mu_{i}, \sigma_{X}^{2}) d\lambda \quad (32) $$

where $\mathbb{P}_i$ is the probability mass function of the Poisson distribution and $f_\Gamma$ is the probability density function of $\Gamma(\mu_i, \sigma_X^2)$. This compound probability distribution is known to be a continuous extension of a negative binomial distribution $\mathbb{N}\mathbb{B}(r = \mu_i^2/\sigma_X^2, p = \mu_i/(\mu_i + \sigma_X^2))$ with the following expression:

$$ \mathcal{L}_{\text{X-ray},i} = \frac{\Gamma(k_i + r)}{k_i! \Gamma(r)} p^r (1 - p)^k \quad (33) $$

The expected value of $\mathbb{N}\mathbb{B}$ is $\mu_i$, and the variance is $\mu_i + \sigma_X^2$, which implies an overdispersion of the Poisson model in the case of $\sigma_X > 0$. We fall back to a Poisson modelling when the uncertainty is zero, with the expected value being equal to the variance. We can now obtain our final likelihood by injecting $r$ and $p$ expressions in the previous definitions and assuming that each bin is statistically independent of the others:

$$ \mathcal{L}_{\text{X-ray}} = \prod_i \frac{\Gamma\left(k_i + \frac{\mu_i^2}{\sigma_X^2}\right)}{k_i! \Gamma\left(\frac{\mu_i^2}{\sigma_X^2}\right)} \left(\frac{\mu_i}{\mu_i + \sigma_X^2}\right)^{\frac{\mu_i^2}{\sigma_X^2}} \left(\frac{\sigma_X^2}{\mu_i + \sigma_X^2}\right)^{k_i} \quad (34) $$
the source and the blank-sky background and measuring where the signal is going below two times the noise (2σ). Thus, we use a square mask of 1250 kpc centred on the position of the centroid of the BCG light distribution.

To constrain the gas distribution, we associate the dPIE potentials with the X-ray surface brightness. We start by obtaining, $\int \rho_{gas}^2$, along the LOS with the analytical formulae computed by Bonamigo et al. (2017). To switch from the previous quantity to the definition of the surface brightness outlined in equation 1, we need to transform $\rho_{gas}^2$, with $\rho_{gas}$ the gas mass density, into $n_e n_p$ for a fully ionized gas, which is done in the following way:

$$\rho_{gas} = \mu(n_e + n_p) \qquad (27)$$

$$\Rightarrow n_e n_p = \frac{\rho^2}{\mu^2 (1 + \text{nhc}) (1 + \frac{1}{\text{nhc}})} \qquad (28)$$

With $\mu$, the mean molecular weight per particle in a fully ionized gas, and nhc, the conversion factor from $n_p$ to $n_e$. For both parameters, we assume the value tabulated in Asplund et al. (2009). The total count number is obtained by multiplying $\int n_e n_p$ by the exposure time, the effective area of the telescope, and the cooling function $\Lambda(Z, T, E)$. $\Lambda(Z, T, E)$ used in this analysis has been computed with a binning of the FoV made with an SN threshold of 10 as detailed in Section 3.1, but we analyse its possible bias in Appendix E with SN thresholds of 6 and 14. As we fit the count number on the broad-band of Chandra science band, we estimate the effective area of the CCD pixels at different energies. In particular, we use the following approximation:

$$\int t_{\text{exp}} A_{\text{eff}}(E) \Lambda(Z, T, E) dE = \sum_{\text{Band}} t_{\text{exp}} A_{\text{eff}}(E_{\text{Band}}) \qquad (29)$$

$$\times \int_{E_{\text{min,Band}}}^{E_{\text{max,Band}}} \Lambda(Z, T, E) dE \qquad (30)$$

Where Band refers to the soft, medium, and hard energy bands of Chandra, $t_{exp}$, the exposure time, $A_{eff}$, effective area of the CCD pixels. $E_{Band}$ is the energy value associated with each of these bands for the computation of the exposure map, and $E_{min/max, Band}$ are the associated bounds. We finally add the blank-sky background to obtain the count model.

## 5.2 Multiply imaged systems

Regarding strong lensing constraints, we build our set of multiply imaged systems from the one put in place by the HFF modelling teams. In particular, we use the Gold sample (see section 4.2 from Lagattuta et al. 2019, for the explanation of all HFF labels), which is a set of 50 images over 19 systems all spectroscopically confirmed. We complete it with other additional spectroscopically confirmed systems presented in Caminha et al. (2016), which brings the total to 61 multiple images from 22 systems. Notably, we have three more images and one more system than previous works, which is due to different selection of multiple image since the objects are present in both sets.

Thanks to our new reduction of the MUSE data cube, we are able to measure the spectroscopic redshift of some objects for the first time, which allows us to also identify a new pair of multiple images. One object (system 303.1) is already present in previous spectroscopic catalogues, while the other comes from these new measurements. In addition, we identify a close pair as a multiply imaged system that was previously considered as a single image (SW - 54 in Karman et al. 2015, 2017, and system 304 in this work) and a second pair

from one confirmed image and the MUSE narrow-band data cube (system 305). These new constraints bring our final data set to a total of 67 images for 25 systems that are indicated with red circles in Fig. 1. The spectra and optical images of these new systems are presented in Appendix D as well as their coordinates and redshifts. The method used to find these new systems is detailed in Richard et al. (2021).

## 5.3 Likelihoods

In our optimization process, we simultaneously constrain the mass distribution with the strong lensing and the X-ray surface brightness. We obtain the combined likelihood, $\mathcal{L}_{tot}$, by assuming that both observables are independent, which leads to the following expression:

$$\mathcal{L}_{\text{tot}} = \mathcal{L}_{\text{SL}} \times \mathcal{L}_{\text{X-ray}} \qquad (31)$$

where $\mathcal{L}_{SL}$ and $\mathcal{L}_{X-ray}$ are the likelihoods for the strong lensing only and the X-ray only, respectively. For the latter one, we use a Poisson-Gamma mixture likelihood. For the former, we use a Gaussian likelihood slightly modified in comparison to previous works based on LENSTOOL (Lagattuta et al. 2017; Mahler et al. 2018).

### 5.3.1 Constructing the X-ray likelihood function

As our modelling of the X-ray gas mass solely uses a couple of dPIE potentials, we can only reproduce the large-scale distribution, excluding all microphysical gas phenomena. These small-scale events in the plasma are linked to surface brightness fluctuations (Eckert et al. 2017), which we take into account by adding an uncertainty on the mean, $\lambda_i$, of the Poisson distribution of the observed photon counts in the $i^{th}$ bin. We assume that $\lambda_i$ follows a Gamma distribution, $\Gamma(\mu_i, \sigma_X^2)$, parametrized by its expected value $\mu_i$ and its variance $\sigma_X^2$ instead of the usual shape and scale parameters. $\mu_i$ is our predicted count model excluding the microphysics, whereas $\sigma_X^2$ represents the induced uncertainty assumed to be the same for all bins. As we only intend to model the count model uncertainty statistically, we marginalize on the realization of the Gamma distribution, which leads to the following likelihood in the $i^{th}$ bin to observe $k_i$ count:

$$\mathcal{L}_{\text{X-ray},i} = \int_{0}^{\infty} \mathbb{P}_{i}(X_{\lambda} = k) f_{\Gamma}(\lambda, \mu_{i}, \sigma_{X}^{2}) d\lambda \qquad (32)$$

where $\mathbb{P}_i$ is the probability mass function of the Poisson distribution and $f_\Gamma$ is the probability density function of $\Gamma(\mu_i, \sigma_X^2)$. This compound probability distribution is known to be a continuous extension of a negative binomial distribution $\mathbb{N}\mathbb{B}(r = \mu_i^2/\sigma_X^2, p = \mu_i/(\mu_i + \sigma_X^2))$ with the following expression:

$$\mathcal{L}_{\text{X-ray},i} = \frac{\Gamma(k_i + r)}{k_i! \Gamma(r)} p^r (1 - p)^{k_i} \qquad (33)$$

The expected value of $\mathbb{N}\mathbb{B}$ is $\mu_i$, and the variance is $\mu_i + \sigma_X^2$, which implies an overdispersion of the Poisson model in the case of $\sigma_X > 0$. We fall back to a Poisson modelling when the uncertainty is zero, with the expected value being equal to the variance. We can now obtain our final likelihood by injecting $r$ and $p$ expressions in the previous definitions and assuming that each bin is statistically independent of the others:

$$\mathcal{L}_{\text{X-ray}} = \prod_{i} \frac{\Gamma\left(k_i + \frac{\mu_i^2}{\sigma_X^2}\right)}{k_i! \Gamma\left(\frac{\mu_i^2}{\sigma_X^2}\right)} \left(\frac{\mu_i}{\mu_i + \sigma_X^2}\right)^{\frac{\mu_i^2}{\sigma_X^2}} \left(\frac{\sigma_X^2}{\mu_i + \sigma_X^2}\right)^{k_i} \qquad (34)$$
Preliminary tests show that results obtained with this likelihood are consistent with the ones given with a Poisson likelihood without uncertainties on its expected value. Indeed, the best-fitting count model obtained through an X-ray-only fit has a relative difference of less than 1 per cent on average with a maximum of less than 10 per cent.

### 5.3.2 Constructing the strong lensing likelihood function

The expression of $\mathcal{L}_{SL}$ is the following (Lagattuta et al. 2017; Mahler et al. 2018):

$$ \mathcal{L}_{SL} = \sum_{j}^{N_{\text{sys}}} \frac{1}{\prod_{i}^{N_{\text{im},j}} \sigma_{ij} \sqrt{2\pi}} \exp(-\chi_j^2/2) \quad (35) $$

with $N_{\text{sys}}$ the total number of multiply imaged systems, $N_{\text{im},j}$, the number of images in the $j^{th}$ system, and $\sigma_{ij}$, the observational error on the $i^{th}$ images of the $j^{th}$ system. $\chi_j^2$ is the $\chi^2$ statistics associated with the $j^{th}$ system and read as follows:

$$ \chi_j^2 = \sum_i^{N_{\text{im},j}} \frac{\left| \bar{\theta}_{i,j}^{\text{obs}} - \bar{\theta}_{i,j}^{\text{pred}} \right|^2}{\sigma_{ij}^2} \quad (36) $$

where $\bar{\theta}_{i,j}^{\text{obs}}$ is the observed multiple images positions and $\bar{\theta}_{i,j}^{\text{pred}}$ is the predicted ones. The difference between our likelihood and previous works based on the LENSTOOL software is hidden in the $\bar{\theta}_{i,j}^{\text{pred}}$ term and specifically in the position of the multiply imaged system source. Previously, the intrinsic source position used to solve the lens equation was chosen to be the barycentre or amplification weighted barycentre of the sources associated with each image individually. This solution is convenient as it does not add extra calculations to the computationally expensive part of this process, which is solving the lens equation, but it introduces a bias to the modelling workflow. Indeed, choosing a different intrinsic source position that will better suit the multiply imaged systems directly modifies the likelihood value by reducing it drastically in some specific cases, for example, when an image is near its associated critical lines. This difference will also modify the Bayesian criteria used to discriminate models. Thus, we add an extra computational step, which consists of optimizing the intrinsic source position to reduce the associated positional error on each image of the system. We do this operation at a fixed model for each step of the nested sampling run, and we solve the non-linear least-square problem associated with the Levenberg-Marquardt algorithm implemented in the GNU Scientific Library. The starting point of the Newton method is the previously defined barycentre.

## 6 MASS RECONSTRUCION OF ABELL S1063

Multiple models are possible with our new method with different complexity. Hence, we describe how discrimination between models is performed in Section 6.1. The reproduction of the observational constraints and the obtained mass distribution for each cluster component is presented in Sections 6.2 and 6.3, respectively.

### 6.1 Model discrimination

Given that our modelling method relies on X-ray and lensing data information of the cluster, we have to adapt the model discrimination to include this new type of non-homogeneous constraints. We choose

to follow a similar two-step process as described in Beauchesne et al. (2021). The first step is modified to include the X-ray part, while the second one includes some changes related to the physically motivated priors on the cluster member parameters as explained in Section 4.4. The original aim of this approach was to solve convergence issues when adding the B-spline surfaces as the number of free parameters increased by a significant number. Indeed, the parametric part of the model struggles to converge to high-likelihood areas. The first step aims at discriminating between models on the parametric side only, while the second one allows us to define the best number of B-spline basis functions.

#### 6.1.1 Parametric model discrimination

To discriminate between parametric models, we treat separately the gas distribution constrained by the X-rays and the rest of the model that relies solely on lensing data. We use a disjoint process as, in our case, the numerical values of the two likelihoods have three to four orders of magnitude difference in favour of the X-ray one. As the gas distribution only represents a small part of the overall mass budget (we know a priori that the gas mass is roughly of the order of ~10 per cent of the total mass), we assume that by discriminating the lensing-only part without the gas mass will not introduce a significant bias.

For the lensing-only part, we rely on previous studies of AS1063, which mostly agree on a model with two large-scale dPIEs as detailed in Section 4.1 (Bonamigo et al. 2018; Bergamini et al. 2019; Limousin, Beauchesne & Jullo 2022). In the case of less studied galaxy clusters where previous models may not be available, we would have based our choice on Bayesian criteria such as the Bayesian evidence (i.e. marginal likelihood) to select which model best describes the data.

For the gas distribution, we use a different approach as the statistical knowledge of the constraints is well known contrary to the lensing ones. Actually, the observational errors on the position of the multiple images integrate both errors due to physical phenomena such as LOS perturbers and/or systematic limitations due to the parametric approach, for example. Hence, we developed a measure of the goodness of fit adapted to our Poisson-Gamma mixture likelihood with a similar approach as in Kaastra (2017). The only differences are the statistics used as well as a numerical approach of the measure instead of an analytical one. We consider that a gas distribution model is explaining well the data when the likelihood of the observations is at least included in the "5σ"CI of the expected value of the X-ray likelihood for the considered model.

We obtain the distribution of the X-ray likelihood through a Monte Carlo approach, where we sample the Poisson-Gamma mixture model based on the parameter of the best fit. We proceed by adding one dPIE at a time, and we stop at the model with less complexity that satisfies the previous condition (i.e. included in the "5σ"CI of the expected value of the likelihood). In our case, we end up with a model composed of three dPIEs as in Bonamigo et al. (2018), but our selection method assures that we reach the right complexity for the data we have. This is only possible through the use of the Poisson-Gamma mixture, as the same model underfits the data and is outside the "5σ"CI with a Poisson statistic. The bound of the "5σ"CI is at 1.5σ from the obtained best-fitting likelihood, but as we are not considering a Gaussian distribution, the "nσ"CI does not scale linearly with the standard deviation.

When both parametric components are defined, we merge them and proceed to a sampling run that is used as a starting point for
Preliminary tests show that results obtained with this likelihood are consistent with the ones given with a Poisson likelihood without uncertainties on its expected value. Indeed, the best-fitting count model obtained through an X-ray-only fit has a relative difference of less than 1 per cent on average with a maximum of less than 10 per cent.

### 5.3.2 Constructing the strong lensing likelihood function

The expression of $\mathcal{L}_{SL}$ is the following (Lagattuta et al. 2017; Mahler et al. 2018):

$$ \mathcal{L}_{SL} = \sum_{j}^{N_{sys}} \frac{1}{\prod_{i}^{N_{im,j}} \sigma_{ij} \sqrt{2\pi}} \exp(-\chi_j^2/2) \quad (35) $$

with $N_{sys}$ the total number of multiply imaged systems, $N_{im,j}$, the number of images in the $j^{th}$ system, and $\sigma_{ij}$, the observational error on the $i^{th}$ images of the $j^{th}$ system. $\chi_j^2$ is the $\chi^2$ statistics associated with the $j^{th}$ system and read as follows:

$$ \chi_j^2 = \sum_i^{N_{im,j}} \frac{\left\| \vec{\theta}_{i,j}^{\text{obs}} - \vec{\theta}_{i,j}^{\text{pred}} \right\|^2}{\sigma_{ij}^2} \quad (36) $$

where $\vec{\theta}_{i,j}^{\text{obs}}$ is the observed multiple images positions and $\vec{\theta}_{i,j}^{\text{pred}}$ is the predicted ones. The difference between our likelihood and previous works based on the LENSTOOL software is hidden in the $\vec{\theta}_{i,j}^{\text{pred}}$ term and specifically in the position of the multiply imaged system source. Previously, the intrinsic source position used to solve the lens equation was chosen to be the barycentre or amplification weighted barycentre of the sources associated with each image individually. This solution is convenient as it does not add extra calculations to the computationally expensive part of this process, which is solving the lens equation, but it introduces a bias to the modelling workflow. Indeed, choosing a different intrinsic source position that will better suit the multiply imaged systems directly modifies the likelihood value by reducing it drastically in some specific cases, for example, when an image is near its associated critical lines. This difference will also modify the Bayesian criteria used to discriminate models. Thus, we add an extra computational step, which consists of optimizing the intrinsic source position to reduce the associated positional error on each image of the system. We do this operation at a fixed model for each step of the nested sampling run, and we solve the non-linear least-square problem associated with the Levenberg-Marquardt algorithm implemented in the GNU Scientific Library. The starting point of the Newton method is the previously defined barycentre.

## 6 MASS RECONSTRUCION OF ABELL S1063

Multiple models are possible with our new method with different complexity. Hence, we describe how discrimination between models is performed in Section 6.1. The reproduction of the observational constraints and the obtained mass distribution for each cluster component is presented in Sections 6.2 and 6.3, respectively.

### 6.1 Model discrimination

Given that our modelling method relies on X-ray and lensing data information of the cluster, we have to adapt the model discrimination to include this new type of non-homogeneous constraints. We choose

to follow a similar two-step process as described in Beauchesne et al. (2021). The first step is modified to include the X-ray part, while the second one includes some changes related to the physically motivated priors on the cluster member parameters as explained in Section 4.4. The original aim of this approach was to solve convergence issues when adding the B-spline surfaces as the number of free parameters increased by a significant number. Indeed, the parametric part of the model struggles to converge to high-likelihood areas. The first step aims at discriminating between models on the parametric side only, while the second one allows us to define the best number of B-spline basis functions.

#### 6.1.1 Parametric model discrimination

To discriminate between parametric models, we treat separately the gas distribution constrained by the X-rays and the rest of the model that relies solely on lensing data. We use a disjoint process as, in our case, the numerical values of the two likelihoods have three to four orders of magnitude difference in favour of the X-ray one. As the gas distribution only represents a small part of the overall mass budget (we know a priori that the gas mass is roughly of the order of ~ 10 per cent of the total mass), we assume that by discriminating the lensing-only part without the gas mass will not introduce a significant bias.

For the lensing-only part, we rely on previous studies of AS1063, which mostly agree on a model with two large-scale dPIEs as detailed in Section 4.1 (Bonamigo et al. 2018; Bergamini et al. 2019; Limousin, Beauchesne & Jullo 2022). In the case of less studied galaxy clusters where previous models may not be available, we would have based our choice on Bayesian criteria such as the Bayesian evidence (i.e. marginal likelihood) to select which model best describes the data.

For the gas distribution, we use a different approach as the statistical knowledge of the constraints is well known contrary to the lensing ones. Actually, the observational errors on the position of the multiple images integrate both errors due to physical phenomena such as LOS perturbers and/or systematic limitations due to the parametric approach, for example. Hence, we developed a measure of the goodness of fit adapted to our Poisson-Gamma mixture likelihood with a similar approach as in Kaastra (2017). The only differences are the statistics used as well as a numerical approach of the measure instead of an analytical one. We consider that a gas distribution model is explaining well the data when the likelihood of the observations is at least included in the "5σ"CI of the expected value of the X-ray likelihood for the considered model.

We obtain the distribution of the X-ray likelihood through a Monte Carlo approach, where we sample the Poisson-Gamma mixture model based on the parameter of the best fit. We proceed by adding one dPIE at a time, and we stop at the model with less complexity that satisfies the previous condition (i.e. included in the "5σ"CI of the expected value of the likelihood). In our case, we end up with a model composed of three dPIEs as in Bonamigo et al. (2018), but our selection method assures that we reach the right complexity for the data we have. This is only possible through the use of the Poisson-Gamma mixture, as the same model underfits the data and is outside the "5σ"CI with a Poisson statistic. The bound of the "5σ"CI is at 1.5σ from the obtained best-fitting likelihood, but as we are not considering a Gaussian distribution, the "nσ"CI does not scale linearly with the standard deviation.

When both parametric components are defined, we merge them and proceed to a sampling run that is used as a starting point for
<figure><img src="image_8.png" /><figcaption>**Figure 5.** From left to right: Maps of the observed X-ray counts, the counts model from our best-fitting model, and the residuals between the two previous maps (i.e. observation minus model). In both the observation and the model, the X-ray point sources have been masked.</figcaption></figure>

the next step. This run, and the following, are performed with the dynamic nested sampling algorithm MLFRIENDS (Buchner et al. 2014; Buchner 2019) implemented in the python package ULTRANEST<fnref n="8" /> (Buchner 2021).

### 6.1.2 Perturbative model discrimination

All variants of perturbative modelling are created in the same way as described in Section 4.4. The discrimination is done with the Bayesian evidence, $log E$ (i.e. marginal likelihood), provided by the nested sampling engine. We select the model that gives the best $log E$ in the case of a parabolic-like shape, or the first modelling reaching the $log E$ plateau as seen in Beauchesne et al. (2021). For all the runs, including a B-spline perturbation, we use an identical observational error estimate for all multiple images as 0.2 arcsec.

We try models with a mesh of 3 × 3, 4 × 4, 5 × 5, and 6 × 6 of B-spline basis. The associated $log E$ shows that B-spline modellings reach a plateau starting at 4 × 4 with $log E = -33309.71$. 5 × 5 and 6 × 6 have values close to the previous one with a $log E$ of -33310.04 and -33304.65, respectively. Models with 3 × 3 mesh of B-spline have a significantly lower marginal log-likelihood with a value of -33331.95. As in Beauchesne et al. (2021), we consider that a model is better than the other if it is more supported by $log E$ at a difference equivalent to $5\sigma$ (i.e. a difference of 12.5 between $log E$ estimations). Thus, the B-spline variant with a mesh of 4 × 4 function is the best according to this selection procedure. The associated best-fitting model has a Root Mean Squared (RMS) error of 0.35 arcsec. Therefore, we construct a new sampling run with this error as the observational one to account for the model systematic uncertainty.

## 6.2 Reproduction of the observational constraints

Regarding the reproduction of the lensing constraints, our modelling reaches an RMS of 0.42 arcsec for the best-fitting model. In comparison to models by B19 and G22, we obtain a better fit of the constraints as these previous works had RMSs of 0.55 and 0.60 arcsec, respectively. The main explanation for this difference is the inclusion of the B-spline surfaces which increase the flexibility of the mass reconstruction. Indeed, by excluding the B-spline surfaces, the RMS of our model is 0.70 arcsec. Note that we include more multiple images than B19 and G22. In particular, the triple system found by Vanzella et al. (2016; i.e. system 11 in our ID system), has been hard to model accurately (Caminha et al. 2016). Our new likelihood scheme has been able to improve the modelling with the inclusion of

system 11, as the barycentre approximation for the source position can struggle with non-linearity due to the presence of critical lines near the images. Our result is also close to the RMS of 0.36 arcsec reported by Limousin, Beauchesne & Jullo (2022) that also use a surface of B-spline but with a mesh of 5 × 5 basis functions, and the same set of constraints as B19 and G22. Moving away from models produced with Lenstool, Raney, Keeton & Brennan (2020) report an RMS of 0.34 arcsec with the inclusion of LOS galaxies. Adding these galaxies does not seem to particularly improve the overall fit as detailed in their study, in particular, system 11 mentioned above that was removed due to potential bias from a LOS galaxy. Previous studies (D’Aloisio, Natarajan & Shapiro 2014; Gilman, Birrer & Treu 2020) tend to favour that LOS structure does not perturb mass models at discernable levels given the current data sets.

Regarding the RMS of our model without B-spline (i.e. 0.70 arcsec), it is higher than models that include system 11 made with a similar method as the Cluster As TelescopeS team<fnref n="9" /> (PIs: Natarajan & Kneib) as part of the Frontier Field lensing modelling effort. Indeed, the RMS obtained for this modelling is of 0.49 arcsec, but it includes an external shear component with a high shear amplitude (i.e. $\gamma = 0.1039$ for the best-fitting model) that seems unlikely to be the sole product of some massive clumps in the cluster outskirts. Removing that component leads to a significant increase in the error of the predicted position of the multiple images, reaching values similar to 0.80 arcsec. This increase is particularly important for system 11 with a jump from 0.83 arcsec to 1.97 arcsec. However, one of our new multiply imaged systems (e.g. system 304 in Appendix D) broke this trend and prevented the external shear from having such undue influence on system 11, and thus on the global RMS. It is partially a reason for the non-inclusion of that extra component, as it did not enhance the model with or without B-spline.

As we are constraining our model with X-ray data in addition to the lensing information, we have to estimate the ability of our modelling to recover the X-ray surface brightness map. The quality of the reconstruction of the X-ray observations can be seen in Fig. 5, which shows, from the left to the right, the observed photon counts with point sources masked, the count model, and the residual. Thanks to the use of multiple dPIE, our gas model has been able to capture the main asymmetry observed in the count map and provide a good fit of the data. However, a banana-shaped pattern of count underestimation can be observed in pink in the residual map below the cluster centre, which indicates the limitation due to the dPIE scale. Such a pattern, combined with the overestimation above in green, can indicate the presence of gas sloshing from the NE clump to the main one (Patterno-

<footnote n="8">⁸https://johannesbuchner.github.io/UltraNest/</footnote>

<footnote n="9">⁹https://archive.stsci.edu/pub/hlsp/frontier/abells1063/models/cats/</footnote>
<figure><img src="image_8.png" /><figcaption>**Figure 5.** From left to right: Maps of the observed X-ray counts, the counts model from our best-fitting model, and the residuals between the two previous maps (i.e. observation minus model). In both the observation and the model, the X-ray point sources have been masked.</figcaption></figure>

the next step. This run, and the following, are performed with the dynamic nested sampling algorithm MLFRIENDS (Buchner et al. 2014; Buchner 2019) implemented in the python package ULTRANEST<fnref n="8" /> (Buchner 2021).

### 6.1.2 Perturbative model discrimination

All variants of perturbative modelling are created in the same way as described in Section 4.4. The discrimination is done with the Bayesian evidence, $log E$ (i.e. marginal likelihood), provided by the nested sampling engine. We select the model that gives the best $log E$ in the case of a parabolic-like shape, or the first modelling reaching the $log E$ plateau as seen in Beauchesne et al. (2021). For all the runs, including a B-spline perturbation, we use an identical observational error estimate for all multiple images as 0.2 arcsec.

We try models with a mesh of 3 × 3, 4 × 4, 5 × 5, and 6 × 6 of B-spline basis. The associated $log E$ shows that B-spline modellings reach a plateau starting at 4 × 4 with $log E = -33309.71$. 5 × 5 and 6 × 6 have values close to the previous one with a $log E$ of -33310.04 and -33304.65, respectively. Models with 3 × 3 mesh of B-spline have a significantly lower marginal log-likelihood with a value of -33331.95. As in Beauchesne et al. (2021), we consider that a model is better than the other if it is more supported by $log E$ at a difference equivalent to $5\sigma$ (i.e. a difference of 12.5 between $log E$ estimations). Thus, the B-spline variant with a mesh of 4 × 4 function is the best according to this selection procedure. The associated best-fitting model has a Root Mean Squared (RMS) error of 0.35 arcsec. Therefore, we construct a new sampling run with this error as the observational one to account for the model systematic uncertainty.

## 6.2 Reproduction of the observational constraints

Regarding the reproduction of the lensing constraints, our modelling reaches an RMS of 0.42 arcsec for the best-fitting model. In comparison to models by B19 and G22, we obtain a better fit of the constraints as these previous works had RMSs of 0.55 and 0.60 arcsec, respectively. The main explanation for this difference is the inclusion of the B-spline surfaces which increase the flexibility of the mass reconstruction. Indeed, by excluding the B-spline surfaces, the RMS of our model is 0.70 arcsec. Note that we include more multiple images than B19 and G22. In particular, the triple system found by Vanzella et al. (2016; i.e. system 11 in our ID system), has been hard to model accurately (Caminha et al. 2016). Our new likelihood scheme has been able to improve the modelling with the inclusion of

system 11, as the barycentre approximation for the source position can struggle with non-linearity due to the presence of critical lines near the images. Our result is also close to the RMS of 0.36 arcsec reported by Limousin, Beauchesne & Jullo (2022) that also use a surface of B-spline but with a mesh of 5 × 5 basis functions, and the same set of constraints as B19 and G22. Moving away from models produced with Lenstool, Raney, Keeton & Brennan (2020) report an RMS of 0.34 arcsec with the inclusion of LOS galaxies. Adding these galaxies does not seem to particularly improve the overall fit as detailed in their study, in particular, system 11 mentioned above that was removed due to potential bias from a LOS galaxy. Previous studies (D’Aloisio, Natarajan & Shapiro 2014; Gilman, Birrer & Treu 2020) tend to favour that LOS structure does not perturb mass models at discernable levels given the current data sets.

Regarding the RMS of our model without B-spline (i.e. 0.70 arcsec), it is higher than models that include system 11 made with a similar method as the Cluster As TelescopeS team<fnref n="9" /> (PIs: Natarajan & Kneib) as part of the Frontier Field lensing modelling effort. Indeed, the RMS obtained for this modelling is of 0.49 arcsec, but it includes an external shear component with a high shear amplitude (i.e. $\gamma = 0.1039$ for the best-fitting model) that seems unlikely to be the sole product of some massive clumps in the cluster outskirts. Removing that component leads to a significant increase in the error of the predicted position of the multiple images, reaching values similar to 0.80 arcsec. This increase is particularly important for system 11 with a jump from 0.83 arcsec to 1.97 arcsec. However, one of our new multiply imaged systems (e.g. system 304 in Appendix D) broke this trend and prevented the external shear from having such undue influence on system 11, and thus on the global RMS. It is partially a reason for the non-inclusion of that extra component, as it did not enhance the model with or without B-spline.

As we are constraining our model with X-ray data in addition to the lensing information, we have to estimate the ability of our modelling to recover the X-ray surface brightness map. The quality of the reconstruction of the X-ray observations can be seen in Fig. 5, which shows, from the left to the right, the observed photon counts with point sources masked, the count model, and the residual. Thanks to the use of multiple dPIE, our gas model has been able to capture the main asymmetry observed in the count map and provide a good fit of the data. However, a banana-shaped pattern of count underestimation can be observed in pink in the residual map below the cluster centre, which indicates the limitation due to the dPIE scale. Such a pattern, combined with the overestimation above in green, can indicate the presence of gas sloshing from the NE clump to the main one (Pатерно-

<footnote n="8">⁸https://johannesbuchner.github.io/UltraNest/</footnote>

<footnote n="9">⁹https://archive.stsci.edu/pub/hlsp/frontier/abells1063/models/cats/</footnote>
<figure><img src="image_9.png" /><figcaption>**Figure 6.** Map of the absolute relative differences between the X-ray loglikelihood computed with a $\sigma_X = 0.48$ counts and the best-fitting value. The white areas are masked point-like sources.</figcaption></figure>

Mahler et al. 2013). Results from Bonamigo et al. (2018) also present the same artefact in their residuals (see Fig. 1 in their analysis).

A special feature in our approach is the inclusion of an uncertainty on our count model, to account for small-scale fluctuations in surface brightness. In contrast to previous similar studies (Bonamigo et al. 2017, 2018), we choose not to use a Poisson statistic to fit the X-ray photon counts. We instead use a Gamma-Poisson mixture model, where the uncertainty, $\sigma_X$, is explicitly included with the Gamma distribution. During the sampling of the joint likelihood, we obtain an estimation of $\sigma_X$ through the standard deviation of the Gamma statistics, which is estimated to be $0.263_{-0.028}^{+0.027}$ counts. If, theoretically, we obtain information on small-scale phenomena such as the microphysics of the gas, in principle, we could assess the robustness of the error estimation. Our goodness-of-fit procedure is able to give us insights into how well our model can describe the observations. Indeed, this Monte Carlo method samples the possible count realization for a given model. Depending on where the likelihood for the best-fitting model is, we may know how much of the observations are plausible according to it. Our results show that this likelihood is between the limits of the "$3\sigma$"CI and "$5\sigma$"CI on the underfitting side of the previously mentioned distribution. Hence, if the real observations are plausible according to the best-fitting model, they are still unlikely.

However, $\sigma_X$ values have a significative influence on the distribution of the plausible X-ray likelihood values, which consequently changes the position of the best-fitting model likelihood in the different "$n\sigma$"CI. Increasing $\sigma_X$ to values between 0.45 and 0.50 photon counts moves the considered likelihood to be close to the expected value, inside the "$1\sigma$"CI. It shows that our method seems to underestimate $\sigma_X$, in the sense that we are not obtaining the model that is the best at explaining the observations. In particular, Fig. 6 presents the absolute relative differences of the log-likelihood per pixel between the best-fitting $\sigma_X$ estimation, and a value of 0.48 counts. As we can see, the discrepancies between both uncertainties are mainly in the low counts area where micro-physics events from say turbulence are not expected to occur, or at least, with a lower amplitude. This highlights the simplicity of our scheme as the central area, and the outskirt would be better treated differently in agreement with the amplitude of the different processes in the gas. Such improvements are needed to extract information on small-scale

phenomena with satisfactory reliability in the estimation of their amplitude.

## 6.3 The mass distribution

Fig. 7 presents the total mass distribution of AS1063 obtained from our best-fitting model. The black and crimson dashed lines indicate the contours of the gas distribution and the smooth mass components, respectively. The statistics derived from the posterior distribution of dPIEs parameters are shown in Appendix F. The dominant component of the total mass is the DM. We can see clearly that the inferred ellipticity of the DM component is different from that of the gas distribution. The latter is rounder on cluster scales. However, the discrepancy between the two components changes towards the centre of the cluster as both ellipticities become more similar to that of the BCG. We also note that the asymmetry in direction to the NE clump is less prominent in the total mass than in the gas and that the overall mass distribution is mainly uni-modal. In particular, our finding does not support the hypothesis of Gómez et al. (2012) of a major merger event between two similar-sized clusters, like in the bullet cluster (Clowe et al. 2006). More likely, per our model, the NE clump was possibly a smaller object, such as a galaxy group infalling into AS1063. We assess this scenario in more detail with the different products and by-products of our new method in Section 8.2.

Notably, the addition of the B-spline surface deforms the elliptical symmetry of the main halo, as shown by the crimson dashed contours. This extends the mass towards the NE clumps, mimicking the effect of another dPIE as found by previous studies (Caminha et al. 2016; Limousin, Beauchesne & Jullo 2022). The smooth mass components also present other deviations from the main halo symmetry, as shown by the fluctuation of the crimson contours around the main dPIE halo iso-masses. These fluctuations are only limited to the core due to the B-spline amplitudes fading away, as highlighted by the farthest contours from the cluster centre. Thus, it is hard to assess if these variations are linked to the mass distribution patterns on larger scales, such as local galaxy overdensities or mass haloes on the cluster outskirt. Further constraints from weak-lensing analysis from the BUFFALO program will expand the area constrained by lensing and, thus, the surface where the B-spline can be defined.

### 6.3.1 Cluster member masses

Fig. 8 presents the central velocity dispersion [i.e. as in the Lenstool code, $\sigma_0 = \frac{2}{3}\sigma_{0,\text{th}}$ from Limousin, Kneib & Natarajan (2005) definition] of the dPIE potentials as a function of their dPIE total mass from this work (green) and from B19 (orange). Thanks to our joint calibration process between the Faber & Jackson and the fundamental plane, all of the modelled galaxies follow the same slope, providing consistency amongst all the cluster members, contrary to G22. Hence, our implementation has been successful at obtaining a scatter in the $\sigma_0$-mass relation without introducing biases for cluster members modelled with different observables. As expected from the calibration in Section 3.2.2, our slope is steeper than the one obtained by B19, even after the addition of the lensing constraints. The parameter estimations obtained with the joint lensing and X-ray fit are listed in Table 2. As the median values are slightly changed, and the statistical errors have been reduced, the new medians are in agreement at ~1σ with the calibration. Thus, the lensing constraints have been successful in improving the accuracy of the fundamental plane estimation. The total mass contributed by cluster members without the BCG is $1.6 \pm 0.4 \ 10^{12} \text{M}_\odot$ at 600 kpc from the BCG.
<figure><img src="image_9.png" /><figcaption>**Figure 6.** Map of the absolute relative differences between the X-ray loglikelihood computed with a $\sigma_X = 0.48$ counts and the best-fitting value. The white areas are masked point-like sources.</figcaption></figure>

Mahler et al. 2013). Results from Bonamigo et al. (2018) also present the same artefact in their residuals (see Fig. 1 in their analysis).

A special feature in our approach is the inclusion of an uncertainty on our count model, to account for small-scale fluctuations in surface brightness. In contrast to previous similar studies (Bonamigo et al. 2017, 2018), we choose not to use a Poisson statistic to fit the X-ray photon counts. We instead use a Gamma-Poisson mixture model, where the uncertainty, $\sigma_X$, is explicitly included with the Gamma distribution. During the sampling of the joint likelihood, we obtain an estimation of $\sigma_X$ through the standard deviation of the Gamma statistics, which is estimated to be $0.263_{-0.028}^{+0.027}$ counts. If, theoretically, we obtain information on small-scale phenomena such as the microphysics of the gas, in principle, we could assess the robustness of the error estimation. Our goodness-of-fit procedure is able to give us insights into how well our model can describe the observations. Indeed, this Monte Carlo method samples the possible count realization for a given model. Depending on where the likelihood for the best-fitting model is, we may know how much of the observations are plausible according to it. Our results show that this likelihood is between the limits of the "$3\sigma$"CI and "$5\sigma$"CI on the underfitting side of the previously mentioned distribution. Hence, if the real observations are plausible according to the best-fitting model, they are still unlikely.

However, $\sigma_X$ values have a significative influence on the distribution of the plausible X-ray likelihood values, which consequently changes the position of the best-fitting model likelihood in the different "$n\sigma$"CI. Increasing $\sigma_X$ to values between 0.45 and 0.50 photon counts moves the considered likelihood to be close to the expected value, inside the "$1\sigma$"CI. It shows that our method seems to underestimate $\sigma_X$, in the sense that we are not obtaining the model that is the best at explaining the observations. In particular, Fig. 6 presents the absolute relative differences of the log-likelihood per pixel between the best-fitting $\sigma_X$ estimation, and a value of 0.48 counts. As we can see, the discrepancies between both uncertainties are mainly in the low counts area where micro-physics events from say turbulence are not expected to occur, or at least, with a lower amplitude. This highlights the simplicity of our scheme as the central area, and the outskirt would be better treated differently in agreement with the amplitude of the different processes in the gas. Such improvements are needed to extract information on small-scale

phenomena with satisfactory reliability in the estimation of their amplitude.

## 6.3 The mass distribution

Fig. 7 presents the total mass distribution of AS1063 obtained from our best-fitting model. The black and crimson dashed lines indicate the contours of the gas distribution and the smooth mass components, respectively. The statistics derived from the posterior distribution of dPIEs parameters are shown in Appendix F. The dominant component of the total mass is the DM. We can see clearly that the inferred ellipticity of the DM component is different from that of the gas distribution. The latter is rounder on cluster scales. However, the discrepancy between the two components changes towards the centre of the cluster as both ellipticities become more similar to that of the BCG. We also note that the asymmetry in direction to the NE clump is less prominent in the total mass than in the gas and that the overall mass distribution is mainly uni-modal. In particular, our finding does not support the hypothesis of Gómez et al. (2012) of a major merger event between two similar-sized clusters, like in the bullet cluster (Clowe et al. 2006). More likely, per our model, the NE clump was possibly a smaller object, such as a galaxy group infalling into AS1063. We assess this scenario in more detail with the different products and by-products of our new method in Section 8.2.

Notably, the addition of the B-spline surface deforms the elliptical symmetry of the main halo, as shown by the crimson dashed contours. This extends the mass towards the NE clumps, mimicking the effect of another dPIE as found by previous studies (Caminha et al. 2016; Limousin, Beauchesne & Jullo 2022). The smooth mass components also present other deviations from the main halo symmetry, as shown by the fluctuation of the crimson contours around the main dPIE halo iso-masses. These fluctuations are only limited to the core due to the B-spline amplitudes fading away, as highlighted by the farthest contours from the cluster centre. Thus, it is hard to assess if these variations are linked to the mass distribution patterns on larger scales, such as local galaxy overdensities or mass haloes on the cluster outskirt. Further constraints from weak-lensing analysis from the BUFFALO program will expand the area constrained by lensing and, thus, the surface where the B-spline can be defined.

### 6.3.1 Cluster member masses

Fig. 8 presents the central velocity dispersion [i.e. as in the Lenstool code, $\sigma_0 = \frac{2}{3}\sigma_{0,\text{th}}$ from Limousin, Kneib & Natarajan (2005) definition] of the dPIE potentials as a function of their dPIE total mass from this work (green) and from B19 (orange). Thanks to our joint calibration process between the Faber & Jackson and the fundamental plane, all of the modelled galaxies follow the same slope, providing consistency amongst all the cluster members, contrary to G22. Hence, our implementation has been successful at obtaining a scatter in the $\sigma_0$-mass relation without introducing biases for cluster members modelled with different observables. As expected from the calibration in Section 3.2.2, our slope is steeper than the one obtained by B19, even after the addition of the lensing constraints. The parameter estimations obtained with the joint lensing and X-ray fit are listed in Table 2. As the median values are slightly changed, and the statistical errors have been reduced, the new medians are in agreement at ~1$\sigma$ with the calibration. Thus, the lensing constraints have been successful in improving the accuracy of the fundamental plane estimation. The total mass contributed by cluster members without the BCG is $1.6 \pm 0.4 \ 10^{12} \text{M}_\odot$ at 600 kpc from the BCG.
<figure><img src="image_10.png" /><figcaption>**Figure 7.** Map of the projected surface mass density from our best-fitting model where the dashed black and crimson contours represent the surface mass density of the X-ray emitting gas only and the smooth mass component (i.e. all mass components including the BCG but without the other cluster members), respectively.</figcaption></figure>

<figure><img src="image_11.png" /><figcaption>**Figure 8.** Central velocity dispersion of cluster members as a function of their dPIE total masses from our model (green). We have also plotted in orange the model from B19 as comparison. The error bars represent the "1σ"CI of the posterior distribution. The slope obtained here is steeper than the one obtained in B19 (see Section 3.2.2).</figcaption></figure>

### 6.3.2 The gas distribution

To compare the obtained gas distribution with other methods, we compute projected gas mass profiles shown in Fig. 9 for several

<table><caption>**Table 2.** Parameters of the Faber & Jackson scaling relation and the fundamental plane of elliptical galaxies obtained after the joint fit from X-ray and lensing constraints.</caption><thead><tr><td>Parameters</td><td>Median</td><td>σ</td></tr></thead><tbody><tr><td>a</td><td>-0.571</td><td>0.026</td></tr><tr><td>b</td><td>0.302</td><td>0.018</td></tr><tr><td>c</td><td>1.478</td><td>0.075</td></tr><tr><td>σ<sup>*</sup><sub>e</sub></td><td>190.5</td><td>5.3</td></tr></tbody></table>

models. To assess the bias produced by the use of a Gamma-Poisson mixture instead of a Poisson statistic, we perform a fit of the gas using only the latter. As shown by the black and orange profiles, the use of the Poisson statistic has only a minimal impact on the overall fit as both curves are almost indiscernible from each other. Hence, it allows us to have a more satisfying statistical explanation of the considered X-ray constraints and provides the same gas distribution as the more standard Poisson likelihood.

We also compare our mass profiles with the gas included in the B19 model (red solid line) as well as the profile obtained with the PYPROFITT multiscale deprojection method (Eckert et al. 2020, green area). The latter method is in good agreement with our results with a mass underestimate of only 3 per cent, but the discrepancy is much larger with the B19 model. Indeed, their model shows an
<figure><img src="image_10.png" /><figcaption>**Figure 7.** Map of the projected surface mass density from our best-fitting model where the dashed black and crimson contours represent the surface mass density of the X-ray emitting gas only and the smooth mass component (i.e. all mass components including the BCG but without the other cluster members), respectively.</figcaption></figure>

<figure><img src="image_11.png" /><figcaption>**Figure 8.** Central velocity dispersion of cluster members as a function of their dPIE total masses from our model (green). We have also plotted in orange the model from B19 as comparison. The error bars represent the "1σ"CI of the posterior distribution. The slope obtained here is steeper than the one obtained in B19 (see Section 3.2.2).</figcaption></figure>

### 6.3.2 The gas distribution

To compare the obtained gas distribution with other methods, we compute projected gas mass profiles shown in Fig. 9 for several

<table><caption>**Table 2.** Parameters of the Faber & Jackson scaling relation and the fundamental plane of elliptical galaxies obtained after the joint fit from X-ray and lensing constraints.</caption><thead><tr><td>Parameters</td><td>Median</td><td>σ</td></tr></thead><tbody><tr><td>a</td><td>-0.571</td><td>0.026</td></tr><tr><td>b</td><td>0.302</td><td>0.018</td></tr><tr><td>c</td><td>1.478</td><td>0.075</td></tr><tr><td>σ<sup>*</sup><sub>e</sub></td><td>190.5</td><td>5.3</td></tr></tbody></table>

models. To assess the bias produced by the use of a Gamma-Poisson mixture instead of a Poisson statistic, we perform a fit of the gas using only the latter. As shown by the black and orange profiles, the use of the Poisson statistic has only a minimal impact on the overall fit as both curves are almost indiscernible from each other. Hence, it allows us to have a more satisfying statistical explanation of the considered X-ray constraints and provides the same gas distribution as the more standard Poisson likelihood.

We also compare our mass profiles with the gas included in the B19 model (red solid line) as well as the profile obtained with the PYPROFITT multiscale deprojection method (Eckert et al. 2020, green area). The latter method is in good agreement with our results with a mass underestimate of only 3 per cent, but the discrepancy is much larger with the B19 model. Indeed, their model shows an
<figure><img src="image_12.png" /><figcaption>**Figure 9.** Projected gas mass density profiles integrated with a disc of increasing radius. The error represents the $3\sigma$ CI of the posterior distribution except for B19 work, where only the best-fitting model is publicly available. Our X-ray constraints and in B19 end at a radius of ~600 kpc.</figcaption></figure>

overestimation of 28 per cent on average compared to our profiles. This model is the only one that has been optimized without the same set of X-ray constraints and parameters of the plasma emission model. Hence, this discrepancy can be partially explained by a difference when computing the effective area of the telescope, for example. We could not reproduce the analysis as Bonamigo et al. (2017, 2018) because the energy at which this effective area has been computed is not publicly available. Hence, we can not assess if there are other reasons for the disagreement between our estimate and B19.

## 7 TESTING OUR RECONSTRUCTION WITH THE MASS DISTRIBUTION OF A MOCK CLUSTER

To estimate the improvement in the mass reconstruction provided by our new method, we create a mock cluster based on our modelling of AS1063 without B-spline perturbations. We start from the best-fitting model and remove all cluster members that are modelled according to the Faber & Jackson law, assuming that the fundamental plane is a good representation of the actual cluster member masses. We use this mass model to predict X-ray counts with the associated Poisson-Gamma mixture while keeping other quantities fixed, such as the background count or the cooling functions. Regarding the lensing constraints, we use the barycentre of the source positions of the actual multiply imaged systems of AS1063 and predict new multiple images with the mock mass. Thanks to this procedure, the realization

of the lensing constraints is fairly similar to reality. We then add a circular Gaussian error of 0.2 arcsec on the positions of the images to account for line-of-sight perturbers and other systematics.

We create four different models to reproduce the mass of the mock cluster while assessing the improvements due to the inclusion of cluster members as well as in the large-scale distribution of the gas. We consider three models: one without the gas and galaxies with only the Faber & Jackson scaling; one including cluster galaxies with the fundamental plane only and a third including only the X-ray emitting gas. And a fourth model that includes both modelling enhancements for cluster member galaxies and the X-ray emitting gas. To account for the combination of the fundamental plane and Faber & Jackson relations, we split the considered galaxies between them to have the same proportions as in the real model. We created the two fundamental plane groups (i.e. galaxies with available $R_e$ or $\sigma_e$ measurements) by taking the most luminous galaxies from the two initial ones, and leaving the remaining elements to the Faber & Jackson scaling. Indeed, inside the MUSE footprint and excluding contamination between cluster members, the faintest ones tend to be modelled with the latter relation as they do not have enough signal for morphological or spectroscopic measurements.

Fig. 10 shows the relative differences of the 2D mass distribution between the simulation and the best fit of each model. From left to right, we have the reference model (i.e. no gas and uncalibrated Faber & Jackson), only gas model, then the fundamental plane on top of the reference and finally, both enhancements. The best-fitting parameters, as well as statistics from the posterior distributions, are available in Appendix F. If we do not consider the discrepancies around the NE halo that are due to the weaker effect of the strong lensing constraints, differences are mainly below 10 per cent in absolute value in the core. As expected, the model with both improvements performs the best with errors mostly below 2 per cent except for the NE halo and some cluster members that switched from the fundamental plane to the Faber & Jackson relation. But the position of the NE halo is mainly constrained by its priors. The real NE position is included in the global posterior distribution of the two models using the fundamental plane scheme. Regarding the two other models, the relative declination of the NE halo is not recovered, with a maximum bias of ~9 arcsec for the model without any improvements. We note that the imperfect modelling scheme of cluster members does not seem to have affected the overall distribution, but only locally, when the discrepancy between relations is higher. We discuss the differences due to the cluster member scaling relations in more detail in Section 7.2.

As shown by the two models where the gas distribution is not considered separately, an underestimation of the mass is introduced

<figure><img src="image_13.png" /><figcaption>**Figure 10.** Maps of the relative differences between the projected surface mass density of each model tested on the mock cluster and the true one. From left to right, the models tested are the reference one, only the gas addition, only the fundamental plane one, and finally, both of them. All maps have the same orientation, with the north in the upward direction and the east directed at the left.</figcaption></figure>
<figure><img src="image_12.png" /><figcaption>**Figure 9.** Projected gas mass density profiles integrated with a disc of increasing radius. The error represents the $3\sigma$ CI of the posterior distribution except for B19 work, where only the best-fitting model is publicly available. Our X-ray constraints and in B19 end at a radius of ~600 kpc.</figcaption></figure>

overestimation of 28 per cent on average compared to our profiles. This model is the only one that has been optimized without the same set of X-ray constraints and parameters of the plasma emission model. Hence, this discrepancy can be partially explained by a difference when computing the effective area of the telescope, for example. We could not reproduce the analysis as Bonamigo et al. (2017, 2018) because the energy at which this effective area has been computed is not publicly available. Hence, we can not assess if there are other reasons for the disagreement between our estimate and B19.

## 7 TESTING OUR RECONSTRUCTION WITH THE MASS DISTRIBUTION OF A MOCK CLUSTER

To estimate the improvement in the mass reconstruction provided by our new method, we create a mock cluster based on our modelling of AS1063 without B-spline perturbations. We start from the best-fitting model and remove all cluster members that are modelled according to the Faber & Jackson law, assuming that the fundamental plane is a good representation of the actual cluster member masses. We use this mass model to predict X-ray counts with the associated Poisson-Gamma mixture while keeping other quantities fixed, such as the background count or the cooling functions. Regarding the lensing constraints, we use the barycentre of the source positions of the actual multiply imaged systems of AS1063 and predict new multiple images with the mock mass. Thanks to this procedure, the realization

of the lensing constraints is fairly similar to reality. We then add a circular Gaussian error of 0.2 arcsec on the positions of the images to account for line-of-sight perturbers and other systematics.

We create four different models to reproduce the mass of the mock cluster while assessing the improvements due to the inclusion of cluster members as well as in the large-scale distribution of the gas. We consider three models: one without the gas and galaxies with only the Faber & Jackson scaling; one including cluster galaxies with the fundamental plane only and a third including only the X-ray emitting gas. And a fourth model that includes both modelling enhancements for cluster member galaxies and the X-ray emitting gas. To account for the combination of the fundamental plane and Faber & Jackson relations, we split the considered galaxies between them to have the same proportions as in the real model. We created the two fundamental plane groups (i.e. galaxies with available $R_e$ or $\sigma_e$ measurements) by taking the most luminous galaxies from the two initial ones, and leaving the remaining elements to the Faber & Jackson scaling. Indeed, inside the MUSE footprint and excluding contamination between cluster members, the faintest ones tend to be modelled with the latter relation as they do not have enough signal for morphological or spectroscopic measurements.

Fig. 10 shows the relative differences of the 2D mass distribution between the simulation and the best fit of each model. From left to right, we have the reference model (i.e. no gas and uncalibrated Faber & Jackson), only gas model, then the fundamental plane on top of the reference and finally, both enhancements. The best-fitting parameters, as well as statistics from the posterior distributions, are available in Appendix F. If we do not consider the discrepancies around the NE halo that are due to the weaker effect of the strong lensing constraints, differences are mainly below 10 per cent in absolute value in the core. As expected, the model with both improvements performs the best with errors mostly below 2 per cent except for the NE halo and some cluster members that switched from the fundamental plane to the Faber & Jackson relation. But the position of the NE halo is mainly constrained by its priors. The real NE position is included in the global posterior distribution of the two models using the fundamental plane scheme. Regarding the two other models, the relative declination of the NE halo is not recovered, with a maximum bias of ~9 arcsec for the model without any improvements. We note that the imperfect modelling scheme of cluster members does not seem to have affected the overall distribution, but only locally, when the discrepancy between relations is higher. We discuss the differences due to the cluster member scaling relations in more detail in Section 7.2.

As shown by the two models where the gas distribution is not considered separately, an underestimation of the mass is introduced

<figure><img src="image_13.png" /><figcaption>**Figure 10.** Maps of the relative differences between the projected surface mass density of each model tested on the mock cluster and the true one. From left to right, the models tested are the reference one, only the gas addition, only the fundamental plane one, and finally, both of them. All maps have the same orientation, with the north in the upward direction and the east directed at the left.</figcaption></figure>
<figure><img src="image_14.png" /><figcaption>**Figure 11.** Maps of the relative differences between the projected gas mass distribution of the model that include only the gas and both the gas and the fundamental plane with the true distribution. The earlier and the latter are presented in the left- and right-panel panels, respectively.</figcaption></figure>

where the gas distribution is locally modifying the main halo ellipticity. However, the gap is reasonable as it is below 4 per cent, which indicates that models from previous studies are still providing a fairly accurate reproduction of the total mass. However, this could have an impact on analyses relying on an accurate estimation of the magnification, as a few per cent offset in mass can lead to significant differences in this quantity. It is also expected that in the case of a more perturbed cluster where the gas is not following the total mass, this underestimation would be greater thus biasing the resulting analysis.

## 7.1 Gas distribution reproduction

Fig. 11 shows the relative differences between the fitted gas distribution and the input one. As one can see, our method recovers it with excellent precision (≤2 per cent) as the discrepancy is mostly ≤1 per cent for the two models. Regarding the X-ray intrinsic error, our models are able to recover the input values (0.538 count) with good agreement (i.e. included in the "1σ"CI). Indeed, the model without the fundamental plane and the one with are presenting X-ray errors of $0.528_{-0.023}^{+0.022}$ and $0.529_{-0.024}^{+0.024}$ counts, respectively.

Regarding our goodness-of-fit procedure detailed in Section 6.1, the models with and without the fundamental plane show likelihoods for the best-fitting model that are inside the "1σ"CI and the "3σ"CI, respectively, of the expected likelihood distribution according to these count models. To assess with more confidence the previous results, we optimize 10 gas-only models with new realizations of X-ray constraints based on the same initial count model. The best-fitting models are showing 5 overfitting models, that are inside the "1σ"CI, and 5 underfitting models split between the limits of the "1σ"CI and the "3σ"CI. Hence, our two mock results made with a

joint fit are consistent with these gas-only models. As we are here in an idealized case, assuming a larger CI of "5σ" to define if a model reaches a good fit is a fair hypothesis according to the results on a mock cluster. We can also add that our goodness-of-fit score is biased by the non-consideration of the whole model posterior distribution for computation cost mitigation. We only use the best-fitting gas model to sample the expected likelihood distribution as it needs more than $10^5$ sample to provide reliable results on one model only, making percentile analysis amongst the whole distribution costly in terms of memory. Thus, we can expect that using the whole posterior distribution could provide different goodness-of-fit score than currently.

## 7.2 Recovering of the cluster member masses

To assess the efficiency of each model to estimate the mass of cluster members, we look at the recovered $r_{cut}$, $\sigma_0$ and their total mass $M$. For each model in the posterior distribution, we build a distribution from the ratio between the considered parameter and its true value amongst all cluster members. We compute the density associated with the previous distribution with the kernel density estimator implemented in Kumar et al. (2019) with a fixed bandwidth of 0.05. From the results on each model, we obtain the posterior distribution of this density. The plain lines and shaded areas in plots in Fig. 12 show the median and the "1σ" CI of the density distribution for each parameter. As one can see, when using fundamental plane, we recover the $r_{cut}$ and $\sigma_0$ more accurately than for the two other models, that rely on the Faber & Jackson solely. However, there are still offsets of ~20 per cent for the $r_{cut}$ when it is down to only a few per cent for the $\sigma_0$. In particular, this gap in the two parameters ($r_{cut}$ and $\sigma_0$) can be explained by the goodness of fit of the scaling factor, $v$, between the light and the mass distribution. Indeed, the model with gas distribution and the one without have $v$ values of $1.87_{-0.34}^{+0.35}$ and $2.45_{-0.36}^{+0.38}$ for an input one of 2.14, respectively. Hence, the biases on the posterior distribution are really similar to the ones on $r_{cut}$ which is expected due to the fact that $r_{cut} \propto v$ (see Appendix C) but also shows that the effect of $r_{core}^*/r_{cut}^*$ have a minor influence on the final estimation.

Regarding the total mass, $M$, shown on the right plot in Fig. 12, the differences when using of the fundamental plane with respect to the other models are less striking. Indeed, the results from the two models using the plane present similar biases than the models relying on the simpler scaling law. The main differences are in the precision of the estimation, which is lower for the two models using the plane, encompassing with a slightly higher probability density the areas around 1. As the $r_{core}$ estimation can not explain these biases (i.e.

<figure><img src="image_15.png" /><figcaption>**Figure 12.** Density distribution of the ratios between estimated dPIE parameters in the different models and their true values. Densities have been estimated for each model in the considered posterior distribution with the kernel density estimator implemented in Kumar et al. (2019) and a fixed bandwidth of 0.05. The shaded areas and the plain lines represent the "1σ"CI, and the median of the distribution of density. From left to right, the parameters considered are the cut radius, $r_{cut}$, the central velocity dispersion, $\sigma_0$, and the total mass, $M$.</figcaption></figure>