# On the nature of the missing mass of galaxy clusters in MOND: the view from gravitational lensing

Benoit Famaey, Lorenzo Pizzuti, Ippocratis D Saltas

► To cite this version:

Benoit Famaey, Lorenzo Pizzuti, Ippocratis D Saltas. On the nature of the missing mass of galaxy clusters in MOND: the view from gravitational lensing. Physical Review D, 2025, 111 (12), pp.123042. 10.1103/dccw-srks. hal-04748532

HAL Id: hal-04748532

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

Submitted on 15 Oct 2025

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
# On the nature of the missing mass of galaxy clusters in MOND: the view from gravitational lensing

Benoit Famaey

Université de Strasbourg, CNRS,
Observatoire astronomique de Strasbourg,
UMR 7550, F-67000 Strasbourg, France

Lorenzo Pizzuti

Dipartimento di Fisica G. Occhialini,
Università degli Studi di Milano Bicocca,
Piazza della Scienza 3, 20126 Milano, Italy

Ippocratis D. Saltas

CEICO, Institute of Physics, Czech Academy of Sciences

Na Slovance 2, 182 21 Praha 8, Czech Republic

(Dated: April 30, 2025)

Modified Newtonian Dynamics (MOND) has long been known to fail in galaxy clusters, implying a residual missing mass problem for clusters in this context. Here, using mass profiles derived from strong- and weak-lensing shear, as well as magnification data, for a sample of clusters from the CLASH survey, we characterize the density profile of this residual MOND missing mass in the central Mpc of galaxy clusters. In line with results obtained in the literature from the hydrostatic equilibrium of hot gas, we find that an inner constant density core and an outer power-law slope sharper than -3.5 provides a good description within ~ 1 Mpc. We also show that the data in the central parts of clusters can be even better represented by a 'dark mass-follows-gas' profile with an exponential cut-off. Clusters in the sample display a remarkable uniformity for the missing-to-hot-gas density ratio in the inner parts, of order ~10, as well as for the exponential cut-off radius, of order ~ 400 kpc. These lensing results can in principle serve as a crucible for relativistic theories of MOND in galaxy clusters, or for any other tentative hypothesis regarding the nature of the clusters residual missing mass in the MOND context.

## I. INTRODUCTION

Dark matter (DM) plays a central role in the standard ACDM cosmological model. It is considered as a meaningful explanation for a large variety of discrepancies between astronomical measurements and the predictions of General Relativity (GR) based on baryonic matter alone. However, the inference of its exact distribution and properties depends crucially on the validity of GR on all scales considered. Various modifications of gravity (see, e.g., [1]) can change the distribution of DM inferred from observational data. In principle, it is even possible to do fully away with DM in galaxies via an appropriate modification of Newtonian dynamics [MOND, 2–4]. The idea is that Newtonian dynamics, the weak-field limit of GR, would break down at ultra-low accelerations typical of the outskirts of galaxies (below $a_0 \approx 10^{-10}$ m/s$^2$) and that the acceleration would fall as $r^{-1}$ instead of $r^{-2}$ in this regime. Such a modification of gravity would imply, observationally, a strong correlation between the Newtonian gravity generated by the baryons and the total gravitational field in stellar systems.

Actually, in rotationally-supported galaxies, observations indicate that the baryonic surface density of galaxies (and therefore the Newtonian gravity generated by the baryons) plays a central role in setting the inner shape of their observed rotation curves (supposedly set by the profile of their DM halo in the GR context), as if baryons were fully dominating even when they are not: the associated diversity of rotation curve

shapes remains challenging to understand in the standard context [5–8], while it is a natural expectation in the context of MOND, see e.g. [3]. At the same time, the baryonic surface density plays absolutely no role in setting the tight relation [9–11] between the baryonic mass $M_b$ and the asymptotic circular speed $V_f$, known as the baryonic Tully-Fisher relation (BTFR), $M_b \propto V_f^4$. The BTFR seems to be surprisingly valid at very large distances from galaxies, as inferred from weak lensing measurements from the Kilo-Degree Survey (KIDS) around both late-type and early-type galaxies [12]. This striking contrast between the central role of surface density in setting observationally the inner shape of galaxy rotation curves and its absence of role in setting the BTFR is not well understood in the standard context, but it is a prediction of MOND. More generally, the close relation observed between the gravitational acceleration generated by baryons and the total one, known as the Radial Acceleration Relation [13, 14], is a direct prediction of MOND too. All this explains why MOND is still taken seriously today in the context of galaxy dynamics, although it is difficult to reconcile with the success of our current standard model of cosmology. MOND-based modifications of gravity also typically require a new scale in addition of acceleration to pass Solar System constraints [15].

For galaxy clusters, where the intra-cluster gas is the dominant baryonic component, the BTFR predicted by MOND can in principle be broadly translated into a baryonic mass ($M_b$)-temperature ($T$) relation, $M_b \propto T^2$. This slope of 2 is reasonably consistent with the data [16], but the normalization pre-
dicted by MOND is completely off. In other words, MOND as presently formulated, fails in galaxy clusters, where there would remain a residual missing mass problem. This missing mass corresponds to a factor of a few globally, but it can rise to very high values in the central parts of clusters. This need for residual missing matter in MOND has long been known [e.g., 17–25], and the gravitational lensing from the Bullet cluster [26, 27] has additionally shown that this residual DM appears to behave in a collisionless way, namely, like usual DM. It would therefore be tempting to combine MOND with a form of (hot) non-baryonic dark matter, such as sterile neutrinos of a few eV, which would naturally be relevant only on galaxy cluster scales and above. However, simulations of structure formation in this MOND-hot dark matter context have shown that it inevitably leads to a considerable overproduction of massive galaxy clusters [28], which severely reduces the appeal of the idea. Another solution may come from the embedding of MOND in a covariant context. Indeed, the MOND behavior, if fundamentally true, should be the weak-field limit of a relativistic theory; such relativistic theories most often imply additional fields to the usual metric [e.g., 29–31], allowing to reproduce MOND in static configurations, whilst preserving the same relation between the lensing potential and gravitational potential as in GR. In recent versions of such relativistic theories, these additional fields can also play the role of DM on cosmic scales [32]: whether or not such fields could also help explain away the residual missing mass in galaxy clusters within such a theoretical framework is, however, an open question. Preliminary studies in this context have e.g. demonstrated a mass-dependent enhancement of the gravitational boost followed by an oscillation at large radii, quite different from the MOND non-relativistic limit for clusters [33]. A related possibility would be that the $a_0$ scale is itself dependent on other quantities, such as the depth of the potential well which clearly separates galaxy clusters from galaxies [34]; however, this would not a priori explain why the residual missing mass would behave in a collisionless way in colliding clusters. Finally, the residual missing mass in the MOND context could simply indicate that MOND actually predicts additional baryonic mass (in a collisionless form) in the center of galaxy clusters, for instance cold gas clouds [35] in a multiphase intergalactic medium. In the latter case, one could expect that the distribution of the missing mass would more or less follow that of the observed hot gas in galaxy clusters.

Irrespective of the possible solution, or absence thereof, we aim hereafter to check whether the residual missing mass profile of galaxy clusters in a MOND context has some universal behavior, and whether it correlates in some ways with the baryonic content of the clusters. Such a profile could then serve as a crucible for relativistic theories of MOND in galaxy clusters or any other tentative hypothesis regarding this residual missing mass in the MOND context. In Section II, we present the theoretical framework and the methods, before presenting the results based on the analysis of 16 clusters from the Cluster Lensing And Supernova survey with Hubble (CLASH) survey in Section III, and concluding in Section IV. In the Appendix, we analyze the residual mass in MOND using kinematics data and discuss why such data yield results

which appear to be in tension with lensing.

## II. THEORY AND METHODS

### A. Gravitational lensing in MOND

Our starting point is any covariant theory [e.g., 29–32] for which the gravitational lensing potential in the weak-field limit is $\Phi_{\text{lens}} \equiv \Phi - \Psi = 2\Phi$, where $\Phi$ is the gravitational potential governing dynamics, and $\Psi$ is the relativistic potential; in other words, theories for which $\Psi = -\Phi$ in the static weak-field limit.

In the classical limit, for a large set of particles of mass $m_i$ moving in a gravitational field generated by the matter density distribution $\rho = \sum_i m_i \delta(\mathbf{x} - \mathbf{x}_i)$, and described by the potential $\Phi$, we have the following action:

$$S = S_{\text{matter}} + S_{\text{gravity}} \qquad (1)$$

$$S_{\text{matter}} = \int d^3 x dt \left( \frac{1}{2} \dot{\mathbf{x}}^2 - \Phi \right) \rho, \qquad (2)$$

and in the Newtonian limit of GR, $S_{\text{gravity}} = -(1/8\pi G) \int d^3x dt(\nabla\Phi)^2$. There is a wide variety of ways to modify this classical gravitational action to get a MOND-like behaviour at low accelerations [36, 37], sometimes only on certain scales depending on the number of other dimensioned constants than acceleration involved or the number of fields involved [38, 39]. It is common nowadays to rely on the so-called “quasi-linear MOND” formulation. This formulation makes the non-linear character of the theory transparent through the introduction of $\Phi_N$ as an auxiliary potential, and is particularly convenient to implement in numerical simulations [e.g. 40–42]. It is important to note here that the conclusions of the present study, assuming spherical symmetry for galaxy clusters, will be largely independent of that particular choice of theory, apart from the fact that only the MOND acceleration constant should play a role (and should remain a constant). It might well not be the case in some recent relativistic formulations of MOND [see, e.g., 33], but it is nevertheless the zeroth order expectation for any relativistic theory boiling down to MOND in the quasi-static weak field limit. The assumption of spherical symmetry, in and by itself, could slightly bias the results, since most known clusters are not strictly speaking spherically symmetric. In this sense, our analysis is not better (but also not worse) than those made within GR with such an assumption [see, e.g., Ref. 3, for an assessment of the moderate effect of geometry on MOND predictions in the context of flat galaxy disks]. Let us also note that our main interest lies in the central parts of clusters where that approximation is reasonable at leading order, at least for clusters in equilibrium.

For the sake of clarity, we briefly remind the basics of the quasi-linear MOND formulation below. The total action is, as usual, a sum of two parts – the matter and the gravity part, respectively – but now we have:

$$S_{\text{gravity}} = -\frac{1}{8\pi G} \int d^3x dt \left( 2\nabla\Phi \cdot \nabla\Phi_N - a_0^2 Q(y^2) \right), \qquad (3)$$
where $y \equiv |\nabla\Phi_N|/a_0$, and $Q(y^2)$ is in principle a non-linear function, which needs to satisfy certain requirements, as we will discuss below. Varying the matter action with respect to the matter density $\rho$, one recovers Newton's second law $\ddot{x} = -\nabla\Phi$. Therefore, $\Phi$ plays the role of the total gravitational potential (to which matter couples), sourced by the total matter density $\rho$. Variation of the action with respect to the potential $\Phi$ yields the classical Poisson equation for baryonic matter as

$$ \nabla^2 \Phi_N = 4\pi G\rho. \qquad (4) $$

The auxiliary gravitational potential $\Phi_N$ is therefore just the usual Newtonian potential sourced by $\rho$, without any modification. The remaining step is the variation of the gravity action with respect to $\Phi_N$, which yields the modified Poisson equation for the total gravitational field to which matter actually couples,

$$ \nabla^2 \Phi = \nabla \cdot (Q'(y^2)\nabla\Phi_N), \qquad (5) $$

with $Q'(y^2) \equiv \frac{d}{ds^2}Q(y^2)$. The non-linearity of the right-hand side of (5) is necessary to account for the “dark-matter-like” effects at galactic scales. The core idea of MOND is that

$$ \begin{aligned} \nabla\Phi &\rightarrow \nabla\Phi_N \text{ for } \nabla\Phi_N \gg a_0, \text{ and} \\ \nabla\Phi &\rightarrow \sqrt{a_0\nabla\Phi_N} \text{ for } \nabla\Phi_N \ll a_0, \end{aligned} \qquad (6) $$

therefore, the function $Q$ just needs to satisfy

$$ Q(y^2) \rightarrow y^2 \text{ for } y^2 \gg 1 \text{ and } Q(y^2) \rightarrow \frac{4}{3}(y^2)^{3/4} \text{ for } y^2 \ll 1. \qquad (7) $$

In this paper, we will use the phenomenologically successful ‘simple’ function for $Q'$ (see [3] and references therein),

$$ Q'(y^2) = \frac{1}{2} \left( 1 + \left( 1 + 4y^{-1} \right)^{\frac{1}{2}} \right). \qquad (8) $$

The general solution of Poisson equation is then equivalent to equation(6) in the asymptotic regimes, up to a curl field correction, and is precisely equal to it in a spherically symmetric system, which we will assume for galaxy clusters.

In GR, we can relate the lensing potential and its source through

$$ \nabla^2 \Phi_{\text{lens}} \equiv \nabla^2 (\Phi - \Psi) = 8\pi G\rho_{\text{GR}}, \qquad (9) $$

where $\rho_{\text{GR}}$ refers to any type of matter, baryonic or dark. Hence in spherical symmetry, we have for the enclosed mass

$$ M_{\text{lensGR}}(r) \equiv \frac{r^2}{2G} \frac{d}{dr} (\Phi - \Psi) = \frac{r^2}{G} \frac{d\Phi(r)}{dr}. \qquad (10) $$

Lensing observations measure the combined potential $\Phi - \Psi$, which in turn can yield a reconstruction of $M_{\text{lensGR}}$. To find the equivalent source mass for a spherically symmetric system in MOND, which we define as $M(r)$, one can use equation (5),

which implies that $\nabla\Phi = Q'(y^2)\nabla\Phi_N$, and then simply solve for the enclosed mass $M(r)$ in

$$ M_{\text{lensGR}}(r) = \frac{r^2}{G} Q'(y^2) \frac{d\Phi_N(r)}{dr} \equiv Q'(y^2)M(r). \qquad (11) $$

Note that $M(r) < M_{\text{lensGR}}(r)$, or in other words, the enclosed mass needed to produce a given lensing potential is smaller in MOND than in GR. Finally,

$$ M(r) = 4\pi \int dr \, r'^2 \rho(r'), \qquad (12) $$

where $\rho$ is the matter density needed to reproduce the observations in MOND, part of which is certainly of baryonic nature, and the rest is the missing mass we are trying to reconstruct here. Note that the assumption of $\Psi = -\Phi$ for the (unspecified) relativistic MOND theory under consideration implies that observations of the galactic velocity field or the X-ray profile of the intracluster gas, providing measurements of the potential $\Phi$ as a function of distance from the cluster's center, will yield a dynamical mass equivalent to the lensing mass by construction, in MONDian gravity as well as in GR and Newtonian gravity, where $M_{\text{dynN}} = M_{\text{lensGR}}$. Importantly, since we assume spherical symmetry, the results of the present analysis are also valid for the quadratic Lagrangian [43] classical limit of modified gravity MOND, since the curl field correction to equation(6) then vanishes in that framework as well as in quasi-linear MOND.

## B. Missing mass profiles in MOND

In a first step, we will reconstruct the enclosed mass and density profiles, $M(r)$ and $\rho(r)$, needed to produce MONDian effects in accordance with lensing observations in galaxy clusters. In practice, assuming spherical symmetry, we rely in the following on the (parametric) estimate of $M_{\text{lensGR}}(r)$ from lensing observations, and compute the enclosed mass profile needed to reproduce the same lensing signal in MOND, $M(r)$, which can be compared to the observed enclosed baryonic mass profile of the cluster. If MOND would not need any residual missing mass in clusters, the two quantities would be identical.

From the MONDian lensing mass $M(r)$, we then subtract the observed gas mass $M_{\text{gas}}(r)$, the stellar mass of the Brightest Cluster Galaxy (BCG) and companion galaxies, $M^*_{\text{BCG}}(r)$, which provide the dominant baryonic component at $r < 0.05$ Mpc, and the mass of the other galaxies of the cluster $M_*(r)$. We do not consider the mass of the Intra-cluster Light (ICL), i.e. stars in the cluster region which are not associated to galaxies, since the contribution to the total stellar mass is small, and is already partially accounted for, through the stellar mass of the BCG (see e.g. [44, 45]). Then, at each radial bin $r$ we can define the missing mass up to radius $r$, and the associated density, as

$$ \begin{aligned} \delta M(r) &= M(r) - M_{\text{gas}}(r) - M_*(r) - M_{\text{BCG}}^*(r), \\ \rho_{\delta M}(r) &= \frac{1}{4\pi r^2} \frac{d[\delta M(r)]}{dr}, \end{aligned} \qquad (13) $$
where $M(r)$ is the total mass predicted by MOND, $M_{gas}(r)$ is the contribution of the X-ray emitting hot gas of the ICM, $M_*(r)$ is the stellar mass in galaxies (except for the BCG), and $M_{BCG}^*(r)$ is the stellar mass of the BCG $^1$.

It is important to notice here that the density profile above may become negative beyond a given radius, depending on the slope of the mass profile computed from equation (11). We will call that point the turnaround radius ($r_0$), and we will always truncate the profile beyond that point in order to preserve the physical character of the density if it is to be made of matter (see, however, [33]). This feature is related, by construction, to the parametric form of the Newtonian/GR mass profile adopted to adjust the lensing observations, and we do not attempt here to determine whether it is an inescapable consequence of cluster lensing. Studying whether the lensing signal in the outskirts of clusters might be reconciled with the isolated MOND expectation, or whether the influence of the external gravitational field [25] could influence those outskirts, or whether it would indicate a sharper drop as in [33], would need a fully non-parametric approach which should be the topic of further works. We can formally define the turnaround radius $r_0$ as the radius at which

$$ \left. \frac{d[\delta M(r)]}{dr} \right|_{r=r_0} = 0. \qquad (14) $$

We are interested in understanding the scaling of the missing mass/density with radius, and the existence of possible universality features in it. In this regard, we will fit the missing mass with a generic density profile, namely, the “$\alpha-\beta-\gamma$” model (e.g., [47]),

$$ \rho_{\delta M}(r) = \frac{\rho_s}{\left(\frac{r}{r_s}\right)^{\gamma} \left[1 + \left(\frac{r}{r_s}\right)^{\alpha} \right]^{(\beta-\gamma)/\alpha}}, \qquad (15) $$

where $\rho_s$ and $r_s$ represent the characteristic density and scale radius of the profile, $\gamma$ represents the (logarithmic) inner slope, $\beta$ represents the outer slope, and $\alpha$ determines the sharpness of the transition between the inner and outer slopes, which we here fix to 1 and do not fit$^2$. We will first fit all other parameters together for each cluster, and then fix the slopes in the search for possible global scaling relations between the characteristic density, scale radius and other observables. Finally,

we will explore whether a “dark mass-follows-gas” profile – which might be expected if the residual missing mass is made of missing baryons – could equally well or even better represent the data.

### C. Lensing data

Our analysis is based on a sample of 16 clusters from the CLASH survey, reconstructed through strong+weak lensing observations [48, 49]. For a review on weak-lensing of galaxy clusters we refer to [50]. The CLASH sample is made of massive and (mostly) relaxed clusters, typically with an X-ray gas temperature larger than 5 keV and displaying a high degree of dynamical relaxation. No lensing data were used in the initial selection of these clusters to avoid any preferential alignment along the line of sight. A few clusters however do show some signs of departures from symmetric X-ray surface brightness distributions. In [51], strong lensing, weak lensing shear, and magnification data for these clusters were combined by integrating data from Hubble observations and ground-based wide-field multicolor imaging, primarily using Suprime-Cam on the Subaru telescope. [51] fitted the lensing data of the CLASH sample by using a Navarro-Frenk-White profile [52] (NFW hereafter), providing the best fit parameters and the uncertainties, in which the relevant observational uncertainties were all taken care of. The error analysis accounts for cosmic noise (i.e cosmic shear acting as noise), the scatter in the lensing profile due to projection effects, as well as for observational uncertainties [50]. It should be noticed that the latter modelling and fitting was performed in the context of General Relativity. Hereafter, the MCMC chains for the NFW fits are used together with the chains for the Burkert profile [53] and the Hernquist profile [54] fits$^3$.

## III. RESULTS

### A. Residual missing mass profiles

Given the lensing observations, we start from the mass profile, $M_{lensGR}(r)$ for each cluster, and by virtue of equation (11) we reconstruct the profile of the total “baryonic” mass needed in MOND, $M(r)$. To account for possible biases in the parametric modeling of the profiles, we consider three ansatz for the GR lensing mass profile, namely the NFW, the Burkert and the Hernquist profiles, and we then average over them$^4$.

<footnote n="1">¹ The profiles for the gas in our sample were provided by Vincenzo Salzano after request, and are computed following the prescription described in [46]. They were double-checked in the Chandra Data Archive, and are provided in the file ClusterInfo.py at https://zenodo.org/records/15299349. We model the minor baryonic component $M_*(r)$ in our sample as $M_*(r) = M_{gas}(r) \times f_{gal}(r/r_{200})$, where $f_{gal}(x)$ is assumed to be the same for all clusters and normalized on our chosen reference cluster MACS 1206. We have further checked that ignoring the mass profiles of the galaxies beyond the BCG entirely does, in fact, not alter significantly any of the final outcomes of our analysis.</footnote>

<footnote n="2">² This assumption is based on a preliminary analysis in which all parameters were treated as free variables. We observed no significant improvement in the fit when optimizing $\alpha$ alongside the other parameters, with very poor BIC. Furthermore, the best-fit values of $\alpha$ for our cluster sample are centered around $\langle\alpha\rangle = 1.02$. Given these preliminary results, we adopted the simpler $\beta-\gamma$ model with $\alpha$ fixed to 1.</footnote>

<footnote n="3">³ The MCMC chains have been provided by Keiichi Umetsu as a private communication</footnote>

<footnote n="4">⁴ We run also the lensing analysis assuming a generalized NFW (gNFW), characterized by an additional free slope parameter with respect to the NFW. However, the uncertainties in this case are extremely large, especially in the cluster core; moreover the additional parameter does not provide any significant improvement on the Posterior. Thus we discarded this model from the averaging procedure.</footnote>
We finally subtract the various observed baryonic contributions as explained in Section II B, see equation (13), to obtain the MOND residual missing mass profiles $\delta M(r)$. These are plotted for each individual cluster on Figure 1.

We then proceed considering the following cases for parametrizing the missing mass profiles:

* **Case 1 ($\beta - \gamma$):** We fit the missing mass $\delta M$ with the $\alpha-\beta-\gamma$ model of equation (15) up to the turnaround radius of each cluster, fixing $\alpha = 1$ after our initial search. The main result of this fit is a very tight constraint on $\gamma$. Averaging over all clusters, we find an average value of

$$ \langle \gamma \rangle = 0.015, \qquad (16) $$

and a weighted geometric mean $< 10^{-3}$, with very low scatter when accounting for the addition in quadrature of the scatter among the sample and the intrinsic uncertainties of the fit for each cluster. This means that the inner slope is very close to a core (which would correspond to $\gamma = 0$). No relevant features have been found for $\gamma$ in relation to the gas mass or any other cluster property. This fitting procedure, however, does not yield a proper upper limit for the outer slope, for which the one sigma upper limit increases together with the size of the uniform prior. The lower limit however never goes below $\beta = 3.7$, and the best-fit outer slope is typically extremely sharp ($\langle \beta \rangle \sim 6$), in line with the general trends previously reported by, e.g., [25]. Concerning the two other parameters, the scale radius $r_s$ and characteristic density $\rho_s$, it is worth noticing that they exhibit a quite strong correlation, even though with large error bars due to the poor constraint on the upper limit of $\beta$. We will explore this correlation hereafter, once fixing all slopes. The main conclusion of this analysis is that the profiles are typically cored ($\gamma = 0$), and have a sharp outer slope with $\beta > 3.7$.

* **Case 2 (fixing all slopes):** Given the results found above, we experiment with a general profile with a lower number of parameters for the fitting mass profile, keeping $\alpha = 1$, and further fixing $\gamma = 0$ (cored inner profile), while considering separately the cases $\beta = 4, 5, 6, 10$ for the outer slope. In each case, only the characteristic (central) density and scale radius are fitted. This case gives a much more favourable Bayesian Information Criterion (BIC hereafter, see e.g. [55]) than the previous one, due to the removal of one free parameter, and it is equivalent to the fits performed by [25], who found a very slight preference for $\beta = 4$ based on hydrostatic equilibrium of hot gas. In our case, the resulting fits to the data also showed no strong statistical preference for any of the four values above. The BIC, however, gives a slight preference for $\beta = 6$.

In particular, we found an average $\Delta$BIC = 0.42 between $\beta = 4$ and $\beta = 6$, in favor of the latter, with a maximum difference of 0.85, always below 1. For the cases $\beta = 5$ and $\beta = 10$, we obtained an average $\Delta$BIC, of 0.08 and 0.16 with $\beta = 6$, still in favor of the latter,

with a maximum difference of 0.23 and 0.74 respectively. The minimal difference with the results of [25] should therefore not be over-interpreted as it is statistically irrelevant, and could also simply result from clusters spanning a slightly different range of parameters (total gas mass, temperature, etc.). On the other hand, when fixing $\beta = 3$, the average $\Delta$BIC with $\beta = 6$ rises to 1.7, with a maximum value of 2.75, and is therefore significant, in accordance with our findings from Case 1. The values of the fitted parameters for each individual cluster when fixing $\beta = 4, 5, 6, 10$ are all given in Table I. We have also included, as a vertical tick on the abscissa, the fitted scale radius for the $\beta = 6$ case in Figure 1.

In Figure 2 we show, for all cases $\beta = 4, 5, 6$ and 10, the central density versus the scale radius, and how the scale radius is correlated with the observed gas mass $M_{\text{gas}}$ (computed within the radius where the enclosed density in GR – hence effective density in MOND – is 200 times the critical density of the universe at the cluster redshift). The latter correlation is particularly significant if we exclude the lowest gas-mass cluster from the sample (MS 2137). Concerning the correlation between the scale radius and the (log) characteristic density, already identified in Case 1, the slope of the relation steepens with increasing $\beta$. In particular, assuming a linear-log relation of the type $r_s = a\log(\rho_s) + b$, $a = -0.11 \pm 0.03, -0.38 \pm 0.04, -0.51 \pm 0.07, -1.07 \pm 0.14$ for $\beta = 4, 5, 6, 10$, respectively. The physical interpretation of this anti-correlation is that the total amount of residual missing mass always remains, within an order of magnitude, within a few $\sim 10^{14}$ solar masses, but can be more or less concentrated, hence with a higher characteristic density when the scale radius is smaller. Moreover, the overall value of $r_s$ tends to grow with $\beta$: as $\beta$ increases, a larger $r_s$ is needed to compensate the faster decrease of the density profile. Conversely, the overall value of $\rho_s$ decreases with $\beta$.

We also search the literature for signs of mild disequilibria, and find that half (8/16) of the sample displays such signs. Such disequilibria should not in principle affect the lensing mass: indeed, such clusters, displayed as black star symbols on Figure 2, do not appear to follow different scaling relations than the other half of the sample.

As a final visual illustration of the parametrized fits performed hereinabove, we provide in Fig. 3, as an example, a comparison of the best-fit profile for the missing mass in the cluster Abell 209, for the four values of $\beta$ considered here. This shows how close these different fits are to each other, and nicely illustrates how the scale radius increases with $\beta$.

* **Case 3 ($\beta$ as a function of $\gamma$):** Finally, in order to check that the above results, especially on the central constant density core, are fully independent from the assumption of a fixed $\beta$, we now consider $\gamma$ as a free parameter with
<figure><img src="image_2.png" /><figcaption>FIG. 1: MOND residual missing mass profiles $M(r)$ for all the galaxy clusters in the sample (black solid lines) together with their 68% C.L. (gray-shaded areas). The red solid lines indicate the gas mass profile following a beta-model (see the file ClusterInfo.py at https://zenodo.org/records/15299349). We overplot as blue-dashed curves the best fit from the ‘dark mass follows gas’ case of equation (17). The vertical blue ticks on the abscissa denote, on the other hand, the best fit scale radius $r_s$ with $\beta = 6$, from Sect. III A, Case 2.</figcaption></figure>

$\gamma \in [0, 0.6]$, and a linear dependence of $\beta$ on $\gamma$ according to $\beta = 4 + 3 \times \gamma$ or $\beta = 8 - 5 \times \gamma$. Fitting these two profiles to all clusters, we confirm that $\gamma \to 0$ is favored for all clusters, in both cases. These profiles actually give a BIC significantly worse than the models with $\alpha = 1, \gamma = 0$ and a fixed $\beta$ between 4 and 10 (due to the presence of an additional free parameter). The values of $\beta$ close to 4 and 8 found from those fits come from $\gamma$ being close to zero and from the fixed relation between $\beta$ and $\gamma$. As discussed earlier, for fixed values of the parameters, there is not statistically meaningful difference in the quality of the fits between $\beta = 4$ and $\beta = 5, 6$ or 10. No relevant features have been found for $\gamma$ in relation to the gas mass.

## B. Does the missing mass follow the gas density ?

So far, we have re-constructed the missing mass $\delta M$ in a systematic way, without caring as to what its underlying nature could be. We found, in accordance with the literature, that the generic missing mass profile is close to following a cored density profile in the central parts followed by a sharp decline with a power-law slope sharper than 3.5 in the outskirts ($r \ge 0.5$ Mpc). If this missing mass would be baryonic, for instance in the form of cold gas clouds, then one might perhaps expect a close correlation with the distribution of the hot gas in the central parts of the cluster. Visually, it indeed appears in Figure 1 that $\delta M$ seems to roughly scale with that of the enclosed gas mass in the central parts.

To check for this possible behaviour in a more quantitative
<figure><img src="image_3.png" /></figure>

$\beta = 5, \gamma = 0$

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

$\beta = 6, \gamma = 0$

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

$\beta = 10, \gamma = 0$

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

FIG. 2: Plots of the parameter $r_s$ vs the observed gas mass (left) and vs Log $\rho_s$; from top to bottom: for Case 2, with $\beta = 4, 5, 6, 10$, as discussed in Section III A. The error-bars represent the 68% confidence intervals, while the blue dotted lines represent the best-fit linear regression of the $r_s$ vs Log $\rho_s$ relation. Black stars identify the clusters with signs of a disturbed state, while clusters at equilibrium are denoted by grey triangles. For more details on the analysis, see Section III.
<table><caption>TABLE I: Parameters of the MOND missing mass density profiles for the CLASH sample, adopting profiles with $\alpha = 1$, $\gamma = 0$, and different values of $\beta$. The characteristic density $\rho_s$ is in units of $M_{\odot}/\text{Mpc}^3$ and the scale radius $r_s$ in Mpc. Column two, three: $\beta = 4$. Column four, five: $\beta = 5$. Column six, seven: $\beta = 6$. Column eight, nine: $\beta = 10$. Constraints are given within $1\sigma$.</caption><thead><tr><th rowspan="2">Cluster</th><th colspan="2">β = 4</th><th colspan="2">β = 5</th><th colspan="2">β = 6</th><th colspan="2">β = 10</th></tr><tr><th>Log ρ<sub>s</sub></th><th>r<sub>s</sub></th><th>Log ρ<sub>s</sub></th><th>r<sub>s</sub></th><th>Log ρ<sub>s</sub></th><th>r<sub>s</sub></th><th>Log ρ<sub>s</sub></th><th>r<sub>s</sub></th></tr></thead><tbody><tr><td>MACS 0329</td><td>16.48<sup>+0.16</sup><sub>-0.13</sub></td><td>0.14<sup>+0.01</sup><sub>-0.03</sub></td><td>16.39 ± 0.14</td><td>0.22<sup>+0.03</sup><sub>-0.04</sub></td><td>16.31<sup>+0.13</sup><sub>-0.13</sub></td><td>0.30<sup>+0.03</sup><sub>-0.05</sub></td><td>16.24<sup>+0.13</sup><sub>-0.12</sub></td><td>0.61<sup>+0.07</sup><sub>-0.09</sub></td></tr><tr><td>MACS 1115</td><td>15.60<sup>+0.14</sup><sub>-0.12</sub></td><td>0.35<sup>+0.04</sup><sub>-0.06</sub></td><td>15.52<sup>+0.13</sup><sub>-0.12</sub></td><td>0.53<sup>+0.06</sup><sub>-0.08</sub></td><td>15.49<sup>+0.08</sup><sub>-0.11</sub></td><td>0.69<sup>+0.08</sup><sub>-0.07</sub></td><td>15.39<sup>+0.08</sup><sub>-0.14</sub></td><td>1.44<sup>+0.19</sup><sub>-0.15</sub></td></tr><tr><td>MS 2137</td><td>15.41<sup>+0.15</sup><sub>-0.15</sub></td><td>0.38<sup>+0.05</sup><sub>-0.08</sub></td><td>15.33<sup>+0.15</sup><sub>-0.13</sub></td><td>0.57<sup>+0.07</sup><sub>-0.11</sub></td><td>15.32<sup>+0.11</sup><sub>-0.18</sub></td><td>0.73<sup>+0.13</sup><sub>-0.12</sub></td><td>15.26<sup>+0.07</sup><sub>-0.16</sub></td><td>1.47<sup>+0.24</sup><sub>-0.15</sub></td></tr><tr><td>RXJ1347</td><td>15.79<sup>+0.12</sup><sub>-0.10</sub></td><td>0.41<sup>+0.04</sup><sub>-0.06</sub></td><td>15.71<sup>+0.12</sup><sub>-0.11</sub></td><td>0.61<sup>+0.06</sup><sub>-0.09</sub></td><td>15.68<sup>+0.09</sup><sub>-0.13</sub></td><td>0.81<sup>+0.10</sup><sub>-0.09</sub></td><td>15.63<sup>+0.10</sup><sub>-0.13</sub></td><td>1.60<sup>+0.22</sup><sub>-0.17</sub></td></tr><tr><td>RXJ2129</td><td>16.24 ± 0.17</td><td>0.14<sup>+0.02</sup><sub>-0.03</sub></td><td>16.15 ± 0.15</td><td>0.21<sup>+0.03</sup><sub>-0.04</sub></td><td>16.06<sup>+0.18</sup><sub>-0.15</sub></td><td>0.29<sup>+0.04</sup><sub>-0.06</sub></td><td>16.05<sup>+0.11</sup><sub>-0.16</sub></td><td>0.55<sup>+0.09</sup><sub>-0.07</sub></td></tr><tr><td>A209</td><td>15.50 ± 0.11</td><td>0.36<sup>+0.04</sup><sub>-0.05</sub></td><td>15.45 ± 0.11</td><td>0.55<sup>+0.06</sup><sub>-0.07</sub></td><td>15.38<sup>+0.15</sup><sub>-0.19</sub></td><td>0.72<sup>+0.19</sup><sub>-0.13</sub></td><td>15.36<sup>+0.08</sup><sub>-0.13</sub></td><td>1.38<sup>+0.15</sup><sub>-0.15</sub></td></tr><tr><td>A611</td><td>15.98<sup>+0.14</sup><sub>-0.12</sub></td><td>0.25<sup>+0.03</sup><sub>-0.04</sub></td><td>15.90<sup>+0.15</sup><sub>-0.11</sub></td><td>0.37<sup>+0.03</sup><sub>-0.06</sub></td><td>15.82<sup>+0.14</sup><sub>-0.12</sub></td><td>0.52<sup>+0.06</sup><sub>-0.08</sub></td><td>15.77<sup>+0.10</sup><sub>-0.12</sub></td><td>1.01<sup>+0.12</sup><sub>-0.10</sub></td></tr><tr><td>A2261</td><td>15.86<sup>+0.11</sup><sub>-0.12</sub></td><td>0.32 ± 0.04</td><td>15.76<sup>+0.10</sup><sub>-0.10</sub></td><td>0.48<sup>+0.05</sup><sub>-0.06</sub></td><td>15.70<sup>+0.09</sup><sub>-0.11</sub></td><td>0.66<sup>+0.08</sup><sub>-0.07</sub></td><td>15.64<sup>+0.08</sup><sub>-0.10</sub></td><td>1.31<sup>+0.13</sup><sub>-0.12</sub></td></tr><tr><td>MACS0416</td><td>15.82<sup>+0.14</sup><sub>-0.13</sub></td><td>0.23<sup>+0.03</sup><sub>-0.04</sub></td><td>15.76<sup>+0.12</sup><sub>-0.13</sub></td><td>0.34<sup>+0.04</sup><sub>-0.05</sub></td><td>15.70<sup>+0.11</sup><sub>-0.11</sub></td><td>0.46<sup>+0.05</sup><sub>-0.06</sub></td><td>15.58<sup>+0.07</sup><sub>-0.12</sub></td><td>0.98<sup>+0.13</sup><sub>-0.08</sub></td></tr><tr><td>MACS0429</td><td>16.11<sup>+0.18</sup><sub>-0.17</sub></td><td>0.19<sup>+0.03</sup><sub>-0.04</sub></td><td>16.04<sup>+0.16</sup><sub>-0.15</sub></td><td>0.28<sup>+0.04</sup><sub>-0.05</sub></td><td>16.00<sup>+0.17</sup><sub>-0.16</sub></td><td>0.38<sup>+0.06</sup><sub>-0.08</sub></td><td>15.97<sup>+0.13</sup><sub>-0.16</sub></td><td>0.72<sup>+0.11</sup><sub>-0.11</sub></td></tr><tr><td>MACS0647</td><td>16.12<sup>+0.14</sup><sub>-0.15</sub></td><td>0.21<sup>+0.03</sup><sub>-0.04</sub></td><td>16.01<sup>+0.15</sup><sub>-0.15</sub></td><td>0.34<sup>+0.04</sup><sub>-0.06</sub></td><td>15.93<sup>+0.16</sup><sub>-0.16</sub></td><td>0.47<sup>+0.07</sup><sub>-0.08</sub></td><td>15.88<sup>+0.14</sup><sub>-0.12</sub></td><td>0.94<sup>+0.10</sup><sub>-0.13</sub></td></tr><tr><td>MACS0744</td><td>16.01<sup>+0.13</sup><sub>-0.11</sub></td><td>0.27<sup>+0.03</sup><sub>-0.04</sub></td><td>15.92 ± 0.12</td><td>0.41<sup>+0.05</sup><sub>-0.06</sub></td><td>15.89<sup>+0.11</sup><sub>-0.11</sub></td><td>0.54<sup>+0.07</sup><sub>-0.07</sub></td><td>15.78<sup>+0.09</sup><sub>-0.17</sub></td><td>1.14<sup>+0.22</sup><sub>-0.14</sub></td></tr><tr><td>MACS1149</td><td>15.62 ± 0.13</td><td>0.39<sup>+0.05</sup><sub>-0.06</sub></td><td>15.53<sup>+0.11</sup><sub>-0.12</sub></td><td>0.59<sup>+0.07</sup><sub>-0.09</sub></td><td>15.51<sup>+0.10</sup><sub>-0.10</sub></td><td>0.76<sup>+0.08</sup><sub>-0.08</sub></td><td>15.51<sup>+0.05</sup><sub>-0.16</sub></td><td>1.43<sup>+0.25</sup><sub>-0.09</sub></td></tr><tr><td>MACS1206</td><td>15.96<sup>+0.17</sup><sub>-0.14</sub></td><td>0.26<sup>+0.03</sup><sub>-0.05</sub></td><td>15.87<sup>+0.14</sup><sub>-0.13</sub></td><td>0.40<sup>+0.05</sup><sub>-0.07</sub></td><td>15.85<sup>+0.13</sup><sub>-0.15</sub></td><td>0.52<sup>+0.08</sup><sub>-0.08</sub></td><td>15.77<sup>+0.09</sup><sub>-0.14</sub></td><td>1.07<sup>+0.14</sup><sub>-0.11</sub></td></tr><tr><td>MACS1720</td><td>16.05<sup>+0.13</sup><sub>-0.14</sub></td><td>0.23<sup>+0.03</sup><sub>-0.04</sub></td><td>15.94 ± 0.13</td><td>0.36<sup>+0.04</sup><sub>-0.06</sub></td><td>15.86<sup>+0.14</sup><sub>-0.12</sub></td><td>0.50<sup>+0.06</sup><sub>-0.08</sub></td><td>15.78<sup>+0.08</sup><sub>-0.15</sub></td><td>1.02<sup>+0.16</sup><sub>-0.09</sub></td></tr><tr><td>RXJ2248</td><td>15.94<sup>+0.13</sup><sub>-0.17</sub></td><td>0.24<sup>+0.04</sup><sub>-0.04</sub></td><td>15.88<sup>+0.20</sup><sub>-0.15</sub></td><td>0.36<sup>+0.04</sup><sub>-0.08</sub></td><td>15.90<sup>+0.16</sup><sub>-0.13</sub></td><td>0.45<sup>+0.05</sup><sub>-0.08</sub></td><td>15.81<sup>+0.07</sup><sub>-0.16</sub></td><td>0.93<sup>+0.13</sup><sub>-0.08</sub></td></tr></tbody></table>

<table><caption>TABLE II: Parameters of the CLASH clusters: name (column 1), temperature (column 2, from [56]), total gas mass (column 3), total MOND missing mass (column 4), missing-to-hot-gas density ratio (column 5), exponential cut-off radius (column 6), and flag for disequilibrium features (column 7, * symbol if signs of disequilibrium exist).</caption><thead><tr><th>Cluster</th><th>kT(keV)</th><th>M<sub>gas</sub>(10<sup>14</sup>M<sub>☉</sub>)</th><th>δM(r<sub>0</sub>)(10<sup>14</sup>M<sub>☉</sub>)</th><th>log(η)</th><th>r<sub>0</sub>/λ(Mpc)</th><th>disequilibrium</th></tr></thead><tbody><tr><td>MACS 0329</td><td>8.0 ± 0.5</td><td>1.13</td><td>1.79 ± 0.47</td><td>1.27<sup>+0.08</sup><sub>-0.07</sub></td><td>0.24 ± 0.05</td><td>* [57]</td></tr><tr><td>MACS 1115</td><td>8.0 ± 0.4</td><td>2.03</td><td>3.17 ± 1.09</td><td>0.66 ± 0.06</td><td>1.40<sup>+0.64</sup><sub>-0.53</sub></td><td></td></tr><tr><td>MS 2137</td><td>5.9 ± 0.3</td><td>0.77</td><td>2.73 ± 1.60</td><td>0.70<sup>+0.03</sup><sub>-0.04</sub></td><td>4.14<sup>+3.00</sup><sub>-0.76</sub></td><td>* [58]</td></tr><tr><td>RXJ1347</td><td>15.5 ± 0.6</td><td>5.28</td><td>7.71 ± 2.37</td><td>0.81 ± 0.05</td><td>1.61<sup>+0.77</sup><sub>-0.48</sub></td><td>* [59]</td></tr><tr><td>RXJ2129</td><td>5.8 ± 0.4</td><td>1.03</td><td>0.85 ± 0.27</td><td>1.06<sup>+0.11</sup><sub>-0.10</sub></td><td>0.23 ± 0.07</td><td></td></tr><tr><td>A209</td><td>7.3 ± 0.5</td><td>2.33</td><td>2.88 ± 1.07</td><td>0.91<sup>+0.08</sup><sub>-0.07</sub></td><td>0.60<sup>+0.14</sup><sub>-0.13</sub></td><td>* [60]</td></tr><tr><td>A611</td><td>7.90 ± 0.35</td><td>1.62</td><td>2.99 ± 1.02</td><td>1.12 ± 0.07</td><td>0.48<sup>+0.12</sup><sub>-0.11</sub></td><td></td></tr><tr><td>A2261</td><td>7.6 ± 0.3</td><td>3.32</td><td>4.46 ± 1.10</td><td>1.05<sup>+0.06</sup><sub>-0.05</sub></td><td>0.65 ± 0.14</td><td></td></tr><tr><td>MACS 0416</td><td>7.5 ± 0.8</td><td>1.58</td><td>1.61 ± 0.61</td><td>1.28<sup>+0.12</sup><sub>-0.11</sub></td><td>0.19 ± 0.04</td><td>* [61]</td></tr><tr><td>MACS 0429</td><td>6.00 ± 0.44</td><td>1.46</td><td>1.61 ± 0.62</td><td>1.03 ± 0.08</td><td>0.38<sup>+0.12</sup><sub>-0.10</sub></td><td></td></tr><tr><td>MACS 0647</td><td>13.3 ± 1.8</td><td>1.28</td><td>2.95 ± 1.24</td><td>1.07<sup>+0.11</sup><sub>-0.09</sub></td><td>0.35 ± 0.09</td><td>* [62]</td></tr><tr><td>MACS 0744</td><td>8.9 ± 0.8</td><td>2.22</td><td>4.18 ± 1.45</td><td>1.01<sup>+0.06</sup><sub>-0.07</sub></td><td>0.55<sup>+0.15</sup><sub>-0.11</sub></td><td></td></tr><tr><td>MACS 1149</td><td>8.7 ± 0.9</td><td>4.84</td><td>4.65 ± 1.65</td><td>1.07<sup>+0.09</sup><sub>-0.10</sub></td><td>0.38<sup>+0.09</sup><sub>-0.06</sub></td><td>* [62]</td></tr><tr><td>MACS 1206</td><td>10.8 ± 0.6</td><td>2.97</td><td>3.24 ± 1.03</td><td>0.98<sup>+0.07</sup><sub>-0.08</sub></td><td>0.40<sup>+0.10</sup><sub>-0.07</sub></td><td></td></tr><tr><td>MACS 1720</td><td>6.6 ± 0.4</td><td>1.58</td><td>2.83 ± 0.96</td><td>1.08 ± 0.06</td><td>0.46 ± 0.09</td><td></td></tr><tr><td>RXJ2248</td><td>12.4 ± 0.6</td><td>2.83</td><td>2.55 ± 1.20</td><td>0.81<sup>+0.12</sup><sub>-0.10</sub></td><td>0.39<sup>+0.13</sup><sub>-0.11</sub></td><td>* [63]</td></tr></tbody></table>

manner, we first consider a case similar to those of the previous sections, but where the $\gamma$ parameter, instead of being fixed to zero, is fixed to have the same central slope as the gas profile. The other parameters are fixed to be $\alpha = 1$ (as always) and $\beta = 6$. We find that the average ΔBIC= 0.5 compared to

the $\beta = 6$, $\gamma = 0$ case remains mostly reasonable, with only one case reaching a ΔBIC of 2. To push our investigations further, we then write down an effective model for the missing mass taking into account the hot gas mass profile of the cluster under consideration,
<figure><img src="image_7.png" /><figcaption>FIG. 3: MOND residual missing mass for Abell 209 (black solid line and gray-shaded area), together with the gas mass (red line) and the best fit profiles for Case 2, with $\beta = 4, 5, 6, 10$, as discussed in Section III A. The dashed vertical lines indicate the corresponding best fit scale radii $r_s$.</figcaption></figure>

$$ \rho_{\delta M}(r) = \eta \cdot \rho_{\text{gas}}(r) e^{-\lambda r/r_0}, \qquad (17) $$

where $\eta$ is a constant indicating the ratio of residual missing mass density over hot gas density. Obviously, the gas mass density does not fall down as sharply as the missing mass density in the outskirts, which can be modelled thanks to the exponential cut-off factor depending on the free parameter $\lambda$ which ensures a smooth and sharp decrease of $\rho_{\delta M}(r)$. The relevant physical parameter of the profile fitted here is actually the cut-off radius $r_0/\lambda$. For each individual cluster, the best-fit profile is shown as a blue dashed curve in Figure 1. Interestingly, this model actually provides an overall good BIC compared to the other cases presented above. The only two exceptions are the clusters MS 2137 and RXJ 1347, for which $\Delta$BIC$\approx$ 7 in favour of the case with fixed $\gamma = 0, \beta = 6$. The best-fit parameters for all clusters of the sample with this profile are listed in Table II, together with a flag for clusters displaying signs of disequilibria.

Figure 4 illustrates the remarkable uniformity for the missing-to-hot-gas density ratio and cut-off radius, which are almost constant, with weighted mean $\langle \log \eta \rangle = 0.93$ and geometric weighted mean $\langle r_0/\lambda \rangle = 0.43$ Mpc. The latter parameter is a bit more scattered than the missing-to-hot-gas density ratio, but the largest values are also associated with large error bars, compatible at $1\sigma$ with values closer to the average. On the other hand, one can note a slight tendency for the missing-to-hot-gas density ratio to increase with decreasing gas mass, although this tendency should be confirmed with a larger sample. As in the previous cases, clusters displaying signs of disequilibria do not display a different behaviour than the other ones, indicating that the lensing results are not sensitive to them.

## IV. DISCUSSION AND CONCLUSIONS

We have used parametric mass profiles derived from strong and weak-lensing shear, as well as magnification data, for a sample of 16 clusters from the CLASH survey, in order to compute the equivalent mass profile needed to produce the same gravitational field in MOND. This has allowed us to study possible scalings of the MOND residual missing mass in the inner parts with the observed baryonic distribution in galaxy clusters.

Fitting a very general profile to the residual missing mass profiles, we found that its inner slope is very close to a core whilst the outer slope is sharp with power-law slopes larger than 3.5, in accordance with the results of [25] from the hydrostatic equilibrium of hot gas in the X-COP sample. We also note a strong anti-correlation between the characteristic density of the residual missing mass and its scale radius, whose normalization and slope depends on the outer slope, whilst the scale radius itself tends to correlate with the total hot gas mass of clusters.

Pushing further the study of a possible connection between the distribution of baryons and that of the missing mass, we also consider the possibility of a “dark mass-follows-gas” profile in the central parts, followed by an exponential cut-off. We show that such a model provides a particularly good description of the missing profiles, with a remarkable uniformity of the missing-to-hot-gas density ratio and cut-off radius. We, however, caution that this uniformity is obtained from lensing alone, but that cross-checks from hydrostatic equilibrium of X-ray emitting gas, and from the kinematics of galaxies, should also be performed to reach a firm conclusion on the typical values. In the Appendix, we perform such an analysis based on the kinematics data of the cluster MACS 1206, which confirms that the scaling of the observed gas mass profile and residual missing mass profile works excellently, but we caution that it yields significantly smaller values of the missing-to-hot-gas density ratio in the central parts ($\log \eta = 0.33$) and a larger cut-off radius $r_0/\lambda = 2.12$ Mpc). Clearly, the lensing and kinematics results are in tension. This inconsistency between kinematics and lensing data is inherent to the data, and not to the gravity theory, and most plausibly due to systematic uncertainties in the interpretation of kinematics data (similar to the hydrostatic bias), despite the fact we used the most accurate available dataset on the market. The kinematics results, therefore, appear to be biased with respect to lensing. We also notice that, before identifying the cause of the discrepancy accurately, including more kinematics datasets in our analysis would then not improve the picture, but act as to strengthen the bias.

Since our study relies on a parametric form of the lensing potential, we do not attempt here to study whether the lensing signal is in tension with MOND in the outskirts of the clusters, which would require a fully non-parametric approach that should be the topic of further works. We emphasize though that the lensing reconstruction of the total mass profile of the cluster has been performed for different mass profile parametrizations. It is worth noting that [25] have shown that the external field effect expected in MOND (and
<figure><img src="image_8.png" /><figcaption>FIG. 4: Left and central plots: fitted parameters in the effective missing density model of equation (17) as a function of the gas mass of the cluster sample. Right: $\eta$ (in log) as a function of the cluster temperatures as found in [56]. For a detailed discussion of these results, see Section III. Note the remarkable constancy of both $\eta$ and $r_0/\lambda$ with both gas mass and temperature. As above, black stars identify the clusters with signs of lack of equilibrium.</figcaption></figure>

mostly irrelevant for our present study of the mass profile in the inner parts) could have a significant effect, whilst [33] have shown that relativistic versions of MOND might also induce a sharp drop of the gravitational field in the outskirts. Our results in the inner parts of clusters can serve as a crucible for such relativistic theories of MOND, or for any other tentative hypothesis regarding the nature of the clusters residual missing mass in the MOND context. For instance, the enhancement of the gravitational boost followed by an oscillation at large radii found by [33] could provide a satisfactory explanation only if the boost is confined to the central parts of the cluster and following the hot gas mass profile, which is not in qualitative agreement with the theoretical profiles presented in [33]. It is worth noting that the same CLASH clusters have previously been studied in the context of the Radial Acceleration Relation by [64], who concluded that a relation similar to that found in galaxies could be deduced from lensing data, but with a characteristic acceleration an order of magnitude larger than expected in MOND. This would appear to be in line with a variation of the $a_0$ scale as a function of other quantities such as the depth of the potential well [34]. This also might appear consistent at zeroth order with our results. However this is not precisely equivalent to our present suggestion. Indeed, an $a_0$ rescaling would imply that the residual missing mass would be entirely correlated (in the usual MONDian way) to the total enclosed baryonic mass in the central parts, including the BCG mass. Our results, on the contrary, exhibit an excellent correlation with the distribution of the gas mass alone. In the X-COP clusters of [25], it also appears clearly that the acceleration profile generated by observable baryons does not follow the shape of the acceleration generated by the residual missing mass. It would, however, be most worthwhile to check whether our findings could be corroborated in these X-COP clusters, namely if removing the influence of the BCG would result in a scaling similar to that found in our CLASH sample. This appears, by eye, to be the case but with typically smaller values of the missing-to-hot-gas density ratio (see also

the Appendix hereafter), which might be related to the hydrostatic bias.

In summary, our study uncovers for the first time a close connection between the distribution of X-ray emitting hot gas and that of the residual missing mass of MOND in galaxy clusters, which can serve as a guide for testing and excluding tentative solutions to the problem of MOND in galaxy clusters. We nevertheless caution that, while the suggested connection between the hot gas mass profile and the missing mass profile appears to hold also from a kinematics analysis (see Appendix), the exact values of the missing-to-hot-gas density ratio can be smaller than reported in our lensing analysis.

## ACKNOWLEDGMENTS

The authors wish to warmly thank the anonymous reference as well as Tobias Mistele for constructive feedback and extremely helpful checks. BF acknowledges funding from the European Research Council (ERC) under the European Union's Horizon 2020 research and innovation program (grant agreement No. 834148) and from the Agence Nationale de la Recherche (ANR projects ANR-18-CE31-0006 and ANR-19-CE31-0017). IDS acknowledges funding by the Czech Grant Agency (GAČR) under the grant number 21-16583M. LP acknowledges funding from the Czech Academy of Sciences under the grant number LQ100102101. The statistical analysis was designed and led by BF and IDS, the Python code was developed by LP under the supervision of the former two, the lensing data used in the analysis were provided by K. Umetsu, and the gas-profile data were provided by V. Salzano respectively. A reproduction package for the statistical analysis of this work in Python is available at: https://zenodo.org/records/15299349
[1] L. Lombriser and J. Peñarrubia, How chameleons core dwarfs with cusps, Phys. Rev. D 91, 084022 (2015), arXiv:1407.7862 [astro-ph.CO].

[2] M. Milgrom, A modification of the Newtonian dynamics as a possible alternative to the hidden mass hypothesis., ApJ 270, 365 (1983).

[3] B. Famaey and S. S. McGaugh, Modified Newtonian Dynamics (MOND): Observational Phenomenology and Relativistic Extensions, Living Reviews in Relativity 15, 10 (2012), arXiv:1112.3960 [astro-ph.CO].

[4] B. Famaey and A. Durakovic, Modified Newtonian Dynamics (MOND), arXiv e-prints , arXiv:2501.17006 (2025), arXiv:2501.17006 [astro-ph.GA].

[5] R. Sancisi, The visible matter – dark matter coupling, in Dark Matter in Galaxies, Vol. 220, edited by S. Ryder, D. Pisano, M. Walker, and K. Freeman (2004) p. 233, arXiv:astro-ph/0311348 [astro-ph].

[6] R. A. Swaters, R. Sancisi, T. S. van Albada, and J. M. van der Hulst, The rotation curves shapes of late-type dwarf galaxies, A&A 493, 871 (2009), arXiv:0901.4222 [astro-ph.CO].

[7] K. A. Oman, J. F. Navarro, A. Fattahi, C. S. Frenk, T. Sawala, S. D. M. White, R. Bower, R. A. Crain, M. Furlong, M. Schaller, J. Schaye, and T. Theuns, The unexpected diversity of dwarf galaxy rotation curves, MNRAS 452, 3650 (2015), arXiv:1504.01437 [astro-ph.GA].

[8] A. Ghari, B. Famaey, C. Laporte, and H. Haghi, Dark matter-baryon scaling relations from Einasto halo fits to SPARC galaxy rotation curves, A&A 623, A123 (2019), arXiv:1811.06554 [astro-ph.GA].

[9] S. S. McGaugh, J. M. Schombert, G. D. Bothun, and W. J. G. de Blok, The Baryonic Tully-Fisher Relation, ApJ 533, L99 (2000), arXiv:astro-ph/0003001 [astro-ph].

[10] F. Lelli, S. S. McGaugh, J. M. Schombert, H. Desmond, and H. Katz, The baryonic Tully-Fisher relation for different velocity definitions and implications for galaxy angular momentum, MNRAS 484, 3267 (2019), arXiv:1901.05966 [astro-ph.GA].

[11] H. Desmond, The scatter, residual correlations and curvature of the SPARC baryonic Tully-Fisher relation, MNRAS 472, L35 (2017), arXiv:1706.01017 [astro-ph.GA].

[12] T. Mistele, S. McGaugh, F. Lelli, J. Schombert, and P. Li, Indefinitely Flat Circular Velocities and the Baryonic Tully–Fisher Relation from Weak Lensing, ApJ 969, L3 (2024), arXiv:2406.09685 [astro-ph.GA].

[13] F. Lelli, S. S. McGaugh, J. M. Schombert, and M. S. Pawlowski, One Law to Rule Them All: The Radial Acceleration Relation of Galaxies, ApJ 836, 152 (2017), arXiv:1610.08981 [astro-ph.GA].

[14] R. Stiskalek and H. Desmond, On the fundamentality of the radial acceleration relation for late-type galaxy dynamics, MNRAS 525, 6130 (2023), arXiv:2305.19978 [astro-ph.GA].

[15] H. Desmond, A. Hees, and B. Famaey, On the tension between the radial acceleration relation and Solar system quadrupole in modified gravity MOND, MNRAS 530, 1781 (2024), arXiv:2401.04796 [astro-ph.GA].

[16] R. H. Sanders, A Faber-Jackson relation for clusters of galaxies: implications for modified dynamics., A&A 284, L31 (1994).

[17] L. S. The and S. D. M. White, Modified Newtonian Dynamics and the Coma Cluster, AJ 95, 1642 (1988).

[18] D. Gerbal, F. Durret, M. Lachieze-Rey, and G. Lima-Neto, Analysis of X-ray galaxy clusters in the framework of modified Newtonian dynamics., A&A 262, 395 (1992).

[19] R. H. Sanders, The Virial Discrepancy in Clusters of Galaxies in the Context of Modified Newtonian Dynamics, ApJ 512, L23 (1999), arXiv:astro-ph/9807023 [astro-ph].

[20] R. H. Sanders, Clusters of galaxies with modified Newtonian dynamics, MNRAS 342, 901 (2003), arXiv:astro-ph/0212293 [astro-ph].

[21] E. Pointecouteau and J. Silk, New constraints on modified Newtonian dynamics from galaxy clusters, MNRAS 364, 654 (2005), arXiv:astro-ph/0505017 [astro-ph].

[22] P. Natarajan and H. Zhao, MOND plus classical neutrinos are not enough for cluster lensing, MNRAS 389, 250 (2008), arXiv:0806.3080 [astro-ph].

[23] G. W. Angus, B. Famaey, and D. A. Buote, X-ray group and cluster mass profiles in MOND: unexplained mass on the group scale, MNRAS 387, 1470 (2008), arXiv:0709.0108 [astro-ph].

[24] S. Ettori, V. Ghirardini, D. Eckert, E. Pointecouteau, F. Gastaldello, M. Sereno, M. Gaspari, S. Ghizzardi, M. Roncarelli, and M. Rossetti, Hydrostatic mass profiles in X-COP galaxy clusters, A&A 621, A39 (2019), arXiv:1805.00035 [astro-ph.CO].

[25] R. Kelleher and F. Lelli, Galaxy clusters in Milgromian dynamics: Missing matter, hydrostatic bias, and the external field effect, arXiv e-prints , arXiv:2405.08557 (2024), arXiv:2405.08557 [astro-ph.CO].

[26] D. Clowe, M. Bradač, A. H. Gonzalez, M. Markevitch, S. W. Randall, C. Jones, and D. Zaritsky, A Direct Empirical Proof of the Existence of Dark Matter, ApJ 648, L109 (2006), arXiv:astro-ph/0608407 [astro-ph].

[27] G. W. Angus, H. Y. Shan, H. S. Zhao, and B. Famaey, On the Proof of Dark Matter, the Law of Gravity, and the Mass of Neutrinos, ApJ 654, L13 (2007), arXiv:astro-ph/0609125 [astro-ph].

[28] G. W. Angus, A. Diaferio, B. Famaey, and K. J. van der Heyden, Cosmological simulations in MOND: the cluster scale halo mass function with light sterile neutrinos, MNRAS 436, 202 (2013), arXiv:1309.6094 [astro-ph.CO].

[29] J. D. Bekenstein, Relativistic gravitation theory for the modified Newtonian dynamics paradigm, Phys. Rev. D 70, 083509 (2004), arXiv:astro-ph/0403694 [astro-ph].

[30] R. H. Sanders, A tensor-vector-scalar framework for modified dynamics and cosmic dark matter, MNRAS 363, 459 (2005), arXiv:astro-ph/0502222 [astro-ph].

[31] M. Milgrom, Broader view of bimetric MOND, Phys. Rev. D 106, 084010 (2022), arXiv:2208.10882 [gr-qc].

[32] C. Skordis and T. Złośnik, New Relativistic Theory for Modified Newtonian Dynamics, Phys. Rev. Lett. 127, 161302 (2021), arXiv:2007.00082 [astro-ph.CO].

[33] A. Durakovic and C. Skordis, Towards galaxy cluster models in Aether-Scalar-Tensor theory: isothermal spheres and curiosities, J. Cosmology Astropart. Phys. 2024, 040 (2024), arXiv:2312.00889 [astro-ph.CO].

[34] H. Zhao and B. Famaey, Unifying all mass discrepancies with one effective gravity law?, Phys. Rev. D 86, 067301 (2012), arXiv:1207.6232 [astro-ph.CO].

[35] M. Milgrom, Marriage à-la-MOND: Baryonic dark matter in galaxy clusters and the cooling flow puzzle, New A Rev. 51, 906 (2008), arXiv:0712.4203 [astro-ph].

[36] J. Bekenstein and M. Milgrom, Does the missing mass problem signal the breakdown of Newtonian gravity?, ApJ 286, 7 (1984).

[37] M. Milgrom, Quasi-linear formulation of MOND, MNRAS 403, 886 (2010), arXiv:0911.5464 [astro-ph.CO].
[38] M. Milgrom, Generalizations of quasilinear MOND, Phys. Rev. D **108**, 084005 (2023), arXiv:2305.01589 [astro-ph.GA].

[39] M. Milgrom, Tripotential MOND theories, Phys. Rev. D **108**, 063009 (2023), arXiv:2305.19986 [gr-qc].

[40] F. Lüghausen, B. Famaey, and P. Kroupa, Phantom of RAMSES (POR): A new Milgromian dynamicsN-body code, Canadian Journal of Physics **93**, 232 (2015), arXiv:1405.5963 [astro-ph.GA].

[41] F. Renaud, B. Famaey, and P. Kroupa, Star formation triggered by galaxy interactions in modified gravity, MNRAS **463**, 3637 (2016), arXiv:1609.04407 [astro-ph.GA].

[42] S. T. Nagesh, I. Banik, I. Thies, P. Kroupa, B. Famaey, N. Wittenburg, R. Parziale, and M. Haslbauer, The Phantom of RAMSES user guide for galaxy simulations using Milgromian and Newtonian gravity, Canadian Journal of Physics **99**, 607 (2021), arXiv:2101.11011 [astro-ph.IM].

[43] J. Bekenstein and M. Milgrom, Does the missing mass problem signal the breakdown of Newtonian gravity?, ApJ **286**, 7 (1984).

[44] E. Contini and Q. Gu, On the mass distribution of the intracluster light in galaxy groups and clusters, The Astrophysical Journal **901**, 128 (2020).

[45] K. E. Furnell, C. A. Collins, L. S. Kelvin, I. K. Baldry, P. A. James, M. Manolopoulou, R. G. Mann, P. A. Giles, A. Bermeo, M. Hilton, R. Wilkinson, A. K. Romer, C. Vergara, S. Bhargava, J. P. Stott, J. Mayers, and P. Viana, The growth of intracluster light in XCS-HSC galaxy clusters from 0.1 < z < 0.5, MNRAS **502**, 2419 (2021), arXiv:2101.01644 [astro-ph.GA].

[46] E. Laudato, V. Salzano, and K. Umetsu, Multi-component DHOST analysis in galaxy clusters, arXiv e-prints, arXiv:2110.11019 (2021), arXiv:2110.11019 [astro-ph.CO].

[47] H. Zhao, Analytical models for galactic nuclei, MNRAS **278**, 488 (1996), arXiv:astro-ph/9509122 [astro-ph].

[48] K. Umetsu, E. Medezinski, M. Nonino, J. Merten, M. Postman, M. Meneghetti, M. Donahue, N. Czakon, A. Molino, S. Seitz, D. Gruen, D. Lemze, I. Balestra, N. Benítez, A. Biviano, T. Broadhurst, H. Ford, C. Grillo, A. Koekemoer, P. Melchior, A. Mercurio, J. Moustakas, P. Rosati, and A. Zitrin, CLASH: Weak-lensing Shear-and-magnification Analysis of 20 Galaxy Clusters, ApJ **795**, 163 (2014), arXiv:1404.1375 [astro-ph.CO].

[49] A. Zitrin, A. Fabris, J. Merten, P. Melchior, M. Meneghetti, A. Koekemoer, D. Coe, M. Maturi, M. Bartelmann, M. Postman, K. Umetsu, G. Seidel, I. Sendra, T. Broadhurst, I. Balestra, A. Biviano, C. Grillo, A. Mercurio, M. Nonino, P. Rosati, L. Bradley, M. Carrasco, M. Donahue, H. Ford, B. L. Frye, and J. Moustakas, Hubble Space Telescope Combined Strong and Weak Lensing Analysis of the CLASH Sample: Mass and Magnification Models and Systematic Uncertainties, ApJ **801**, 44 (2015), arXiv:1411.1414 [astro-ph.CO].

[50] K. Umetsu, Cluster-galaxy weak lensing, A&A Rev. **28**, 7 (2020), arXiv:2007.00506 [astro-ph.CO].

[51] K. Umetsu, A. Zitrin, D. Gruen, J. Merten, M. Donahue, and M. Postman, CLASH: Joint Analysis of Strong-lensing, Weak-lensing Shear, and Magnification Data for 20 Galaxy Clusters, ApJ **821**, 116 (2016), arXiv:1507.04385.

[52] J. F. Navarro, C. S. Frenk, and S. D. M. White, A Universal Density Profile from Hierarchical Clustering, ApJ **490**, 493 (1997), arXiv:astro-ph/9611107.

[53] A. Burkert, The Structure of Dark Matter Halos in Dwarf Galaxies, ApJ **447**, L25 (1995), astro-ph/9504041.

[54] L. Hernquist, An analytical model for spherical galaxies and bulges, ApJ **356**, 359 (1990).

[55] G. A. Mamon, A. Cava, A. Biviano, A. Moretti, B. Poggianti, and D. Bettoni, Structural and dynamical modeling of WINGS clusters. II. The orbital anisotropies of elliptical, spiral, and lenticular galaxies, A&A **631**, A131 (2019), arXiv:1901.06393 [astro-ph.GA].

[56] M. Postman, D. Coe, N. Benítez, L. Bradley, T. Broadhurst, M. Donahue, H. Ford, O. Graur, et al., The Cluster Lensing and Supernova Survey with Hubble: An Overview, ApJS **199**, 25 (2012), arXiv:1106.3328 [astro-ph.CO].

[57] M. Girardi, W. Boschin, A. Mercurio, N. Nocerino, M. Nonino, P. Rosati, A. Biviano, R. Demarco, C. Grillo, B. Sartoris, P. Tozzi, and E. Vanzella, Clash-vlt: Galaxy cluster macs j0329-0211 and its surroundings using galaxies as kinematic tracers, Astronomy and Astrophysics **692**, A175 (2024).

[58] I.-N. T. Chiu and S. M. Molnar, Testing hydrostatic equilibrium in galaxy cluster ms 2137, The Astrophysical Journal **756**, 1 (2012).

[59] R. E. Johnson, J. Zuhone, C. Jones, W. R. Forman, and M. Markevitch, Sloshing gas in the core of the most luminous galaxy cluster rxj1347.5-1145, The Astrophysical Journal **751**, 95 (2012).

[60] W.-C. Feng, H. Yu, H.-H. Zhao, X.-L. Hou, S.-M. Jia, C.-K. Li, Y.-L. Cheng, P. Tozzi, M. Sun, and Y. Chen, Weak Merging Scenario of CLASH Cluster A209, Research in Astronomy and Astrophysics **24**, 045016 (2024).

[61] C. A. Raney, C. R. Keeton, and D. T. Zimmerman, A new framework for understanding systematic errors in cluster lens modelling - I. Selection and treatment of cluster member galaxies, MNRAS **508**, 5587 (2021).

[62] G. Giovannini, M. Cau, A. Bonafede, H. Ebeling, L. Feretti, M. Girardi, M. Gitti, F. Govoni, A. Ignesti, M. Murgia, G. B. Taylor, and V. Vacca, Diffuse radio sources in a statistically complete sample of high-redshift galaxy clusters, A&A **640**, A108 (2020), arXiv:2006.08494 [astro-ph.CO].

[63] A. Mercurio, P. Rosati, A. Biviano, M. Annunziatella, M. Girardi, B. Sartoris, M. Nonino, M. Brescia, G. Riccio, C. Grillo, I. Balestra, G. B. Caminha, G. De Lucia, R. Gobat, S. Seitz, P. Tozzi, M. Scodeggio, E. Vanzella, G. Angora, P. Bergamini, S. Borgani, R. Demarco, M. Meneghetti, V. Strazzullo, L. Tortorelli, K. Umetsu, A. Fritz, D. Gruen, D. Kelson, M. Lombardi, C. Maier, M. Postman, G. Rodighiero, and B. Ziegler, CLASH-VLT: Abell S1063. Cluster assembly history and spectroscopic catalogue, A&A **656**, A147 (2021), arXiv:2109.03305 [astro-ph.GA].

[64] Y. Tian, K. Umetsu, C.-M. Ko, M. Donahue, and I. N. Chiu, The Radial Acceleration Relation in CLASH Galaxy Clusters, ApJ **896**, 70 (2020), arXiv:2001.08340 [astro-ph.CO].

[65] P. Rosati, I. Balestra, C. Grillo, A. Mercurio, M. Nonino, A. Biviano, M. Girardi, E. Vanzella, and Clash-VLT Team, CLASH-VLT: A VIMOS Large Programme to Map the Dark Matter Mass Distribution in Galaxy Clusters and Probe Distant Lensed Galaxies, The Messenger **158**, 48 (2014).

[66] G. B. Caminha, C. Grillo, P. Rosati, M. Meneghetti, A. Mercurio, S. Ettori, I. Balestra, A. Biviano, K. Umetsu, E. Vanzella, M. Annunziatella, M. Bonamigo, C. Delgado-Correal, M. Girardi, M. Lombardi, M. Nonino, B. Sartoris, P. Tozzi, M. Bartelmann, L. Bradley, K. I. Caputi, D. Coe, H. Ford, A. Fritz, R. Gobat, M. Postman, S. Seitz, and A. Zitrin, Mass distribution in the core of MACS J1206. Robust modeling from an exceptionally large sample of central multiple images, A&A **607**, A93 (2017), arXiv:1707.00690 [astro-ph.GA].

[67] A. Biviano, L. Pizzuti, A. Mercurio, B. Sartoris, P. Rosati, S. Ettori, M. Girardi, C. Grillo, G. B. Caminha, and M. Nonino, CLASH-VLT: The Inner Slope of the MACS J1206.2-
0847 Dark Matter Density Profile, ApJ 958, 148 (2023),
arXiv:2307.06804 [astro-ph.CO].

[68] M. Annunziatella, A. Biviano, A. Mercurio, M. Nonino, P. Rosati, I. Balestra, V. Presotto, M. Girardi, R. Gobat, C. Grillo, E. Medezinski, D. Kelson, M. Postman, M. Brescia, B. Sartoris, R. Demarco, A. Koekemoer, D. Lemze, L. Bradley, and B. Ziegler, Clash-vlt: The stellar mass function and stellar mass density profile of the z=0.44 cluster of galaxies macs j1206.2-0847, Astronomy and Astrophysics (2014).

[69] B. et al., in prep., (2023).

[70] L. Pizzuti, I. D. Saltas, and L. Amendola, mg-mamposst: a code to test modifications of gravity with internal kinematics and lensing analyses of galaxy clusters, MN-RAS 506, 595 (2021), https://academic.oup.com/mnras/article-pdf/506/1/595/38884757/stab1727.pdf.

[71] L. Pizzuti, I. D. Saltas, A. Biviano, G. Mamon, and L. Amendola, Mg-mamposst, a fortran code to test gravity at galaxy-cluster scales, Journal of Open Source Software 8, 4800 (2023).

[72] O. Tiret, F. Combes, G. W. Angus, B. Famaey, and H. S. Zhao, Velocity dispersion around ellipticals in MOND, A&A 476, L1 (2007), arXiv:0710.4070 [astro-ph].

[73] J. Einasto, On the Construction of a Composite Model for the Galaxy and on the Determination of the System of Galactic Parameters, Trudy Astrofizicheskogo Instituta Alma-Ata 5, 87 (1965).

## Appendix A: Kinematics of the cluster MACS 1206

As a complementary investigation to the above lensing analysis, we have compared the results with those obtained from kinematics data in the cluster MACS 1206. As we also commented briefly in the main text, these results appear to be in tension with those derived from lensing. This is independent of the gravity theory, and likely due to inherent systematics despite the good accuracy of the MACS 1206 kinematics. Therefore, caution is needed in their interpretation and comparison with the lensing analysis. The MACS 1206 cluster has been extensively studied in the spectroscopic follow-up campaign of CLASH with the VIMOS spectrograph at Very Large Telescope (CLASH-VLT [65]), complemented with ob-

servation from the MUSE integrated field instrument for the cluster core [66]. Data about the stellar velocity dispersion of the BCG itself have also been collected (see [67] and references therein<fnref n="5" />). The galaxy masses are taken from [68]. The reconstruction of the total dynamical mass follows the method presented in [69], where an upgraded version of the MG-MAMPOSS code [70, 71] has been applied to fit the member galaxies velocity field, along with the velocity dispersion of the BCG. Note that in kinematics' mass determination of clusters, one has to account for an additional unknown, namely the velocity anisotropy profile $\beta = 1 - (\sigma_{\theta}^2 + \sigma_{\phi}^2)/(2\sigma_r^2)$, where $\sigma_{r,\theta,\phi}$ are the velocity dispersion components in the radial, tangential and azimuthal direction, respectively. The velocity anisotropy profile $\beta(r)$ is modeled using a “generalized Tiret” [gT, 72] and a “generalized Osipkov-Merritt” (gOM) models [55], which account for a broad range of orbits phenomenology in clusters. These are applied to four different mass profile ansatz, the generalized NFW (gNFW) model, the Einasto profile [73] with the characteristic exponent $m$ fixed to $m = 3$, as well as the Burkert and Hernquist models already mentioned in Section III A. As done for the lensing analyses, we then average over all the mass and anisotropy profiles in the kinematics analysis, compute the equivalent MOND mass profile $M(r)$, and subtract all the baryonic components (i.e. the BCG, gas and galaxy profiles) to obtain $\delta M(r)$, according to equation (13). We find a total residual missing mass of $(3.10 \pm 0.63) \times 10^{14} \text{M}_\odot$, only slightly lower than for the lensing analysis, but significantly lower in the most central regions. As already pointed out, this inconsistency is likely due to systematic uncertainties inherent to kinematics data, leading to biased results compared to lensing. Another reason for the discrepancy might however be linked to the fact that the kinematics analysis excludes from the phase-space a background group of galaxies located in the projection very near the cluster core, which is perhaps affecting in the lensing mass reconstruction. Importantly, the model of equation (17) is still a very good match to the data, but the main difference with the lensing analysis is that $\log \eta = 0.33_{-0.05}^{+0.04}$, much lower than in the lensing analysis. To compensate for a similar total amount of residual missing mass, the value of $r_0/\lambda = 2.12 \pm 0.45$ Mpc is on the other hand much larger than what was found by the lensing analysis.

<footnote n="5">⁵ The public redshift catalog, references and strong lensing models for this clusters are available at the CLASH-VLT website</footnote>