# CHEX-MATE: CLUster Multi-Probes in Three Dimensions (CLUMP-3D), I. Gas Analysis Method using X-ray and Sunyaev-Zel'dovich Effect Data

Junhan Kim, Jack Sayers, Mauro Sereno, Iacopo Bartalucci, Loris Chappuis, Sabrina de Grandi, Federico de Luca, Marco de Petris, Megan E Donahue, Dominique Eckert, et al.

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

To cite this version:

Junhan Kim, Jack Sayers, Mauro Sereno, Iacopo Bartalucci, Loris Chappuis, et al. CHEX-MATE: CLUster Multi-Probes in Three Dimensions (CLUMP-3D), I. Gas Analysis Method using X-ray and Sunyaev-Zel'dovich Effect Data. Astronomy & Astrophysics - A&A, 2024, 686, pp.A97. 10.1051/0004-6361/202347399. hal-04171366

HAL Id: hal-04171366

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

Submitted on 31 Oct 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_2.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# CHEX-MATE: CLUster Multi-Probes in Three Dimensions
(CLUMP-3D)

## I. Gas analysis method using X-ray and Sunyaev–Zel’dovich effect data

Junhan Kim<fnref n="1" />,<fnref n="2" />, Jack Sayers<fnref n="1" />, Mauro Sereno³,⁴, Iacopo Bartalucci⁵, Loris Chappuis⁶, Sabrina De Grandi⁷, Federico De Luca⁸,⁹, Marco De Petris¹⁰, Megan E. Donahue¹¹, Dominique Eckert⁶, Stefano Ettori³,⁴,⁵, Massimo Gaspari²⁴,¹², Fabio Gastaldello⁵, Raphael Gavazzi¹³,¹⁴, Adriana Gavidia<fnref n="1" />, Simona Ghizzardi⁵, Asif Iqbal¹⁵, Scott T. Kay¹⁶, Lorenzo Lovisari⁵,¹⁷, Ben J. Maughan¹⁸, Pasquale Mazzotta⁸,⁹, Nobuhiro Okabe¹⁹,²⁰,²¹, Etienne Pointecouteau²², Gabriel W. Pratt¹⁵, Mariachiara Rossetti⁵, and Keiichi Umetsu²³

*(Affiliations can be found after the references)*

Received 7 July 2023 / Accepted 6 March 2024

### ABSTRACT

Galaxy clusters are the products of structure formation through myriad physical processes that affect their growth and evolution throughout cosmic history. As a result, the matter distribution within galaxy clusters, or their shape, is influenced by cosmology and astrophysical processes, in particular the accretion of new material due to gravity. We introduce an analysis method for investigating the three-dimensional triaxial shapes of galaxy clusters from the Cluster HERITAGE project with XMM-Newton – Mass Assembly and Thermodynamics at the Endpoint of structure formation (CHEX-MATE). In this paper, the first in a CHEX-MATE triaxial analysis series, we focus on utilizing X-ray data from XMM-Newton and Sunyaev–Zel’dovich (SZ) effect maps from Planck and the Atacama Cosmology Telescope to obtain a three-dimensional triaxial description of the intracluster medium (ICM) gas. We present the forward modeling formalism of our technique, which projects a triaxial ellipsoidal model for the gas density and pressure, to be compared directly with the observed two-dimensional distributions in X-rays and the SZ effect. A Markov chain Monte Carlo is used to estimate the posterior distributions of the model parameters. Using mock X-ray and SZ observations of a smooth model, we demonstrate that the method can reliably recover the true parameter values. In addition, we applied the analysis to reconstruct the gas shape from the observed data of one CHEX-MATE galaxy cluster, PSZ2 G313.33+61.13 (Abell 1689), to illustrate the technique. The inferred parameters are in agreement with previous analyses for the cluster, and our results indicate that the geometrical properties, including the axial ratios of the ICM distribution, are constrained to within a few percent. With a much better precision than previous studies, we thus further establish that Abell 1689 is significantly elongated along the line of sight, resulting in its exceptional gravitational lensing properties.

**Key words.** galaxies: clusters: general – galaxies: clusters: intracluster medium – galaxies: clusters: individual: Abell 1689 – cosmology: observations – X-rays: galaxies – X-rays: galaxies: clusters

## 1. Introduction

Galaxy clusters are useful probes of structure formation, astrophysical processes such as shocks and feedback from active galactic nuclei, and cosmology (Davis et al. 1985; Voit 2005; Allen et al. 2011; Kravtsov & Borgani 2012; Markevitch & Vikhlinin 2007; McNamara & Nulsen 2007). For instance, they are fundamental to the science goals of numerous ongoing and upcoming large survey projects, including those carried out by eROSITA (extended ROentgen Survey with an Imaging Telescope Array; Predehl et al. 2021), Euclid (Euclid Collaboration 2019), and the Rubin Observatory (Ivezić et al. 2019). In order to maximize the scientific reach of such programs, particularly with regard to cosmological parameter constraints, it is crucial to accurately characterize the ensemble average physical properties of galaxy clusters along with the intrinsic scatter relative to these averages (e.g., Lima & Hu 2005; Zhan & Tyson 2018; Euclid Collaboration 2019). One such example are the scaling relations used to connect global galaxy cluster observables to the underlying halo mass (Rozo & Rykoff 2014; Mantz et al. 2016; Pratt et al. 2019). While these scaling relations are generally sensitive to a range of astrophysical processes (e.g., Ansarifard et al. 2020), some

observables, including the gravitational weak lensing measurements often used to determine the absolute mass, have deviations from average relations that are dominated by projection effects related to asphericity and orientation (Meneghetti et al. 2010; Becker & Kravtsov 2011).

The Cluster HERITAGE project with XMM-Newton – Mass Assembly and Thermodynamics at the Endpoint of structure formation (CHEX-MATE; CHEX-MATE Collaboration 2021)<fnref n="1" /> is an effort to provide a more accurate understanding of the population of galaxy clusters at low-z and with high masses, particularly in the context of cosmology and mass calibration, including the shape of their matter distributions and the effects of the baryonic physics on their global properties. The project is based on a 3 Ms XMM-Newton program that observed 118 galaxy clusters, from two equal-sized subsamples selected from the Planck all-sky Sunyaev–Zel’dovich (SZ)<fnref n="2" /> effect survey. The CHEX-MATE Tier-1 and Tier-2 samples each include 61 galaxy clusters with four overlapping clusters and represent a volume-limited sample ($0.05 < z \le 0.2$) in the local universe and a

<footnote n="1">¹ http://xmm-heritage.oas.inaf.it/</footnote>

<footnote n="2">² Throughout this paper, we use the abbreviation “SZ” to represent the thermal Sunyaev–Zel’dovich effect.</footnote>
mass-limited sample ($M_{500} \ge 7.25 \times 10^{14} M_{\odot}$)<fnref n="3" /> of the most massive objects in the universe, respectively. The X-ray observing program has recently been completed, and initial results from the analyses of these data along with publicly available SZ data have already been published (Campitiello et al. 2022; Oppizzi et al. 2023; Bartalucci et al. 2023).

We utilized triaxial modeling techniques (e.g., Limousin et al. 2013) to investigate the three-dimensional mass distribution within the CHEX-MATE galaxy clusters to infer their intrinsic properties. This approach was motivated by two factors: (1) three-dimensional triaxial shapes provide a better approximation of galaxy clusters than spherical models, and the parameters obtained from such an analysis, such as mass, have lower levels of bias and intrinsic scatter (Becker & Kravtsov 2011; Khatri & Gaspari 2016). (2) A correlation between the triaxial shape of the dark matter (DM) halo and its formation history has been established in simulations (e.g., Ho et al. 2006; Lau et al. 2021; Stapelberg et al. 2022), suggesting that triaxial shape measurements can be a powerful probe of cosmology independent of other techniques currently in use. For instance, some lensing-based shape measurements have found good agreement with ΛCDM predictions (Oguri et al. 2010; Chiu et al. 2018), while a recent multi-probe triaxial analysis suggests a ≈2σ discrepancy between the observed and predicted minor to major axial ratio distributions (Sereno et al. 2018). This discrepancy could indicate that clusters formed more recently than predicted. Alternatively, elevated merger rates (Despali et al. 2017), a reduced influence of baryons on the DM (Suto et al. 2017), or an enhanced feedback (Kazantzidis et al. 2004) could also explain the observed cluster shapes. CHEX-MATE offers a uniform selection of galaxy clusters with consistent measurements of the intracluster medium (ICM) density and temperature. This clean, well-characterized selection with a large sample size (~80 clusters excluding major mergers; Campitiello et al. 2022) will enable a robust cosmological measurement of the triaxial shape distribution.

For our analysis, we adopted the CLUster Multi-Probes in Three Dimensions (CLUMP-3D; Sereno et al. 2017, 2018; Chiu et al. 2018; Sayers et al. 2021) project and implemented significant updates to the modeling package. CLUMP-3D incorporates multiwavelength data from X-ray observations (surface brightness and spectroscopic temperature), millimeter-wave observations (SZ surface brightness), and optical observations (gravitational lensing), which are the projected observables. Then, it assumes triaxial distributions of the ICM gas and matter density profiles. Taking advantage of the different dependences of the X-ray and SZ signals on the gas density and temperature, these signals are used to probe the line-of-sight extent of the ICM, and gravitational lensing data are used to probe the projected matter distribution. In particular, the X-ray emission observed from the ICM is proportional to the line-of-sight integral of the squared electron density ($n_e$) multiplied by the X-ray cooling function, $\Lambda$, represented as $S_X \propto \int n_e^2 \Lambda d\ell$. Meanwhile, the detected SZ signal is proportional to the line-of-sight integral of the product of electron density and temperature ($T_e$), denoted as $B_{SZ} \propto \int n_e T_e d\ell$. Given that the ICM temperature ($T_X$) can be spectroscopically measured using X-ray observations, the line-of-sight elongation ($\Delta l$) can subsequently be determined through the combination of these three measurements as

$\Delta l \sim (B_{SZ}^2 \Lambda)/(S_X T_X^2)$. Assuming a co-alignment of the triaxial axes of the ICM and DM distributions, while still allowing for different axial ratios for the two quantities, our multi-probe analysis can thus constrain the three-dimensional shapes of galaxy clusters. CLUMP-3D was introduced in Sereno et al. (2017), where the authors inferred the triaxial matter and gas distribution of the galaxy cluster MACS J1206.2-0847. The technique built upon similar methods developed to constrain cluster morphology (e.g., Sereno & Umetsu 2011; Sereno et al. 2012). It was then applied to measure the shapes of the Cluster Lensing And Supernova survey with Hubble (CLASH<fnref n="4" />; Postman et al. 2012) clusters, to probe the ensemble-averaged three-dimensional geometry (Sereno et al. 2018; Chiu et al. 2018) as well as the radial profile of the nonthermal pressure fraction (Sayers et al. 2021). These results demonstrate the potential of the three-dimensional triaxial shape measurement technique, but they were relatively imprecise due to the sample size, data quality, and systematics related to cluster selection. Thus, the much larger CHEX-MATE galaxy cluster sample, with a well-understood selection function and more uniform and higher-quality X-ray data, will provide improved statistics and more robust constraints on the shape measurements.

In this paper we demonstrate several improvements to the original CLUMP-3D formalism while modeling the ICM distributions observed by XMM-Newton and Planck and ground-based SZ data from the Atacama Cosmology Telescope (ACT). As detailed below, we performed a fully two-dimensional analysis of the X-ray temperature (Lovisari et al. 2024) and SZ data, whereas the original CLUMP-3D only treated the X-ray surface brightness (SB) in two dimensions while using one-dimensional azimuthally averaged profiles of both the X-ray spectroscopic temperatures and the SZ effect data. In addition, we now model the ICM gas density and pressure instead of its density and temperature. This allows us to fit the data with fewer parameters, thus accelerating the model fitting process. Additionally, we fully rewrote the code in Python to facilitate the future public release of the package. In Sect. 2 we summarize the triaxial analysis formalism and describe the model fitting method. In Sect. 3 we introduce the X-ray and SZ data from our program and apply the technique to a CHEX-MATE galaxy cluster. In a subsequent paper, we will include gravitational lensing constraints in a manner that also builds upon, and improves, the existing CLUMP-3D technique. With these X-ray, SZ effect, and gravitational lensing data, we will be able to model the triaxial distributions of both the ICM and the DM. Throughout this study, we adopt a ΛCDM cosmology characterized by $H_0 = 70 \text{ km s}^{-1} \text{ Mpc}^{-1}$, $\Omega_m = 0.3$, and $\Omega_\Lambda = 0.7$. $E(z)$ represents the ratio of the Hubble constant at redshift, $z$, to its present value, $H_0$, and $h_{70} = H_0/(100 \text{ s}^{-1} \text{ Mpc}^{-1})/0.7$.

## 2. Triaxial analysis: Formalism and the model fit

While the mathematical description of a triaxial geometry for astronomical objects and their physical profiles has been introduced in previous studies (e.g., Stark 1977; Binggeli 1980; Binney 1985; Oguri et al. 2003; De Filippis et al. 2005; Corless & King 2007; Sereno et al. 2010, 2012, 2018), these works lack consistency in their notation. To prevent confusion, we present our mathematical formalism for triaxial modeling in this section. We then describe our model fitting procedure and the implementation of the fitting algorithm in our software package.

<footnote n="3">³ The parameter $M_{500}$ denotes the mass enclosed within a radius ($R_{500}$) where the mean overdensity is 500 times the critical density at a specific redshift, and we used the $M_{500}$ and $R_{500}$ values from Planck Collaboration XXVII (2016).</footnote>

<footnote n="4">⁴ https://www.stsci.edu/~postman/CLASH/</footnote>
As this paper focuses on the analysis method of studying ICM distributions, we do not include the gravitational lensing data in our fits. Future works in this series will expand our formalism to include total mass density profiles constrained by gravitational lensing measurements. For instance, in the case of a Navarro–Frenk–White (NFW; Navarro et al. 1997) profile the gravitational lensing analysis requires two additional parameters – total mass ($M_{200}$) and concentration ($c_{200}$) – assuming that the gas and matter distributions are co-aligned along the ellipsoidal axes. This assumption is well supported by a two-dimensional weak-lensing and X-ray analysis of 20 high-mass galaxy clusters (Umetsu et al. 2018), as well as by cosmological hydrodynamical simulations (Okabe et al. 2018).

## 2.1. Geometry and projection

To connect the intrinsic cluster geometry to the projected properties observed in the plane of the sky, we assume a triaxial ellipsoidal model for the gas distribution, where the thermodynamic profiles of the ICM are represented as a function of $\zeta$, the ellipsoidal radius. In the intrinsic coordinate system of the ellipsoid ($x_1, x_2, x_3$), it is defined as:

$$ \zeta^2 = \frac{x_1^2}{q_1^2} + \frac{x_2^2}{q_2^2} + x_3^2, \qquad (1) $$

where $q_1$ and $q_2$ are minor-to-major and intermediate-to-major axial ratios, respectively ($0 < q_1 \le q_2 \le 1$). Given a semi-major axis of the ellipsoid $l_s$, the volume of the ellipsoid is $(4\pi/3)l_s^3q_1q_2$. The ellipsoid becomes a prolate shape if $q_1 = q_2 \le 1$ and an oblate shape if $q_1 \le q_2 = 1$.

Figure 1 illustrates the geometry of the ellipsoid and the involved coordinate systems. It is essential to note that the axes defining the ICM model may not align with the observer's frame. To relate the ellipsoid's intrinsic coordinate system ($x_1^{\text{int}}$, $x_2^{\text{int}}$, $x_3^{\text{int}}$) to the observer's coordinate frame ($x_1^{\text{obs}}$, $x_2^{\text{obs}}$, $x_3^{\text{obs}}$), we employ three Euler angles. These angles describe the relationship between the two coordinate systems: (1) the angle between $x_1^{\text{int}}$, aligned with the major axis of the ellipsoid, and $x_3^{\text{obs}}$, which lies along the observer's line of sight ($\theta$), (2) the angle between $x_1^{\text{int}}$ and the line of nodes ($\varphi$), and (3) the angle between $x_1^{\text{obs}}$ and the line of nodes ($\psi$). The line of nodes is the intersection of the $x_1^{\text{int}}-x_2^{\text{int}}$ plane and the $x_1^{\text{obs}}-x_2^{\text{obs}}$ plane, and it is aligned with the vector $x_3^{\text{int}} \times x_3^{\text{obs}}$.

We can derive the geometric properties of the projected ellipse from the intrinsic parameters of the ellipsoid when it is projected onto the plane from any direction. These properties encompass the semimajor axis of the projected ellipse $l_p$, its ellipticity $\epsilon$, the orientation of the ellipse in the plane of the sky $\theta_\epsilon$, and the elongation parameter $e_\parallel$. The projected profiles are expressed as a function of $\xi$, the elliptical radius of the ellipse in the plane of the sky.

The ellipticity of the projected ellipse ($\epsilon$) is

$$ \epsilon = 1 - q_p, \qquad (2) $$

where $q_p$ is the minor-to-major axial ratio of the observed projected isophote ($q_p \le 1$), which is the inverse of $e_p$ used in Sereno et al. (2012), and

$$ q_p = \sqrt{\frac{j + l - \sqrt{(j - l)^2 + 4k^2}}{j + l + \sqrt{(j - l)^2 + 4k^2}}}, \qquad (3) $$

<figure><img src="image_3.png" /><figcaption>Fig. 1. Triaxial ellipsoid model and coordinate systems used in the triaxial analysis. The intrinsic coordinate system of the ellipsoid is denoted by dotted gray arrows ($x_1^{\text{int}}$, $x_2^{\text{int}}$, and $x_3^{\text{int}}$), where $x_1^{\text{int}}$ represents the major axis. The black arrows ($x_1^{\text{obs}}$, $x_2^{\text{obs}}$, and $x_3^{\text{obs}}$) correspond to the observer's coordinate system, where $x_1^{\text{obs}}$ is aligned with the observer's line of sight. In other words, an observer views the ellipsoid in the $-x_3^{\text{obs}}$ direction. The three Euler angles ($\theta$, $\varphi$, $\psi$) characterize the intrinsic coordinate system of the ellipsoid in relation to the observer's coordinate system. The blue line represents the line of nodes, which is the intersection of the $x_1^{\text{int}}-x_2^{\text{int}}$ plane and the $x_1^{\text{obs}}-x_2^{\text{obs}}$ plane, and it is aligned with the vector $x_3^{\text{int}} \times x_3^{\text{obs}}$. The red ellipse denotes the projection of the ellipsoid on the sky plane, with $l_p$ representing its semimajor axis. The dashed black line on the ellipse shows the projected major axis of the ellipsoid on the sky plane. The green ellipse is the projection of the ellipsoid onto the plane that is perpendicular to the sky plane, and $l_{\text{los}}$ is the half size of the ellipse along the observer line of sight. See also Figs. 2 and 3 in Sereno et al. (2012).</figcaption></figure>

where

$$ \begin{aligned} j &= \frac{1}{2} \left[ \left( \frac{1}{q_1^2} + \frac{1}{q_2^2} \right) - \frac{\sin^2 \theta \cos^2 \psi (q_1^2 + q_2^2 - 2)}{q_1^2 q_2^2} \right. \\ & \qquad \left. + \left( \frac{1}{q_1^2} - \frac{1}{q_2^2} \right) \left\{ \cos 2\varphi (\cos^2 \theta \cos^2 \psi - \sin^2 \psi) - \cos \theta \sin 2\varphi \sin 2\psi \right\} \right], \\ k &= \frac{1}{4q_1^2 q_2^2} \left[ 2 \cos \theta (q_1^2 - q_2^2) \cos 2\psi \sin 2\varphi \right. \\ & \qquad \left. + \left\{ \sin^2 \theta (q_1^2 + q_2^2 - 2) + (1 + \cos^2 \theta) (q_1^2 - q_2^2) \cos 2\varphi \right\} \sin 2\psi \right], \\ l &= \frac{1}{2} \left[ \left( \frac{1}{q_1^2} + \frac{1}{q_2^2} \right) - \frac{\sin^2 \theta \sin^2 \psi (q_1^2 + q_2^2 - 2)}{q_1^2 q_2^2} \right. \\ & \qquad \left. + \left( \frac{1}{q_1^2} - \frac{1}{q_2^2} \right) \left\{ \cos 2\varphi (\cos^2 \theta \sin^2 \psi - \cos^2 \psi) + \cos \theta \sin 2\varphi \sin 2\psi \right\} \right]. \end{aligned} \qquad (4) $$

It is worth noting that the expressions of $j$, $k$, and $l$ in Stark (1977) and Binggeli (1980) differ from those presented above, as they assumed $\psi = 0$, using only two angles to align the major ellipsoidal axis with the observer's line of sight. However, a coordinate transformation requiring $\psi$ is necessary to align the remaining axes.

The orientation angle in the plane of the sky of the projected ellipse is

$$ \theta_{\epsilon} = \tan^{-1} \left( \frac{l - j + \sqrt{(j - l)^2 + 4k^2}}{2k} \right), \qquad (5) $$
and the elongation parameter of the ellipsoid is

$$e_{\parallel} \equiv \frac{l_{\text{los}}}{l_{\text{p}}} = \sqrt{\frac{q_{\text{p}}}{q_{1}q_{2}}} f^{-3/4}, \qquad (6)$$

where

$$f = \sin^2 \theta \left[ \left( \frac{\sin \varphi}{q_1} \right)^2 + \left( \frac{\cos \varphi}{q_2} \right)^2 \right] + \cos^2 \theta. \qquad (7)$$

The elongation parameter, $e_{\parallel}$, represents the ratio of the size of the ellipsoid along the observer's line of sight to the major axis of the projected ellipse in the sky plane, providing a measure of the three-dimensional geometry of the triaxial ellipsoid model of the ICM. In the gas analysis presented in Sereno et al. (2012), the orientation angle ($\theta_e$; represented as $\epsilon$ by Sereno et al. 2012) was determined from the X-ray map, while the elongation parameter (represented as $e_{\Delta}$, which is equivalent to $1/e_{\parallel}$) was estimated from the combined X-ray and SZ analysis. Later, Sereno et al. (2017) simultaneously constrained the individual Euler angles by treating the axial ratios and three angles as free parameters.

Then, the semimajor axis of the projected ellipse becomes

$$l_{\text{p}} = \frac{l_{\text{s}}}{e_{\parallel} \sqrt{f}} \qquad (8)$$

$$= l_{\text{s}} \sqrt{\frac{q_{1}q_{2}}{q_{\text{p}}}} f^{1/4}, \qquad (9)$$

and the projected length scales $l_s$ and $l_{los}$ are related by the elongation parameter, that is,

$$l_{\text{los}} = l_{\text{s}} / \sqrt{f}. \qquad (10)$$

In the plane of the sky, an elliptical radius $\xi$ becomes

$$\xi^2 = \left( x_1^2 + \frac{x_2^2}{q_p^2} \right) \left( \frac{l_s}{l_p} \right)^2 \qquad (11)$$

(Sereno et al. 2010)<fnref n="5" />.

Finally, three-dimensional volume density can be projected onto the sky plane by utilizing the geometric parameters

$$F_{2\text{D}}(\xi; l_{\text{p}}, p_{\text{i}}) = \frac{2}{\sqrt{f}} \int_{\xi}^{\infty} F_{3\text{D}}(\zeta; l_{\text{s}}, p_{\text{i}}) \frac{\zeta}{\sqrt{\zeta^2 - \xi^2}} d\zeta, \qquad (12)$$

$$F_{2\text{D}}(x_{\xi}; l_{\text{p}}, p_{\text{i}}) = 2l_{\text{p}}e_{\parallel} \int_{x_{\xi}}^{\infty} F_{3\text{D}}(x_{\xi}; l_{\text{s}}, p_{\text{i}}) \frac{x_{\xi}}{\sqrt{x_{\xi}^2 - x_{\xi}^2}} dx_{\xi}, \qquad (13)$$

where $x_{\xi} = \zeta/l_s$, $x_{\xi} = \xi/l_p$, and $p_i$ are the parameters describing the intrinsic density profile (Stark 1977; Sereno 2007; Sereno et al. 2010). Using this projection, we calculated the SZ and X-ray maps on the sky plane from the three-dimensional ellipsoidal distribution of the ICM profiles and fit the model to the data. We describe the analytic profiles ($F_{3D}$) for the physical quantities related to the direct observables ($F_{2D}$) in the next section.

## 2.2. Electron density and pressure profiles

We used smooth analytic functions of the electron density and pressure profiles to describe the thermodynamics and spatial distribution of the ICM, and then used these functions to compute observable quantities, such as the SZ effect map, the X-ray SB map, and the X-ray temperature map. The model lacks the ability to effectively constrain small-scale structures that deviate from its assumptions. However, the three-dimensional description of the profiles provides a better approximation compared to spherical models. After accounting for instrumental effects, such as the point spread function (PSF), these model maps are then compared to the observed data. The original CLUMP-3D package, as detailed in Sereno et al. (2017), instead assumed smooth analytic functions for the gas density and temperature (Vikhlinin et al. 2006; Baldi et al. 2012). However, because the presence (or not) of a cool core alters the overall shape of the temperature profile (e.g., Pratt et al. 2007), the analytic function needs to be sufficiently flexible to allow for either a decrease or increase in temperature at small radii. Pressure profiles are more regular in their global shape (e.g., Arnaud et al. 2010), and therefore a simpler function with fewer free parameters can be used to describe the ICM. Thus, our overall model can be more easily constrained than the one used by Sereno et al. (2017). Table 1 lists the model parameters used in our gas analysis, including the geometric parameters described in the previous section.

The electron density profile is described as

$$n_e(\zeta) = n_0 \left(\frac{\zeta}{\zeta_c}\right)^{-\eta_e} \left[1 + \left(\frac{\zeta}{\zeta_c}\right)^2\right]^{-3\beta_e/2 + \eta_c/2} \left[1 + \left(\frac{\zeta}{\zeta_t}\right)^3\right]^{-\gamma_e/3}, \qquad (14)$$

where $n_0$ is the central electron density, $\zeta_c$ is the core radius, and $\zeta_t$ is the tidal radius ($\zeta_t > \zeta_c$). ($\beta_e$, $\eta_e$, $\gamma_e$) represent the power law exponent of the electron density distribution for the intermediate, inner, and external slope of the profile, respectively (Vikhlinin et al. 2006; Ettori et al. 2009). The electron pressure profile is modeled using a generalized NFW (gNFW) profile (Navarro et al. 1996; Nagai et al. 2007; Arnaud et al. 2010). It is described as

$$\frac{P_e(x)}{P_{500}} = \frac{P_0}{(c_{500}x)^{\gamma_p} [1 + (c_{500}x)^{\alpha_p}]^{(\beta_p - \gamma_p)/\alpha_p}}, \qquad (15)$$

where $x = \zeta/R_{500}$, ($\gamma_p$, $\alpha_p$, $\beta_p$) describes the power law exponent for the central ($r \ll r_s$), intermediate ($r \sim r_s = R_{500}/c_{500}$), and outer ($r \gg r_s$) regions, and the characteristic pressure is

$$P_{500} = 1.65 \times 10^{-3} E(z)^{8/3} \times \left[ \frac{M_{500}}{3 \times 10^{14} h_{70}^{-1} M_{\odot}} \right]^{2/3} h_{70}^2 \text{ keV cm}^{-3}. \qquad (16)$$

The expressions for $P_{500}$ provided in Nagai et al. (2007) and Arnaud et al. (2010) represent the gas pressure and the electron pressure, respectively. We opted to use the electron pressure formulation. In order to convert the electron pressure, $P_e$, into gas pressure, it is necessary to incorporate both the mean molecular weight and the mean molecular weight per free electron into the calculations. As noted by Nagai et al. (2007), strong degeneracies between the pressure profile parameters generally prevent meaningful constraints when all are varied (see also Battaglia et al. 2012). For our baseline fits, we thus fixed the values of $c_{500}$ and $\gamma_p$ to 1.4 and 0.3 as in Sayers et al. (2023). In addition, because $\beta_p$ characterizes the pressure profile in the outer regions, it may not be well constrained depending on the

<footnote n="5">⁵ Assuming that the ellipse is expressed as $(x_1^2 + x_2^2)/(q_p^2) = 1$, $q_p$ is the minor-to-major axial ratio ($b/a$), and the elliptical radius, which is the corresponding major axis length, becomes $\sqrt{x_1^2 + x_2^2/(q_p^2)}$ because $x_1^2 + x_2^2/(q_p^2) = a^2$.</footnote>
<table><caption>**Table 1.** Gas model parameters.</caption><thead><tr><th>Parameter</th><th>Units</th><th>Description</th><th>Default prior</th></tr></thead><tbody><tr><td colspan="4" style="text-align: center;">Geometrical parameters of a triaxial ellipsoid (Eqs. (1) and (4))</td></tr><tr><td>q<sub>ICM,1</sub></td><td></td><td>Minor-to-major axial ratio of the ICM distribution</td><td>U(0, 1)</td></tr><tr><td>q<sub>ICM,2</sub></td><td></td><td>Intermediate-to-major axial ratio of the ICM distribution</td><td>U(q<sub>ICM,1</sub>, 1)</td></tr><tr><td>cos θ</td><td></td><td>Cosine of the inclination angle of the ellipsoid major axis</td><td>U(0, 1)</td></tr><tr><td>φ</td><td>deg</td><td>Second Euler angle</td><td>U(-π/2, π/2)</td></tr><tr><td>ψ</td><td>deg</td><td>Third Euler angle</td><td>U(-π/2, π/2)</td></tr><tr><td colspan="4" style="text-align: center;">Electron density profile (Eq. (24))</td></tr><tr><td>n<sub>0</sub></td><td>cm<sup>-3</sup></td><td>Central scale density of the distribution of electrons</td><td>U(10<sup>-6</sup>, 10)</td></tr><tr><td>ζ<sub>c</sub></td><td>kpc</td><td>Ellipsoidal core radius of the gas distribution</td><td>U(0, 10<sup>3</sup>)</td></tr><tr><td>ζ<sub>t</sub></td><td>Mpc</td><td>Ellipsoidal truncation radius of the gas distribution (ζ<sub>t</sub> > ζ<sub>c</sub>)</td><td>U(ζ<sub>c</sub>/10<sup>3</sup>, 3)</td></tr><tr><td>β<sub>e</sub></td><td></td><td>Slope of the gas distribution (in the intermediate region)</td><td>U(0, 3)</td></tr><tr><td>η<sub>e</sub></td><td></td><td>Slope of the gas distribution (inner)</td><td>U(0, 1)</td></tr><tr><td>γ<sub>e</sub></td><td></td><td>Slope of the gas distribution (outer)</td><td>U(0, 5)</td></tr><tr><td colspan="4" style="text-align: center;">Gas pressure profile (Eq. (15))</td></tr><tr><td>P<sub>0</sub></td><td></td><td>Normalization for the gNFW pressure profile</td><td>U(0, 10<sup>2</sup>)</td></tr><tr><td>c<sub>500</sub></td><td></td><td>Pressure profile concentration (r ~ r<sub>s</sub> = R<sub>500</sub>/c<sub>500</sub>)</td><td>δ(1.4)</td></tr><tr><td>γ<sub>p</sub></td><td></td><td>Slope parameter for central region (r ≪ r<sub>s</sub>)</td><td>δ(0.3)</td></tr><tr><td>α<sub>p</sub></td><td></td><td>Slope parameter for intermediate region (r ~ r<sub>s</sub>)</td><td>U(0, 5)</td></tr><tr><td>β<sub>p</sub></td><td></td><td>Slope parameter for outer region (r ≫ r<sub>s</sub>)</td><td>U(0, 15)<sup>(a)</sup></td></tr></tbody></table>

**Notes.** We consider five geometric parameters ($q_{ICM,1}$, $q_{ICM,2}$, $\theta$, $\varphi$, $\psi$), six electron density parameters ($n_0$, $\zeta_c$, $\zeta_t$, $\beta_e$, $\eta_e$, $\gamma_e$), and five gas pressure parameters ($P_0$, $c_{500}$, $\gamma_p$, $\alpha_p$, $\beta_p$). For the geometric and electron pressure profile parameters, we primarily adopt the priors in Sereno et al. (2018). We also assign delta priors to $c_{500}$ (=1.4), and $\gamma_p$ (=0.3) as a default, resulting in 14 free parameters. In the default prior column, $\mathcal{U}$ refers to a uniform prior and $\delta$ refers to a delta function that fixes the parameter for the model fit. <sup>(a)</sup>For the cluster PSZ2 G313+61.13, to which we applied the model fit in this paper (Sect. 3), we employed a delta prior (Eq. (17)) because we limited the map size to be within $R_{500}$, which results in very little sensitivity to $\beta_p$ (Sayers et al. 2023).

map size chosen for the actual fit. For the demonstration of our approach using actual CHEX-MATE data in Sect. 3, we restrict the map size of the X-ray and SZ observational data to within $R_{500}$ to mask out potential spurious signal at large radii that do not originate from a target cluster, and therefore an external constraint on the value of $\beta_p$ is required. In such cases, we used a value that depends on the mass and redshift, given by

$$ \beta_p = 5.495 \left( \frac{M_{500}}{10^{15} M_{\odot}} \right)^{0.15} (1+z)^{0.02}. \quad (17) $$

This relation is derived from a combined X-ray and SZ analysis of galaxy clusters with a redshift range of $0.05 \le z \le 0.60$ and mass range of $4 \times 10^{14} \le M_{500} \le 30 \times 10^{14} M_{\odot}$ (Sayers et al. 2023). This fit is thus valid for the mass and redshift ranges of the CHEX-MATE clusters, with Tier-1 covering $0.05 < z < 0.2$ and $2 \times 10^{14} < M_{500} < 9 \times 10^{14} M_{\odot}$, and Tier-2 encompassing $z < 0.6$ and $M_{500} > 7.25 \times 10^{14} M_{\odot}$.

## 2.3. Sunyaev–Zel'dovich effect and X-ray observables

In this section, we summarize the observables associated with the SZ effect and the X-ray emissivity, and explain their relationship to the electron density and pressure profiles introduced earlier. The SZ effect is characterized by the Compton-y parameter, which is proportional to the integrated electron pressure along the line of sight.

$$ y \equiv \frac{\sigma_T}{m_e c^2} \int_{\parallel} P_e dl = \frac{\sigma_T k_B}{m_e c^2} \int_{\parallel} n_e T_e dl, \quad (18) $$

where $\sigma_T$ is the Thomson cross-section, $k_B$ is the Boltzmann constant, $n_e$ is the electron number density, and $T_e$ is the electron temperature. The X-ray observations are primarily sensitive

to the SB of the ICM,

$$ \mathrm{SB} = \frac{1}{4\pi(1+z)^3} \int_{\parallel} n_e^2 \Lambda_{\mathrm{eff}}(T_e, Z) dl \quad (19) $$

(Reese et al. 2010), due to thermal Bremsstrahlung, where the cooling function $\Lambda_{\mathrm{eff}}(T_e, Z)$ quantifies the thermal radiation emitted from a fully ionized plasma due to collisions, taking into account the relative abundance of each chemical element. It can be calculated using software such as XSPEC (Arnaud 1996). We used a precalculated table and interpolated the value in the temperature ($T_e$)-metallicity ($Z$) space during the model computation. To calculate the emissivity, the instrument response within the chosen energy band [0.7–1.2] keV and the Galactic hydrogen column density must be taken into account, as explained in Bartalucci et al. (2023), which describes the details of the data analysis used to produce the SB maps. With our software, we performed the calculation using the Python package pyproffit<sup>6</sup> (Eckert et al. 2020).

The XMM-Newton data can also be used to derive projected temperature maps of ICM via spectroscopic fits (Lovisari et al. 2024). Within our model, we approximated this spectroscopic temperature based on the formalism of Mazzotta et al. (2004) as follows:

$$ T_{\mathrm{sp}} = \frac{\int W T_e dV}{\int W dV} \mathrm{keV}; \quad W = \frac{n_e^2}{T_e^{3/4}}, \quad (20) $$

which is valid for Bremsstrahlung ($T_e \ge 3$ keV).

The SZ and X-ray observables (Eqs. (18) and (19)) are modeled as projections of the three-dimensional profiles parameterized by the ellipsoidal radius $\zeta$ (or $x_\zeta$). The three-dimensional

<footnote n="6"><sup>6</sup> https://pyproffit.readthedocs.io/en/latest/intro.html</footnote>
volume density of the models, F<sub>3D</sub>(x<sub>ζ</sub>; l<sub>s</sub>, p<sub>i</sub>), can be written analytically, and the two-dimensional maps are calculated following Eq. (12). The model Compton-y parameter is

$$y_{\text{model}}(x_{\zeta}; l_p, p_i) = (2l_p e_{\parallel}) \left(\frac{\sigma_{\text{T}}}{m_e c^2}\right) \int_{x_{\zeta}}^{\infty} P_e(x_{\zeta}) \frac{x_{\zeta}}{\sqrt{x_{\zeta}^2 - x_{\zeta}^2}} dx_{\zeta}, \quad (21)$$

where

$$P_e(x_{\zeta}) = \frac{P_0 P_{500}}{\left(c_{500} x_{\zeta} \frac{l_s}{R_{500}}\right)^{\gamma_p} \left[1 + \left(c_{500} x_{\zeta} \frac{l_s}{R_{500}}\right)^{\alpha_p}\right]^{\left(\beta_p - \gamma_p\right)/\alpha_p}} \text{keV cm}^{-3}. \quad (22)$$

This integration can be computationally expensive, depending on the size of the map. To expedite the calculation, we created a linearly spaced sample of the (normalized) elliptical radius x<sub>ζ</sub> and interpolated the integration results while generating a model map. We applied the same technique in the X-ray observable calculation. Lastly, we convolved the model map with the appropriate PSF shape (e.g., a 7' full width at half maximum Gaussian map in the case of Planck and a 1.6' full width at half maximum in the case of ACT; see Fig. 2).

Similarly, the X-ray SB (Eq. (19)) model becomes

$$\text{SB}_{\text{model}}(x_{\zeta}; l_p, p_i) = (2l_p e_{\parallel}) \frac{1}{4\pi(1+z)^3} \times \int_{x_{\zeta}}^{\infty} n_e^2(x_{\zeta}) \Lambda_{\text{eff}}(T_e(x_{\zeta}), Z(x_{\zeta})) \frac{x_{\zeta}}{\sqrt{x_{\zeta}^2 - x_{\zeta}^2}} dx_{\zeta}, \quad (23)$$

where

$$n_e(x_{\zeta}) = n_0 \left(x_{\zeta} \frac{l_s}{\zeta_c}\right)^{-\eta_c} \left[1 + \left(x_{\zeta} \frac{l_s}{\zeta_c}\right)^2\right]^{-3\beta_s/2 + \eta_c/2} \left[1 + \left(x_{\zeta} \frac{l_s}{\zeta_t}\right)^3\right]^{-\gamma_c/3}, \quad (24)$$

and the electron temperature is

$$T_e(x_{\zeta}) = \frac{P_e(x_{\zeta})}{n_e(x_{\zeta}) k_B}. \quad (25)$$

We used a radius-dependent metallicity profile Z(x<sub>ζ</sub>) obtained from the X-COP galaxy cluster samples (Ghizzardi et al. 2021) for calculating the cooling function.

Upon generating the model, instrumental responses are incorporated to facilitate a direct comparison between the model and the data. For the XMM-Newton X-ray maps, the sky background in the [0.7–1.2] keV band (2.49 × 10<sup>-4</sup> cts s<sup>-1</sup> arcmin<sup>-2</sup>; Bartalucci et al. 2023) is considered. Specifically, we adopted the sky and particle background measured by the European Photon Imaging Camera (EPIC; Strüder et al. 2001; Turner et al. 2001) M2 CCD in the [0.5–2] keV band and converted it for the [0.7–1.2] keV band. After adding the sky background, the vignetting is applied. Subsequently, the resulting map is convolved with a Gaussian profile to account for the PSF. The nominal PSF of XMM-Newton can be closely represented using a Gaussian function with a 6" full width at half maximum (FWHM)<sup>7</sup>. However, the actual FWHM of the PSF is

dependent on the angle relative to the optical axis, and combining images from different cameras could potentially deteriorate the final PSF. Therefore, we followed the convention of Bartalucci et al. (2023) and assumed the Gaussian has a FWHM of 10". The line-of-sight integration of the observed quantities described above was performed to a depth of 10 Mpc in radius.

To summarize, the observational data used in our analysis includes two-dimensional images of the SZ signal, X-ray SB, and X-ray temperature. Then we used our triaxial model to generate analogous images based on the model parameters delineated in Table 1. The observed and model-generated images can then be directly compared to facilitate our fitting process, and the method employed for this fitting procedure is elaborated upon in the following section.

## 2.4. Fitting formalism

The χ<sup>2</sup> statistic is used to define the likelihood of the model. We used emcee (Foreman-Mackey et al. 2013), a Python-based affine-invariant ensemble Markov chain Monte Carlo (MCMC; Goodman & Weare 2010) package, for the model fitting process. By performing MCMC sampling (Hogg & Foreman-Mackey 2018), we determined the posterior distribution of the parameters that describe the triaxial model. When conducting a model fit with the data, we occasionally needed to adjust the scale parameter of the stretch move within the affine-invariant ensemble sampling algorithm implemented in the package to enhance performance (Huijser et al. 2015).

We define the χ<sup>2</sup> functions for our analysis below; they are based on two-dimensional maps of the SZ and X-ray data rather than the original one-dimensional radial profiles used in the CLUMP-3D method presented in Sereno et al. (2017). The χ<sup>2</sup> function for the two-dimensional SZ map is

$$\chi^2_{\text{SZ}} = \sum_{i,j=1}^{N_s} [y_i - \hat{y}_i] (C_{\text{SZ}}^{-1})_{ij} [y_j - \hat{y}_j], \quad (26)$$

where ŷ<sub>i</sub> is the model Compton-y within a pixel, and y<sub>i</sub> is the observed value. To deal with the correlated noise in the SZ data, we used the inverse of the uncertainty covariance matrix (C<sup>-1</sup><sub>sz</sub>). Similarly, the χ<sup>2</sup> function for the X-ray temperature map becomes

$$\chi^2_{\text{T}} = \sum_{i=1}^{N_{\text{T}}} \left( \frac{T_{\text{sp},i} - \hat{T}_{\text{sp},i}}{\delta T_{\text{sp},i}} \right)^2, \quad (27)$$

where T<sub>sp,i</sub> is the model spectroscopic temperature within a pixel, and T<sub>sp,i</sub> is the observed value with uncertainty δT<sub>sp,i</sub>.

For the X-ray SB, we employed a dual approach. We used a two-dimensional model fit within the circular region that encloses 80% of the emission and one-dimensional analysis for the outside region where the background and the source emission is comparable and signal-to-noise ratio is relatively low. In the exterior region, we computed azimuthal medians in annular bins to mitigate biases in measuring X-ray SB caused by gas clumping, as suggested by Eckert et al. (2015). While our current analysis solely uses the two-dimensional map of X-ray temperature, in future work we intend to implement an approach that is fully consistent with our treatment of the X-ray SB to also mitigate local deviations from homogeneity in the X-ray temperature data (Lovisari et al. 2024). Then, the combined likelihood becomes

$$\chi^2_{\text{SB}} = \chi^2_{\text{SB,1D}} + \chi^2_{\text{SB,2D}} \quad (28)$$

<footnote n="7">7 https://xmm-tools.cosmos.esa.int/external/xmm_user_support/documentation/uhb/onaxisraypsf.html</footnote>
where

$$ \chi_{\text{SB,1D}}^2 = \sum_{i=1}^{N_{\text{SB,1D}}} \left( \frac{S_{\text{X,1D},i} - \hat{S}_{\text{X,1D},i}}{\delta S_{\text{X,1D},i}} \right)^2, \qquad (29) $$

and

$$ \chi_{\text{SB,2D}}^2 = \sum_{i=1}^{N_{\text{SB,2D}}} \left( \frac{S_{\text{X,2D},i} - \hat{S}_{\text{X,2D},i}}{\delta S_{\text{X,2D},i}} \right)^2 . \qquad (30) $$

Here $\hat{S}_{\text{X},i}$ is the model SB, and $S_{\text{X},i}$ and $\delta S_{\text{X},i}$ are obtained from the observational data. We employed SB measurements and the corresponding error for our two-dimensional analysis here, assuming Gaussian statistics. This should be a valid assumption, as we define regions with sufficiently large photon counts (i.e., $\ge 20$). However, the formally correct approach is to use the Cash statistic, which accounts for Poisson fluctuations in the photon counts (Cash 1979). Fits using the Cash statistic for photon counting in the low count regime will be explored in future works.

Finally, the total $\chi^2$ statistic becomes

$$ \chi_{\text{X+SZ}}^2 = \chi_{\text{SZ}}^2 + \chi_{\text{T}}^2 + \chi_{\text{SB}}^2, \qquad (31) $$

and the MCMC is used to sample $\chi^2_{\text{X+SZ}}$ within the parameter space near the best fit.

## 2.5. Parameter estimation with mock data

To validate the accuracy of our model fitting algorithm, we conducted a full analysis using mock observations of a galaxy cluster described by our model from known input parameter values. Using the model parameters outlined in Table 1, we generate model SZ, X-ray SB and temperature maps, incorporating the instrument PSF response. For this test, we generated a wide range of mock datasets by varying the inclination angle ($\cos\theta$) and minor-to-major axial ratio of the ICM distribution ($q_{\text{ICM,1}}$), using the reference model as a basis. The reference mock cluster has the following characteristics: $z = 0.18$, $M_{500} = 8.8 \times 10^{14} M_{\odot}$, $R_{500} = 7.4'$. Additionally, we used the following values for the geometric configuration and electron density and pressure parameters, with ($q_{\text{ICM,1}}$, $q_{\text{ICM,2}}$, $\cos\theta$, $\phi$, $\psi$) = (0.6, 0.75, 0.8, -25, 60), ($n_0$, $\zeta_c$, $\zeta_t$, $\beta_c$, $\eta_c$, $\gamma_c$) = (0.002, 175, 1.5, 0.6, 0.3, 1.8), $P_0$, $\alpha_p$ = (10.0, 1.0). In this case, $\epsilon_\parallel = 1.02$.

The model maps generated with the input parameters of the reference model are presented in Fig. 2. For each pixel based on its coordinates within the observed map, we calculated the observables projected onto the two-dimensional sky plane (Sect. 2.3). Then, instrumental effects, including the PSF response, are applied. As we discuss in the next section, our baseline analysis of the observed data uses the combined Planck and ACT SZ effect map (Aiola et al. 2020), and we assumed a PSF with a FWHM of 1.6'. To ensure adequate angular sampling of the PSF, we required a maximum pixel size equal to the FWHM divided by 3.

In addition, we incorporated noise into each mock observation. Using the error maps for the observed data, we randomly sampled Gaussian noise distributions for the SZ, X-ray SB, and X-ray temperature maps, respectively. Figure 3 shows the posterior distribution of the parameters from our fit to this mock observation. The posterior distributions indicate that we can accurately recover most of the varied parameter values within the expected deviations due to noise fluctuations. Thus, our

<figure><img src="image_4.png" /><figcaption>Fig. 2. Projected PSF-convolved SZ model map (top) and PSF convolved X-ray SB (middle) and temperature (bottom) maps. The contours for the models are overlaid to improve the visual representation of the maps. The simulated PSFs used for the SZ and X-ray maps correspond to 1.6' and 10'' FWHM, respectively. To ensure accurate PSF convolution, we binned the pixels in the maps such that the FWHM of each instrument's PSF is covered by at least three pixels.</figcaption></figure>

fitting methodology is able to reliably determine the underlying shape and thermodynamics of the observed mock galaxy cluster.

The use of both SZ and X-ray data in our analysis allows us to measure the three-dimensional geometry of the ICM distribution by constraining the elongation parameter (Sect. 2.1), since the two observational probes redundantly measure the thermodynamic properties of the gas along the line of sight. However, it should be noted that there may be degeneracies in determining cluster shape through this multi-probe approach depending on the relative orientation of the geometry, especially in inferring the geometric parameters of the three-dimensional structure, as discussed in Sereno (2007). These degeneracies can cause bias in the recovered shape parameters along with multi-modality in the posterior distributions.

The model fitting technique has been validated using the reference mock dataset, as shown in Fig. 3. Additionally, we present
<figure><img src="image_5.png" /><figcaption>**Fig. 3.** Posterior distributions estimated from our MCMC for a mock observation generated from a smooth model. The vertical green lines in each plot indicate the input parameters used to generate the mock observation maps, while the vertical red lines represent the median value from the accepted MCMC samples. The values displayed above each histogram show the median of the distribution, along with the $1\sigma$ (68.3%) credible region, which is indicated by dashed vertical lines in every plot. Additionally, the solid black line in the two-dimensional distributions encloses the 68% credible region for the parameter pairs. As highlighted by Vikhlinin et al. (2006) and Nagai et al. (2007), there is a correlation between the model parameters related to the ICM radial profiles, and individual parameter values exhibit a degeneracy. However, our objective is to ascertain smooth analytic functions that accurately represent the electron density and pressure profiles, thereby providing a comprehensive description of the ICM thermodynamics.</figcaption></figure>

the results of parameter estimation across a broader range of parameters. This analysis particularly concentrates on geomet- ric parameters that define the cluster’s shape and its projected appearance in the sky. In Table 2 and Fig. 4, we provide con- straints for the geometrical parameters of a triaxial ellipsoid ($q_{ICM,1}$, $q_{ICM,2}$, and $\cos\theta$, $e_\parallel$). These constraints are derived from various model fits, each with different values of $q_{ICM,1}$ and $\cos\theta$, while keeping the other parameters constant. We used the model parameters from the previous fit (Fig. 3) as a reference for these analyses.

In general, the input model parameters are successfully recovered within the test parameter ranges, although the uncer- tainties in the inclination angle increase as the major axis of the triaxial ellipsoid approaches a perpendicular orientation rel- ative to the observer’s line of sight. Most of the input values fall within the 68% confidence interval obtained from our fits, while the remaining values are included within the 95% confi- dence interval. Considering the expected fluctuations caused by the random noise added to the input maps, this suggests that any potential bias in our fitting method is minimal compared
<table><caption>**Table 2.** Model parameters and their values estimated using mock data.</caption><thead><tr><th colspan="8">Model parameters (input/triaxial fit results)</th></tr><tr><th colspan="2">q<sub>ICM,1</sub></th><th colspan="2">q<sub>ICM,2</sub></th><th colspan="2">cos θ</th><th colspan="2">e<sub>||</sub></th></tr><tr><th>Input</th><th>Fitted</th><th>Input</th><th>Fitted</th><th>Input</th><th>Fitted</th><th>Input</th><th>Fitted</th></tr></thead><tbody><tr><td>0.60</td><td>0.60<sup>+0.02</sup><sub>-0.05</sub></td><td>0.75</td><td>0.69<sup>+0.12</sup><sub>-0.05</sub></td><td>0.0</td><td>0.14<sup>+0.11</sup><sub>-0.10</sub></td><td>0.71</td><td>0.66<sup>+0.14</sup><sub>-0.07</sub></td></tr><tr><td>0.60</td><td>0.61<sup>+0.02</sup><sub>-0.04</sub></td><td>0.75</td><td>0.86<sup>+0.10</sup><sub>-0.15</sub></td><td>0.2</td><td>0.17<sup>+0.18</sup><sub>-0.13</sub></td><td>0.72</td><td>0.83<sup>+0.11</sup><sub>-0.13</sub></td></tr><tr><td>0.60</td><td>0.61<sup>+0.04</sup><sub>-0.04</sub></td><td>0.75</td><td>0.75<sup>+0.09</sup><sub>-0.08</sub></td><td>0.4</td><td>0.34<sup>+0.16</sup><sub>-0.19</sub></td><td>0.77</td><td>0.75<sup>+0.09</sup><sub>-0.06</sub></td></tr><tr><td>0.60</td><td>0.62<sup>+0.03</sup><sub>-0.04</sub></td><td>0.75</td><td>0.79<sup>+0.11</sup><sub>-0.09</sub></td><td>0.6</td><td>0.52<sup>+0.14</sup><sub>-0.15</sub></td><td>0.86</td><td>0.85<sup>+0.08</sup><sub>-0.10</sub></td></tr><tr><td>0.60</td><td>0.60<sup>+0.02</sup><sub>-0.04</sub></td><td>0.75</td><td>0.65<sup>+0.05</sup><sub>-0.04</sub></td><td>0.8</td><td>0.75<sup>+0.06</sup><sub>-0.05</sub></td><td>1.02</td><td>0.96<sup>+0.09</sup><sub>-0.06</sub></td></tr><tr><td>0.60</td><td>0.59<sup>+0.03</sup><sub>-0.04</sub></td><td>0.75</td><td>0.72<sup>+0.04</sup><sub>-0.06</sub></td><td>1.0</td><td>0.95<sup>+0.04</sup><sub>-0.05</sub></td><td>1.33</td><td>1.28<sup>+0.07</sup><sub>-0.09</sub></td></tr><tr><td>0.15</td><td>0.13<sup>+0.03</sup><sub>-0.02</sub></td><td>0.75</td><td>0.80<sup>+0.09</sup><sub>-0.10</sub></td><td>0.8</td><td>0.75<sup>+0.05</sup><sub>-0.05</sub></td><td>0.59</td><td>0.45<sup>+0.13</sup><sub>-0.07</sub></td></tr><tr><td>0.30</td><td>0.29<sup>+0.03</sup><sub>-0.02</sub></td><td>0.75</td><td>0.71<sup>+0.07</sup><sub>-0.09</sub></td><td>0.8</td><td>0.81<sup>+0.06</sup><sub>-0.08</sub></td><td>0.87</td><td>0.87<sup>+0.12</sup><sub>-0.14</sub></td></tr><tr><td>0.45</td><td>0.44<sup>+0.03</sup><sub>-0.04</sub></td><td>0.75</td><td>0.68<sup>+0.09</sup><sub>-0.08</sub></td><td>0.8</td><td>0.80<sup>+0.08</sup><sub>-0.09</sub></td><td>0.98</td><td>0.96<sup>+0.11</sup><sub>-0.11</sub></td></tr><tr><td>0.60</td><td>0.60<sup>+0.02</sup><sub>-0.04</sub></td><td>0.75</td><td>0.65<sup>+0.05</sup><sub>-0.04</sub></td><td>0.8</td><td>0.75<sup>+0.06</sup><sub>-0.05</sub></td><td>1.02</td><td>0.96<sup>+0.09</sup><sub>-0.06</sub></td></tr></tbody></table>

to the measurement noise. In a subsequent paper, we will perform a more detailed exploration spanning the range of masses, redshifts, and data quality in the CHEX-MATE sample. In addition, we will fit mock observations of simulated clusters to determine the impact of deviations from triaxial symmetry due to, for example, sub-structures or mergers.

### 3. Application to CHEX-MATE data

In this section, we introduce the X-ray and SZ data collected from our program. We applied the triaxial analysis technique to analyze a CHEX-MATE galaxy cluster PSZ2 G313.33+61.13 (Abell 1689), and the cluster serves as an illustrative example to demonstrate the method.

#### 3.1. Data

Table 3 summarizes the SZ and X-ray data from CHEX-MATE available for our multiwavelength analysis of the ICM distribution. The foundation of our analysis is the 3 Ms XMM-Newton observing program CHEX-MATE (CHEX-MATE Collaboration 2021), from which we obtained two-dimensional X-ray SB and temperature maps produced using the Voronoi tessellation method (Cappellari & Copin 2003; Diehl & Statler 2006). The details of the image production are reported in Bartalucci et al. (2023) and Lovisari et al. (2024), and here we report briefly the main analysis steps.

The XMM-Newton observations of the clusters were obtained using the EPIC instrument (Strüder et al. 2001; Turner et al. 2001). To create the X-ray SB map, photon-count images in the [0.7–1.2] keV range were extracted from the data acquired using the MOS1, MOS2, and pn cameras on the instrument. The energy band was selected to optimize the contrast between the emission from the source and the background (Ettori et al. 2010). The images from all three cameras were combined to maximize the statistical significance while accounting for the energy-band responses. Additionally, the X-ray maps are instrumental-background subtracted and corrected for the exposure. Point sources are removed from the analysis (Ghirardini et al. 2019) by masking them with circular regions that appear as empty circular regions in the X-ray maps in Fig. 5.

<figure><img src="image_6.png" /><figcaption>**Fig. 4.** Summary of the model parameters inferred from the mock data, as detailed in Table 2. The summary includes comparisons of $q_{ICM,1}$ to $q_{ICM,12}$ (top) and $q_{ICM,2}$ to $\cos\theta$ (bottom), along with their estimated uncertainties. The crosses in the figures indicate the input model parameters.</figcaption></figure>

Furthermore, they are spatially binned to have at least 20 counts per bin using the Voronoi technique. X-ray temperature maps (Lovisari et al. 2024) were prepared in a similar manner for the data obtained in the [0.3–7] keV band, with background modeling (Lovisari & Reiprich 2019) and spectral fitting performed. The fitting procedure to ascertain the temperature was done utilizing XSPEC (Arnaud 1996), which was employed to minimize the modified Cash statistics (Cash 1979) with the assumption of Asplund et al. (2009) metallicity. Subsequently, Voronoi-binned
<table><caption>**Table 3.** X-ray and SZ observation data and instruments used for the analysis of PSZ2 G313.33+61.13.</caption><thead><tr><th>Wavelength</th><th>Type</th><th>Instrument</th><th>Reference</th></tr></thead><tbody><tr><td>X-ray</td><td>Surface brightness (SB)</td><td>XMM-Newton</td><td>Bartalucci et al. (2023)</td></tr><tr><td>X-ray</td><td>Temperature</td><td>XMM-Newton</td><td>Lovisari et al. (2024)</td></tr><tr><td>mm-wave</td><td>SZ y-map</td><td>Planck <sup>(a)</sup></td><td>Planck Collaboration XXII (2016)</td></tr><tr><td>mm-wave</td><td>SZ y-map</td><td>ACT (ACTPol) <sup>(b)</sup></td><td>Madhavacheril et al. (2020)</td></tr></tbody></table>

Notes. The publicly released datasets of the surveys are available in the links below. <sup>(a)</sup>https://heasarc.gsfc.nasa.gov/W3Browse/all/plancksz2.html. <sup>(b)</sup>https://github.com/ACTCollaboration/DR4_DR5_Notebooks.

maps were generated to achieve a high signal-to-noise ratio (~30) for each cell.

Planck SZ maps are available for all of the CHEX-MATE galaxy clusters by definition (Planck Collaboration Int. V 2013; Pointecouteau et al. 2021). From these data we generated a custom y-map using the Modified Internal Linear Component Algorithm (MILCA; Hurier et al. 2013) with an improved angular resolution of 7' FWHM compared to the one publicly released by Planck with an angular resolution of 10' FWHM (Planck Collaboration XXII 2016). Also, ground-based SZ observations from cosmic microwave background surveys, including the ACT and the South Pole Telescope (SPT; Bleem et al. 2022)<fnref n="8" />, as well as the Caltech Submillimeter Observatory (CSO) Bolocam galaxy cluster archive<fnref n="9" /> (Sayers et al. 2013a), provide higher angular resolution data for a subset of CHEX-MATE clusters. Some of these ground-based data are currently publicly accessible, while others are slated for future release.

In this demonstration paper, we make use of the ACT SZ component-separated maps. The recent data release 4 (DR4) from the ACT provides component-separated maps, one of which is the SZ (Aiola et al. 2020; Madhavacheril et al. 2020). These maps were generated by analyzing data from a 2100 sq. deg area of the sky, captured using the ACTPol receiver (Henderson et al. 2016) at 98 and 150 GHz. These data offer more than four times finer angular resolution compared to the Planck map. Then, the maps were jointly analyzed and combined with Planck data. Rather than using the noise estimate provided with these data, which is quantified as a two-dimensional power spectral density, we instead followed an approach based on the recent analysis of similar joint ACT and Planck maps in Pointecouteau et al. (2021). Specifically, we randomly sampled 10 000 maps, ensuring that their size aligns with that of the input SZ data, in the corresponding ACT region (for instance, the region designated as “BN” for the cluster Abell 1689 analyzed in the next section). Then, we computed the covariance using these maps to estimate the noise covariance matrix. The resulting noise rms for the y-map is approximately ~9 × 10⁻⁶ per 0.5' square pixel, and the diagonal elements of the noise covariance matrix are shown along with the y-map in Fig. 5.

## 3.2. PSZ2 G313.33+61.13 (Abell 1689)

Using the datasets described above, we demonstrate our fitting method for PSZ2 G313.33+61.13 (Abell 1689), which is a Tier-2 cluster in the CHEX-MATE sample located at z = 0.1832 with a Planck SZ estimated mass of $M_{500} = 8.77 \times 10^{14} M_{\odot}$. We note the lensing mass measurement of the cluster is ~70% higher than

the Planck hydrostatic mass estimate (see Umetsu et al. 2015). We conducted a triaxial fit aligning the model center with the X-ray peak (Bartalucci et al. 2023). For a morphologically regular cluster, like Abell 1689, any deviations or offsets between the SZ and X-ray measurements are expected to have minimal impact on the overall model fit. The Planck + ACT SZ y-maps, along with the XMM-Newton X-ray SB and temperature maps, are shown in Fig. 5. Maps of the rms noise for each observable are also included, and indicate that the cluster is imaged at high signal to noise.

In the XMM-Newton maps presented, we masked bright point sources and excluded these areas from our X-ray analysis. As outlined in Ghirardini et al. (2019) and Bartalucci et al. (2023), the X-ray images were prepared to ensure a consistent cosmic X-ray background flux across the entire field of view during the point source subtraction process. Consequently, X-ray sources that fall below the masking threshold are not expected to have a significant impact on the model fit.

For the SZ data, we identified radio and dusty point sources using the NRAO VLA Sky Survey (NVSS; Condon et al. 1998) at 1.4 GHz and the Herschel SPIRE Point Source Catalog. Our search was centered on the region with $r \le R_{500}$ around Abell 1689 (RA 197.860 deg, Dec -1.336 deg). Within this area, we found five NVSS radio sources, with the brightest emitting 59.6 mJy. To infer their flux densities at 150 GHz, we identified counterparts in the OVRO/BIMA 30 GHz catalog (Coble et al. 2007) and calculated the spectral index $\alpha$ (flux density $S \propto \nu^{\alpha}$), following the methodology described in Sayers et al. (2013b). For dusty sources, we performed a similar analysis and fitted a gray-body spectrum ($S \propto \nu^{\beta}B_{\nu}(\nu, T)$) to 14 sources that were identified in the 500, 350, and 250 µm bands of the Herschel catalogs (Sayers et al. 2013c). These radio and dusty sources' flux densities, when extrapolated to 150 GHz, were all found to be below ~2 mJy, which is the typical rms noise level of the ACT maps (Madhavacheril et al. 2020). Specifically, the flux densities of the five radio sources are below 0.16 mJy with the mean of 0.06 mJy, and the 14 dusty sources have an average flux density of 1.04 mJy with a standard deviation of 0.3 mJy, with the brightest at 1.63 mJy. Consequently, given that none of the point sources are found to be brighter than the map noise rms, we determined that point source contamination in the SZ map is negligible, allowing us to use the data without additional masking.

This particular cluster was chosen for our demonstration because its triaxial shape has been well studied in the literature (Morandi et al. 2011; Sereno & Umetsu 2011; Sereno et al. 2012; Umetsu et al. 2015). For example, Sereno et al. (2012) performed a gas-only analysis using radial profiles of the X-ray and SZ observations from Chandra, XMM-Newton, and WMAP, along with various ground-based SZ facilities, and constrained the shape and orientation of the cluster's triaxial model with $q_{ICM,1} = 0.70 \pm 0.15$, $q_{ICM,2} = 0.81 \pm 0.16$, and $\cos\theta = 0.70 \pm 0.29$. A subsequent study by Umetsu et al. (2015) presented a

<footnote n="8">⁸ https://pole.uchicago.edu/public/data/sptsz_ymap/</footnote>

<footnote n="9">⁹ https://irsa.ipac.caltech.edu/data/Planck/release_2/ancillary-data/bolocam/bolocam.html</footnote>
<figure><img src="image_7.png" /></figure>

**Fig. 5.** SZ and X-ray maps of the CHEX-MATE galaxy cluster PSZ2 G313.33+61.13 (Abell 1689). The ACT Compton-y map (top left) and its error map (bottom left), the X-ray SB map (top middle) and its error map (bottom middle), and the X-ray temperature map (top right) and its error map (bottom right) are shown. The ACT SZ map is one of the component-separated map products that were produced using the internal linear combination method and combined with data from *Planck* (i.e., this is a joint map from ACT + *Planck* data; Madhavacheril et al. 2020). The X-ray maps are the data products from Bartalucci et al. (2023) and Lovisari et al. (2024). Bright point sources in the X-ray SB maps are indicated by white circles and were masked in the analysis, and the same point-source regions were also removed from the spectral analysis to obtain the temperatures. The regions are excluded from the model fit. The X-ray SB maps were binned using a Voronoi algorithm to ensure an adequate number of photon counts per bin, with smaller bin sizes used in the central region where the count rates are higher. For the temperature maps, Voronoi binning was similarly applied using a fixed signal-to-noise ratio of 30 instead of a fixed number of counts, ensuring a roughly uniform statistical uncertainty per bin. In the X-ray and SZ maps, red circles indicate the two-dimensional map regions included in our analysis. We incorporated a circular region with a radius of *r* = *R*<sub>500</sub> around the galaxy cluster center for the SZ and X-ray data by applying a circular mask to the maps, and the radius is shown as red circles. For this particular cluster, *R*<sub>500</sub> is equal to 7.42' (~1.37 Mpc). In the X-ray SB and its corresponding error maps, dashed red circles represent the region that encompasses 80% of the emission in the SB map, which is where the two-dimensional (inner region) and one-dimensional (outer region) analyses are separated, and it is located at *r* = 2.58'. We will explore and implement a comparable approach that integrates both two- and one-dimensional analysis techniques for temperature data, as described by Lovisari et al. (2024), in a forthcoming analysis.

combined multiwavelength analysis that included lensing data,
with the inferred ICM distribution being $q_{\text{ICM},1} = 0.60 \pm 0.14$,
$q_{\text{ICM},2} = 0.70 \pm 0.16$. Their derived value of $\cos\theta$, obtained
from the combined lensing and X-ray/SZ analysis, was found
to be $0.93 \pm 0.06$. The large $\cos\theta$ suggests that the major axis
of the triaxial ellipsoid ($\lambda_3^{\text{int}}$ in Fig. 1) is closely aligned with the
observer's line of sight.

Figure 6 shows the posterior of the model parameters that
describe our triaxial fit of PSZ2 G313.33+61.13, using the data
from *Planck*, ACT, and *XMM-Newton*. We find axial ratios of
$q_{\text{ICM},1} = 0.65 \pm 0.02$ and $q_{\text{ICM},2} = 0.79 \pm 0.02$. These values are
consistent with previous results, but an order of magnitude more
precise (Table 4). Our fits indicate the major axis of Abell 1689 is
almost perfectly aligned with the line of sight, with $\cos\theta \ge 0.96$
at 90% confidence. While previous works also indicated such
an alignment, a much wider range of orientations were allowed
in those fits. We note that our analysis only includes statistical
uncertainties on the fit, and the uncertainty due to data cali-
bration is not taken into account here. Also, as the elongation
parameter (Eq. (6)), which is the ratio of the size of the ellipsoid
along the observed line of sight to the major axis of the projected

ellipse in the sky plane, quantifies the three-dimensional geom-
etry of the triaxial ellipsoid model of the ICM. We thus present
constraints on $e_\parallel$ rather than on $\varphi$ and $\psi$. The inferred $e_\parallel$ is well
constrained in the fit to a value of $1.24 \pm 0.03$ and is consistent
with the gas analysis result of Sereno et al. (2012), who found
$e_\Delta = 0.66 \pm 0.21$, which corresponds to $1.15 \le e_\parallel \le 2.22$.
Figure 7 shows the reconstructed SZ, X-ray SB and temper-
ature maps of PSZ2 G313.33+61.13, incorporating the instru-
ment response, generated using the recovered parameters from
Fig. 6. The difference map, which is created by taking the input
data and subtracting the reconstructed model from it, reveals that
the majority of the pixels exhibit relative errors that are spread
within a range of $\pm 4\sigma$ (Fig. 6). The residuals for the SZ, X-ray
SB, and X-ray temperatures are distributed around zero. Their
respective standard deviations are equivalent to $1.5\sigma$, $0.6\sigma$, and
$1.1\sigma$ when fitted by a Gaussian.

For comparison, we performed an additional X-ray + SZ fit
using only the *Planck* SZ data, without incorporating the ground-
based ACT data. We obtain posteriors that significantly devi-
ate from our baseline fit with ACT data. We attributed this
to the coarse angular resolution of *Planck*, which prevents it
<figure><img src="image_8.png" /><figcaption>**Fig. 6.** Posterior distribution of the model fit parameters for the galaxy cluster PSZ2 G313.33+61.13 obtained using SZ data from *Planck* and ACT, as well as X-ray data from *XMM-Newton*. The vertical red lines indicate the median value from the accepted MCMC samples, with values displayed along with their 68% credible regions above each histogram. Instead of the Euler angles $\phi$ and $\psi$, we present $e_\parallel$, which is a function of five geometric parameters of a triaxial ellipsoid (Eq. (6)).</figcaption></figure>

<table><caption>**Table 4.** Parameters describing the triaxial geometry of PSZ2 G313.33+61.13 (Abell 1689).</caption><thead><tr><td>$q_{ICM,1}$</td><td>$q_{ICM,2}$</td><td>$\cos \theta$</td><td>$e_\parallel$</td><td>Reference</td></tr></thead><tbody><tr><td>0.70 ± 0.15</td><td>0.81 ± 0.16</td><td>0.70 ± 0.29</td><td>1.68 ± 0.53</td><td>Sereno et al. (2012)</td></tr><tr><td>0.60 ± 0.14</td><td>0.70 ± 0.16</td><td>0.93 ± 0.06</td><td>1.16 ± 0.10</td><td>Umetsu et al. (2015)</td></tr><tr><td>0.65 ± 0.02</td><td>0.79 ± 0.02</td><td>≥0.96 <sup>(a)</sup></td><td>1.24 ± 0.03</td><td>This work</td></tr></tbody></table>

**Notes.** <sup>(a)</sup> At 90% confidence.

from resolving morphological features given the angular size of Abell 1689 at $z = 0.1832$. To test this, we generated two sets of mock observations using the recovered parameters from our baseline fit to the observed data from both *Planck* and ACT (along with *XMM-Newton*). One mock was based on the proper-

ties of the y-map from *Planck* + ACT, while the other mimicked the y-map with only *Planck* SZ data, including the appropriate noise and PSF shape for each case. Our fit to the mock mul- tiwavelength data with the *Planck* + ACT y-map yields recov- ered parameters closely aligned with the input model, suggesting
<figure><img src="image_9.png" /><figcaption>**Fig. 7.** Reconstructed SZ and X-ray models of PSZ2 G313.33+61.13 generated using the recovered parameters from Fig. 6 (top). The difference between the observational data and the reconstructed model map above, in units of pixel-based error (middle). The histogram of the distribution of the relative error in the middle panels (bottom). The X-ray SB histogram takes into account both the residuals in the inner two-dimensional region, which includes 80% of the emissions observed in the data, and the outer map region where we implemented one-dimensional analysis using azimuthal medians (see Fig. 5). In all cases, the residuals are distributed within ±4σ level compared to the error. When the relative errors of the SZ, X-ray SB, and X-ray temperature are modeled with a Gaussian fit, their standard deviations align with 1.5σ, 0.6σ, and 1.1σ respectively. Comparatively, 63.8%, 64.4%, and 100% of the data from SZ, X-ray SB, and X-ray temperature maps respectively exhibit pixel-based signal-to-noise ratios exceeding 1.5.</figcaption></figure>

these data can accurately recover the input ICM shape. In contrast, the second mock observation based on the *Planck*-only y-map produces a set of parameters significantly deviating from the input. This suggests that the SZ data from *Planck* alone are insufficient to reliably fit our triaxial model, at least for a galaxy cluster with this specific shape at this specific redshift. This confirms that our fit to observed data using the *Planck*-only y-map are likely biased. In a subsequent paper we will explore this issue in more detail, to better understand which types of galaxy clusters can (or cannot) be reliably reconstructed with the data available for CHEX-MATE.

Furthermore, in order to evaluate how the much higher overall signal to noise of the X-ray SB compared to the SZ and X-ray temperature impacts the results, we carried out an additional fit using the reduced $\chi^2$ for each of the three observables in order to weight them equally in the fit. The results of this fit indicate that there is only a minimal shift in the values of the derived geometric parameters based on this equal weighting of

the three observables. Specifically, in the reduced $\chi^2$ fit, $q_{ICM,1}$ has a value of 0.70 ± 0.04, $q_{ICM,2}$ is 0.78 ± 0.05, and $e_\|$ stands at 1.22 ± 0.07. We also attempted to account for fluctuations in the calibration uncertainty, which can be especially important for the temperature profile (e.g., Schellenberger et al. 2015; Wallbank et al. 2022). We conducted model fits by introducing an additional ~10% uncertainty of the temperature, but observed little changes in the parameters, with posteriors displaying similar levels of variation.

Table 4 presents a comparison of axial ratios, the (cosine of) inclination angle, and elongation parameter from our study. These metrics show notable improvement when compared to the gas analysis results reported by Sereno et al. (2012). This marked improvement is attributable to multiple factors: our use of deeper new *XMM-Newton* data not available to Sereno et al. (2012); our use of an *XMM-Newton* SB image rather than a shallower *Chandra* SB image; our use of much higher quality SZ data from *Planck* and ACT rather than from WMAP and
SZA/OVRO/BIMA; and our improved analysis formalism making use of fully two-dimensional images for all of the observables rather than a projected elliptically averaged profile of X-ray SB and temperature along with a single-aperture photometric measurement of the SZ signal. While a direct comparison between the current and previous ICM analysis of Abell 1689 is not straightforward due to data availability and differences in the ICM model, we performed an additional test. Instead of using the full two-dimensional images from the SZ and X-ray data, we processed the input data into one-dimensional projected profiles, which are more analogous to the data used in previous studies, and ran a model fit. The analysis yielded inferred geometrical parameters: $q_{ICM,1} = 0.61^{+0.14}_{-0.12}$, $q_{ICM,2} = 0.86^{+0.10}_{-0.14}$, $\cos\theta = 0.78^{+0.15}_{-0.39}$. These results align with the uncertainty levels found in the previous analysis, and suggest that much of the improvement in parameter constraints is related to the incorporation of higher quality, fully two-dimensional data available for our analysis.

As we will illustrate in subsequent studies, the derived geometric parameters of the ICM distribution, such as the elongation that quantifies the three-dimensional geometry, can be applied in conjunction with gravitational lensing measurements. For these fits, we will work under the assumption that the triaxial axes of the ICM and DM are co-aligned, but with axial ratios that are allowed to vary. The lensing analysis becomes crucial for discerning the triaxial shapes of DM, circumventing the need to rely on hydrostatic equilibrium or simulation-based corrections. Consequently, a comprehensive multi-probe analysis facilitates a characterization of the total matter distribution, which is essential for precise lensing-based mass calibrations (Sereno et al. 2018), along with allowing for a determination of the distribution of nonthermal pressure support (Sayers et al. 2021).

## 4. Conclusions

We have improved a multi-probe analysis package to fit the three-dimensional ellipsoidal shapes of CHEX-MATE galaxy clusters. This package builds upon CLUMP-3D (Sereno et al. 2017), which has been employed to analyze the triaxial shapes of CLASH clusters (Sereno et al. 2018; Chiu et al. 2018; Sayers et al. 2021). Specifically, we made the following improvements: (1) we model two-dimensional distributions of the SZ and X-ray temperature data, in contrast to the one-dimensional azimuthally averaged profiles in these quantities used by Sereno et al. (2017), (2) we parametrize electron density and pressure rather than density and temperature, reducing the number of parameters and speeding up the fit, and (3) we have ported the code to Python to facilitate a future public release. For the two-dimensional map analyses, we have added the capability to include publicly available SZ data from ground-based cosmic microwave background surveys such as that conducted with ACT, in addition to the default Planck SZ maps.

We verified the triaxial analysis method through mock data analysis and applied it to the actual CHEX-MATE galaxy cluster, PSZ2 G313.33+61.13 (Abell 1689). The analysis effectively constrains the model geometry, in particular at the few percent level for the axial ratios. Our results are consistent with previous analyses of Abell 1689 available in the literature. Specifically, we find axial ratios of $q_{ICM,1} = 0.65 \pm 0.02$ and $q_{ICM,2} = 0.79 \pm 0.02$ and elongation parameter $e_\parallel = 1.24 \pm 0.03$. Compared to the similar gas-only analysis using X-ray and SZ data presented in Sereno et al. (2012), the axial ratios and elongation parameters in our study demonstrate a substantial improvement, with uncertainties an order of magnitude lower. Our results indicate

that Abell 1689 has axial ratios typical of what is expected for the general population of galaxy clusters (Jing & Suto 2002; Lau et al. 2011), but a remarkably close alignment between the major axis and the line of sight. This alignment has resulted in exceptional lensing properties for Abell 1689, such as an abundance of strong lensing features (e.g., Broadhurst et al. 2005; Limousin et al. 2007), one of the largest Einstein radii ever observed (47", Coe et al. 2010), and an extremely large concentration of mass when fitted to a spherically symmetric model ($c_{vir} = 12.8^{+3.1}_{-2.4}$ or $c_{200} = 10.2^{+2.6}_{-1.9}$, Umetsu et al. 2011; Umetsu 2020). We thus conclude that there is nothing unusual about the triaxial shape of Abell 1689, other than its orientation. In addition, the estimated axial ratios of the cluster yield a triaxiality parameter $t = 0.66$ (Franx et al. 1991). While the incorporation of lensing data is necessary for a direct quantitative comparison with DM axial ratios, the calculated $t$ classifies this halo as being close to the “prolate” population that comprises ~80% of the total cluster fraction in the DM-only simulations (Vega-Ferrero et al. 2017). The integration of lensing data for a comprehensive multiwavelength analysis, as well as the public release of the software and data products, will be addressed in subsequent papers of this series.

**Acknowledgements.** J.K. and J.S. were supported by NASA Astrophysics Data Analysis Program (ADAP) Grant 80NSSC21K1571. J.K. is supported by a Robert A. Millikan Fellowship from the California Institute of Technology (Caltech). J.K. is supported by Korea Advanced Institute of Science & Technology (KAIST) research fund for new faculty settlement. M.S. acknowledges financial contribution from contract ASI-INAF n.2017-14-H.0 and from contract INAF mainstream project 1.05.01.86.10. M.E.D. acknowledges partial support from the NASA ADAP, primary award to SAO with a subaward to MSU, SV9-89010. S.E., F.G., and M.R. acknowledge the financial contribution from the contracts Prin-MUR 2022 supported by Next Generation EU (M4.C2.1.1, n. 20227RNLY3 The concordance cosmological model: stress-tests with galaxy clusters), ASI-INAF Athena 2019-27-HH.0, “Attività di Studio per la comunità scientifica di Astrofisica delle Alte Energie e Fisica Astroparticellare” (Accordo Attuativo ASI-INAF n. 2017-14-H.0), and from the European Union’s Horizon 2020 Programme under the AHEAD2020 project (grant agreement n. 871158). This research was supported by the International Space Science Institute (ISSI) in Bern, through ISSI International Team project #565 (Multi-Wavelength Studies of the Culmination of Structure Formation in the Universe). A.I., E.P., and G.W.P. acknowledge support from CNES, the French space agency. K.U. acknowledges support from the National Science and Technology Council of Taiwan (grant 109-2112-M-001-018-MY3) and from the Academia Sinica (grants AS-IA-107-M01 and AS-IA-112-M04). B.J.M. acknowledges support from STFC grant ST/V000454/1. L.L. acknowledges financial contribution from the INAF grant 1.05.12.04.01. M.G. acknowledges partial support from the ERC Consolidator Grant BlackHoleWeather (101086804).

## References

Aiola, S., Calabrese, E., Maurin, L., et al. 2020, JCAP, 2020, 047

Allen, S. W., Evrard, A. E., & Mantz, A. B. 2011, ARA&A, 49, 409

Ansarifard, S., Rasia, E., Biffi, V., et al. 2020, A&A, 634, A113

Arnaud, K. A. 1996, in Astronomical Data Analysis Software and Systems V, eds. G. H. Jacoby, & J. Barnes, ASP Conf. Ser., 101, 17

Arnaud, M., Pratt, G. W., Piffaretti, R., et al. 2010, A&A, 517, A92

Asplund, M., Grevesse, N., Sauval, A. J., & Scott, P. 2009, ARA&A, 47, 481

Baldi, A., Ettori, S., Molendi, S., & Gastaldello, F. 2012, A&A, 545, A41

Bartalucci, I., Molendi, S., Rasia, E., et al. 2023, A&A, 674, A179

Battaglia, N., Bond, J. R., Pfrommer, C., & Sievers, J. L. 2012, ApJ, 758, 75

Becker, M. R., & Kravtsov, A. V. 2011, ApJ, 740, 25

Binggeli, B. 1980, A&A, 82, 289

Binney, J. 1985, MNRAS, 212, 767

Bleem, L. E., Crawford, T. M., Ansarinejad, B., et al. 2022, ApJS, 258, 36

Broadhurst, T., Benítez, N., Coe, D., et al. 2005, ApJ, 621, 53

Campitiello, M. G., Ettori, S., Lovisari, L., et al. 2022, A&A, 665, A117

Cappellari, M., & Copin, Y. 2003, MNRAS, 342, 345

Cash, W. 1979, ApJ, 228, 939

CHEX-MATE Collaboration (Arnaud, M., et al.) 2021, A&A, 650, A104

Chiu, I. N., Umetsu, K., Sereno, M., et al. 2018, ApJ, 860, 126

Coble, K., Bonamente, M., Carlstrom, J. E., et al. 2007, AJ, 134, 897
Coe, D., Benítez, N., Broadhurst, T., & Moustakas, L. A. 2010, ApJ, 723, 1678

Condon, J. J., Cotton, W. D., Greisen, E. W., et al. 1998, AJ, 115, 1693

Corless, V. L., & King, L. J. 2007, MNRAS, 380, 149

Davis, M., Efstathiou, G., Frenk, C. S., & White, S. D. M. 1985, ApJ, 292, 371

De Filippis, E., Sereno, M., Bautz, M. W., & Longo, G. 2005, ApJ, 625, 108

Despali, G., Giocoli, C., Bonamigo, M., Limousin, M., & Tormen, G. 2017, MNRAS, 466, 181

Diehl, S., & Statler, T. S. 2006, MNRAS, 368, 497

Eckert, D., Roncarelli, M., Ettori, S., et al. 2015, MNRAS, 447, 2198

Eckert, D., Finoguenov, A., Ghirardini, V., et al. 2020, Open J. Astrophys., 3, 12

Ettori, S., Morandi, A., Tozzi, P., et al. 2009, A&A, 501, 61

Ettori, S., Gastaldello, F., Leccardi, A., et al. 2010, A&A, 524, A68

Euclid Collaboration (Adam, R., et al.) 2019, A&A, 627, A23

Foreman-Mackey, D., Hogg, D. W., Lang, D., & Goodman, J. 2013, PASP, 125, 306

Franx, M., Illingworth, G., & de Zeeuw, T. 1991, ApJ, 383, 112

Ghirardini, V., Eckert, D., Ettori, S., et al. 2019, A&A, 621, A41

Ghizzardi, S., Molendi, S., van der Burg, R., et al. 2021, A&A, 646, A92

Goodman, J., & Weare, J. 2010, Commun. Appl. Math. Comput. Sci., 5, 65

Henderson, S. W., Allison, R., Austermann, J., et al. 2016, J. Low Temp. Phys., 184, 772

Ho, S., Bahcall, N., & Bode, P. 2006, ApJ, 647, 8

Hogg, D. W., & Foreman-Mackey, D. 2018, ApJS, 236, 11

Huijser, D., Goodman, J., Brewer, B. J., et al. 2015, arXiv e-prints [arXiv:1509.02230]

Hurier, G., Macías-Pérez, J. F., & Hildebrandt, S. 2013, A&A, 558, A118

Ivezić, Ž., Kahn, S. M., Tyson, J. A., et al. 2019, ApJ, 873, 111

Jing, Y. P., & Suto, Y. 2002, ApJ, 574, 538

Kazantzidis, S., Kravtsov, A. V., Zentner, A. R., et al. 2004, ApJ, 611, L73

Khatri, R., & Gaspari, M. 2016, MNRAS, 463, 655

Kravtsov, A. V., & Borgani, S. 2012, ARA&A, 50, 353

Lau, E. T., Nagai, D., Kravtsov, A. V., & Zentner, A. R. 2011, ApJ, 734, 93

Lau, E. T., Hearin, A. P., Nagai, D., & Cappelluti, N. 2021, MNRAS, 500, 1029

Lima, M., & Hu, W. 2005, Phys. Rev. D, 72, 043006

Limousin, M., Richard, J., Jullo, E., et al. 2007, ApJ, 668, 643

Limousin, M., Morandi, A., Sereno, M., et al. 2013, Space Sci. Rev., 177, 155

Lovisari, L., & Reiprich, T. H. 2019, MNRAS, 483, 540

Lovisari, L., Ettori, S., Rasia, E., et al. 2024, A&A, 682, A45

Madhavacheril, M. S., Hill, J. C., Næss, S., et al. 2020, Phys. Rev. D, 102, 023534

Mantz, A. B., Allen, S. W., Morris, R. G., et al. 2016, MNRAS, 463, 3582

Markevitch, M., & Vikhlinin, A. 2007, Phys. Rep., 443, 1

Mazzotta, P., Rasia, E., Moscardini, L., & Tormen, G. 2004, MNRAS, 354, 10

McNamara, B. R., & Nulsen, P. E. J. 2007, ARA&A, 45, 117

Meneghetti, M., Rasia, E., Merten, J., et al. 2010, A&A, 514, A93

Morandi, A., Limousin, M., Rephaeli, Y., et al. 2011, MNRAS, 416, 2567

Nagai, D., Kravtsov, A. V., & Vikhlinin, A. 2007, ApJ, 668, 1

Navarro, J. F., Frenk, C. S., & White, S. D. M. 1996, ApJ, 462, 563

Navarro, J. F., Frenk, C. S., & White, S. D. M. 1997, ApJ, 490, 493

Oguri, M., Lee, J., & Suto, Y. 2003, ApJ, 599, 7

Oguri, M., Takada, M., Okabe, N., & Smith, G. P. 2010, MNRAS, 405, 2215

Okabe, T., Nishimichi, T., Oguri, M., et al. 2018, MNRAS, 478, 1141

Oppizzi, F., De Luca, F., Bourdin, H., et al. 2023, A&A, 672, A156

Planck Collaboration XXII. 2016, A&A, 594, A22

Planck Collaboration XXVII. 2016, A&A, 594, A27

Planck Collaboration Int. V. 2013, A&A, 550, A131

Pointecouteau, E., Santiago-Bautista, I., Douspis, M., et al. 2021, A&A, 651, A73

Postman, M., Coe, D., Benítez, N., et al. 2012, ApJS, 199, 25

Pratt, G. W., Böhringer, H., Croston, J. H., et al. 2007, A&A, 461, 71

Pratt, G. W., Arnaud, M., Biviano, A., et al. 2019, Space Sci. Rev., 215, 25

Predehl, P., Andritschke, R., Arefiev, V., et al. 2021, A&A, 647, A1

Reese, E. D., Kawahara, H., Kitayama, T., et al. 2010, ApJ, 721, 653

Rozo, E., & Rykoff, E. S. 2014, ApJ, 783, 80

Sayers, J., Czakon, N. G., Mantz, A., et al. 2013a, ApJ, 768, 177

Sayers, J., Mroczkowski, T., Czakon, N. G., et al. 2013b, ApJ, 764, 152

Sayers, J., Mroczkowski, T., Zemcov, M., et al. 2013c, ApJ, 778, 52

Sayers, J., Sereno, M., Ettori, S., et al. 2021, MNRAS, 505, 4338

Sayers, J., Mantz, A. B., Rasia, E., et al. 2023, ApJ, 944, 221

Schellenberger, G., Reiprich, T. H., Lovisari, L., Nevalainen, J., & David, L. 2015, A&A, 575, A30

Sereno, M. 2007, MNRAS, 380, 1207

Sereno, M., & Umetsu, K. 2011, MNRAS, 416, 3187

Sereno, M., Jetzer, P., & Lubini, M. 2010, MNRAS, 403, 2077

Sereno, M., Ettori, S., & Baldi, A. 2012, MNRAS, 419, 2646

Sereno, M., Ettori, S., Meneghetti, M., et al. 2017, MNRAS, 467, 3801

Sereno, M., Umetsu, K., Ettori, S., et al. 2018, ApJ, 860, L4

Stapelberg, S., Tchernin, C., Hug, D., Lau, E. T., & Bartelmann, M. 2022, A&A, 663, A17

Stark, A. A. 1977, ApJ, 213, 368

Strüder, L., Briel, U., Dennerl, K., et al. 2001, A&A, 365, L18

Suto, D., Peirani, S., Dubois, Y., et al. 2017, PASJ, 69, 14

Turner, M. J. L., Abbey, A., Arnaud, M., et al. 2001, A&A, 365, L27

Umetsu, K. 2020, A&ARv, 28, 7

Umetsu, K., Broadhurst, T., Zitrin, A., Medezinski, E., & Hsu, L.-Y. 2011, ApJ, 729, 127

Umetsu, K., Sereno, M., Medezinski, E., et al. 2015, ApJ, 806, 207

Umetsu, K., Sereno, M., Tam, S.-I., et al. 2018, ApJ, 860, 104

Vega-Ferrero, J., Yepes, G., & Gottlöber, S. 2017, MNRAS, 467, 3226

Vikhlinin, A., Kravtsov, A., Forman, W., et al. 2006, ApJ, 640, 691

Voit, G. M. 2005, Rev. Mod. Phys., 77, 207

Wallbank, A. N., Maughan, B. J., Gastaldello, F., Potter, C., & Wik, D. R. 2022, MNRAS, 517, 5594

Zhan, H., & Tyson, J. A. 2018, Rep. Progr. Phys., 81, 066901

1 California Institute of Technology, 1200 E. California Blvd., MC 367-17, Pasadena, CA 91125, USA

2 Department of Physics, Korea Advanced Institute of Science and Technology (KAIST), 291 Daehak-ro, Yuseong-gu, Daejeon 34141, Republic of Korea
e-mail: junhan@kaist.ac.kr

3 INAF, Osservatorio di Astrofisica e Scienza dello Spazio, Via Piero Gobetti 93/3, 40129 Bologna, Italy

4 INFN, Sezione di Bologna, Viale Berti Pichat 6/2, 40127 Bologna, Italy

5 INAF, Istituto di Astrofisica Spaziale e Fisica Cosmica di Milano, Via A. Corti 12, 20133 Milano, Italy

6 Department of Astronomy, University of Geneva, Ch. d'Ecogia 16, 1290 Versoix, Switzerland

7 INAF – Osservatorio Astronomico di Brera, Via E. Bianchi 46, 23807 Merate, LC, Italy

8 Università degli Studi di Roma “Tor Vergata”, Via della Ricerca Scientifica 1, 00133 Roma, Italy

9 INFN, Sezione di Roma “Tor Vergata”, Via della Ricerca Scientifica 1, 00133 Roma, Italy

10 Dipartimento di Fisica, Sapienza Università di Roma, Piazzale Aldo Moro 5, 00185 Roma, Italy

11 Department of Physics and Astronomy, Michigan State University, 567 Wilson Road, East Lansing, MI 48864, USA

12 Department of Astrophysical Sciences, Princeton University, Princeton, NJ 08544, USA

13 Laboratoire d'Astrophysique de Marseille, Aix-Marseille Univ., CNRS, CNES, Marseille, France

14 Institut d'Astrophysique de Paris, UMR 7095, CNRS & Sorbonne Université, 98 bis Boulevard Arago, 75014 Paris, France

15 Université Paris-Saclay, Université Paris Cité, CEA, CNRS, AIM, 91191 Gif-sur-Yvette, France

16 Jodrell Bank Centre for Astrophysics, Department of Physics and Astronomy, The University of Manchester, Manchester M13 9PL, UK

17 Center for Astrophysics | Harvard & Smithsonian, 60 Garden Street, Cambridge, MA 02138, USA

18 HH Wills Physics Laboratory, University of Bristol, Bristol BS8 1TL, UK

19 Physics Program, Graduate School of Advanced Science and Engineering, Hiroshima University, Hiroshima 739-8526, Japan

20 Hiroshima Astrophysical Science Center, Hiroshima University, 1-3-1 Kagamiyama, Higashi-Hiroshima, Hiroshima 739-8526, Japan

21 Core Research for Energetic Universe, Hiroshima University, 1-3-1, Kagamiyama, Higashi-Hiroshima, Hiroshima 739-8526, Japan

22 IRAP, Université de Toulouse, CNRS, CNES, UT3-UPS, Toulouse, France

23 Academia Sinica Institute of Astronomy and Astrophysics (ASIAA), No. 1, Section 4, Roosevelt Road, Taipei 10617, Taiwan

24 Department of Physics, Informatics and Mathematics, University of Modena and Reggio Emilia, 41125 Modena, Italy