EDP Sciences logo
Subscriber Authentication Point
Search
Menu
Home All issues Volume 688 (August 2024) A&A, 688 (2024) A78 Full HTML
Open Access
Issue
A&A
Volume 688, August 2024
Article Number A78
Number of page(s) 20
Section Cosmology (including clusters of galaxies)
DOI https://doi.org/10.1051/0004-6361/202449968
Published online 06 August 2024

© The Authors 2024

Licence Creative CommonsOpen Access article, published by EDP Sciences, under the terms of the Creative Commons Attribution License (https://creativecommons.org/licenses/by/4.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

This article is published in open access under the Subscribe to Open model. Subscribe to A&A to support open access publication.

1. Introduction

The nature of dark matter (DM) is one of the most puzzling mysteries in physics. A major alternative to non-baryonic DM is Milgromian dynamics, also known as modified Newtonian dynamics (MOND), which posits a departure from the standard non-relativistic laws of gravity and/or inertia below a characteristic acceleration scale of a0 ≃ 1.2 × 10−10 m s−1 (Milgrom 1983a,b,c). The MOND paradigm has been very successful on galaxy scales, reproducing the dynamical properties of different types of galaxies without the need for DM (see Famaey & McGaugh 2012; Milgrom 2014; Banik & Zhao 2022 for reviews). However, the main motivation of MOND is not merely to eliminate the need for DM, but rather to find a natural explanation for the dynamical regularities and scaling laws that galaxies obey (e.g., Lelli et al. 2017; McGaugh et al. 2020), some of which were indeed predicted by MOND before being observed (see McGaugh 2020; Lelli 2022 for reviews).

Moving beyond galaxy scales, the MOND paradigm is able to explain the dynamical properties of galaxy groups without the need for DM (Milgrom 2002, 2018, 2019). In addition, basic MOND cosmological models predicted the existence of massive galaxies at z ≃ 10 (Sanders 1998, 2008) well in advance of their discovery by the James Webb Space Telescope; these still represent a challange for the ΛCDM cosmological model (see Boylan-Kolchin 2023 and references therein). However, the development of a fully fledged MOND cosmology is still in progress and is concurrent with the construction of a proper relativistic extension of MOND. A recent example is the aether scalar tensor (AeST) theory proposed by Skordis & Złośnik (2021), which is able to reproduce the power spectrum of the cosmic microwave background (CMB), the linear mass power spectrum, and the equal speed of gravitational waves and electromagnetic waves.

A long-standing challenge for MOND appears at the intermediate scales of galaxy clusters (The & White 1988; Gerbal et al. 1992; Sanders 1994, 1999, 2003; Aguirre et al. 2001; Pointecouteau & Silk 2005; Takahashi & Chiba 2007; Natarajan & Zhao 2008; Angus et al. 2008; Ettori et al. 2019; Tian et al. 2020, 2021; Li et al. 2023). As non-relativistic systems with internal accelerations of around a0, galaxy clusters should behave as predicted by the basic MOND tenets. Instead, it has been clear for decades that MOND requires more mass in galaxy clusters than is observed, although the deficit is drastically reduced with respect to that implied by Newtonian dynamics. This missing mass cannot be contained inside individual galaxies but must be embedded within the intracluster medium (ICM), which largely dominates the known baryonic mass budget of galaxy clusters.

Sanders (2003, 2007) proposed that the MOND ‘missing mass’ may consist of standard active neutrinos with masses of ∼2 eV, which were recently ruled out by the KATRIN experiment (Aker et al. 2022a). Other proposed solutions to the MOND galaxy cluster problem include: (1) undetected baryons, such as compact clouds of cold gas (Milgrom 2008), (2) sterile neutrinos (Angus et al. 2008, 2010), (3) extended MOND theories in which the MOND acceleration scale varies with the depth of the gravitational potential (Zhao & Famaey 2012; Hodson & Zhao 2017), and (4) the effect of additional fields in relativistic extensions of MOND (Durakovic & Skordis 2024). Interestingly, cold gas clouds and/or sterile neutrinos could explain the Bullet Cluster in MOND similarly to active neutrinos (Angus et al. 2007) because the missing mass would consist of non-collisional particles that become displaced from the hot gas during the clusters collision. Numerical calculations of colliding clusters in extended MOND theories have not yet been carried out.

In the present paper, we present a study of five galaxy clusters from the XMM-Newton cluster outskirts project (X-COP, Ghirardini et al. 2019). The aims of this study are to constrain the spatial distribution of the MOND missing mass and to gain new insights into its possible nature. The X-COP sample represents the state of the art in terms of cluster data, providing high-quality X-ray imaging and spectroscopy for studying ICM properties, as well as optical imaging and spectroscopy that can be used to study the galaxy properties. In Sect. 2 we describe the cluster dataset, while in Sect. 3 we present our fitting methodology. In Sect. 4, the main fitting results are described. In Sect. 5 we discuss current observational constraints on the nature of the MOND missing mass, while in Sect. 6 we provide a short summary of our findings.

2. Dataset

2.1. The cluster sample

We consider galaxy clusters from the X-COP project (Ghirardini et al. 2019), which measured the cluster gravitational field out to large radii combining X-ray data from XMM-Newton with Sunyaev–Zeldovich (SZ) data from Planck (Planck Collaboration XXIX 2014). The X-COP clusters were selected to have (1) a signal-to-noise ratio (S/N) of higher than 12 in the SZ effect from Planck, (2) an expected halo masses of M500 > 3 × 1014 M⊙, (3) a redshift in the range of 0.04 < z < 0.1, (4) an apparent angular size of θ500 > 10 arcmin, and (5) a hydrogen column density of NH < 1021 cm−2 along the line of sight. These criteria yielded 15 clusters; however, 3 clusters were excluded due to their complex morphologies, leaving a sample of 12 clusters (Ghirardini et al. 2019). Of these 12 clusters, five objects (A644, A1795, A2029, A2142, A2319) had additional measurements of the stellar mass distribution in the brightest cluster galaxy (BCG) and the surrounding satellite galaxies (Eckert et al. 2022). In this work, we focus on these five galaxy clusters with full baryonic information. Basic information on these galaxy clusters is given in Table 1.

Table 1.

Properties of the cluster sample.

Importantly, the clusters A644 and A2319 show signatures of merger activity, so the hot gas may potentially be out of dynamical equilibrium. A2319 is a well-studied cluster undergoing a major merger (Molendi et al. 1999; O’Hara et al. 2004; Govoni et al. 2004; Farnsworth et al. 2013; Yan et al. 2014; Storm et al. 2015). A644 displays various indications of past and/or ongoing merger activity (Bauer & Sarazin 2000; Buote et al. 2005; Fusco-Femiano et al. 2009). We keep these two clusters in our sample for the sake of comparison, but we stress that they may not be appropriate to test a dynamical theory such as MOND.

2.2. Observed gravitational field

Assuming that the hot gas is a spherically symmetric system in hydrostatic equilibrium, the pressure gradient exactly balances the gravitational attraction, so that

g obs ( r ) = − 1 ρ gas ( r ) d P gas ( r ) d r , $$ \begin{aligned} g_{\rm obs}(r) = -\frac{1}{\rho _{\rm gas}(r)} \frac{\mathrm{d}P_{\rm gas}(r)}{\mathrm{d}r}, \end{aligned} $$(1)

where gobs is the observed gravitational acceleration, Pgas is the gas pressure, and ρgas is the gas volume density. X-ray imaging provides the emissivity ϵν that is proportional to n e 2 ∝ ρ gas 2 $ n_{\rm e}^2 \propto \rho_{\rm gas}^2 $, where ne is the electron density of the ionised gas. To infer ρgas, one typically assumes no gas clumping and spherical geometry in order to deproject line-of-sight integrated quantities into 3D radial quantities. In addition, X-ray spectroscopy provides the gas temperature profile Tgas, again after appropriate emission-weighted deprojections (see, e.g., Eckert et al. 2022). Assuming the ideal gas law, we have

P gas ( r ) = k B κ m p ρ gas ( r ) T gas ( r ) , $$ \begin{aligned} P_{\rm gas}(r) = \frac{k_{\rm B}}{\kappa m_{\rm p}} \rho _{\rm gas}(r) T_{\rm gas}(r), \end{aligned} $$(2)

where kB is Boltzmann’s constant, mp is the proton mass, and κ is the mean atomic weight of the gas that depends on its chemical composition (κ = 0.6 for fully ionised gas with solar abundances). Equation (2) allows measuring the gas pressure profile and calculating its radial derivative to infer gobs(r) from Eq. (1). For numerical reasons, it is convenient to describe ρgas(r) and Tgas(r) via the sum of N analytic functions with a set of free parameters α, so that the derivative dPgas/dr can be calculated analytically once the values of α are specified (see Eckert et al. 2022).

At large cluster radii, the gas pressure can also be measured from the SZ effect through a unit-less Comptonization parameter y, which is given by

y = σ T m e c 2 ∫ l P gas ( l ) d l , $$ \begin{aligned} { y} = \frac{\sigma _{\rm T}}{m_{\rm e} c^2} \int _l P_{\rm gas}(l) \mathrm{d}l, \end{aligned} $$(3)

where σT is the Thomson cross section, me the mass of the electron, and c is the speed of light. Basically, the SZ effect provides the gas pressure Pgas integrated along the line of sight l, so the 3D pressure profile can be determined by deprojecting the y-profiles assuming spherical symmetry. In this work, we use the radial profiles of ρgas, Tgas, Pgas, and gobs computed by Eckert et al. (2022) using the so-called ‘non-parametric log-normal mixture reconstruction’. This is the most empirical, model-independent method to compute gobs(r).

In Newtonian dynamics, one equates gobs to the Newtonian gravitational field given by

g N ( r ) = G N M tot ( < r ) r 2 , $$ \begin{aligned} g_{\rm N}(r) = \dfrac{G_{\rm N} M_{\rm tot}( < r)}{r^2}, \end{aligned} $$(4)

where GN is Newton’s constant and Mtot is the ‘total’ dynamical mass, which is generally much larger than the visible baryonic mass Mbar, implying large amounts of DM. The total mass profile can then be measured combining Eqs. (1), (2), and (4).

In Milgromian dynamics, it is most natural to work in terms of accelerations, so that

g MOND μ ( g MOND a 0 ) = g N , bar , $$ \begin{aligned} g_{\rm MOND} \mu \left(\frac{g_{\rm MOND}}{a_0}\right) = g_{\rm N, bar}, \end{aligned} $$(5)

where gN, bar is the Newtonian baryonic gravitational field (discussed in the next section), gMOND is the expected MOND acceleration to be equated to gobs, and μ(x) is the MOND interpolating function with x = gMOND/a0. The functional form of μ(x) is not specified by the general MOND paradigm, but its asymptotic limits must be

μ ( x ≫ 1 ) ≈ 1 and μ ( x ≪ 1 ) ≈ x . $$ \begin{aligned} \mu (x \gg 1) \approx 1 \quad \mathrm{and}\quad \mu (x \ll 1) \approx x. \end{aligned} $$(6)

In particular, the low acceleration limit can be derived assuming scale invariance after the empirical normalisation of a0 (Milgrom 2009). Equation (5) is strictly valid in MOND modified gravity theories for isolated, spherically symmetric systems (Bekenstein & Milgrom 1984; Milgrom 2010, 2023). Equation (5) is also valid in MOND modified inertia theories in the case of isolated systems with purely circular orbits (Milgrom 1994), which is clearly not the case for the random gas motions in the ICM. In this work, therefore, we are primarily testing MOND modified gravity theories, albeit we expect Eq. (5) to provide the correct order of magnitude also in modified inertia theories (for isolated systems).

2.3. Newtonian baryonic gravitational field

The Newtonian gravitational field sourced by baryons (gN, bar) is given by the contributions from the ICM (gN, ICM) and from cluster galaxies (gN, gal). The contributions from globular clusters and intracluster light (made of free-floating stars inside the galaxy cluster) are negligible. At large cluster radii, the majority (∼90%) of the baryonic contribution is from the ICM, while a minor fraction (∼10%) is from gas and stars inside galaxies. At small radii, however, the BCG contribution can be important.

In this work, we use the values of gN, bar = gN, ICM + gN, gal provided by Eckert et al. (2022). The values of gN, ICM are computed using ρgas from X-ray imaging. In particular, the observed emissivity profile ϵ ν (r)∝ ρ gas 2 (r) $ \epsilon_\nu(r)\propto\rho_{\rm gas}^2(r) $ is described as the sum of several King functions, for which the correspondence between projected 2D profiles and intrinsic 3D ones can be written analytically, so the deprojection is trivial (see Eckert et al. 2020). The radial profile of gN, gal, instead, consists of two different components: the central BCG (gN, BCG) and satellite galaxies (gN, sat). Their derivation is briefly summarised in the following.

The mass distribution of the BCG is measured using r-band imaging and assuming a radially constant mass-to-light ratio Υ*. The value of Υ* was dynamically estimated by Loubser et al. (2020) by modelling the observed stellar kinematics. In particular, the observed stellar kinematics was fitted considering – in addition to the stellar contribution – a supermassive black hole and a DM halo, whose mass was determined using weak lensing data. Ideally, in a MOND context, one should redo the stellar kinematic fits without a DM halo using Milgromian dynamics, but the final values of Υ* are expected to be similar within the uncertainties because the DM contribution is small in the central parts of BCGs (in MONDian terminology, the galaxy is in the Newtonian regime at a ≫ a0). Moreover, the uncertainties on gN, BCG due to the DM halo parameters are included in the error budget (see Eckert et al. 2022 for details). For simplicity, therefore, we use the values computed by Loubser et al. (2020).

The mass distribution of satellite galaxies was computed using u-, g-, r- and i- band imaging; the values of Υ* were estimated from spectral energy distribution fitting with stellar population models (see van der Burg et al. 2015 for details). The gas content inside satellite galaxies is neglected, but it is expected to be very small (smaller than the uncertainties on Υ*) in passive cluster galaxies dominated by old stellar populations.

The Newtonian baryonic gravitational field is given by

g N , bar = g N , ICM + g N , BCG + g N , sat $$ \begin{aligned} g_{\rm N, bar} = g_{\rm N, ICM} + g_{\rm N, BCG} + g_{\rm N, sat} \end{aligned} $$(7)

Any additional ‘missing mass’ component could be trivially added to Eq. (7). Since we are assuming spherical symmetry, the Newtonian enclosed mass profile and the corresponding circular-velocity curve of a component ‘i’ are given by

M i ( < r ) = g N , i ( r ) · r 2 G N and V N , i ( r ) = g N , i ( r ) · r . $$ \begin{aligned} M_{\rm i}( < r) = \frac{g_{\rm N, i}(r) \cdot r^2}{G_{\rm N}} \quad \mathrm {and} \quad V_{\rm N,i}(r) = \sqrt{g_{\rm N, i}(r) \cdot r}. \end{aligned} $$(8)

3. Methodology

3.1. MOND theoretical formalism

The best approach to constrain the distribution of the missing matter in MOND is to assume a parametric density profile and directly fit it to the data, as is traditionally done for DM halos in Newtonian dynamics (e.g., Li et al. 2020). This method is more robust to uncertainties than the subtraction approach used by Eckert et al. (2022), which we revisit in Appendix A.

For practical reasons, it is convenient to rewrite the MOND law for isolated systems (Eq. (5)) as follows:

g MOND , iso = ν ( g N a 0 ) g N , $$ \begin{aligned} g_{\rm MOND, iso} = \nu \left(\frac{g_{\rm N}}{a_0}\right) g_{\rm N}, \end{aligned} $$(9)

where ν(y)⋅μ(x) = 1 with x = gMOND/a0 and y = gN/a0. Equation (9) is mathematically equivalent to the radial acceleration relation of galaxies (Lelli et al. 2017), but it is important to stress that the former is a MOND prediction that applies in specific situations, while the latter is an empirical description of the observed dynamics of galaxies. For example, Eq. (9) does not apply in cases where the so-called external field effect (EFE) is relevant.

The EFE is a characteristic prediction of MOND due to the violation of the strong equivalence principle (SEP). In general relativity, the SEP dictates that the internal dynamics of a gravitational system are independent of any external gravitational field in which the system is embedded (apart from possible tidal forces). MOND breaks the SEP due to its non-linearity (Bekenstein & Milgrom 1984), while preserving the weak equivalence principle (the universality of free fall). Thus, in MOND, the internal dynamics of a gravitational system can depend on its location in space and time due to the EFE.

In general, accounting for the EFE requires complex numerical 3D computations that solve the MOND modified Poisson equation (e.g., Chae & Milgrom 2022). To a first order approximation, the EFE can be analytically calculated considering a 1D solution in which the internal and external accelerations are summed in modulus as if their vectors always have the same direction. This gives the formula (Famaey & McGaugh 2012):

g MOND , EFE = g N ν ( g N + g N , e a 0 ) + g N , e [ ν ( g N + g N , e a 0 ) − ν ( g N , e a 0 ) ] $$ \begin{aligned} g_{\rm MOND, EFE} = g_{\rm N} \nu \left( \dfrac{g_{\rm N} + g_{\rm N,e}}{a_0} \right) + g_{\rm N, e} \left[ \nu \left( \dfrac{g_{\rm N} + g_{\rm N,e}}{a_0}\right) - \nu \left( \dfrac{g_{\rm N,e}}{a_0}\right)\right] \end{aligned} $$(10)

where gN, e is the Newtonian external field from the large-scale distribution of baryonic mass in the Universe.

In MOND modified gravity theories, Eq. (10) is an approximated formula that maximises the EFE because gN and gN, e are summed in modulus rather than in a vectorial way. Equations (9) and (10), therefore, provide two extreme MOND scenarios: no EFE and maximal EFE. Any other EFE implementation is expected to give intermediate results. A possible exception is represented by MOND modified inertia theories in which the EFE is effectively given by several times the instantaneous value of gN, e (Milgrom 2022), so Eq. (10) may possibly underestimate the EFE in those specific theories.

In Eqs. (9) and (10), we set

g N ( r ) = g N , mm ( r ; M mm , tot , r s ) + Υ bar g N , bar ( r ) , $$ \begin{aligned} g_{\rm N}(r) = g_{\rm N, mm}(r; M_{\rm mm, tot}, r_{\rm s}) + \Upsilon _{\rm bar} g_{\rm N, bar}(r), \end{aligned} $$(11)

where gN, mm is the Newtonian gravitational field of the missing mass component, which depends on two free parameters (Mmm, rs) as we describe in Sect. 3.3. The quantity Υbar, instead, is a nuisance parameter that scales up or down the baryonic component, similarly to a baryonic mass-to-light ratio. The purpose of Υbar is to account for systematic uncertainties in gbar (or equivalently in the total baryonic mass), which may be due to uncertainties in absolute X-ray flux calibration, deprojection methodology, gas clumpiness, assumed mean atomic weight, and so on. We impose a lognormal prior on log10(Υbar) in a Bayesian context, centered at 0 with a standard deviation of 0.1 dex, corresponding to average systematic uncertainties in the measured baryonic mass of ∼25%.

3.2. Bayesian fitting formalism

The fitting parameters are determined using a Markov-chain-Monte-Carlo (MCMC) method in Bayesian statistics. We define the likelihood ℒ = exp(−0.5χ2) with

χ 2 = ∑ k N [ g obs − g MOND ( p ) ] 2 δ g obs 2 , $$ \begin{aligned} \chi ^2 = \sum _{k}^{N} \frac{[g_{\rm obs} - g_{\rm MOND}(\boldsymbol{p})]^{2}}{\delta ^{2}_{g_{\rm obs}}}, \end{aligned} $$(12)

where gobs is the observed acceleration (Eq. (1)) at radius Rk, δgobs is the associated error, and gMOND is given by Eq. (9) or Eq. (10). The model acceleration depends on the fitting parameters p = {Mmm, rs, Υbar} plus gN, e in the EFE case. The posterior probability distributions of p are mapped using emcee (Foreman-Mackey et al. 2013). The MCMC chains are initialised with 200 walkers. We run 1000 burn-in iterations, then the sampler is run for another 2000 iterations. The emcee parameter a, which controls the size of the stretch move, is set equal to 2. This generally gives acceptance fractions larger than 50%. In all cases, the MCMC posterior probability distributions are well-behaved and show a single peak (see Appendix C), indicating that the fitting parameters are well determined. The maximum-likelihood values and their 68% confidence intervals (1σ errors) are given in Table 2.

Table 2.

Estimated parameters from MCMC fits (see Sect. 3.1 for details).

At large cluster radii (r ≳ 1 Mpc), an important systematic uncertainty is the so-called ‘hydrostatic bias’, which we discuss in detail in Sect. 5.4. In short, the outer gas may not be in hydrostatic equilibrium, so Eq. (1) may not provide a robust measurement of gobs. We aim to estimate the possible amount of hydrostatic bias in MOND, but its effect is degenerate with the strength of the EFE (see also Appendix A). To have upper and lower limits on hydrostatic bias in MOND, we proceed as follows. In the no-EFE case (Eq. (9)) we fit only the data at r < 1 Mpc, so that we obtain the maximum amount of hydrostatic bias needed in MOND. In the EFE case (Eq. (10)) we fit all data points, so that we estimate the minimum amount of hydrostatic bias needed in MOND. The fiducial radius of 1 Mpc is empirically inferred using the basic subtraction approach described in Appendix A.

3.3. Density profiles for the missing mass

We explored various density profiles for the missing mass and found good results with cored profiles of the following type:

ρ mm ( r ) = ρ 0 ( 1 + r r s ) n $$ \begin{aligned} \rho _{\rm mm}(r) = \frac{\rho _0}{(1+\frac{r}{r_{\rm s}})^n} \end{aligned} $$(13)

where ρ0 is the central core density, rs is a scale radius, and n is the outer slope of the density profile. We choose density profiles with a denominator of the type (1 + ℛ)n rather than (1 + ℛn) for the sake of simplicity (where ℛ = r/rs). In the former case, indeed, the enclosed mass can be expressed by simple analytic formulas for integer values of n (see, e.g., Eq. (14)), which is not true in the latter case for most values of n. For a fixed value of n, we do not expect major differences between these two types of cored profiles. Notably, a cored profile is expected in the case that the missing mass is made of massive neutrinos (Sanders 2007).

We also tested truncated spheres for which ρmm(r) = 0 at r > rt, where rt is a truncation radius. We explored different behaviours of ρmm(r) at r < rt, such as ρmm(r) = const, ρmm(r)∝r, and ρmm(r)∝1/r, but in all cases we found that a truncated sphere did not improve the fits with respect to Eq. (13).

In Newtonian dynamics, to obtain a flat velocity curve, the mass density profile must have an outer slope n ≃ 2 (e.g., an isothermal sphere). A value of n = 2 gives an enclosed mass profile that linearly increases with r and unphysically diverges to infinity at large radii. In ΛCDM cosmology, DM halos are predicted to have n ≃ 2 only at intermediate radii, while they tend to n ≃ 3 at large radii (e.g., the NFW profile), so the enclosed DM mass logarithmically diverges to infinity. For n > 3, instead, the enclosed mass converges to a physical, finite value.

We explored various values of n (from 2 to 6) and found that n = 4 provides the best fits (but we note that the Bayesian Information Criterion provides very small differences for fits with n ≥ 4). An outer slope n = 4 is particularly interesting because it corresponds to that expected for a MOND isothermal sphere (Milgrom 1984). For n = 4, the total enclosed mass is

M mm ( < r ) = ∫ 0 r 4 π r ′ 2 ρ 0 ( 1 + r ′ r s ) 4 d r ′ = 4 π ρ 0 r s 3 3 ( r r s ) 3 ( 1 + r r s ) 3 . $$ \begin{aligned} M_{\rm mm}( < r) = \int _0 ^{r} 4 \pi r^{\prime }^2 \frac{\rho _0}{(1+\frac{r^\prime }{r_{\rm s}})^4} \mathrm{d}r^\prime = \frac{4 \pi \rho _0 r_{\rm s}^3}{3} \frac{\big (\frac{r}{r_{\rm s}}\big )^3}{\big (1+\frac{r}{r_{\rm s}}\big )^3}. \end{aligned} $$(14)

In the limit of r → ∞, the mass converges to M mm,tot =4/3π ρ 0 r s 3 $ M_{\rm mm, tot} = 4/3 \pi \rho_0 r_{\rm s}^3 $, so it is convenient to rewrite Eq. (14) as

M mm ( R ) = M mm , tot R 3 ( 1 + R ) 3 . $$ \begin{aligned} M_{\rm mm}(\mathcal{R} ) = M_{\rm mm, tot}\frac{\mathcal{R} ^3}{(1+\mathcal{R} )^3}. \end{aligned} $$(15)

with the dimensionless variable ℛ = r/rs. The Newtonian gravitational field of the missing mass component (gmm) is then given by Eq. (4) and depends on two fitting parameters (Mmm, tot, rs).

4. Results

4.1. MOND fits with no EFE

Figure 1 shows the MCMC results for the isolated, no-EFE case (Eq. (9)). The fits are performed in acceleration space, but we also show the corresponding enclosed mass and circular velocity profiles from Eq. (8). We stress that Mobs and MMOND are fictitious dynamical concepts: the former represents the Newtonian mass inferred from the observed acceleration (Eq. (1)) while the latter is that inferred from the MOND-predicted acceleration. In fact, both Mobs and MMOND do not correspond to the sum of Mbar and Mmm (the observable masses of the physical components) because of the MOND boost. On the other hand, Vobs and VMOND represent the ‘measurable’ velocities that a test particle on a circular orbit should display to be in equilibrium with the cluster gravitational field.

thumbnail Fig. 1.

MCMC fit results for the isolated case (Eq. (9)). From left to right: acceleration profiles, enclosed mass profiles, circular-velocity profiles, and relative residual difference. In the first three panels, observed data (black points) are fitted with a mass model (black line) that includes the contribution of visible baryons (red line) and missing mass (blue line). We fit the data up to 1 Mpc (black points), where the hydrostatic bias should play a minor role, and then extrapolate the model to larger radii (grey line) to estimate the maximum hydrostatic bias in MOND (rightmost panels).

In general, the observations are well-fitted at r < 1 Mpc. The rightmost panels in Fig. 1 show the relative residuals in percentage, given by

δ = 100 · g M − g obs g obs . $$ \begin{aligned} \delta = 100 \cdot \frac{g_{\rm M} - g_{\rm obs}}{g_{\rm obs}}. \end{aligned} $$(16)

Considering mass residuals rather than acceleration residuals would give the same results (cf. with Eq. (8)). In general, the residuals at r < 1 Mpc are within ±20%, demonstrating that the fits are acceptable considering usual astronomical accuracy. In fact, acceleration residuals of about 10–20% may be due to mere deviations from spherical symmetry rather than hydrostatic bias. The residuals at r > 1 Mpc, where the MOND fit is extrapolated, quantify the maximum amount of hydrostatic bias in the case of no EFE. Such a maximum hydrostatic bias ranges from about 10% to 100%, apart for A2319 that is known to be involved in a major merger (see Sect. 2.1).

In some clusters (A2029 and A2142), the observed profiles show ‘bumps and wiggles’ that may not necessarily trace the equilibrium gravitational potential (i.e., a static high-density mass shell) but may rather be driven by local deviations from the hydrostatic equilibrium or some other idiosyncrasies in the observations. Our best-fit models generally pass through those bumps and wiggles, so the final parameters should not be heavily affected by them.

4.2. MOND fits with EFE

Figure 2 shows the MCMC results for the EFE case (Eq. (10)). In the EFE case we consider all data, even those at r > 1 Mpc, because we aim to quantify the minimum amount of hydrostatic bias needed in MOND when a maximal external field is taken into account. The observations are well fitted at all radii apart from the outermost ones. In the EFE case, the hydrostatic bias generally ranges from 10% to 40% up to R ≃ 2 Mpc. This level of hydrostatic bias is sensible and comparable to that inferred in a ΛCDM context, as we discuss in Sect. 5.4. For the very outermost points beyond 2 Mpc, the hydrostatic bias may go up to 50–60% but these points are more uncertain due to the lower S/N ratio of the SZ effect.

thumbnail Fig. 2.

Same as Fig. 1 but for the non-isolated case with a basic modelling of the EFE (Eq. (10)). In this case we fit all data points, even those at r > 1 Mpc, because we aim to estimate the minimum amonut of hydrostatic bias in MOND (rightmost panel).

For well-behaved clusters without known merging signatures (A1795, A2029, and A2142), the values of gN, e are around 0.001 − 0.002a0 and consistent with those inferred from the large-scale distribution of baryonic mass in the nearby Universe (see Figs. 7 and 8 in Chae et al. 2021), as well as those inferred from galaxy rotation-curve fits (Chae et al. 2021; Chae & Milgrom 2022). For clusters with known merger signatures (A644, A2319), the values of gN, e increase to 0.1 − 0.3a0. These high values may indicate either a real EFE signal due to the merging subcluster, or be a ‘compensating’ fitting artifact due to out-of-equilibrium dynamics. In any case, it is remarkable that known merging clusters clearly stand out in the sample in terms of gN, e, suggesting that this quantity is not merely an additional free parameter but could rather carry proper physical information.

The ratio of missing-to-baryonic mass (Mmm/Mbar) as a function of radius is shown in Fig. 3. The shape of these profiles explicitly show that the missing matter is more centrally concentrated than the ICM, as suggested by the early work of Sanders (1999, 2003). When the EFE is taken into account, the values of Mmm/Mbar are systematically higher than for the isolated case because the MOND boost is decreased, so more missing mass is required (see also Fig. A.1). For non-merging clusters (A1795, A2029, and A2142), Mmm/Mbar reaches a peak value of about 1 − 5 at r ≃ 200 − 300 kpc and decreases to about 0.4 − 1.5 at r ≳ 1 Mpc, indicating that the total amount of missing mass is comparable to or smaller than the ICM mass. For merging clusters (A644 and A2319), the ratio Mmm/Mbar display systematically higher values at all radii, which may again be driven by non-equilibrium dynamics rather than by actual missing mass. In any case, it is remarkable that known galaxy mergers (again) stand out in terms of MOND-inferred physical properties.

thumbnail Fig. 3.

Missing-to-baryonic mass ratio in galaxy clusters as a function of radius in the EFE case (solid line) and isolated case (dashed line). The high values of A644 and A2319 should be interpreted with caution because these two galaxy clusters display merger signatures.

4.3. Comparison with galaxy kinematics data

Another approach to measure the total mass profiles of galaxy clusters (so their acceleration profiles) is to model the observed kinematics of its galaxy members. Recently, Li et al. (2023) studied a sample of 16 galaxy clusters from HIFLUGCS (Tian et al. 2021), which includes the non-merging clusters in our sample (A1795, A2029, and A2142). In particular, Li et al. (2023) solved the spherical Jeans equations parametrising the galaxy velocity anisotropy and using two ‘virial shape parameters’, which ameliorate the well-known mass-anisotropy degeneracy.

Figure 4 compares the acceleration profiles from X-COP using hydrostatic equilibrium with those from Li et al. (2023) using galaxy kinematics. For A1795 and A2142, the two methods generally agree, but at large radii (r ≳ 1 Mpc) galaxy kinematics give systematically higher accelerations than the hydrostatic equilibrium, pointing to hydrostatic bias (see also Sect. 5.4). Remarkably, the MOND model reproduces the galaxy kinematics data better than the X-COP data at large radii, despite it was fitted to the latter dataset, not the former one. This fact testifies the predictive power of MOND even on galaxy cluster scales.

thumbnail Fig. 4.

Comparison between acceleration profiles from hydrostatic equilibrium (black dots; from X-COP) and Jeans modelling of galaxy kinematics (red dots; from Li et al. 2023). The black line is the MOND EFE model fitted to the X-COP data (same as in Fig. 2). Remarkably, for A1795 and A2142, the MOND model at large radii predicts the behaviour of the stellar kinematic data. For A2029 the comparison between the two datasets is poor, probably due to incompleteness in the galaxy kinematic data and/or unusual galaxy velocity anisotropy (see Li et al. 2023).

For A2029, the two methods significantly disagree. We cannot tell which one of the two profiles (if any) is the most reliable. However, we note that A2029 has been extensively discussed in Li et al. (2023) because it clearly stands out from the rest of their sample. Firstly, A2029 is the only galaxy cluster in Li et al. (2023) for which radially varying incompleteness may be a concern (see their Appendix A). Secondly, A2029 is the only galaxy clusters in Li et al. (2023) that show a strong radial variation in the galaxy velocity anisotropy, going from −0.4 at small radii to 0.7 at large radii (see their Fig. 4). This behaviour in anisotropy is quite unusual and cast some doubts on the enclosed mass profile, given the mass-anisotropy degeneracy. Future investigations of A2029 using both hydrostatic equilibrium and galaxy kinematics may share new light on this cluster.

5. Discussion

It has been known for decades that MOND can greatly reduce the need for DM in galaxy clusters but cannot entirely eradicate it. While the global mass discrepancies (Mtot/Mbar) are typically about 5–10 in Newtonian dynamics, they are reduced to a factor of about 2 in MOND (Sanders 1999). In this work, we studied the spatial distribution of the missing matter required in MOND. Physically acceptable solutions for the missing mass profile can be trivially found as long as there is some modest level of hydrostatic bias at large radii (r > 1 Mpc). The situation is summarised in Fig. 5, which shows the location of our galaxy clusters on the RAR defined by disc galaxies (Lelli et al. 2017). After a sensible missing mass component is included, galaxy clusters lie on the RAR within the observed scatter. At low acceleration (large radii), the cluster data systematically deviate below the average RAR, which may be due to hydrostatic bias and/or the EFE (e.g., Chae et al. 2020, 2021).

thumbnail Fig. 5.

Locations of galaxy clusters on the radial acceleration relation after accounting for a missing mass component. The top panel corresponds to the isolated case (Eq. (9)), while the bottom panel to the EFE case (Eq. (10)). Grey points show disc galaxies from the SPARC database (Lelli et al. 2016), while coloured points show the five galaxy clusters in our sample. The solid line is the MOND prediction for isolated systems; the dashed line is the line of unity (Newtonian prediction with no DM). In the top panel, open symbols indicate data at R > 1 Mpc, where the mass models with no EFE have been extrapolated.

If the missing mass required by MOND is a real physical entity (rather than a more complex modification of the gravitational law), then its nature needs to be identified. In the following, we discuss the viability and the existing constraints on (1) undetected ‘missing baryons’ (Sect. 5.1), (2) standard active neutrinos (Sect. 5.2), and (3) sterile neutrinos (Sect. 5.3). Finally, we discuss the long-standing issue of hydrostatic bias in both ΛCDM and MOND contexts (Sect. 5.4).

5.1. The nature of the missing mass: Undetected baryons

The ΛCDM cosmological model has a well-studied ‘missing baryons’ problem: the total amount of baryons observed in galaxies, galaxy groups, and galaxy clusters makes up only ∼18% of the cosmic baryon density (Ωbar) expected from ΛCDM fits to the CMB and from big bang nucleosynthesis (BBN) calculations (e.g., Fukugita & Peebles 2004). The standard solution is that the vast majority of baryons reside in the diffuse warm-hot intergalactic medium (WHIM) that exist in between galaxies and galaxy clusters (Shull et al. 2012; Macquart et al. 2020), but whose precise amount is still subject to significant uncertainties.

In a MOND cosmology, one may expect a similar missing baryon problem because the physics of BBN primarily depends on the baryon-to-photon ratio Ωbar/Ωγ and to a minor degree on the details of the expansion history, so on the underlying gravitational law1. In this case, the amount of missing baryons that one would need inside galaxy clusters is of the order of only 1–4% Ωbar, representing a small fraction of the cosmologically available missing baryons (82% Ωbar). In such a MOND cosmology, the bulk of the baryons should still reside in the WHIM.

Milgrom (2008) proposed that the missing baryons in galaxy clusters could be made of dense clouds of cold gas. A multi-phase, multi-temperature ICM is a sensible possibility to consider. For example, the intergalactic medium (IGM) around isolated galaxies is known to be multi-phase, showing evidence for cold gas clouds with T ≃ 104 K embedded in a hot diffuse plasma with T ≃ 106 − 107 K (e.g., Afruni et al. 2019, 2021, 2022). State-of-the-art H I surveys of the nearest galaxy clusters reach a point-source mass limit of 0.5 − 1.0 × 106 M⊙ at distances of 15 − 20 Mpc (Kleiner et al. 2023; Boselli et al. 2023). Then, if the gas clouds proposed by Milgrom (2008) are made of atomic hydrogen, they must be less massive than 5 × 105 M⊙. The upcoming square kilometre array (SKA) will improve the H I mass limits by at least one order of magnitude (e.g., Blyth et al. 2015). If the cold gas clouds are made of molecular hydrogen (H2), instead, they would be more difficult to detect because they may share the same metallicity of the ICM (∼0.3 Z⊙), so their CO emission would be weak due to variations in the CO-to-H2 conversion factor (e.g., Bolatto et al. 2013). In addition, there should be a mechanism that prevents the molecular clouds to collapse and form stars (otherwise they would emit UV radiation), so H I clouds may be more appealing than H2 ones.

Milgrom (2008) noted that the kinetic energy of these hypothetical clouds is about ten times larger than the thermal energy of the ICM within 200–300 kpc, providing a substantial energy reservoir. Then, if the gas clouds interact at a sufficient rate (via cloud-cloud collisions or cloud-ICM dynamical friction), their kinetic energy can be converted into thermal energy, heating up the ICM and solving the long-standing ‘cooling flow’ problem in galaxy clusters (Fabian 1994). Interestingly, this type of heating source would be smoothly distributed in the cluster core and steady with time, in contrast to feedback from active galactic nuclei that is generally proposed to solve the cooling flow problem in a ΛCDM context (Sijacki et al. 2007). Assuming that the kinetic heating rate from collisions balances the cooling rate of the ICM, Milgrom (2008) finds that the area covering factor of the clouds fA is about 10−3 − 10−4 depending on the ratio Mmm/Mbar in the cluster core. The volume filling factor fV of the clouds must be even smaller than this because fV ≃ fArcl/(2rs), where rc is the typical cloud radius and rs is the scale radius of the missing mass component (Table 2), so rcl/(2rs)≪1. Intriguingly, cold gas clouds with similarly low volume filling factors have been inferred to exist inside the giant Lyα-emitting nebulae of ionised gas surrounding massive quasar-host galaxies at z ≃ 3 (Pezzulli & Cantalupo 2019), which are presumably the progenitors of local BCGs at the center of galaxy clusters.

We propose an additional argument to infer the cloud properties: assuming that they are pressure confined by the hot ICM in the same way as the H I high-velocity clouds around the Milky Way are thought to be pressure confined by its hot corona (e.g., Spitzer 1956). Pressure equilibrium therefore implies that

n cold n hot = T hot T cold , $$ \begin{aligned} \dfrac{n_{\rm cold}}{n_{\rm hot}} = \dfrac{T_{\rm hot}}{T_{\rm cold}}, \end{aligned} $$(17)

where Tcold and ncold are the unknown temperatures and hydrogen number density of the cold gas clouds, while Thot ≃ 108 K and nhot ≃ 10−3 cm−3 are typical values for X-ray-emitting gas. If T ≃ 103 − 104 K as for typical H I gas, then ncold ≃ 10 − 100 cm−3. If T ≃ 10 − 100 K as for typical H2 gas, then ncold ≃ 103 − 104 cm−3. If we consider T ≃ 104 K and a conservative H I upper limit of Mc < 105 M⊙, we expect a cloud radius Rc < 50 pc. Clouds with lower masses and/or lower temperatures (such as molecular clouds) would have even smaller radii. These small sizes are consistent with the low values of fA inferred by Milgrom (2008) and make these clouds very difficult to detect in both emission (because they would not properly ‘fill the beam’ of existing radio/mm interferometers at cluster distances) and absorption (because they would have a very small chance to align with a bright background source such as a quasar).

To explore the possible connection between missing mass and hot gas, we plotted Mmm, tot and rs versus several observed ICM properties. We found a tentative correlation with the median value of the 3D temperature T3D (see Fig. 6). Formally, these relations are statistically significant, having high values of the Pearson’s, Spearman’s and Kendall’s correlation parameters (given in Table 3). However, we stress that these values are based on only five galaxy clusters, two of which may be out of dynamical equilibrium, possibly providing unrealiable values for Mmm and rs. A similar Mmm − T3D relation was reported by Angus et al. (2008) for a larger cluster sample, but using a different methodology and definitions of Mmm and T3D. If these correlations are confirmed by future studies, they could point to a real connection between the properties of the missing matter and those of the ICM, possibly favouring the missing baryons hypothesis.

thumbnail Fig. 6.

Relations between the 3D temperature of the ICM and the properties of the missing mass component: Total mass Mmm (leftmost panels) and scale radius rs (rightmost panels) for both the isolated and EFE case. In all cases there are statistically significant correlations but we warn that they are based only on five clusters, two of which (A644 and A2319) may be out of dynamical equilibrium.

Table 3.

Results of correlation tests between the 3D temperature of the ICM and the properties of the missing mass component (see Fig. 6).

5.2. The nature of the missing mass: Active neutrinos

Neutrino oscillation experiments provide clear evidence for neutrino masses and flavour mixing. These experiments do not measure the actual neutrino masses mν, but indicate that the largest mass difference between active neutrinos is Δmν ≃ 0.05 eV (e.g., Gonzalez-Garcia & Nir 2003). Since the number density of neutrinos produced in the early Universe is expected to be similar to that of photons, there must be a cosmic neutrino fluid with

Ω ν h 2 = 1 94 ∑ i N ν m ν , i eV $$ \begin{aligned} \Omega _{\rm \nu }h^2 = \dfrac{1}{94} \sum _i^{N_\nu } \dfrac{m_{\nu , \mathrm i}}\mathrm{eV} \end{aligned} $$(18)

where the sum is over Nν neutrino types (Sanders 2003). If mν ≃ Δmν, then Ων ≃ 10−3 and neutrinos have no significant cosmological mass density, so they cannot contribute to the mass of any bound system, as one expects in the ΛCDM context. The Planck data combined with baryonic acustic oscillations, indeed, imply that ∑imν, i < 0.12 eV in a flat ΛCDM cosmology (Planck Collaboration VI 2020) otherwise Ων would be too large and leave too little room to ΩCDM, reducing the power on small scales. In a MOND context, these constraints on ∑imν, i do not apply, so one can consider the case where mν ≫ Δmν and the masses of the three active neutrinos are nearly equal. Then, cosmic neutrinos would behave like hot dark matter (HDM), undergoing gravitational instability and collapse on spatial scales that depend on their mass (e.g., Sanders 2003).

Sanders (2003, 2007) found that standard active neutrinos with mν ≃ 2 eV could provide the missing matter required by MOND on galaxy cluster scales (R ≃ 1 Mpc), while they would be too light (too hot) to provide significant gravitational contributions on galaxy scales (R ≃ 100 kpc), preserving the MOND successes in that regime. On the contrary, Takahashi & Chiba (2007), Angus et al. (2008) and Natarajan & Zhao (2008) found that active neutrinos with mν ≃ 2 eV cannot entirely explain the MOND missing mass in the center of galaxy clusters, especially in low-mass and low-temperature ones. In any case, the latest results of the KATRIN experiment provide an upper limit of mν < 0.8 eV at 90% confidence level for the electron anti-neutrino mass (Aker et al. 2022a), ruling out active neutrinos with mν ≃ 2 eV. Hereafter, we briefly recap the basic neutrino argument and revisit the issue in light of our new measurements.

Cosmological neutrinos are created with a maximum phase-space density of (2πℏ)−3 per type (including antineutrinos) which is conserved during their subsequent evolution (Tremaine & Gunn 1979). During the formation of a collapsed object out of the mixed neutrino-baryon fluid, the two fluids are expected to attain the same velocity dispersion (temperature) via violent relaxation, so one can derive a relation between the maximum final density of the neutrino component and the velocity dispersion (temperature) of the system (cf. Sanders 2003, 2007):

ρ ν g cm − 3 ≤ 10 − 28 ( T keV ) 1.5 ∑ i N ν ( m ν , i eV ) 4 , $$ \begin{aligned} \dfrac{\rho _{\nu }}\mathrm{g\,cm^{-3}} \le 10^{-28} \left( \dfrac{T}{\mathrm{keV}} \right)^{1.5} \sum _i^{N_\nu } \left(\dfrac{m_{\nu , \mathrm i}}\mathrm{eV}\right)^4, \end{aligned} $$(19)

or equivalently

ρ ν M ⊙ pc − 3 ≤ ( 1.5 × 10 − 6 ) ( T keV ) 1.5 ∑ i N ν ( m ν , i eV ) 4 . $$ \begin{aligned} \dfrac{\rho _{\nu }}{M_\odot \,\mathrm{pc}^{-3}} \le (1.5\times 10^{-6}) \left( \dfrac{T}{\mathrm{keV}} \right)^{1.5} \sum _i^{N_\nu } \left(\dfrac{m_{\nu , \mathrm i}}\mathrm{eV}\right)^4. \end{aligned} $$(20)

If we consider three active neutrinos with the same mass mν, a, the mean T3D of the X-COP galaxy clusters (Table 1), and our measurements of ρ0 (Table 2), we find mν, a > 2 − 5 eV that is surely ruled out by the KATRIN upper limit. Thus, it is clear that active neutrinos cannot form the missing mass needed by MOND in galaxy clusters, confirming the findings of Takahashi & Chiba (2007), Angus et al. (2008) and Natarajan & Zhao (2008).

5.3. The nature of the missing mass: Sterile neutrinos

Another proposal is that of sterile neutrinos (Angus et al. 2008): hypothetical right-handed Fermions that are neutral under both weak interactions (‘sterile’) and electromagnetic ones (‘neutrinos’). Sterile neutrinos are empirically motivated by observed anomalies in neutrino oscillations and can be trivially added to the standard model of particle physics (e.g., Böser et al. 2020; Dasgupta & Kopp 2021). Assuming that MOND exactly converges to General Relativity (GR) in the early Universe, Angus (2009) showed that a single type of sterile neutrino can fully replace CDM in fitting the CMB data from WMAP, providing Ωνsh2 = 0.117 and mν,s ≃ 11 eV (in the case of thermal production for which Eq. (18) holds). The CMB fit of Angus (2009) fostered the investigation of a MOND cosmology supplemented by 11-eV sterile neutrinos, dubbed νHDM, which lead to some promising results but also severe challenges (Angus & Diaferio 2011; Katz et al. 2013; Angus et al. 2014; Haslbauer et al. 2020; Asencio et al. 2021; Wittenburg et al. 2023).

Importantly, the most recent CMB data from Planck strongly constrain the equation of state of a generic ‘dark fluid’ during recombination, showing no significant deviation from a dust fluid (CDM) and leaving little room for HDM (Thomas et al. 2016; Kopp et al. 2018; Ilić et al. 2021). In addition, relativistic extensions of MOND may significantly deviate from GR in the early Universe (Skordis & Złośnik 2021), so the constraints on neutrino masses (active and sterile) are expected to be different and need to be recomputed accounting for CMB lensing.

Our measurements of ρ0 and the Gunn-Tremaine limit (Eq. (19)) imply that mν, s > 10 ± 3 eV, considering the average mass from our five galaxy clusters and its standard deviation. This lower limit is consistent with mν, s ≃ 11 eV. Similar results were found by Angus et al. (2010) using self-consistent models of semi-degenerate neutrinos in equilibrium with the clusters’ gravitational potential. To date, the KATRIN experiment indicates that sterile neutrinos with mν,s ≃ 11 eV remain a viable possibility only for small values of active-to-sterile neutrino mixing (Aker et al. 2022b), which actually is the regime where the basic argument of Angus (2009) holds. Stronger constraints on sterile neutrinos will be available in future KATRIN runs, so we will know whether a mass of 11 eV is viable in the coming years.

5.4. Hydrostatic bias: Insights from both MOND and ΛCDM

The issue of hydrostatic bias has been amply studied in the standard cosmological context (e.g., Smith et al. 2016; Eckert et al. 2016; Henson et al. 2017; Angelinelli et al. 2020; Barnes et al. 2021). Clusters’ mass profiles are obtained under the assumption that the hot gas is in hydrostatic equilibrium with the gravitational potential, but this assumption may break down for a variety of reasons. For example, there may be non-thermal gravitational support that is unaccounted for by Eqs. (1) and (2). Such non-thermal support may be due to gas turbulence, bulk gas motions, or even gas rotation (Bartalesi et al. 2024). In addition, the ICM is not expected to be completely smooth and uniform, but to contain gas clumps that can alter the inferred gas density ρgas (e.g., Roncarelli et al. 2013; Eckert et al. 2015; Towler et al. 2023). Indeed, the observed X-ray emissivity ϵ is proportional to ρ gas 2 $ \rho_{\rm gas}^2 $, so small variations in ρgas can lead to large variations in ϵ. Finally, some galaxy clusters may be undergoing mergers in their outskirts, so they are not entirely relaxed structures in dynamical equilibrium, as it is probably the case for A644 and A2319 in our sample (see Sect. 2.1). These possible effects, among others, may cause the assumption of hydrostatic equilibrium to be invalid. Importantly, these effects are not unique to ΛCDM, nor to MOND, so they can play a role in both contexts.

In a ΛCDM context, hydrostatic bias is often invoked to solve two related cosmological problems: (1) discrepancies in the σ8 parameter (the normalisation of the power spectrum of density fluctuations) obtained from cluster counts with respect to that obtained from CMB fits (Planck Collaboration XX 2014), and (2) discrepancies between the baryonic fraction measured in galaxy clusters (the ratio between observed baryonic mass and halo virial mass) with the cosmic baryonic fraction Ωbar/Ωm from CMB fits (e.g., Eckert et al. 2019; Li et al. 2023; Wicker et al. 2023). The amount of hydrostatic bias has been measured using a variety of methods, leading to contradictory results, but typically range between 10 and 50%. In a MOND context, the level of hydrostatic bias may be different from that required in ΛCDM, but we actually find it to be comparable (10–40%) with a clear radial dependence (Fig. 2).

In the specific case of the X-COP cluster sample, Eckert et al. (2022) investigated the issue of hydrostatic bias by comparing masses from hydrostatic equilibrium with those from weak lensing. They concluded that hydrostatic bias is important at radii beyond R500, which is defined as the radius where the enclosed ‘dynamical mass’ density is 500 times the critical density of the Universe. Such a radius approximately corresponds to 1 Mpc (see Table 2 in Eckert et al. 2022), similar to the radius at which our MOND fits imply an increasing level of hydrostatic bias. For A1795 and A2142, the effect of hydrostatic bias is confirmed by comparing the acceleration profiles from X-COP with those from galaxy kinematics (Fig. 4), as also pointed out by Li et al. (2023). For A2029, instead, the situation remains unclear. In any case, a sensible level of hydrostatic bias allows for a MOND missing mass component with a physical mass density profile that converges to a finite total mass, so the concerns raised in Eckert et al. (2022) might simply be due to hydrostatic bias.

6. Conclusions

We built mass models of galaxy clusters in Milgromian dynamics (MOND). We focused on five clusters from the X-COP sample (Ghirardini et al. 2018) for which both high-quality baryonic (stars in galaxies and hot gas in the ICM) and dynamical information (from hydrostatic equilibrium and the SZ effect) are available. This sample contains two merging clusters (A644 and A2319), which are studied only for the sake of comparison with clusters without known merger signatures (A1795, A2029, A2142). Our results can be summarised as follows:

  1. We confirm the well-known result that galaxy clusters require additional ‘missing matter’ in MOND. Using a basic subtraction approach, the mass profiles of the missing matter decline at a rate of r > 1 Mpc and are therefore unphysical. This effect largely vanishes considering the MOND EFE and/or sensible levels of hydrostatic bias at large radii.

  2. Using a Bayesian MCMC approach, we fit the acceleration profiles of the clusters, adding a missing mass component. We find good results using a density profile with an inner core and an outer slope of −4, which gives a finite total mass for the missing matter component.

  3. MOND fits without the EFE imply a maximum amount of hydrostatic bias at R > 1 Mpc of between 10%–100%, whereas MOND fits considering the EFE reduce the implied amount of hydrostatic bias at R > 1 Mpc to 10%–40%. The required external field strengths (gNe/a0 ≃ 10−3) are consistent with those expected from the large-scale baryonic mass distribution (Chae et al. 2021) except for the two merging clusters, which may be out of dynamical equilibrium.

  4. For non-merging clusters, the missing-to-visible mass ratio (Mmm/Mbar) is about 1 − 5 at R ≃ 200 − 300 kpc and decreases to 0.4 − 1.1 at large radii, indicating that the total amount of missing mass is smaller than or comparable to the ICM mass. For merging clusters, the values of Mmm/Mbar are systematically higher but may be driven by out-of-equilibrium dynamics rather than actual missing mass.

In conclusion, galaxy clusters do not seem to be an insurmountable challenge for MOND as long as there is a sensible extra component with similar mass to the hot gas. Such missing matter may be baryonic, such as pressure-confined dense clouds of cold gas, or may require a minimal extension of the standard model of particle physics, such as sterile neutrinos with mν, s ≳ 10 eV. In comparison, in a ΛCDM cosmology, non-baryonic DM needs to be at least five times more abundant than the visible baryonic matter and cannot consist of light particles. We find tentative evidence for a possible correlation between the properties of the missing matter (mass and scale radius) and the temperature of the hot gas, suggesting that the missing matter may be a real physical entity related to the ICM properties. This tentative correlation could favour a purely baryonic interpretation of the missing matter, but we stress that it is based on only five galaxy clusters, two of which are possibly out of dynamical equilibrium. Clearly, a larger cluster sample with full baryonic information is needed to explore possible statistical correlations in a more comprehensive way and to shed new light on the nature of the missing mass in MOND.

1

More generally, any cosmological model with standard BBN and the usual interpretation of the CMB radiation (providing Ωγ) is expected to display a missing baryons problem.

Acknowledgments

We are grateful to Dominique Eckert and Stefano Ettori for providing the X-COP data, and to Pengfei Li for providing the galaxy kinematics data. We also thank Andrea Biviano, Harry Desmond, Stacy McGaugh, Moti Milgrom, James Schombert, Constantinos Skordis, and Paolo Tozzi for precious comments and suggestions about this work.

References

  1. Afruni, A., Fraternali, F., & Pezzulli, G. 2019, A&A, 625, A11 [Google Scholar]
  2. Afruni, A., Fraternali, F., & Pezzulli, G. 2021, MNRAS, 501, 5575 [NASA ADS] [Google Scholar]
  3. Afruni, A., Pezzulli, G., & Fraternali, F. 2022, MNRAS, 509, 4849 [Google Scholar]
  4. Aguirre, A., Schaye, J., & Quataert, E. 2001, ApJ, 561, 550 [NASA ADS] [CrossRef] [Google Scholar]
  5. Aker, M., Balzer, M., Batzler, D., et al. 2022a, J. Phys. G Nucl. Phys., 49, 100501 [NASA ADS] [CrossRef] [Google Scholar]
  6. Aker, M., Batzler, D., Beglarian, A., et al. 2022b, Phys. Rev. D, 105, 072004 [CrossRef] [Google Scholar]
  7. Angelinelli, M., Vazza, F., Giocoli, C., et al. 2020, MNRAS, 495, 864 [NASA ADS] [CrossRef] [Google Scholar]
  8. Angus, G. W. 2009, MNRAS, 394, 527 [NASA ADS] [CrossRef] [Google Scholar]
  9. Angus, G. W., & Diaferio, A. 2011, MNRAS, 417, 941 [NASA ADS] [CrossRef] [Google Scholar]
  10. Angus, G. W., Shan, H. Y., Zhao, H. S., & Famaey, B. 2007, ApJ, 654, L13 [CrossRef] [Google Scholar]
  11. Angus, G. W., Famaey, B., & Buote, D. A. 2008, MNRAS, 387, 1470 [CrossRef] [Google Scholar]
  12. Angus, G. W., Famaey, B., & Diaferio, A. 2010, MNRAS, 402, 395 [NASA ADS] [CrossRef] [Google Scholar]
  13. Angus, G. W., Diaferio, A., Famaey, B., Gentile, G., & Heyden, K. J. V. D. 2014, J. Cosmol. Astropart. Phys., 2014, 079 [CrossRef] [Google Scholar]
  14. Asencio, E., Banik, I., & Kroupa, P. 2021, MNRAS, 500, 5249 [Google Scholar]
  15. Banik, I., & Zhao, H. 2022, Symmetry, 14, 1331 [NASA ADS] [CrossRef] [Google Scholar]
  16. Barnes, D. J., Vogelsberger, M., Pearce, F. A., et al. 2021, MNRAS, 506, 2533 [NASA ADS] [CrossRef] [Google Scholar]
  17. Bartalesi, T., Ettori, S., & Nipoti, C. 2024, A&A, 682, A31 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  18. Bauer, F., & Sarazin, C. L. 2000, ApJ, 530, 222 [NASA ADS] [CrossRef] [Google Scholar]
  19. Bekenstein, J., & Milgrom, M. 1984, ApJ, 286, 7 [NASA ADS] [CrossRef] [Google Scholar]
  20. Blyth, S., van der Hulst, T. M., Verheijen, M. A. W., et al. 2015, in Advancing Astrophysics with the Square Kilometre Array (AASKA14), 128 [CrossRef] [Google Scholar]
  21. Bolatto, A. D., Wolfire, M., & Leroy, A. K. 2013, ARA&A, 51, 207 [CrossRef] [Google Scholar]
  22. Boselli, A., Serra, P., de Gasperin, F., et al. 2023, A&A, 676, A92 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  23. Böser, S., Buck, C., Giunti, C., et al. 2020, Progr. Part. Nucl. Phys., 111, 103736 [CrossRef] [Google Scholar]
  24. Boylan-Kolchin, M. 2023, Nat. Astron., 7, 731 [NASA ADS] [CrossRef] [Google Scholar]
  25. Buote, D. A., Humphrey, P. J., & Stocke, J. T. 2005, ApJ, 630, 750 [NASA ADS] [CrossRef] [Google Scholar]
  26. Chae, K.-H., & Milgrom, M. 2022, ApJ, 928, 24 [NASA ADS] [CrossRef] [Google Scholar]
  27. Chae, K.-H., Lelli, F., Desmond, H., et al. 2020, ApJ, 904, 51 [NASA ADS] [CrossRef] [Google Scholar]
  28. Chae, K.-H., Desmond, H., Lelli, F., McGaugh, S. S., & Schombert, J. M. 2021, ApJ, 921, 104 [NASA ADS] [CrossRef] [Google Scholar]
  29. Dasgupta, B., & Kopp, J. 2021, Phys. Rep., 928, 1 [NASA ADS] [CrossRef] [Google Scholar]
  30. Durakovic, A., & Skordis, C. 2024, JCAP, 2024, 040 [CrossRef] [Google Scholar]
  31. Eckert, D., Roncarelli, M., Ettori, S., et al. 2015, MNRAS, 447, 2198 [Google Scholar]
  32. Eckert, D., Ettori, S., Coupon, J., et al. 2016, A&A, 592, A12 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  33. Eckert, D., Ghirardini, V., Ettori, S., et al. 2019, A&A, 621, A40 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  34. Eckert, D., Finoguenov, A., Ghirardini, V., et al. 2020, Open J. Astrophys., 3, 12 [NASA ADS] [CrossRef] [Google Scholar]
  35. Eckert, D., Ettori, S., Pointecouteau, E., van der Burg, R. F. J., & Loubser, S. I. 2022, A&A, 662, A123 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  36. Ettori, S., Ghirardini, V., Eckert, D., et al. 2019, A&A, 621, A39 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  37. Fabian, A. C. 1994, ARA&A, 32, 277 [Google Scholar]
  38. Famaey, B., & Binney, J. 2005, MNRAS, 363, 603 [NASA ADS] [CrossRef] [Google Scholar]
  39. Famaey, B., & McGaugh, S. S. 2012, Liv. Rev. Relat., 15, 10 [Google Scholar]
  40. Farnsworth, D., Rudnick, L., Brown, S., & Brunetti, G. 2013, ApJ, 779, 189 [Google Scholar]
  41. Foreman-Mackey, D., Hogg, D. W., Lang, D., & Goodman, J. 2013, PASP, 125, 306 [Google Scholar]
  42. Fukugita, M., & Peebles, P. J. E. 2004, ApJ, 616, 643 [Google Scholar]
  43. Fusco-Femiano, R., Cavaliere, A., & Lapi, A. 2009, ApJ, 705, 1019 [NASA ADS] [CrossRef] [Google Scholar]
  44. Gerbal, D., Durret, F., Lachieze-Rey, M., & Lima-Neto, G. 1992, A&A, 262, 395 [NASA ADS] [Google Scholar]
  45. Ghirardini, V., Ettori, S., Eckert, D., et al. 2018, A&A, 614, A7 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  46. Ghirardini, V., Eckert, D., Ettori, S., et al. 2019, A&A, 621, A41 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  47. Gonzalez-Garcia, M. C., & Nir, Y. 2003, Rev. Mod. Phys., 75, 345 [NASA ADS] [CrossRef] [Google Scholar]
  48. Govoni, F., Markevitch, M., Vikhlinin, A., et al. 2004, ApJ, 605, 695 [Google Scholar]
  49. Haslbauer, M., Banik, I., & Kroupa, P. 2020, MNRAS, 499, 2845 [Google Scholar]
  50. Henson, M. A., Barnes, D. J., Kay, S. T., McCarthy, I. G., & Schaye, J. 2017, MNRAS, 465, 3361 [NASA ADS] [CrossRef] [Google Scholar]
  51. Hodson, A. O., & Zhao, H. 2017, A&A, 598, A127 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  52. Ilić, S., Kopp, M., Skordis, C., & Thomas, D. B. 2021, Phys. Rev. D, 104, 043520 [Google Scholar]
  53. Katz, H., McGaugh, S., Teuben, P., & Angus, G. W. 2013, ApJ, 772, 10 [NASA ADS] [CrossRef] [Google Scholar]
  54. Kleiner, D., Serra, P., Maccagni, F. M., et al. 2023, A&A, 675, A108 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  55. Kopp, M., Skordis, C., Thomas, D. B., & Ilić, S. 2018, Phys. Rev. Lett., 120, 221102 [NASA ADS] [CrossRef] [Google Scholar]
  56. Lelli, F. 2022, Nat. Astron., 6, 35 [NASA ADS] [CrossRef] [Google Scholar]
  57. Lelli, F., McGaugh, S. S., & Schombert, J. M. 2016, AJ, 152, 157 [Google Scholar]
  58. Lelli, F., McGaugh, S. S., Schombert, J. M., & Pawlowski, M. S. 2017, ApJ, 836, 152 [Google Scholar]
  59. Li, P., Lelli, F., McGaugh, S., & Schombert, J. 2020, ApJS, 247, 31 [Google Scholar]
  60. Li, P., Tian, Y., Júlio, M. P., et al. 2023, A&A, 677, A24 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  61. Loubser, S. I., Babul, A., Hoekstra, H., et al. 2020, MNRAS, 496, 1857 [Google Scholar]
  62. Macquart, J. P., Prochaska, J. X., McQuinn, M., et al. 2020, Nature, 581, 391 [Google Scholar]
  63. McGaugh, S. 2020, Galaxies, 8, 35 [NASA ADS] [CrossRef] [Google Scholar]
  64. McGaugh, S., Lelli, F., Li, P., & Schombert, J. 2020, in Galactic Dynamics in the Era of Large Surveys, eds. M. Valluri, & J. A. Sellwood, 353, 144 [NASA ADS] [Google Scholar]
  65. Milgrom, M. 1983a, ApJ, 270, 371 [Google Scholar]
  66. Milgrom, M. 1983b, ApJ, 270, 384 [Google Scholar]
  67. Milgrom, M. 1983c, ApJ, 270, 365 [Google Scholar]
  68. Milgrom, M. 1984, ApJ, 287, 571 [CrossRef] [Google Scholar]
  69. Milgrom, M. 1994, Ann. Phys., 229, 384 [NASA ADS] [CrossRef] [Google Scholar]
  70. Milgrom, M. 2002, ApJ, 577, L75 [NASA ADS] [CrossRef] [Google Scholar]
  71. Milgrom, M. 2008, Nat. Astron Rev., 51, 906 [NASA ADS] [CrossRef] [Google Scholar]
  72. Milgrom, M. 2009, ApJ, 698, 1630 [NASA ADS] [CrossRef] [Google Scholar]
  73. Milgrom, M. 2010, MNRAS, 403, 886 [NASA ADS] [CrossRef] [Google Scholar]
  74. Milgrom, M. 2014, Scholarpedia, 9, 31410 [NASA ADS] [CrossRef] [Google Scholar]
  75. Milgrom, M. 2018, Phys. Rev. D, 98, 104036 [NASA ADS] [CrossRef] [Google Scholar]
  76. Milgrom, M. 2019, Phys. Rev. D, 99, 044041 [NASA ADS] [CrossRef] [Google Scholar]
  77. Milgrom, M. 2022, Phys. Rev. D, 106, 064060 [NASA ADS] [CrossRef] [Google Scholar]
  78. Milgrom, M. 2023, Phys. Rev. D, 108, 063009 [NASA ADS] [CrossRef] [Google Scholar]
  79. Molendi, S., De Grandi, S., Fusco-Femiano, R., et al. 1999, ApJ, 525, L73 [NASA ADS] [CrossRef] [Google Scholar]
  80. Natarajan, P., & Zhao, H. 2008, MNRAS, 389, 250 [NASA ADS] [CrossRef] [Google Scholar]
  81. O’Hara, T. B., Mohr, J. J., & Guerrero, M. A. 2004, ApJ, 604, 604 [CrossRef] [Google Scholar]
  82. Pezzulli, G., & Cantalupo, S. 2019, MNRAS, 486, 1489 [NASA ADS] [CrossRef] [Google Scholar]
  83. Planck Collaboration XX. 2014, A&A, 571, A20 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  84. Planck Collaboration XXIX. 2014, A&A, 571, A29 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  85. Planck Collaboration VI. 2020, A&A, 641, A6 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  86. Pointecouteau, E., & Silk, J. 2005, MNRAS, 364, 654 [NASA ADS] [CrossRef] [Google Scholar]
  87. Roncarelli, M., Ettori, S., Borgani, S., et al. 2013, MNRAS, 432, 3030 [NASA ADS] [CrossRef] [Google Scholar]
  88. Sanders, R. H. 1994, A&A, 284, L31 [NASA ADS] [Google Scholar]
  89. Sanders, R. H. 1998, MNRAS, 296, 1009 [NASA ADS] [CrossRef] [Google Scholar]
  90. Sanders, R. H. 1999, ApJ, 512, L23 [NASA ADS] [CrossRef] [Google Scholar]
  91. Sanders, R. H. 2003, MNRAS, 342, 901 [NASA ADS] [CrossRef] [Google Scholar]
  92. Sanders, R. H. 2007, MNRAS, 380, 331 [NASA ADS] [CrossRef] [Google Scholar]
  93. Sanders, R. H. 2008, MNRAS, 386, 1588 [NASA ADS] [CrossRef] [Google Scholar]
  94. Shull, J. M., Smith, B. D., & Danforth, C. W. 2012, ApJ, 759, 23 [NASA ADS] [CrossRef] [Google Scholar]
  95. Sijacki, D., Springel, V., Di Matteo, T., & Hernquist, L. 2007, MNRAS, 380, 877 [NASA ADS] [CrossRef] [Google Scholar]
  96. Skordis, C., & Złośnik, T. 2021, Phys. Rev. Lett., 127, 161302 [NASA ADS] [CrossRef] [Google Scholar]
  97. Smith, G. P., Mazzotta, P., Okabe, N., et al. 2016, MNRAS, 456, L74 [Google Scholar]
  98. Spitzer, L., Jr 1956, ApJ, 124, 20 [NASA ADS] [CrossRef] [Google Scholar]
  99. Storm, E., Jeltema, T. E., & Rudnick, L. 2015, MNRAS, 448, 2495 [NASA ADS] [CrossRef] [Google Scholar]
  100. Takahashi, R., & Chiba, T. 2007, ApJ, 671, 45 [NASA ADS] [CrossRef] [Google Scholar]
  101. The, L. S., & White, S. D. M. 1988, AJ, 95, 1642 [NASA ADS] [CrossRef] [Google Scholar]
  102. Thomas, D. B., Kopp, M., & Skordis, C. 2016, ApJ, 830, 155 [NASA ADS] [CrossRef] [Google Scholar]
  103. Tian, Y., Umetsu, K., Ko, C.-M., Donahue, M., & Chiu, I. N. 2020, ApJ, 896, 70 [NASA ADS] [CrossRef] [Google Scholar]
  104. Tian, Y., Yu, P.-C., Li, P., McGaugh, S. S., & Ko, C.-M. 2021, ApJ, 910, 56 [NASA ADS] [CrossRef] [Google Scholar]
  105. Towler, I., Kay, S. T., & Altamura, E. 2023, MNRAS, 520, 5845 [NASA ADS] [CrossRef] [Google Scholar]
  106. Tremaine, S., & Gunn, J. E. 1979, Phys. Rev. Lett., 42, 407 [NASA ADS] [CrossRef] [Google Scholar]
  107. van der Burg, R. F. J., Hoekstra, H., Muzzin, A., et al. 2015, A&A, 577, A19 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  108. Wicker, R., Douspis, M., Salvati, L., & Aghanim, N. 2023, A&A, 674, A48 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  109. Wittenburg, N., Kroupa, P., Banik, I., Candlish, G., & Samaras, N. 2023, MNRAS, 523, 453 [NASA ADS] [CrossRef] [Google Scholar]
  110. Yan, P.-F., Yuan, Q.-R., Zhang, L., & Zhou, X. 2014, AJ, 147, 106 [NASA ADS] [CrossRef] [Google Scholar]
  111. Zhao, H., & Famaey, B. 2012, Phys. Rev. D, 86, 067301 [NASA ADS] [CrossRef] [Google Scholar]

Appendix A: MOND missing mass profiles from a basic subtraction approach

In this section, we estimate the MOND missing mass profile using a basic subtraction approach, similarly to Eckert et al. (2022). Differently from Eckert et al. (2022), however, we use the non-parametric log-normal mixture reconstruction of gobs(r) rather than the parametric Einasto reconstruction. In addition, we investigate the EFE that was neglected by Eckert et al. (2022).

In the isolated case, we consider Eq. 5 with

g N ( r ) = g N , bar ( r ) + g N , mm ( r ) . $$ \begin{aligned} g_{\rm N}(r) = g_{\rm N, bar}(r) + g_{\rm N, mm}(r). \end{aligned} $$(A.1)

Combining Eq. 5 with Eq. 8, we can then infer the radial profile of the missing mass component:

M mm ( < r ) = r 2 G N [ g obs μ ( g obs a 0 ) − g N , bar ] . $$ \begin{aligned} M_{\rm mm}( < r) = \dfrac{r^2}{G_{\rm N}}\left[ g_{\rm obs} \mu \Big (\frac{g_{\rm obs}}{a_0}\Big ) - g_{\rm N, bar}\right]. \end{aligned} $$(A.2)

The resulting mass profiles are shown in Fig. A.1 using the so-called ‘simple’ μ function (Famaey & Binney 2005). Different interpolations function give similar results. The mass profiles are unphysical because they start to decrease after ∼1 Mpc and even become negative in some cases. Similar results were obtained by Eckert et al. (2022). Taking the data at face value, these results would imply that MOND is not a viable theory. There are, however, three important caveats: (1) data at large radii comes from the SZ effect and there could be unknown systematics, (2) Eq. 5 neglects the cosmic EFE from the large-scale matter distribution, and (3) hydrostatic bias may be important beyond 1 Mpc.

thumbnail Fig. A.1.

Enclosed mass profiles for the MOND missing matter using a simplistic subtraction approach (see Sect. A). Dots show data from X-ray observations only, while triangles show those adding the SZ effect. For each cluster, the left panel shows the results for the isolated case with no EFE, while the right panels consider a basic implementation of the EFE from the large-scale structure using MOND external field strengths of 10% (red), 30% (blue), 50% (orange) and 100% (green) a0. In most cases, considering a relatively strong EFE leads to physical (non-decreasing) enclosed mass profiles; exceptions are A644 and A2319 that have clear signatures of merger activity and may be out of dynamical equilibrium.

Regarding the first caveat, Fig. A.1 shows that X-ray data and SZ data give consistent results in the radial range where both data are available, so we suspect that possible systematics between the two different methods play a minor role.

Regarding the second and third caveat, both hydrostatic bias and EFE may flatten the missing mass profiles by decreasing the value of gobs, so their effects are degenerate. The EFE can be approximately taken into account using Eq. 59 of Famaey & McGaugh (2012). Rearranging that equation, we derive a new equation for Mmm:

M mm ( < r ) = r 2 G N × { g obs μ ( g obs + g M , e a 0 ) + g M , e [ μ ( g obs + g M , e a 0 ) − μ ( g M , e a 0 ) ] − g N , bar } , $$ \begin{aligned} \begin{split}&M_{\rm mm}( < r) = \dfrac{r^2}{G_{\rm N}} \times \\&\left\{ g_{\rm obs} \mu \Bigg (\frac{g_{\rm obs}+g_{\rm M, e}}{a_0}\Bigg ) + g_{\rm M, e} \Bigg [\mu \Bigg (\frac{g_{\rm obs}+g_{\rm M, e}}{a_0}\Bigg ) - \mu \Bigg (\frac{g_{\rm M, e}}{a_0}\Bigg ) \Bigg ] - g_{\rm N, bar} \right\} , \end{split} \end{aligned} $$(A.3)

where gM, e is the mean MOND gravitational field due to large-scale mass distribution in the Universe. Eq. A.3 is mathematically equivalent to the EFE implementation in Sect. 3 (see Chae et al. 2020, 2021). However, we stress that gN, e in Eq. 10 is a very different concept than gM, e in Eq. A.3. The former one is the Newtonian external field from the large-scale distribution of baryons, so it can be approximately measured by summing the individual Newtonian contributions of galaxies and galaxy clusters in the nearby Universe (Chae et al. 2021). The latter one is the MOND external gravitational field, which is a difficult quantity to measure because it requires non-linear cosmological calculations of the large-scale structure in MOND. In addition, we stress that these 1D EFE formulas are approximated because they neglect the vectorial nature of the gravitational fields.

We repeated the calculation of Mmm using Eq. A.3, assuming the simple interpolation function and fixing gM, e to be 10%, 30%, 50% and 100% the value of a0. The results are shown in Fig. A.1. In three out of five cases (A1795, A2029, A2142), the EFE can flatten out the missing mass profiles, making them physical. The required external field strength is between 0.3-0.5 a0, which is relatively high but still within the realm of physical reality. Higher values of gM, e make the mass profiles increasing with radius, implying that the missing mass has not converged to a finite value, but there is a saturation effect in Mmm for gM, e ≳ a0. Clearly, these basic calculations do not consider hydrostatic bias, so the values of gM, e are as large as needed.

For A644 and A2319, there is no value of gM, e that can avoid the decreasing profiles for the missing mass component. These two clusters are known to show signatures of merger activity (see Sect. 2.1), so the hot gas at large radii is likely out of dynamical equilibrium. If so, unphysical mass profiles would be a natural outcome of a falsifiable dynamical theory. Actually, they would testify the ability of MOND in identifying merging clusters.

Appendix B: Baryonic scaling

For completeness, we present MOND fits where the observed baryonic contribution is scaled using a single free parameter Υbar, as it is traditionally done in rotation-curve analyses of disc galaxies. In the context of galaxy clusters, these models assume that the missing matter has the same distribution as the observed matter. We consider the isolated MOND case (Eq. 9) and the simple interpolation function. The fit results are shown in Figure B.1 and the best-fit values of Υbar are provided in Table B.1.

thumbnail Fig. B.1.

MCMC fit results for the baryonic-scaling case in isolation. Panels are the same as in Fig. 2.

Table B.1.

Best-fit values of Υbar for the MOND baryonic scaling.

The best-fit values of Υbar suggests that there must be about 3-5 times more matter than observed to reproduce the observations in MOND, but this applies only if the missing matter has the same distribution as the observed matter. The fact that the acceleration profiles are not well fit, instead, indicates that this cannot be the case: the missing mass must have a different distribution that the visible baryons, as we explicitly modeled in Sect. 3.2. In that case, the total amount of missing mass is about 0.4-1 time the amount of visible mass.

Appendix C: Corner plots

Figures C.1 and C.2 show ‘corner plots’ from MCMC fits in the isolated and EFE cases, respectively. For each cluster, the various panels show the posterior probability distribution of pairs of fitting parameters, and the marginalised probability distribution of each fitting parameter. In the inner panels, individual MCMC samples outside the 2σ confidence region are shown with black dots, while binned MCMC samples inside the 2σ confidence region are shown by a greyscale; black contours correspond to the 1σ and 2σ confidence regions; the red squares and red solid lines show median values. In the outer panels (histograms), solid and dashed lines correspond, respectively, to the median and ±1σ values, which are reported at the top of the panel. In general, the posterior probability distributions are well-behaved and show clear peaks, indicating that the fitting quantities and their uncertainties are well measured.

thumbnail Fig. C.1.

Posterior probability distributions from MCMC fits in the isolated case. See Appendix C for details.

thumbnail Fig. C.2.

Posterior probability distributions from MCMC fits in the EFE case. See Appendix C for details.

All Tables

Table 1.

Properties of the cluster sample.

In the text
Table 2.

Estimated parameters from MCMC fits (see Sect. 3.1 for details).

In the text
Table 3.

Results of correlation tests between the 3D temperature of the ICM and the properties of the missing mass component (see Fig. 6).

In the text
Table B.1.

Best-fit values of Υbar for the MOND baryonic scaling.

In the text

All Figures

thumbnail Fig. 1.

MCMC fit results for the isolated case (Eq. (9)). From left to right: acceleration profiles, enclosed mass profiles, circular-velocity profiles, and relative residual difference. In the first three panels, observed data (black points) are fitted with a mass model (black line) that includes the contribution of visible baryons (red line) and missing mass (blue line). We fit the data up to 1 Mpc (black points), where the hydrostatic bias should play a minor role, and then extrapolate the model to larger radii (grey line) to estimate the maximum hydrostatic bias in MOND (rightmost panels).
In the text
thumbnail Fig. 2.

Same as Fig. 1 but for the non-isolated case with a basic modelling of the EFE (Eq. (10)). In this case we fit all data points, even those at r > 1 Mpc, because we aim to estimate the minimum amonut of hydrostatic bias in MOND (rightmost panel).
In the text
thumbnail Fig. 3.

Missing-to-baryonic mass ratio in galaxy clusters as a function of radius in the EFE case (solid line) and isolated case (dashed line). The high values of A644 and A2319 should be interpreted with caution because these two galaxy clusters display merger signatures.
In the text
thumbnail Fig. 4.

Comparison between acceleration profiles from hydrostatic equilibrium (black dots; from X-COP) and Jeans modelling of galaxy kinematics (red dots; from Li et al. 2023). The black line is the MOND EFE model fitted to the X-COP data (same as in Fig. 2). Remarkably, for A1795 and A2142, the MOND model at large radii predicts the behaviour of the stellar kinematic data. For A2029 the comparison between the two datasets is poor, probably due to incompleteness in the galaxy kinematic data and/or unusual galaxy velocity anisotropy (see Li et al. 2023).
In the text
thumbnail Fig. 5.

Locations of galaxy clusters on the radial acceleration relation after accounting for a missing mass component. The top panel corresponds to the isolated case (Eq. (9)), while the bottom panel to the EFE case (Eq. (10)). Grey points show disc galaxies from the SPARC database (Lelli et al. 2016), while coloured points show the five galaxy clusters in our sample. The solid line is the MOND prediction for isolated systems; the dashed line is the line of unity (Newtonian prediction with no DM). In the top panel, open symbols indicate data at R > 1 Mpc, where the mass models with no EFE have been extrapolated.
In the text
thumbnail Fig. 6.

Relations between the 3D temperature of the ICM and the properties of the missing mass component: Total mass Mmm (leftmost panels) and scale radius rs (rightmost panels) for both the isolated and EFE case. In all cases there are statistically significant correlations but we warn that they are based only on five clusters, two of which (A644 and A2319) may be out of dynamical equilibrium.
In the text
thumbnail Fig. A.1.

Enclosed mass profiles for the MOND missing matter using a simplistic subtraction approach (see Sect. A). Dots show data from X-ray observations only, while triangles show those adding the SZ effect. For each cluster, the left panel shows the results for the isolated case with no EFE, while the right panels consider a basic implementation of the EFE from the large-scale structure using MOND external field strengths of 10% (red), 30% (blue), 50% (orange) and 100% (green) a0. In most cases, considering a relatively strong EFE leads to physical (non-decreasing) enclosed mass profiles; exceptions are A644 and A2319 that have clear signatures of merger activity and may be out of dynamical equilibrium.
In the text
thumbnail Fig. B.1.

MCMC fit results for the baryonic-scaling case in isolation. Panels are the same as in Fig. 2.
In the text
thumbnail Fig. C.1.

Posterior probability distributions from MCMC fits in the isolated case. See Appendix C for details.
In the text
thumbnail Fig. C.2.

Posterior probability distributions from MCMC fits in the EFE case. See Appendix C for details.
In the text

Current usage metrics show cumulative count of Article Views (full-text article views including HTML views, PDF and ePub downloads, according to the available data) and Abstracts Views on Vision4Press platform.

Data correspond to usage on the plateform after 2015. The current usage metrics is available 48-96 hours after online publication and is updated daily on week days.

Initial download of the metrics may take a while.