Issue 
A&A
Volume 658, February 2022



Article Number  A26  
Number of page(s)  29  
Section  Extragalactic astronomy  
DOI  https://doi.org/10.1051/00046361/202142060  
Published online  27 January 2022 
Probing the radial acceleration relation and the strong equivalence principle with the Coma cluster ultradiffuse galaxies
^{1}
Université de Strasbourg, CNRS UMR 7550, Observatoire Astronomique de Strasbourg, 67000 Strasbourg, France
email: jonathan.freundlich@astro.unistra.fr
^{2}
European Southern Observatory, KarlSchwarzschildStr. 2, 85748 Garching, Germany
Received:
20
August
2021
Accepted:
22
September
2021
The tight radial acceleration relation (RAR) obeyed by rotationally supported disk galaxies is one of the most successful a priori predictions of the modified Newtonian dynamics (MOND) paradigm on galaxy scales. Another important consequence of MOND as a classical modification of gravity is that the strong equivalence principle (SEP) – which requires the dynamics of a small, freefalling, selfgravitating system not to depend on the external gravitational field in which it is embedded – should be broken. Multiple tentative detections of this socalled external field effect (EFE) of MOND have been made in the past, but the systems that should be most sensitive to it are galaxies with low internal gravitational accelerations residing in galaxy clusters within a strong external field. Here, we show that ultradiffuse galaxies (UDGs) in the Coma cluster do lie on the RAR, and that their velocity dispersion profiles are in full agreement with isolated MOND predictions, especially when including some degree of radial anisotropy. However, including a breaking of the SEP via the EFE seriously deteriorates this agreement. We discuss various possibilities to explain this within the context of MOND, including a combination of tidal heating and higher baryonic masses. We also speculate that our results could mean that the EFE is screened in cluster UDGs. The fact that this would happen precisely within galaxy clusters, where classical MOND fails, could be especially relevant to the nature of the residual MOND missing mass in clusters of galaxies.
Key words: gravitation / dark matter / galaxies: evolution / galaxies: clusters: general / galaxies: clusters: individual: Coma / galaxies: kinematics and dynamics
© J. Freundlich et al. 2022
Open 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.
1. Introduction
Alternatives to general relativity (GR) have been the subject of active research over the past decades (e.g. Clifton et al. 2012, and references therein), the reason being that GR alone cannot explain the observed dynamics on scales ranging from galaxies to the largest scales in the observable Universe without additional degrees of freedom. These are denoted dark matter (DM) and dark energy, and their actual nature is probably the most pressing issue for modern physics. A great deal of research into modified gravity concentrates on the dark energy question, which is related to the cosmological constant problem; however, it has also been suggested since the early 1980s (Milgrom 1983a,b,c, 1984; Bekenstein & Milgrom 1984) that the phenomena attributed to DM might also be, at least partly, related to new gravitational degrees of freedom rather than to new particles in the matter sector as generally assumed.
This approach, known as modified Newtonian dynamics (MOND; see e.g. Sanders & McGaugh 2002; Famaey & McGaugh 2012; Milgrom 2014, for reviews) postulates that the gravitational acceleration g approaches when the Newtonian gravitational acceleration g_{N} falls below a characteristic acceleration scale a_{0} ≈ 10^{−10} m s^{−2} but remains Newtonian above this threshold. This allows one to directly predict the dynamics of galaxies from their baryonic mass distribution alone. This empirical modification of the gravitational law was initially proposed (Milgrom 1983a,b,c) to solve the missing mass problem in the high surface brightness galaxies known at the time, and especially their asymptotically flat circular velocity curves (e.g., Bosma 1978; Rubin et al. 1978; Faber & Gallagher 1979). It is particularly intriguing that this simple recipe has survived almost 40 years of scrutiny at galactic scales, as it has been able to predict the dynamics of a wide variety of galaxies (e.g. Begeman et al. 1991; Sanders 1996, 2019; McGaugh & de Blok 1998; Sanders & Verheijen 1998; de Blok & McGaugh 1998; Sanders & Noordermeer 2007; Gentile et al. 2007a, 2011; Swaters et al. 2010; Famaey & McGaugh 2012; Milgrom 2012; McGaugh & Milgrom 2013a,b), including low surface brightness and dwarf galaxies where internal accelerations can be well below a_{0}, such that the MOND acceleration should a priori deviate significantly from the Newtonian acceleration. This was a core prediction of the original MOND papers and one of its most intriguing successes. MOND can also provide possible answers to various other puzzles in galaxy dynamics, such as the prevalence of bulgeless discs (Combes 2014) and fast bars (Tiret & Combes 2007, 2008; Roshan et al. 2021), or the detailed kinematics of polar ring galaxies (Lüghausen et al. 2013). Generally speaking, there now appears to be a clear and direct connection between the baryonic mass distribution and the rotation curve in most disc galaxies, known as the radial acceleration relation (RAR; McGaugh et al. 2016; Lelli et al. 2017), and this empirical relation is actually indistinguishable from the original MOND prescription (Li et al. 2018). Evaluating whether the RAR holds for all types of galaxies and in all environments is thus of high importance when assessing the viability of MOND as an alternative to particle DM in galaxies.
It is important to note that it was also originally predicted that the nonlinearity of the MOND acceleration would typically lead to a violation of the strong equivalence principle of GR, according to which the internal dynamics of a selfgravitating system embedded in a constant gravitational field should not depend on the external field strength. Within MOND, systems embedded in an external field stronger than their internal one should experience an ‘external field effect’ (EFE; Milgrom 1983b, 2014; Bekenstein & Milgrom 1984; Famaey & McGaugh 2012; McGaugh & Milgrom 2013a,b; Wu & Kroupa 2015; Haghi et al. 2019a) whose consequence is notably that the deviations from Newtonian dynamics are suppressed if the external field is strong enough, and in particular if it is larger than a_{0}. Its influence can be important for the stability and secular evolution of galaxies even when it is weak (Banik et al. 2020), and it can create interesting features such as asymmetric tidal tails of globular clusters (Thomas et al. 2018). The EFE is for instance an observational necessity to allow the dynamics of wide binary stars to remain consistent with MOND (Pittordis & Sutherland 2019; Banik 2019, although see also Hernandez et al. 2022). Because of this EFE, a rotationally supported (pressure supported) system in isolation is expected to have a higher rotational velocity (velocity dispersion) than the same system around a massive host (e.g. Wu et al. 2007; Gentile et al. 2007b; McGaugh & Milgrom 2013a,b; Pawlowski & McGaugh 2014; Pawlowski et al. 2015; McGaugh 2016; Hees et al. 2016; Haghi et al. 2016; Müller et al. 2019; Chae et al. 2020, 2021). In particular, the latter should not follow the RAR, which is contrary to the case of the more isolated systems that should lie on the RAR. This breaking of the strong equivalence principle should be an unmistakable sign of MOND, and it is therefore important to test it for galaxies with internal gravitational accelerations lower than the external field in which they are embedded: this is the focus of the present work.
The need for DM within GR is of course not limited to galaxies. Expanding MOND predictions to the cosmological regime requires a relativistic framework for the paradigm. In order to retain the success of the standard ΛCDM cosmological model on large scales, some hybrid models have been proposed, for instance, where GR is retained but gravity is effectively modified in galaxies through some exotic properties of DM itself, such as in dipolar DM (Blanchet & Le Tiec 2009; Bernard & Blanchet 2015; Blanchet & Heisenberg 2015) or superfluid DM (Khoury 2015; Berezhiani & Khoury 2015; Berezhiani et al. 2018, 2019). More traditional relativistic MOND theories rely on a multifield framework (typically with a scalar and a vector field in addition to the metric), as originally proposed by Bekenstein (2004), but adapted to pass the most recent constraints from gravitational waves (Skordis & Złośnik 2019). It has recently been shown, as a proof of concept, how the angular power spectrum of the cosmic microwave background (CMB) could be reproduced in such a framework (Skordis & Złośnik 2021); the scalar field, which gives rise to the MOND behaviour in the quasistatic limit, also plays the role of DM in the timedependent cosmological regime, thereby providing an analogue to cosmological DM for the CMB. However, the real challenge for such an approach, and for MOND in general, is to explain the mass discrepancy in galaxy clusters. It has indeed long been known that applying the MOND recipe to galaxy clusters yields a residual missing mass problem in these objects (e.g., Sanders 1999, 2003; Pointecouteau & Silk 2005; Natarajan & Zhao 2008; Angus et al. 2008). This is essentially because, contrary to the case of galaxies, there is observationally a need for DM even where the observed acceleration is larger than a_{0}, meaning that the MOND prescription is not enough to explain the observed discrepancy. In the central parts of clusters, the ratio of MOND dynamical mass to observed baryonic mass can reach a value of 10. This cluster missing mass problem extends to giant ellipticals residing at the centre of clusters (Bílek et al. 2019a). It is also clear that this residual missing mass must be collisionless (Clowe et al. 2006; Angus et al. 2007), and it has hence been proposed that it could be made of cold, dense molecular gas clouds (Milgrom 2008) or some form of hot dark matter (HDM) such as sterile neutrinos, which would not condense on galaxy scales (Angus et al. 2010; Haslbauer et al. 2020). In such cases, the residual missing mass should be an important gravitational source contributing to the EFE acting on galaxies residing in clusters. On the other hand, if the residual MOND missing mass problem would itself be a gravitational phenomenon, it would then not necessarily contribute to the EFE as a source. Therefore, studying the dynamics of galaxies residing in galaxy clusters, and in particular whether the EFE can be detected there, should provide powerful constraints for relativistic modelbuilding in the MOND context and also illuminate our understanding of scaling relations with environment in the cold dark matter (CDM) paradigm. Galaxies with a very low internal gravity, hence ultradiffuse ones, are best suited for such a study.
Ultradiffuse galaxies (UDGs) are low central surface brightness (μ_{g, 0} > 24 mag arcsec^{−2}) objects with optical luminosities typical of dwarf galaxies (L ∼ 10^{7} − 10^{8} L_{⊙}) but effective radii comparable to that of the Milky Way (r_{eff} > 1.5 kpc). Such galaxies have been identified for several decades (e.g. Fosbury et al. 1978; Sandage & Binggeli 1984; Impey et al. 1988; Karachentsev et al. 2000) but have recently undergone a revival of interest as deep imaging observations revealed their ubiquity both in groups and clusters (e.g. van Dokkum et al. 2015a,b; Janowiecki et al. 2015; Koda et al. 2015; Mihos et al. 2015, 2017; Muñoz et al. 2015; MartínezDelgado et al. 2016; van der Burg et al. 2016; Yagi et al. 2016; Merritt et al. 2016; Janssens et al. 2017; Shi et al. 2017; Venhola et al. 2017; Müller et al. 2018; Zaritsky et al. 2019; Habas et al. 2020) and in the field (e.g. Leisman et al. 2017; Román & Trujillo 2017; Prole et al. 2019). Different scenarios have been proposed to explain their formation, invoking them as either (i) failed Milky Waylike galaxies that lost their gas after forming their first stars (van Dokkum et al. 2015a,b, 2016; Yozin & Bekki 2015; Peng & Lim 2016; Lim et al. 2018; Pandya et al. 2018; Villaume et al. 2021), possibly due to denser group and cluster environments; (ii) the highspin tail of the dwarf galaxy population (Amorisco & Loeb 2016; Rong et al. 2017; Shi et al. 2017); (iii) tidal debris from mergers or tidally disrupted dwarfs (Beasley & Trujillo 2016; Greco et al. 2018; Toloba et al. 2018; Jiang et al. 2019; Carleton et al. 2019); or (iv) galaxies whose spatial extent results from dynamical heating due to stellar feedback (Di Cintio et al. 2017; Chan et al. 2018; Jiang et al. 2019; Freundlich et al. 2020a,b). Within the CDM model of structure formation, their particularly low surface brightness allows us to probe the galaxyhalo connection in the lowmass regime, with an ongoing debate regarding their DM content (Beasley et al. 2016; Peng & Lim 2016; van Dokkum et al. 2016, 2017, 2018, 2019a; Martin et al. 2018; Lim et al. 2018; Wasserman et al. 2018; Nusser 2018, 2019; Hayashi & Inoue 2018; Dutta Chowdhury et al. 2019, 2020; Laporte et al. 2019; Danieli et al. 2019; Emsellem et al. 2019; Fensch et al. 2019; Trujillo et al. 2019; Haslbauer et al. 2019; Gannon et al. 2020, 2021; Müller et al. 2021). The possible diversity of their DM and gas content (Janowiecki et al. 2015; Papastergis et al. 2017; Leisman et al. 2017; Shi et al. 2017; Karunakaran et al. 2020) may indicate that more than one of the formation scenarios coexist (cf. also Toloba et al. 2018; RuizLara et al. 2018; FerréMateu et al. 2018; Chan et al. 2018; Jiang et al. 2019; Forbes et al. 2020; Sales et al. 2020).
Ultradiffuse galaxies in clusters provide a testing ground for MOND and the EFE given the singularly low internal accelerations stemming from their low surface brightness and the strong external field. The small velocity dispersion observed in the two group UDGs NGC 1052DF2 and NGC 1052DF4, inferring dynamical masses close to their stellar masses, was initially interpreted as a challenge for MOND (van Dokkum et al. 2018, 2019a). Indeed, the dynamical effect attributed to DM in the CDM model, and to a modification of the gravitational law within MOND in isolation, would be absent. However, taking the EFE into account removes or significantly lessens the tension (Famaey et al. 2018; Kroupa et al. 2018; Müller et al. 2019; Haghi et al. 2019a). On the other hand, the large velocity dispersion of the Coma cluster UDG DF44 (van Dokkum et al. 2016, 2019b) and its relative agreement with the isolated MOND prediction without EFE has been used to place constraints on its distance from the cluster centre within MOND, or on a potential need for an additional baryonic mass (Bílek et al. 2019b; Haghi et al. 2019b). Different approaches have been used to take the EFE into account in this context: Kroupa et al. (2018) and Haghi et al. (2019a) used fitting functions for the onedimensional lineofsight velocity dispersion in an external field stemming from MOND Nbody simulations by Haghi et al. (2009), Famaey et al. (2018) and Müller et al. (2019) used the onedimensional analytical expression for the acceleration field in the presence of an external field from Famaey & McGaugh (2012), Eq. (59), together with the Wolf et al. (2010) relation for the lineofsight velocity dispersion. Bílek et al. (2019b) and Haghi et al. (2019b) did not quantitatively assess the EFE on the velocity dispersion for DF44, since the data required as little of it as possible. We propose examining the question of the EFE in this galaxy more quantitatively and expanding the study to a larger sample of UDGs.
In this article, we compare observed stellar velocity dispersion measurements of 11 Coma cluster UDGs with MOND modified gravity predictions with and without EFE. We also propose a new analytic formula in quasilinear MOND (QUMOND, Milgrom 2010), assuming a spherical selfgravitating system embedded in a constant external field, which is tested against numerical computations and compared with existing formulae. Throughout this work, we adopt a distance of 100 Mpc, a distance modulus (m − M)_{0} = 35, and central coordinates RA: 12:59:48.75 and Dec: +27:58:50.9 (J2000) (e.g., Carter et al. 2008; Alabi et al. 2020) for the Coma cluster. In Sect. 2, we present the sample of Coma cluster UDGs with stellar velocity dispersion measurements. In Sect. 3, we model their dynamics within MOND in isolation. In Sect. 4, we present our new QUMOND formula for the EFE and model these UDGs including this effect. We discuss the implications of the confrontation between the MOND predictions and the observations in Sect. 5, and we conclude in Sect. 6.
2. A sample of Coma cluster UDGs
The Coma cluster is one of the nearest and moststudied rich and dense cluster environments. van Dokkum et al. (2015a) notably coined the term ultradiffuse galaxies when they uncovered its large population of low surface brightness objects with the Dragonfly Telephoto Array (Abraham & van Dokkum 2014), and the first spatially resolved kinematics of a UDG was obtained for a Coma cluster UDG: DF44 (van Dokkum et al. 2019b).
2.1. DF44
The velocity and velocity dispersion profiles of DF44 were obtained with Keck Cosmic Web Imager (KCWI, Morrissey et al. 2012, 2018), an integral field unit (IFU) spectrograph on the Keck II telescope optimised for low surface brightness objects. These profiles were studied in the context of Newtonian gravity assuming CDM by van Dokkum et al. (2019b), in the context of MOND by Bílek et al. (2019b) and Haghi et al. (2019b). Wasserman et al. (2019) further considered Newtonian gravity with fuzzy dark matter, an exotic form of DM (cf. e.g. Hu et al. 2000). The relatively high velocity dispersion profile, with a luminosityweighted km s^{−1} at the halflight radius, is interpreted within the CDM model as indicative of a galaxy dominated by DM even at its centre (van Dokkum et al. 2019b). Bílek et al. (2019b) and Haghi et al. (2019b) showed that this velocity dispersion profile fits well with the MOND prediction in isolation. These two studies used the van Dokkum et al. (2015a, 2016, 2017, 2019b) Iband luminosity L_{I} = 3.0 ± 0.6 × 10^{8} L_{⊙} and projected halflight radius R_{e} = 4.7 kpc of DF44 to infer the gravitational acceleration within MOND and, from it, the lineofsight velocity dispersion profile through Jeans modelling. Bílek et al. (2019b) assumed a spherical Sérsic profile given the projected Sérsic index n = 0.94, fixed the stellar masstolight ratio to 1.3, and considered different anisotropy parameters. Haghi et al. (2019b) assumed a spherical Plummer profile and left both the dynamical masstolight ratio and the anisotropy parameter free. Bílek et al. (2019b) interpreted the agreement between the velocity dispersion profile with the MOND prediction in isolation as an indication that DF44 is either far from the cluster centre (the projected distance is 1.8 Mpc), undergoing disruption, or embedded in a dark baryonic halo. It was also noted that some degree of tangential anisotropy allowed a better fit of the outer data points, which is also true in the CDM context.
2.2. Ten additional UDGs
While van Dokkum et al. (2015a) initially reported 47 UDGs in the Coma cluster using the Dragonfly Telephoto Array, Koda et al. (2015) and Yagi et al. (2016) increased the number to 854 using Subaru archival data, and Alabi et al. (2020) discovered a further 29 of them. Amongst all Coma cluster UDGs, about 30 have been spectroscopically confirmed (Kadowaki et al. 2017; Alabi et al. 2018; RuizLara et al. 2018; Chilingarian et al. 2019). Most lack stellar velocity and velocity dispersion measurements, which are particularly challenging to obtain given their low surface brightness, such that internal dynamics and DM content within the CDM model are generally poorly constrained. Chilingarian et al. (2019) obtained spatially resolved stellar velocity profiles and velocity dispersions for nine UDGs in the Coma cluster, seven of them being taken from the Yagi et al. (2016) catalogue, using the Binospec spectrograph on the 6.5 mm MMT telescope (Fabricant et al. 2019). They modelled them within the CDM model, obtaining DM fractions between 50 and 90% within the halflight radius. We note that amongst the nine UDGs for which Chilingarian et al. (2019) obtained spatially resolved velocity profiles, only three have spatially resolved velocity dispersion profiles. Together with DF44, this thus yields ten Coma cluster UDGs with spatially resolved stellar velocity profiles and velocity dispersion measurements, four of them having spatially resolved velocity dispersion profiles. We note that one galaxy amongst them (J130005.40+275333.0) displays a 20 km s^{−1} rotation, and three others (J125848.94+281037.1, J125904.06+281422.4, and J130028.34+274820.5) show hints of rotation (< 10 km s^{−1}), always at least twice as low as the typical level of the lineofsight velocity dispersion.
To this sample, we further add DFX1, a UDG whose central velocity dispersion is reported in van Dokkum et al. (2017) together with its Iband Sérsic index, halflight radius, and magnitude. We do not include UDGs for which the velocity dispersion has been measured only from the globular cluster population, notably as some disagreement between stellar and globular cluster measurements has been reported (cf. for example Müller et al. 2020 versus Forbes et al. 2021).
2.3. Sample properties
In the present article, we model the 11 Coma cluster UDGs with velocity dispersion measurements within MOND, namely DF44 and the ten UDGs mentioned above. DF44 was already modelled within MOND by Bílek et al. (2019b) and Haghi et al. (2019a), but the EFE was not quantitatively assessed, which justifies its inclusion in our current study.
Figure 1 shows the Coma cluster as seen by Dragonfly in the r and g bands together with the positions of the UDGs at stake. For DF44 and DFX1, we retain the Iband luminosity, projected halflight radius, and Sérsic index from van Dokkum et al. (2015a, 2016, 2017). van Dokkum et al. (2019b) obtained stellar masstolight (M/L) ratios in the 1−1.5 range. We assume here M/L = 1.3 as Bílek et al. (2019b). For the Chilingarian et al. (2019) UDGs, we used the bestfit Rband luminosities, projected halflight radii, Sérsic indices from the Yagi et al. (2016) catalogue when available, and the stellar M/L ratios obtained by Chilingarian et al. (2019) by fitting the full Binospec spectrum with simple stellar population models. We did not leave the stellar masstolight ratio as a free parameter, so our results are direct predictions, not fits. For the two UDGs identified by Chilingarian et al. (2019) that are not in the Yagi et al. (2016) catalogue, we used the Rband luminosities and projected halflight radii obtained by Chilingarian et al. (2019), and we assumed Sérsic indices n = 1. We note that while van Dokkum et al. (2019b) assumed a distance of 100 Mpc to the Coma cluster, Chilingarian et al. (2019) used 99 Mpc, but we neglected the 1% difference. Table 1 summarises the properties of the sample UDGs, including the luminosityweighted effective stellar velocity dispersion σ_{eff}.
Fig. 1. Spatial distribution of the Coma cluster UDGs modelled in this article projected on a colour image of the cluster created from the Dragonfly g and r images (Abraham & van Dokkum 2014; van Dokkum et al. 2015a). The adopted centre of the Coma cluster is indicated by the white cross. DF44 and DFX1, which are the furthest away from the centre, are indicated in red, the Chilingarian et al. (2019) UDGs are displayed in orange. 
Sample UDGs and their properties.
The sample distribution of radial velocities (which can be found in van Dokkum et al. 2015b, 2017; Chilingarian et al. 2019) is centred at 6654 km s^{−1} with a dispersion of 922 km s^{−1}. This is in line with typical Coma cluster galaxies (e.g. Colless & Dunn 1996, with galaxies centred at 6853 km s^{−1} and a dispersion of 1082 km s^{−1}), such that there is a priori no reason to believe the sample kinematics are biased.
3. MOND models in isolation
3.1. Radial acceleration relation
Since the radial acceleration relation (RAR; McGaugh et al. 2016; Milgrom 2016; Lelli et al. 2017; Li et al. 2018) between the radial dynamical acceleration at the halflight radius, inferred from the velocity dispersion, and the radial gravitational acceleration predicted by the observed baryonic distribution within this radius, is a core prediction of MOND in isolation, we first check whether the Coma cluster UDGs do fall on this relation or not. Within MOND, the tight scatter of this relation linking the dynamical and baryonic masses of galaxies over a wide range of galaxy types, sizes, and morphologies is indeed a consequence of the modification of the gravitational law.
3.1.1. Observed versus inferred acceleration
From the effective stellar velocity dispersion σ_{eff} within the deprojected halflight radius r_{1/2} ≈ (4/3)R_{e}, where R_{e} is the halflight radius on the sky, we derived the Newtonian dynamical mass within r_{1/2} using the Wolf et al. (2010) formula:
(cf. also Kretschmer et al. 2021 for a discussion on the proportionality factor based on cosmological zoomin simulations), and we estimated the corresponding ‘observed’ gravitational acceleration:
This gravitational acceleration can be compared to that due to the baryons stemming from the stellar distribution M(r):
We assume that the uncertainties on g_{obs} are dominated by those on σ_{eff}, i.e., dg_{obs}/g_{obs} = 2dσ_{eff}/σ_{eff}, and that the uncertainties on g_{bar} are a combination of those on the M/L ratio and a conservative 30% uncertainty on the luminosity.
We note that this approach is obviously less precise than the detailed Jeans modelling that we performed afterwards, but it gives a first indication as to whether the Coma UDGs agree with the isolated MOND predictions or not.
3.1.2. UDG spherical mass model
To describe the threedimensional stellar mass distribution of each UDG, we assumed a uniform stellar masstolight M/L ratio and deprojected the twodimensional Sérsic light profile using the semianalytical approximation proposed by Lima Neto et al. (1999):
where n is the Sérsic index, a = R_{e}exp[0.1789 − 0.6950 × n − nlnn], and p = 1 − 0.6097/n + 0.05563/n^{2}. These expressions are given in their Eqs. (13) and (19), and in Appendix A of Márquez et al. (2000) for the updated value of the 1/n^{2} term of p, which was already used in Bílek et al. (2019b). The twodimensional Sérsic density profile indeed lacks an actual analytical deprojection in three dimensions, and the Lima Neto et al. (1999) approximation is accurate within 5% for n between 0.5 and 10 and radii between 10^{−2} and 10^{3} R_{e}. The corresponding enclosed mass profile is
where is the incomplete gamma function, such that ρ_{0} = M_{0}/[4πna^{3}Γ(3n − np)], given the total mass M_{0}.
3.1.3. Resulting RAR
Table 2 indicates the values of g_{obs} and g_{bar} for the sample UDGs; Fig. 2 displays them on the radial acceleration relation obtained by McGaugh et al. (2016) for local spiral galaxies from the Spitzer Photometry and Accurate Rotation Curves (SPARC) sample (Lelli et al. 2016a) and for local dwarf spheroidals (dSphs) by Lelli et al. (2017), including satellites of the Milky Way and of Andromeda. As for the latter, the UDG data points are generally above but in reasonable agreement with the following empirical fitting function:
Fig. 2. Radial acceleration relation for the Coma cluster UDGs (magenta error bars) together with that obtained for SPARC local spiral galaxies by McGaugh et al. (2016, blue colour scale, red squares, and dashed lines for the binned mean and standard deviation) and for dSphs of the Local Group (Lelli et al. 2017, light grey error bars for their least reliable data and dark grey error bars for their better quality sample). The plain black line shows the fit using Eq. (6) obtained by McGaugh et al. (2016) on the SPARC sample of spiral galaxies. The Coma cluster UDGs fall within the scatter obtained for the dSphs and in reasonable agreement with the empirical function of McGaugh et al. (2016), which reproduces MOND phenomenology. 
Observed and baryonic accelerations at the halflight radius for the radial acceleration relation.
with a_{0} = 1.20 × 10^{−10} m s^{−2} obtained for the local spirals and shown as a plain black line (McGaugh 2008; McGaugh et al. 2016; Lelli et al. 2017). This function reproduces the MOND phenomenology in isolation: g_{obs} ≈ g_{bar} when g_{bar} ≫ a_{0} and when g_{bar} ≪ a_{0}. The fact that the Coma cluster UDGs fall on the same radial acceleration relation as spirals and dSphs not only shows that their dynamical and baryonic masses are similarly correlated (cf. also McGaugh et al. 2000, 2021; McGaugh 2005, 2012; Lelli et al. 2016b, 2019 on the baryonic Tully–Fisher relation), but it also hints towards the UDGs being well described by MOND in isolation. We now turn to a more detailed Jeans modelling of the velocity dispersion profiles to confirm this rough first hint.
3.2. MOND velocity dispersion in isolation
We compare the observed velocity dispersion measurements with the profiles expected within MOND, assuming that the galaxies are sufficiently far away from the centre of the Coma cluster and from other galaxies to neglect the EFE, that they can be approximated as spherical, and that M/L is uniform for each galaxy. We followed similar steps to those of Bílek et al. (2019b) to determine the expected velocity dispersion profiles in isolation, which are summarised below. We highlight that the predicted profiles are obtained without any free fitting parameters; they depend only on the stellar distribution, approximated by a Sérsic sphere, the fixed M/L, and the assumed anisotropy. The parameters describing the stellar distribution and M/L (fixed at its value from the literature) are indicated in Table 1. We neglect any gaseous component.
3.2.1. Lineofsight velocity dispersion
The velocity dispersion of a spherical collisionless system in equilibrium is related to the gravitational acceleration g through the Jeans equation (Binney & Tremaine 2008, Eq. (4.215)):
where ρ is the density, σ_{r} the radial velocity dispersion, and the anisotropy parameter, with σ_{t} being the tangential velocity dispersion. The lineofsight velocity dispersion σ_{los}, which is the observable quantity, can be expressed as a function of σ_{r} by projecting the velocity ellipsoid along the line of sight, with
(Binney & Mamon 1982), where Σ is the surface density. Given g, it is thus possible to recover σ_{los} by injecting Eq. (7) into Eq. (8). The appendix of Mamon & Łokas (2005) provides formulas expressing σ_{los} as a single integral for specific anisotropy profiles, including the case of a uniform β, on which we focused in this work.
3.2.2. MOND phenomenology
In spherical symmetry, MOND modified gravity reduces to the following relation, sometimes dubbed Milgrom’s relation:
where ν is an interpolating function such that ν(y)→1 for y ≫ 1 and ν(y)→y^{−1/2} for y ≪ 1, and g_{N} is the Newtonian gravitational field deduced from the baryonic mass distribution (e.g. Famaey & McGaugh 2012). Although yν(y) should be monotonically increasing for g and g_{N} to be equivocally determined, MOND allows some freedom in the shape of the interpolating function (Milgrom 1983b). In the following, we adopt the function
corresponding to Eq. (6), which seems to be favoured by galaxy rotation curves (Famaey & McGaugh 2012). However, Hees et al. (2016) noted that this function is excluded in the strong acceleration regime by Solar System data. Therefore, we also considered the more general family of functions (e.g. Hees et al. 2016):
for which α = 0.5 gives the previously adopted interpolating function, and for which α = 2 safely passes the Solar System constraints. We also consider the socalled simple interpolating function
which yields good fits in the intermediate and weak gravity regime of galaxies (Famaey & Binney 2005; Zhao & Famaey 2006; Sanders & Noordermeer 2007; Gentile et al. 2011). We note that the adopted interpolating function (Eq. (10) or Eq. (11) with α = 0.5) is very similar to the simple interpolating function in the intermediate and weak gravity regime, as can be seen in Fig. 3.
Fig. 3. Comparison between the three MOND interpolating functions considered in this work, namely ν from Eq. (10), the simple interpolating function from Eq. (12), and ν_{2} from Eq. (11) with α = 2. The fiducial ν is very similar to the simple interpolating function. 
As in Sect. 3.1, we assumed a fixed and uniform M/L and deprojected the twodimensional Sérsic light profile using the semianalytical approximation proposed by Lima Neto et al. (1999) into the spherical threedimensional mass profile M(r) given by Eq. (5). Assuming the UDG to be isolated, the Newtonian acceleration is simply g_{N}(r) = GM(r)/r^{2}, from which we can derive the MOND acceleration g using Eq. (9) and the lineofsight velocity dispersion using the Mamon & Łokas (2005) formulas.
3.2.3. Resulting velocity dispersion profiles
Figure 4 compares the observed velocity dispersion measurements with the MOND predictions in isolation given the stellar distribution, as outlined in the previous sections. For each UDG, we consider uniform anisotropy parameters β = 0 (isotropic, in black), β = +0.5 (radially biased, in magenta), and β = −0.5 (tangentially biased, in blue). We also indicate the Newtonian prediction given the stellar distribution without dark matter (dotted lines). The interpolating function adopted was that of Eq. (10), namely ν_{0.5} in the terminology of Eq. (11), but we also considered ν_{2} and the simple interpolating function of Eq. (12), yielding similar results because these galaxies are all in the weakacceleration regime where all interpolating functions are similar. Given that there is no free parameter involved on each curve, the figure shows a remarkably good agreement between the observed velocity dispersion measurements and the isolated MOND models. The only parameter that has been changed from one curve to the other is the anisotropy, and we note that the agreement is even better for a slight radial anisotropy, β = +0.5. Such a radial anisotropy may inform us about the formation of these galaxies in MOND (cf. Sect. 5.1). As expected, the Newtonian models without dark matter largely underestimate the velocity dispersions, for example by more than 3σ in the case of DF44, as was already reported by Bílek et al. (2019b). We noted that four galaxies display some – or a hint of – rotation (cf. Sect. 2.2). Taking this into account in detail would largely complicate the present analysis, but the magnitude of the correction would typically only be of the order of 1 km s^{−1} or less. In other words, given the level of rotation hinted at in these four galaxies, the same galaxies without rotation would typically see an increase of their velocity dispersion data by about 1 km s^{−1}, which would obviously not change our conclusions.
Fig. 4. Comparison between the measured lineofsight velocity dispersion of the sample UDGs (grey error bars and/or shaded area), the MOND model in isolation (plain lines), and the Newtonian model (dotted lines). Black corresponds to a uniform anisotropy parameter β = 0 (isotropic), magenta to β = +0.5 (radial), and blue to β = −0.5 (tangential). Data points stem from van Dokkum et al. (2019b) for DF44, van Dokkum et al. (2017) for DFX1, and from Chilingarian et al. (2019) for the other UDGs. The shaded area corresponds to the effective velocity dispersion or, when not available (in the case of DFX1), central stellar velocity dispersion. 
Fig. 4. continued. 
4. MOND models with external field effect
4.1. Mass and UDG distributions of the Coma cluster
To assess the MOND velocity dispersion with EFE for the Coma cluster UDGs, one first needs to estimate the external gravitational field they inhabit, stemming from the cluster mass distribution. We derived the mass profile M_{C}(R) of the Coma cluster within MOND as in Sanders (2003) from the hydrostatic equilibrium of the Xrayemitting hot gas, assumed to be isothermal. Using the simple interpolation function (Eq. (12)) instead of that used by Sanders (2003) in their Eq. (8), since Fig. 3 shows how closely the simple interpolation function follows our fiducial interpolation function (Eq. (10)), the mass profile yields
where
is the gravitational acceleration,
is the Newtonian mass profile, and the intracluster Xrayemitting hot gas density is described by a βmodel:
with an exponent β_{C} = 0.71, a temperature kT_{C} = 8.6 keV, a mean atomic weight μ = 0.61, a characteristic radius r_{C} = 276 kpc, a central density ρ_{C0} = 9.0 × 10^{4} M_{⊙} kpc^{−3}, and m_{p} being the proton mass (Reiprich 2001). Figure 5 shows the corresponding mass profile. The external Newtonian acceleration at a distance R from the centre of the Coma cluster is g_{Ne}(R) = − GM_{C}(R)/R^{2}.
Fig. 5. Coma cluster cumulative dynamical mass derived by Sanders (2003) from the hydrostatic equilibrium of the Xrayemitting hot gas in Newtonian dynamics (dotted black line) and in MOND (plain black line). While the Newtonian dynamical mass (M_{NC}) continues to increase above 1 Mpc, the MOND mass (M_{C}) converges towards 5.6 × 10^{14} M_{⊙}. 
UDGs are observed at given projected distances, but their actual distances from the cluster centre may be significantly larger. van der Burg et al. (2016) shows that the projected radial distribution of UDGs in clusters is generally steeper than NFW at the outskirts, flatter than NFW towards the centre, and on average well fit with an Einasto (1965) profile:
with concentration c_{Ein} ≡ R_{200}/r_{−2} = 1.83 and exponent α_{Ein} = 0.92. Assuming that Coma cluster UDGs follow this distribution enables one to retrieve the most probable distance and the average distance of a UDG from the cluster centre given its projected distance. The probability for an object at projected distance d to be between R and R + dR given the threedimensional distribution ρ_{Ein} is
from which the mode and the average distance can be derived. We note that the average distance d_{mean} is always larger than the mode, both being larger than the projected distance d. The resulting average distances are indicated in Table 1. We used R_{200} = 2.94 Mpc (Kubo et al. 2007).
In Fig. 6, we show the projected distance distribution of the sample UDGs (magenta) together with that of the parent Coma cluster UDG population from the Yagi et al. (2016) catalogue (black) and that corresponding to the van der Burg et al. (2016) Einasto distribution (dashed black), which was obtained for a sample of eight nearby galaxy clusters excluding Coma. The Einasto distribution provides a reasonable fit to the Yagi et al. (2016) distribution, hence validating our estimate of the average distance d_{mean}. For comparison, the figure also shows the distribution that would result from having no UDGs below 10 Mpc and a UDG density decreasing as 1/r^{3} beyond (dotted blue). The increase of this latter curve is dominated by the increase of the area spanned by the linearly spaced radial bins. A uniform distribution within 10 Mpc would thus yield a similar curve in the regime shown in the plot (d < 2.5 Mpc). The discrepancy between the two black curves and the dotted blue curve together with the van der Burg et al. (2016) Einasto distribution generally advocates for the UDG distribution to peak towards the centre rather than being pushed beyond 10 Mpc or remaining constant.
Fig. 6. Distribution of the projected distance of the sample UDGs (plain red) and of the parent Yagi et al. (2016) Coma cluster UDG catalogue (plain black) compared to the average Einasto profile derived by van der Burg et al. (2016, dashed black) and a distribution with no UDGs within 10 Mpc (dotted blue). A constant UDG density within 10 Mpc would closely follow the latter, given the xaxis range. The van der Burg et al. (2016) Einasto distribution provides a reasonable fit to the Coma cluster UDGs; a distribution without UDGs within 10 Mpc or with a constant UDG density are disfavoured. 
4.2. An analytic formula for the EFE in QUMOND
In Newtonian dynamics, the internal dynamics of a selfgravitating system embedded in a constant gravitational field does not depend on the external field strength because of linearity. This is not the case in MOND, where nonNewtonian effects appear where the absolute values of the internal acceleration g and of the external acceleration g_{e} are both less than a_{0}. In MOND, a system with g_{e} < g < a_{0} presents a deep MOND behaviour, a system with g < a_{0} < g_{e} is Newtonian, and a system with g < g_{e} < a_{0} is Newtonian with a renormalised gravitational constant, depending on g_{e}. Famaey & McGaugh (2012, Eqs. (59) and (60)) provide a formula expressing the MOND internal gravitational field g as a function of the internal Newtonian field g_{N} and the external field g_{e} (or its Newtonian counterpart g_{Ne}) in one dimension, i.e., when g and g_{e} are aligned, namely, for the case of QUMOND:
It turns out that simulations within the Bekenstein & Milgrom (1984) version of MOND tend to agree well with this estimate (Haghi et al. 2019a), which yields a rather strong effect of the external field on the internal dynamics. However, in the QUMOND version (Milgrom 2010), this onedimensional formula would a priori appear to be an upper limit on the impact of the EFE.
In the context of QUMOND, we hereafter propose a simple approximate analytical formula for spherical systems by averaging the typical effect over the sphere. We validate the formula using a numerical Poisson solver in QUMOND in Appendix A. We consider a spherical system, characterised by a internal gravitational field g embedded in a constant external field g_{e}. Following Eq. (9) in the QUMOND formalism, the external field can be expressed in terms of its Newtonian counterpart g_{Ne} through
while the total field yields
where g_{N} is the Newtonian internal field. Hence, the internal gravitational field
whose intensity depends on the position on the sphere. Using the standard spherical coordinates such that and with , as shown in Fig. 7, we have
Fig. 7. Spherical selfgravitating system embedded in a constant external field. 
and
We now average this expression over the sphere of a given radius r and approximate ⟨ν(x)⟩ by ν(⟨x⟩). Integrating over the sphere yields an average
This norm yields when g_{N} ≥ g_{Ne} and when g_{N} ≤ g_{Ne}. Approximating
we obtain the average radial and tangential internal gravitational fields
and
The expression for ⟨g_{r}⟩ yields
from which we can retrieve the MOND behaviour when g_{Ne} < g_{N} < a_{0} and the Newtonian behaviour with a renormalised gravitational constant when g_{N} < g_{Ne} < a_{0}. As already noted, for instance by Banik & Zhao (2018), an unusual consequence of the MOND EFE is that the gravitational attraction due to a given mass may not always be directed towards it. The expression for ⟨g_{r}⟩ can be injected in the Jeans equation to retrieve the average velocity dispersion at any radius.
Equations (22) and (23) with cos θ = 0 further enable us to express the internal gravitational field in the case where g and g_{e} are perpendicular:
The second term in the square root of Eq. (30) is small compared to the first one, such that
Given that the form of Eq. (27), ⟨g_{r}⟩ is close to g_{⊥}: the average radial acceleration field approaches the perpendicular one. Thus, either Eq. (29) or Eq. (31) can be used to retrieve the internal acceleration field with EFE in QUMOND.
Figure 8 shows the behaviour of the different EFE formulae as a function of the Newtonian external field, for a given typical internal Newtonian acceleration g_{N} = a_{0}/100 (cf. Fig. A.1). The different formulas transition from the deep MOND regime (g = ν(g_{N}/a_{0})g_{N}, highlighted as a magenta horizontal line) to the Newtonian regime (g = g_{N}). Since g_{∥} assumes the fields to be aligned in the same direction, it systematically yields a lower gravitational acceleration than ⟨g_{r}⟩≈g_{⊥}. Therefore, using our QUMOND formula minimises the EFE compared to g_{∥}, which might be closer to estimates in the Bekenstein & Milgrom (1984) formulation (Haghi et al. 2019a).
Fig. 8. MOND gravitational acceleration g stemming from the different analytic formulae as a function of the Newtonian external field g_{Ne}, for a given internal Newtonian acceleration g_{N} = a_{0}/100. The magenta data points correspond to the radial acceleration field where g_{N} = a_{0}/100 in the QUMOND numerical integrations plotted in Fig. A.1, which consider a Plummer sphere at distances d = 0.5, 1, 2, 5, 10, and 20 Mpc (from right to left) from a point mass representing the Coma cluster. We refer the reader to Fig. A.1 for the values of g_{Ne} corresponding to each distance d within the Coma cluster. The grey line corresponds to a numerical integration of the average radial acceleration of Eq. (24). 
We refer the reader to Appendix A for a validation of our formula in the context of a numerical Poisson solver in QUMOND, as developed by Oria et al. (2021). Points corresponding to the acceleration where g_{N} = a_{0}/100 in the QUMOND numerical integrations shown in Fig. 8 (in magenta) make it possible to visualise the success of Eq. (29) to approximate the QUMOND acceleration with EFE throughout the transition from the deepMOND, lowg_{Ne} regime to the Newtonian, highg_{Ne} regime. We also compute the mean radial acceleration of Eq. (24) numerically (grey line), i.e., integrating over θ numerically, and we find that it overlaps perfectly with the data points from the Poisson solver. The expressions for ⟨g_{r}⟩≈g_{⊥} only very slightly underestimate the EFE, and thus provide a very good estimate of the acceleration field and of the resulting velocity dispersion in QUMOND. All other formulae will generate a stronger EFE.
4.3. Application to the Coma cluster UDGs
We applied Eq. (29) and the Mamon & Łokas (2005) formulas to derive the velocity dispersion profiles predicted in MOND with EFE for the Coma cluster UDGs of the sample. The internal Newtonian acceleration g_{N} is obtained as in Sect. 3 from the stellar distribution, assuming a Sérsic sphere, a uniform M/L (assumed to be 1.3 for DF44 and DFX1, measured for the other galaxies), and no gas. The external Newtonian acceleration g_{Ne} is derived for fiducial distances from the centre of the Coma cluster using the cluster mass distribution obtained by Sanders (2003) and explained in detail in Sect. 4.1.
Figure 9 compares the observed velocity dispersion measurements of each Coma cluster UDG to the MOND predictions with EFE for different distances from the cluster centre, namely the average distance d_{mean} inferred from the Einasto UDG distribution of van der Burg et al. (2016) (cf. Sect. 4.1), 5 Mpc, and 10 Mpc. The EFE pushes the MOND prediction towards the Newtonian case, significantly altering the agreement with the data compared to the isolated MOND prediction: the prediction with EFE is systematically below the measured velocity dispersion for all 11 galaxies of the sample, from about 1σ below to more than 3σ. In addition to ⟨g_{r}⟩ from Eq. (29), we also considered g_{∥} from Eq. (60) of Famaey & McGaugh (2012), recalled in the current Eq. (19), and g_{⊥} from Eq. (30). While g_{⊥} yields similar profiles as ⟨g_{r}⟩, g_{∥} pushes the velocity dispersion profile towards the Newtonian case even further. This result confirms the tension already noted by Bílek et al. (2019b) and Haghi et al. (2019a) for DF44 between observations and the EFE, expanding the sample to 11 Coma cluster UDGs. In particular, Fig. 6 has demonstrated that the projected distribution of the considered UDGs would be expected to be very different if they were all further than 10 Mpc from the cluster centre. As in Sect. 3.2.3, we note that taking into account the rotation hints for four UDGs from Chilingarian et al. (2019) would typically only cause differences at the level of a km s^{−1} in the lineofsight dispersions. This would not change our conclusions, and it would actually only go into the direction of (slightly) increasing the tension between the data and the MOND models with EFE. We now discuss the possible implications of this tension.
Fig. 9. Comparison between the measured lineofsight velocity dispersion of the sample UDGs (grey error bars and/or shaded area) and the MOND prediction at different distances from the cluster centre. Black lines correspond to a uniform anisotropy parameter β = 0 (isotropic), decreasing magenta lines to β = +0.5 (radially biased). The plain and dotted light grey (red) lines recall the predicted isotropic (radiallybiased) velocity dispersion profiles in the isolated MOND and Newtonian cases, respectively (cf. Fig. 4). Data points stem from van Dokkum et al. (2019b) for DF44, van Dokkum et al. (2017) for DFX1, and from Chilingarian et al. (2019) for the other UDGs. The shaded area corresponds to the effective velocity dispersion or, when not available (in the case of DFX1), the central stellar velocity dispersion. 
Fig. 9. continued. 
5. Discussion
Assuming sphericity, uniform M/L ratios, a mass model for the Coma cluster, and an accurate analytical QUMOND formula for the EFE, we find through Jeans modelling that the observed velocity dispersion profiles of all 11 Coma cluster UDGs with velocity dispersion measurements are described well by MOND in isolation but are at odds with the MOND prediction with EFE. We hence do not see any evidence for a violation of the strong equivalence principle in Coma cluster UDGs, contrarily to, for instance, Chae et al. (2020, 2021), for disc galaxies in the field. Our work extends that of Bílek et al. (2019b) and Haghi et al. (2019a), which is limited to DF44 and makes the result all the more compelling. We recall that the MOND predictions do not involve any free parameter.
We hereafter state different possible interpretations for the tension between the measurements and the MOND prediction with EFE, constraining either MOND itself or the formation and evolution of UDGs within this theory: (1) The observed UDGs are further away from the cluster centre than they seem, have fallen inside the cluster relatively recently, and/or are disrupted by tides in the cluster environment. (2) They have higher stellar masstolight ratios than assumed here. (3) They are surrounded by additional baryonic dark matter halos. (4) The EFE varies from one galaxy to another depending on its individual history. (5) The characteristic acceleration scale of MOND varies with the environment; it is higher in clusters. (6) The cluster environment shuts down the EFE within the parent relativistic theory of MOND. Alternatively, MOND being an effective dark matter scaling relation of course also remains a serious possibility; in that context, the fact that cluster UDGs obey the same scaling relation as field spirals, despite their very different environments and likely different formation scenarios, is still particularly intriguing, irrespectively of the underlying theoretical framework.
5.1. Survivor bias
5.1.1. Further away from the cluster centre
The fact that the predictions within MOND in isolation are compatible with the observed stellar velocity dispersion measurements for all 11 Coma cluster UDGs with existing such measurements (Fig. 4), without any free parameter, while the MOND predictions with EFE are far away from the data (Fig. 9), might at first suggest that the UDGs may be further away from the cluster centre than they seem from their projected distance. If they were sufficiently close to the cluster centre for a sufficiently long time and at equilibrium, we would expect their velocity dispersions to follow the predictions with EFE, which are close to the Newtonian values and much lower than the measured ones. As shown in Fig. 9, it would require distances d ≳ 10 Mpc for the EFE to be negligible and the predicted velocity dispersions to approach the high measured values.
As already noted by Bílek et al. (2019b) for DF44, the wide field map of UDG candidates around the Coma cluster by Zaritsky et al. (2019) is compatible with UDGs occurring up to at least 15 Mpc from the cluster centre, with overdensities inside the ∼3 Mpc virial radius and along largescale structure features such as cosmic filaments (cf. their Fig. 11). We note that van Dokkum et al. (2015a) did not detect any UDGs within a projected distance of 300 kpc from the centre of Coma, which was used to put an approximate constraint on the binding mass of UDGs, but the number of objects towards the centre, the extended light profiles of the brighter galaxies, and the intracluster light may have prevented such detections. In particular, van der Burg et al. (2016) did find UDGs within 300 kpc of the centre.
It was also argued that the difference in lineofsight velocity observed by van Dokkum et al. (2015b) between DF44 and the cluster centre (Δv ≈ −800 km s^{−1}) may indicate that this UDG is in the background of the cluster and hence significantly more distant from the centre than its 1.8 Mpc projected distance. This assumed the UDG to be on a relatively radial orbit towards the cluster centre, which would notably happen if it were falling along a cosmic filament. However, we note that the Newtonian rotation velocity at the projected distance of DF44 is about 1100 km s^{−1} given the mass of the Coma cluster (5.5 × 10^{14} M_{⊙}), such that the velocity difference does not necessarily imply a large distance from the centre.
Figure 6 compared the observed projected distance distribution of Coma cluster UDGs from the Yagi et al. (2016) catalogue (black histogram) with the expectation from a distribution with no UDGs within 10 Mpc and a UDG density decreasing as 1/r^{3} beyond (dotted blue). The observed UDG distribution is obviously very different, since it decreases beyond ∼1 Mpc and is relatively well fit by an Einasto profile such as that obtained by van der Burg et al. (2016). Actually, the figure indicates that the Yagi et al. (2016) UDG sample is even more centrally concentrated than the van der Burg et al. (2016) Einasto profile. Unless the considered UDGs were not chosen randomly, we therefore conclude that the centrally depleted UDG distribution required to neglect the EFE is unlikely.
5.1.2. Recent infall onto the cluster
As can be seen in Figs. 4 and 9, a radially biased anisotropy β = +0.5 often enables us to increase the velocity dispersion towards the centre and to make it closer to the observed measurements. Two Chilingarian et al. (2019) galaxies (J130026.26+272735.2 and J130038.63+272835.3) further have decreasing velocity dispersion profiles consistent with a radially biased anisotropy. A radiallybiased anisotropy indeed leads to a decreasing velocity dispersion profile, while a tangentially biased anisotropy leads to an increasing velocity dispersion profile near the centre. DF44 has a generally rising velocity dispersion profile better fit by a tangentially biased anisotropy, as discussed in Bílek et al. (2019b) and Haghi et al. (2019b); however, van Dokkum et al. (2019b) reported a nonGaussian line profile for DF44 with a positive h_{4} GaussHermite moment, which instead favours a radial rather than tangential anisotropy. This positive h_{4} moment indeed reflects the presence of wings in the velocity distribution, which may result from radial orbits (e.g. van der Marel & Franx 1993). Such wings may, however, also stem from deviations from spherical symmetry.
Radial orbits would be expected if the UDGs are on radial orbits themselves, which could provide an explanation for the absence of EFE. While the MOND EFE may contribute to puffup UDGs as they fall in the cluster environment, it also means that their internal dynamical time increases (Milgrom 2015), such that their observed velocity dispersion at a given time may not reflect the EFE at this time but at an earlier, larger distance from the cluster centre. The facts that UDGs may be flowing along cosmic filaments (Zaritsky et al. 2019) and that several UDGs around DF44 share similar radial velocities (van Dokkum et al. 2019b) are compatible with them falling into the Coma cluster on radial orbits for the first time. The cluster orbital time at 1 Mpc of its centre is ∼3 Gyr within MOND, while the orbital times at halfmass radius within the UDGs lie between 0.4 and 1.0 Gyr in the deep MOND regime (between 1.3 and 3.3 Gyr in the Newtonian, strong EFE regime) given their deprojected Sérsic mass profiles. If the difference may result in a higher velocity dispersion than expected from the EFE at the location of the UDG, it may not be sufficient to explain the fact that measured velocity dispersions are close to the isolated MOND predictions for all sample UDGs. This will need to be confirmed by detailed simulations. In any case, failing to explain their observed velocity dispersions in this way would not mean that observed UDGs are not significantly biased towards objects that recently fell inside the cluster.
5.1.3. Disrupted by tides
Both a centrally depleted UDG distribution and observed UDGs that only recently fell onto the cluster could be explained by enhanced tides in MOND and rapid tidal disruption: observations would be biased towards notyet disrupted objects, as in the survivor bias. The MOND EFE is indeed expected to decrease the gravitational support of satellite galaxies and hence make them more prone to tidal effects and tidal disruption (e.g. Brada & Milgrom 2000). If the velocity dispersion of the UDG was set before entering the cluster environment in a region where the EFE was negligible, the UDG may become unbound and dissolve when the gravitational support decreases because of the EFE. Bílek et al. (2019b) accordingly interpreted the observed rising velocity dispersion profile of DF44 as being due to an extra acceleration by tidal forces, which could also explain its relative elongation (with an axis ratio b/a = 0.69) observed by van Dokkum et al. (2019b). van der Burg et al. (2016) further reported a significant central deficit of UDGs compared to compact dwarfs in a sample of eight nearby clusters, possibly because of strong dynamical interactions between galaxies in the central cluster region and tidal destruction due to the cluster potential. Signs of tidal interactions have been directly observed for some UDGs; for example, Mihos et al. (2015), Wittmann et al. (2017), and Bennet et al. (2018) found UDGs associated with tidal streams, Merritt et al. (2016), Toloba et al. (2016), Venhola et al. (2017), and Lim et al. (2020) reported UDGs with elongated or distorted optical morphologies, Scott et al. (2021) observed atomic gas morphology and kinematics of a UDG that are consistent with a recent tidal interaction, and Montes et al. (2020) claimed that the small velocity dispersion of the peculiar UDG NGC 1052DF4 results from tidal disruption. For the current UDG sample, the significant elongations reported by Chilingarian et al. (2019) in many of the galaxies (b/a being respectively 0.69, 0.62, 0.65, 0.56, 0.50, 0.44, and 0.42 for DF44, DFX1, J125848.94+281037.1, J125904.06+281422.4, J125929.89+274303.0, J130026.26+272735.2, and J130038.63+272835.3) and the Sshaped image fitting residuals they obtain for J130028.34+274820.5 may indicate tidal interactions.
To validate the possibility that these Coma cluster UDGs are subjected to tides, we estimate their tidal susceptibility:
which is the ratio between their halfmass radius r_{1/2} ≈ (4/3)R_{e} and the Roche lobe radius r_{2}. The Roche lobe radius r_{2} is derived from the radius r_{1} of the inner Lagrange point using Zhao (2005) and Zhao & Tian (2006), r_{1} being obtained by numerically solving the equation describing the equilibrium between the UDG internal gravitational force, taking the EFE and the external tidal force into account, as explained in Appendix B. Figure 10 shows η as a function of the distance R from the centre of the cluster. As highlighted by the plain red data points, the tidal susceptibility at the average distance d_{mean} (cf. Sect. 4.1) is in the 0.5−1 range for all UDGs, so we expect them to be at least partially affected by tides. We further note that the UDGs are most likely not observed at their pericentre, where the tidal susceptibility would be higher, and that the tidal susceptibility at the projected distance is always above one. The figure also includes the tidal susceptibility of the Coma cluster UDGs without EFE, namely using MOND in isolation. In this case, the tidal susceptibility at the average distance (and even at the projected distance) is always below 0.5. From Fig. 10, we conclude that UDGs should be affected by tides beyond their halflight radius when taking the MOND EFE into account, which would bias observations towards survivors and/or objects that recently fell onto the cluster and could even deplete the central region of the cluster (although this latter point does not seem to be supported by the observations: cf. Sect. 5.1.1). Since MOND in isolation provides an approximate estimate of the measured velocity dispersion (Sect. 3), the internal gravitational accelerations of the sample UDGs may be specifically high enough to escape tides, but this requires us to ignore the EFE. Alternatively, they may be outofequilibrium objects whose stars have been heated by tides, but this would beg the question of why their inflated velocity dispersions would precisely match up with the isolated MOND prediction, which would then appear as a coincidence. Again, future detailed simulations could help to answer these questions more precisely.
Fig. 10. Tidal susceptibility η as a function of the distance R from the centre of the Coma cluster for the sample UDGs. The distance is assumed to be larger than the projected distance d, i.e., R > d. Points correspond to the average distance d_{mean} assuming the van der Burg et al. (2016) Einasto UDG distribution (cf. Sect. 4.1): plain magenta for MOND with EFE, open black for MOND in isolation. 
5.2. Higher stellar masses
Since UDGs are rather unusual objects compared to other galaxy classes, that is, with Sérsic indices close to one but no disc morphology and pressure support, standard procedures to derive their stellar M/L ratios (such as that based on the Pégase.HR evolutionary synthesis package by Le Borgne et al. 2004 used by Chilingarian et al. 2019) may not directly apply. In which case, higher stellar M/L ratios than those used here may alleviate the discrepancy between the predicted and measured velocity dispersions. In Fig. 11, we show the bestfit stellar M/L ratio that enabled us to recover the effective stellar velocity dispersion σ_{eff}, assuming MOND with EFE at the average distance d_{mean}, both in the isotropic case (β = 0) and in the radially biased case (β = +0.5). The average bestfit M/L ratio is 11.5 with β = 0, and 7.0 with β = +0.5. Such values could have been expected from Fig. 9, where the predicted velocity dispersion profile with EFE generally falls short of the observed measurements by a factor close to . We note that there is a small increasing trend of the bestfit M/L ratio towards the cluster centre, compensating for the enhanced EFE. A combination of higher M/L ratio values than fiducial and tidal heating (cf. Sect. 5.1.3) could be invoked to explain the high velocity dispersions. In Appendix C, we explicitly display the radially anisotropic predictions with a stellar masstolight ratio M/L = 4 at the average distance for each galaxy, together with the bestfit M/L curves. If stellar M/L ratios up to 4 may be allowed (e.g. Famaey et al. 2018), the particularly high bestfit values may call for a significant amount of baryonic dark matter.
Fig. 11. Bestfit M/L ratio matching the effective velocity dispersion σ_{eff}, assuming MOND with EFE at the average distance d_{mean}. Crosses correspond to the fiducial M/L ratio values from Table 1, plain black circles to the bestfit values with an anisotropy parameter β = 0, and open magenta circles to the bestfit values with β = 0.5. Plain lines correspond to linear fits, dashed and dotted lines to the averages. The average bestfit M/L ratio is 11.5 with β = 0 and 7.0 with β = +0.5. The fits with β = +0.5 are displayed in Fig. C.1. 
5.3. Baryonic dark matter
We recall that although MOND does reduce the mass discrepancy in clusters compared to Newtonian dynamics, clusters within MOND still need to contain an undetected component whose mass is comparable to the observed baryonic mass (e.g. Sanders 1999, 2003; Milgrom 2008, 2015; Angus et al. 2008) and is rather concentrated in the central parts (Angus et al. 2008). If this undetected component were composed of baryons (e.g. molecular hydrogen), this cluster baryonic dark matter would only account for a small fraction (∼15%) of the unobserved missing baryons. Moreover, if it is in collisionless form (such as compact molecular hydrogen clouds), it may notably help explain the ‘bullet cluster’ (Clowe et al. 2006), where the peak of the mass distribution deduced from weak lensing lies away from the observed baryons. Milgrom (2015) argued that the existence of many cluster UDGs may indeed indicate the presence of baryonic dark matter to stabilise them against tidal disruption. Bílek et al. (2019b) already interpreted the rising velocity dispersion of DF44 and its large number of globular clusters as possible clues indicating the presence of a surrounding baryonic dark halo. However, they noted that (i) if the UDG were far away from the cluster centre and experiencing little or no EFE, the baryonic halo surrounding it (if any) should be sufficiently extended to have no strong gravitational influence towards the centre of the galaxy, and (ii) if the UDG were, on the contrary, experiencing a substantial EFE, the baryonic halo should conspire to compensate for the decrease from the isolated MOND prediction. The two options a priori make the hypothesis of a baryonic dark halo either unnecessary or unappealing, especially as they would have to be generalised to all 11 UDGs of the current sample. However, with the (strong) assumption of the baryonic dark matter following the light in all 11 UDGs, a constant dynamicaltostellar mass ratio of about 10 could be enough to explain the observations, as can also be inferred from Fig. 11. The option of cluster baryonic dark matter is obviously difficult to formally exclude, but would nevertheless imply a conspiracy to make it match the isolated MOND prediction for all the galaxies. From a heuristic point of view, this would incidentally decrease the appeal of modified gravity altogether to explain the observed phenomenology in galaxies.
5.4. Modified inertia
While current MOND theories are of the modified gravity type, only involving a modification of the action of the gravitational field, Milgrom (1994, 2002, 2011) has also considered theories of the modified inertia type, where it is the kinetic action (the coupling to matter for relativistic theories) that is modified. For example in the nonrelativistic case, this would leave the Poisson equation intact but modify the equation of motion of a test particle into A(r, t, a_{0}) = − ∇Φ(r(t)) where A is a function of the trajectory of the particle and Φ is the potential. In such theories, the MOND effects would depend on the whole trajectory of a system and not just on its instantaneous state. As explicitly shown in Milgrom (2011), it is actually possible to construct such a theory without EFE, or with an EFE which differs significantly from that of modified gravity theories, in which case the isolated MOND behaviour of UDGs might reflect their different history of formation and trajectories (cf. for instance the field spirals of Chae et al. 2020). However, such a framework may slightly hinder the predictive power of MOND if the EFE varies from one galaxy to another depending on its individual history. If the EFE were absent altogether, this would also be at odds with the observational evidence and actual need for the EFE in certain places; for instance, in wide binaries (Pittordis & Sutherland 2019; Banik 2019, although see also Hernandez et al. 2022), Andromeda dwarf galaxies (McGaugh & Milgrom 2013a,b), the crater ultrafaint dwarf (McGaugh et al. 2016; Caldwell et al. 2017), or the UDGs NGC 1052DF2 and NGC 1052DF4 (Famaey et al. 2018; Kroupa et al. 2018).
5.5. EMOND
Rather than assuming an additional dark matter component to explain the subsiding mass discrepancy in clusters within MOND, it has been suggested that MOND phenomenology may be the limiting case of a more general relation depending, for example, on the depth of the gravitational potential (or of some scalar field in relativistic versions of the theory). In particular, Zhao & Famaey (2012) proposed an ‘extended MOND’ (EMOND) framework where the critical acceleration a_{0} is replaced by an acceleration A_{0}(Φ) increasing with the depth of the potential well to boost the MOND effect in galaxy clusters. This acceleration would vary from A_{0} ≈ 1 × 10^{−10} m s^{−2} inside galaxies and the Solar System to ∼8 × 10^{−10} m s^{−2} in galaxy clusters (cf. Fig. 1 of Zhao & Famaey 2012), but it could reach values as high as ∼100 × 10^{−10} m s^{−2} in their central regions (Hodson & Zhao 2017a,b). A high critical acceleration in clusters means that the gravitational acceleration would deviate from the Newtonian one at higher accelerations, thus increasing the velocity dispersion in cluster UDGs compared to the MOND predictions with EFE (Sect. 4), even when the external field of the cluster is strong. Although using an unresolved and indirect estimate of the velocity dispersion of Coma cluster UDGs (derived from their effective radius and surface brightness using the fundamental manifold of Zaritsky et al. 2006a,b, 2008), Hodson & Zhao (2017b) already studied Coma cluster UDGs in the EMOND framework. By comparing the baryonic mass needed in EMOND to the stellar mass inferred from photometry, they conclude that this framework is consistent with the observations.
To emulate EMOND, which should not harbour any missing mass in galaxy clusters by construction, we hereafter determine the Newtonian external field from the dominant intracluster hot gas mass distribution,
with ρ_{C} given by the βmodel of Eq. (16), and replace the critical acceleration a_{0} by a variable A_{0}. As shown in Fig. 12, the bestfit A_{0} matching the effective velocity dispersion σ_{eff} when including the EFE and assuming the UDGs to be at their average distance d_{mean} is of the order of 10 × a_{0}, as expected by Zhao & Famaey (2012) in clusters, and decreases with d_{mean} in the isotropic β = 0 case, which follows the increasing need to compensate for the EFE closer to the cluster centre and reach the high measured velocity dispersions.
Fig. 12. Bestfit acceleration A_{0} matching the effective velocity dispersion when emulating EMOND at the average distance d_{mean}, with the external field derived from the dominant hot gas distribution. Plain black circles correspond to an anisotropy parameter β = 0 and the open magenta circles to β = 0.5. Plain lines correspond to linear fits and dashed lines to the averages. A cluster acceleration A_{0} = 10 × a_{0} provides reasonably good fits to the velocity dispersion measurements, especially when some radial anisotropy is allowed (cf. also Fig. D.1). 
We note that a slightly higher cluster mass distribution – to take into account galaxies in addition to the intracluster hot gas – would increase the EFE, decrease the predicted velocity dispersion, and thus lead to slightly higher bestfit accelerations. In Fig. D.1, we compare the measured velocity dispersions to the predicted EMOND profiles with A_{0} = 10 × a_{0}, which are obtained without any free parameter. We find that the predictions fit particularly well with the data, especially since we adopted a single fixed value for A_{0}, and notably for DF44 (contrarily to the indirect inference of Hodson & Zhao 2017b). Incidentally, the above case would also roughly correspond (modulo the weaker Newtonian external field from the cluster in the EMOND case) to the case of increasing the baryonic mass with a massfollowslight profile and with a dynamicaltostellar mass ratio of about 10, as previously discussed in Sect. 5.2. We also note that galaxy J130005.40+285333.0, for which the EMOND agreement is less convincing, is the closest to the cluster centre, namely that exposed to the strongest tidal forces and residing the deepest in the cluster potential well.
However, the more general problem with EMOND is to account for the displaced peaks of the lensing convergence maps of colliding clusters such as the bullet cluster. While peaks of convergence displaced from the main baryonic peaks are easy to produce, both in MOND and EMOND, producing them outside the baryonic peaks – and not in between them – has never been demonstrated to be feasible.
5.6. Screening the EFE in galaxy clusters
We note that the apparent absence of EFE happens precisely within galaxy clusters, where classical MOND fails to explain the overall dynamics of the cluster, and we conjecture here that these two facts might possibly be related. In the case of EMOND, this would be explained by an effective increase of the MOND acceleration constant, not by a screening of the EFE itself. But another possibility is that the EFE is severely damped in galaxy clusters.
In a theory such as that of Skordis & Złośnik (2021), the action harbours a freefunction, playing the role of the MOND interpolating function, depending both on the spatial gradient squared of the scalar field ∇φ^{2} (with a 3/2 exponent, characteristic of MOND actions) and on its temporal derivative having a nonzero minimum leading to gravitating ‘dust’. It is this timedependent term that allows to reproduce a reasonable angular power spectrum for the CMB, and one could therefore speculate that it may also give rise to additional gravitating ‘dust’ inside galaxy clusters, to explain the residual missing mass of MOND. However, it is not clear that, if the scalar field is dominated by this dust component inside the cluster itself, it would couple to the scalar field within the UDG in the same way as in the fully quasistatic limit. Therefore, one could imagine that, precisely because the residual missing mass in galaxy clusters would be caused by the same scalar field as that creating the MOND effect inside the UDG, the EFE could be effectively screened within clusters. We note that this is especially relevant for any model that would try to explain the residual MOND missing mass in clusters of galaxies, as such an explanation would not work if the residual missing mass is made of additional hot DM such as light sterile neutrinos.
In this context of EFE screening, one could imagine two possibilities: one where the EFE would be solely produced by the baryonic mass of the cluster, and one where it would be almost fully screened, with the UDG living in a MOND bubble, effectively decorrelated from the dynamics of the cluster itself. We tested the first hypothesis by redoing the analysis in Sect. 4 using only the Coma cluster hot gas mass distribution M_{gas} (Eq. (33)) derived from the βmodel of Eq. (16) as a source of EFE, instead of the mass distribution M_{C}(r) inferred from hydrostatic equilibrium. This mass is about 1 dex below M_{C} at a distance of 1 Mpc but reaches M_{C} at 10 Mpc. As a consequence, the resulting velocity dispersions at distances smaller than 10 Mpc are higher than in Fig. 9, as shown in Fig. E.1. However, the difference is not sufficient to significantly alter our conclusion on the mismatch between the observed velocity dispersions and the predictions with EFE at the average d_{mean}.
This means that, to explain our results with the nominal values of the stellar M/L ratios, the EFE should be almost fully screened for the UDGs residing inside clusters. This is actually also the case in some hybrid versions of MOND such as the superfluid DM theory (Khoury 2015; Berezhiani & Khoury 2015; Berezhiani et al. 2018, 2019). As discussed in detail in Sect. IX.B of Berezhiani et al. (2018), the superfluid core would be rather small in galaxy clusters (of the order of a few hundred kpc at most), and no EFE would be expected for cluster UDGs, which is contrary to the case of satellite galaxies orbiting within the superfluid core of their host, where the EFE would be expected to be similar to the MOND case.
6. Conclusion
Given their singularly low internal accelerations, ultradiffuse galaxies (UDGs) in clusters provide a testing ground for modified Newtonian dynamics (MOND) and its external field effect (EFE) breaking the strong equivalence principle (Sect. 1). In this article, we compared observed velocity dispersion measurements in 11 UDGs of the Coma cluster (Sect. 2) with MOND predictions with and without EFE, hence continuing the work initiated by Bílek et al. (2019b) and Haghi et al. (2019b) on DF44. Our modelling relied on deprojected Sérsic spherical mass models for the UDGs (Sect. 3.1.2), uniform stellar masstolight ratios, a mass model stemming from Xray observations for the Coma cluster (Sect. 4.1), and quasilinear MOND (QUMOND, Milgrom 2010) for the EFE. Our work can be summarised as follows:

We first showed that the Coma cluster UDGs of the sample lie roughly on the same radial acceleration relation (RAR) as spirals in the field, giving a first indication that their dynamical and baryonic masses are similarly correlated (Sect. 3.1).

We then proceeded to compare the velocity dispersion profiles with Jeans modelling in isolated MOND, showing a remarkably good agreement given that there is no free parameter involved. The agreement is even better when introducing a fixed degree of radial anisotropy (Sect. 3.2).

To take into account the EFE, we proposed a simple approximate analytical formula for the internal acceleration of a spherical system embedded in a constant external field in the context of QUMOND (Sect. 4.2, Eq. (29)), which we validated using a numerical Poisson solver in QUMOND (Appendix A).

We used this formula to compare the velocity dispersion measurements with MOND models including the EFE, which pushes the predicted velocity dispersion towards the Newtonian case and away from the observations at most plausible distances from the centre of the cluster (Sect. 4.3). We noted that all other formulae would generate a stronger EFE and thus increase the tension even more.
These results therefore disfavour the violation of the strong equivalence principle expected in MOND for the UDGs of the Coma cluster, and imply that the EFE has to be severely damped or screened within these systems. We then discussed several options to explain this:

We considered the option that observed UDGs are biased towards objects further away from the centre of the cluster than it seems, but deemed it unlikely given their projected distance distribution (Sect. 5.1.1).

We did not totally exclude the possibility that observed UDGs are falling into the Coma cluster on radial orbits for the first time before being disrupted by tides (Sects. 5.1.2 and 5.1.3), but we noted that the difference in dynamical times within the cluster at the location of the UDGs and within the UDGs might not be sufficient to explain the fact that measured velocity dispersions are all consistent with the fully isolated MOND prediction.

We discussed the option of higher stellar masstolight M/L ratios than measured with standard procedures and/or baryonic dark matter located within the UDGs to increase their velocity dispersion and stabilise them against tidal disruption, but noted that it would imply a conspiracy to make it match the isolated MOND prediction for all the sampled UDGs (Sects. 5.2 and 5.3). A relatively high M/L ≃ 4 ratio for all galaxies combined with some tidal heating was nevertheless not excluded.

We also briefly discussed modified inertia, where the MOND effects would depend on the whole trajectory of a system (Sect. 5.4), and in which the EFE could even effectively be absent from most situations. We noted that the latter case of killing the EFE altogether would be in contradiction with the constraints from, for example, wide binaries.

We showed that the extended MOND (EMOND) framework where the critical acceleration a_{0} is replaced by an acceleration A_{0}(Φ) depending on the depth of the gravitational potential does provide a good match to the velocity dispersion measurements with A_{0} ∼ 10 × a_{0} (Sect. 5.5). We noted, however, that EMOND still struggles to explain the apparently collisionless nature of the residual missing mass in galaxy clusters.

Finally, we considered options where the EFE would be screened in galaxy clusters (Sect. 5.6), which is the case within the superfluid DM theory, and could potentially also be the case in theories where the additional degrees of freedom responsible for the MOND behaviour in galaxies would simultaneously be responsible for an additional ‘dust’ component in galaxy clusters (Skordis & Złośnik 2021).
Of course, we note that MOND simply being an effective dark matter scaling relation also remains a serious possibility; however, even in this case, it is still particularly intriguing that UDGs obey the same scaling relations as field spirals despite their very different environments and formation scenarios, irrespectively of the underlying theoretical framework.
Acknowledgments
We thank the anonymous referee for their uplifting report. The authors would like to thank M. Milgrom, H.S. Zhao, E. Asencio, I. Banik, G. Mamon, R. Errani, F. Jiang, D. Zaritsky, F. Lelli and F. Combes for stimulating discussions and comments. B.F. is grateful for enlightening discussions with C. Skordis and A. Durakovic. J.F., B.F., P.A.O. and R.I. acknowledge funding from the European Research Council (ERC) under the European Union’s Horizon 2020 research and innovation programme (grant agreement No. 834148, GreatDigInTheSky). B.F. and R.I. also acknowledge funding from the Agence Nationale de la Recherche (ANR project ANR18CE310006 and ANR19CE310017). M.B. acknowledges support from the Polish National Science Centre under the grant 2017/26/D/ST9/00449 and is thankful for the financial support by Cercle Gutenberg. O.M. is grateful to the Swiss National Science Foundation for financial support.
References
 Abraham, R. G., & van Dokkum, P. G. 2014, PASP, 126, 55 [Google Scholar]
 Alabi, A., FerréMateu, A., Romanowsky, A. J., et al. 2018, MNRAS, 479, 3308 [Google Scholar]
 Alabi, A. B., Romanowsky, A. J., Forbes, D. A., Brodie, J. P., & Okabe, N. 2020, MNRAS, 496, 3182 [NASA ADS] [CrossRef] [Google Scholar]
 Amorisco, N. C., & Loeb, A. 2016, MNRAS, 459, L51 [NASA ADS] [CrossRef] [Google Scholar]
 Angus, G. W., Shan, H. Y., Zhao, H. S., & Famaey, B. 2007, ApJ, 654, L13 [CrossRef] [Google Scholar]
 Angus, G. W., Famaey, B., & Buote, D. A. 2008, MNRAS, 387, 1470 [CrossRef] [Google Scholar]
 Angus, G. W., Famaey, B., & Diaferio, A. 2010, MNRAS, 402, 395 [NASA ADS] [CrossRef] [Google Scholar]
 Banik, I. 2019, MNRAS, 487, 5291 [NASA ADS] [CrossRef] [Google Scholar]
 Banik, I., & Zhao, H. 2018, ScieFed J. Astrophys., 1, 1000008 [Google Scholar]
 Banik, I., Thies, I., Famaey, B., et al. 2020, ApJ, 905, 135 [Google Scholar]
 Beasley, M. A., & Trujillo, I. 2016, ApJ, 830, 23 [CrossRef] [Google Scholar]
 Beasley, M. A., Romanowsky, A. J., Pota, V., et al. 2016, ApJ, 819, L20 [NASA ADS] [CrossRef] [Google Scholar]
 Begeman, K. G., Broeils, A. H., & Sanders, R. H. 1991, MNRAS, 249, 523 [Google Scholar]
 Bekenstein, J. D. 2004, Phys. Rev. D, 70, 083509 [Google Scholar]
 Bekenstein, J., & Milgrom, M. 1984, ApJ, 286, 7 [NASA ADS] [CrossRef] [Google Scholar]
 Bennet, P., Sand, D. J., Zaritsky, D., et al. 2018, ApJ, 866, L11 [NASA ADS] [CrossRef] [Google Scholar]
 Berezhiani, L., & Khoury, J. 2015, Phys. Rev. D, 92, 103510 [Google Scholar]
 Berezhiani, L., Famaey, B., & Khoury, J. 2018, JCAP, 2018, 021 [CrossRef] [Google Scholar]
 Berezhiani, L., Elder, B., & Khoury, J. 2019, JCAP, 2019, 074 [Google Scholar]
 Bernard, L., & Blanchet, L. 2015, Phys. Rev. D, 91, 103536 [NASA ADS] [CrossRef] [Google Scholar]
 Bílek, M., Samurović, S., & Renaud, F. 2019a, A&A, 625, A32 [Google Scholar]
 Bílek, M., Müller, O., & Famaey, B. 2019b, A&A, 627, L1 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Binney, J., & Mamon, G. A. 1982, MNRAS, 200, 361 [Google Scholar]
 Binney, J., & Tremaine, S. 2008, Galactic Dynamics: Second Edition (Princeton: Princeton University Press) [Google Scholar]
 Blanchet, L., & Heisenberg, L. 2015, JCAP, 2015, 026 [CrossRef] [Google Scholar]
 Blanchet, L., & Le Tiec, A. 2009, Phys. Rev. D, 80, 023524 [CrossRef] [Google Scholar]
 Bosma, A. 1978, PhD Thesis, University of Groningen [Google Scholar]
 Brada, R., & Milgrom, M. 2000, ApJ, 541, 556 [NASA ADS] [CrossRef] [Google Scholar]
 Caldwell, N., Walker, M. G., Mateo, M., et al. 2017, ApJ, 839, 20 [NASA ADS] [CrossRef] [Google Scholar]
 Carleton, T., Errani, R., Cooper, M., et al. 2019, MNRAS, 485, 382 [NASA ADS] [CrossRef] [Google Scholar]
 Carter, D., Goudfrooij, P., Mobasher, B., et al. 2008, ApJS, 176, 424 [Google Scholar]
 Chae, K.H., Lelli, F., Desmond, H., et al. 2020, ApJ, 904, 51 [NASA ADS] [CrossRef] [Google Scholar]
 Chae, K. H., Desmond, H., Lelli, F., McGaugh, S. S., & Schombert, J. M. 2021, ApJ, 921, 104 [NASA ADS] [CrossRef] [Google Scholar]
 Chan, T. K., Kereš, D., Wetzel, A., et al. 2018, MNRAS, 478, 906 [NASA ADS] [CrossRef] [Google Scholar]
 Chilingarian, I. V., Afanasiev, A. V., Grishin, K. A., Fabricant, D., & Moran, S. 2019, ApJ, 884, 79 [NASA ADS] [CrossRef] [Google Scholar]
 Clifton, T., Ferreira, P. G., Padilla, A., & Skordis, C. 2012, Phys. Rep., 513, 1 [NASA ADS] [CrossRef] [Google Scholar]
 Clowe, D., Bradač, M., Gonzalez, A. H., et al. 2006, ApJ, 648, L109 [NASA ADS] [CrossRef] [Google Scholar]
 Colless, M., & Dunn, A. M. 1996, ApJ, 458, 435 [Google Scholar]
 Combes, F. 2014, A&A, 571, A82 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Danieli, S., van Dokkum, P., Conroy, C., Abraham, R., & Romanowsky, A. J. 2019, ApJ, 874, L12 [NASA ADS] [CrossRef] [Google Scholar]
 de Blok, W. J. G., & McGaugh, S. S. 1998, ApJ, 508, 132 [NASA ADS] [CrossRef] [Google Scholar]
 Di Cintio, A., Brook, C. B., Dutton, A. A., et al. 2017, MNRAS, 466, L1 [NASA ADS] [CrossRef] [Google Scholar]
 Dutta Chowdhury, D., van den Bosch, F. C., & van Dokkum, P. 2019, ApJ, 877, 133 [NASA ADS] [CrossRef] [Google Scholar]
 Dutta Chowdhury, D., van den Bosch, F. C., & van Dokkum, P. 2020, ApJ, 903, 149 [NASA ADS] [CrossRef] [Google Scholar]
 Einasto, J. 1965, Trudy Astrofizicheskogo Instituta AlmaAta, 5, 87 [NASA ADS] [Google Scholar]
 Emsellem, E., van der Burg, R. F. J., Fensch, J., et al. 2019, A&A, 625, A76 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Faber, S. M., & Gallagher, J. S. 1979, ARA&A, 17, 135 [NASA ADS] [CrossRef] [Google Scholar]
 Fabricant, D., Fata, R., Epps, H., et al. 2019, PASP, 131, 075004 [NASA ADS] [CrossRef] [Google Scholar]
 Famaey, B., & Binney, J. 2005, MNRAS, 363, 603 [NASA ADS] [CrossRef] [Google Scholar]
 Famaey, B., & McGaugh, S. S. 2012, Liv. Rev. Rel., 15, 10 [Google Scholar]
 Famaey, B., McGaugh, S., & Milgrom, M. 2018, MNRAS, 480, 473 [NASA ADS] [CrossRef] [Google Scholar]
 Fensch, J., van der Burg, R. F. J., Jeřábková, T., et al. 2019, A&A, 625, A77 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 FerréMateu, A., Alabi, A., Forbes, D. A., et al. 2018, MNRAS, 479, 4891 [CrossRef] [Google Scholar]
 Forbes, D. A., Alabi, A., Romanowsky, A. J., Brodie, J. P., & Arimoto, N. 2020, MNRAS, 492, 4874 [Google Scholar]
 Forbes, D. A., Gannon, J. S., Romanowsky, A. J., et al. 2021, MNRAS, 500, 1279 [Google Scholar]
 Fosbury, R. A. E., Mebold, U., Goss, W. M., & Dopita, M. A. 1978, MNRAS, 183, 549 [NASA ADS] [CrossRef] [Google Scholar]
 Freundlich, J., Dekel, A., Jiang, F., et al. 2020a, MNRAS, 491, 4523 [NASA ADS] [CrossRef] [Google Scholar]
 Freundlich, J., Jiang, F., Dekel, A., et al. 2020b, MNRAS, 499, 2912 [NASA ADS] [CrossRef] [Google Scholar]
 Gannon, J. S., Forbes, D. A., Romanowsky, A. J., et al. 2020, MNRAS, 495, 2582 [NASA ADS] [CrossRef] [Google Scholar]
 Gannon, J. S., Dullo, B. T., Forbes, D. A., et al. 2021, MNRAS, 502, 3144 [Google Scholar]
 Gentile, G., Salucci, P., Klein, U., & Granato, G. L. 2007a, MNRAS, 375, 199 [CrossRef] [Google Scholar]
 Gentile, G., Famaey, B., Combes, F., et al. 2007b, A&A, 472, L25 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Gentile, G., Famaey, B., & de Blok, W. J. G. 2011, A&A, 527, A76 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Greco, J. P., Greene, J. E., Strauss, M. A., et al. 2018, ApJ, 857, 104 [NASA ADS] [CrossRef] [Google Scholar]
 Habas, R., Marleau, F. R., Duc, P.A., et al. 2020, MNRAS, 491, 1901 [NASA ADS] [Google Scholar]
 Haghi, H., Baumgardt, H., Kroupa, P., et al. 2009, MNRAS, 395, 1549 [NASA ADS] [CrossRef] [Google Scholar]
 Haghi, H., Bazkiaei, A. E., Zonoozi, A. H., & Kroupa, P. 2016, MNRAS, 458, 4172 [NASA ADS] [CrossRef] [Google Scholar]
 Haghi, H., Kroupa, P., Banik, I., et al. 2019a, MNRAS, 487, 2441 [CrossRef] [Google Scholar]
 Haghi, H., Amiri, V., Hasani Zonoozi, A., et al. 2019b, ApJ, 884, L25 [NASA ADS] [CrossRef] [Google Scholar]
 Haslbauer, M., Banik, I., Kroupa, P., & Grishunin, K. 2019, MNRAS, 489, 2634 [NASA ADS] [CrossRef] [Google Scholar]
 Haslbauer, M., Banik, I., & Kroupa, P. 2020, MNRAS, 499, 2845 [Google Scholar]
 Hayashi, K., & Inoue, S. 2018, MNRAS, 481, L59 [NASA ADS] [CrossRef] [Google Scholar]
 Hees, A., Famaey, B., Angus, G. W., & Gentile, G. 2016, MNRAS, 455, 449 [NASA ADS] [CrossRef] [Google Scholar]
 Hernandez, X., Cookson, S., & Cortés, R. A. M. 2022, MNRAS, 509, 2304 [NASA ADS] [Google Scholar]
 Hodson, A. O., & Zhao, H. 2017a, A&A, 598, A127 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Hodson, A. O., & Zhao, H. 2017b, A&A, 607, A109 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Hu, W., Barkana, R., & Gruzinov, A. 2000, Phys. Rev. Lett., 85, 1158 [NASA ADS] [CrossRef] [Google Scholar]
 Impey, C., Bothun, G., & Malin, D. 1988, ApJ, 330, 634 [Google Scholar]
 Janowiecki, S., Leisman, L., Józsa, G., et al. 2015, ApJ, 801, 96 [NASA ADS] [CrossRef] [Google Scholar]
 Janssens, S., Abraham, R., Brodie, J., et al. 2017, ApJ, 839, L17 [NASA ADS] [CrossRef] [Google Scholar]
 Jiang, F., Dekel, A., Freundlich, J., et al. 2019, MNRAS, 487, 5272 [CrossRef] [Google Scholar]
 Kadowaki, J., Zaritsky, D., & Donnerstein, R. L. 2017, ApJ, 838, L21 [Google Scholar]
 Karachentsev, I. D., Karachentseva, V. E., Suchkov, A. A., & Grebel, E. K. 2000, A&AS, 145, 415 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Karunakaran, A., Spekkens, K., Zaritsky, D., et al. 2020, ApJ, 902, 39 [CrossRef] [Google Scholar]
 Khoury, J. 2015, Phys. Rev. D, 91, 024022 [NASA ADS] [CrossRef] [Google Scholar]
 Koda, J., Yagi, M., Yamanoi, H., & Komiyama, Y. 2015, ApJ, 807, L2 [NASA ADS] [CrossRef] [Google Scholar]
 Kretschmer, M., Dekel, A., Freundlich, J., et al. 2021, MNRAS, 503, 5238 [NASA ADS] [CrossRef] [Google Scholar]
 Kroupa, P., Haghi, H., Javanmardi, B., et al. 2018, Nature, 561, E4 [CrossRef] [Google Scholar]
 Kubo, J. M., Stebbins, A., Annis, J., et al. 2007, ApJ, 671, 1466 [NASA ADS] [CrossRef] [Google Scholar]
 Laporte, C. F. P., Agnello, A., & Navarro, J. F. 2019, MNRAS, 484, 245 [NASA ADS] [CrossRef] [Google Scholar]
 Le Borgne, D., RoccaVolmerange, B., Prugniel, P., et al. 2004, A&A, 425, 881 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Leisman, L., Haynes, M. P., Janowiecki, S., et al. 2017, ApJ, 842, 133 [NASA ADS] [CrossRef] [Google Scholar]
 Lelli, F., McGaugh, S. S., & Schombert, J. M. 2016a, AJ, 152, 157 [Google Scholar]
 Lelli, F., McGaugh, S. S., & Schombert, J. M. 2016b, ApJ, 816, L14 [Google Scholar]
 Lelli, F., McGaugh, S. S., Schombert, J. M., & Pawlowski, M. S. 2017, ApJ, 836, 152 [Google Scholar]
 Lelli, F., McGaugh, S. S., Schombert, J. M., Desmond, H., & Katz, H. 2019, MNRAS, 484, 3267 [NASA ADS] [CrossRef] [Google Scholar]
 Li, P., Lelli, F., McGaugh, S., & Schombert, J. 2018, A&A, 615, A3 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Lim, S., Peng, E. W., Côté, P., et al. 2018, ApJ, 862, 82 [Google Scholar]
 Lim, S., Côté, P., Peng, E. W., et al. 2020, ApJ, 899, 69 [CrossRef] [Google Scholar]
 Lima Neto, G. B., Gerbal, D., & Márquez, I. 1999, MNRAS, 309, 481 [Google Scholar]
 Lüghausen, F., Famaey, B., Kroupa, P., et al. 2013, MNRAS, 432, 2846 [CrossRef] [Google Scholar]
 Lüghausen, F., Famaey, B., & Kroupa, P. 2015, Can. J. Phys., 93, 232 [CrossRef] [Google Scholar]
 Mamon, G. A., & Łokas, E. L. 2005, MNRAS, 363, 705 [Google Scholar]
 Márquez, I., Lima Neto, G. B., Capelato, H., Durret, F., & Gerbal, D. 2000, A&A, 353, 873 [NASA ADS] [Google Scholar]
 Martin, N. F., Collins, M. L. M., Longeard, N., & Tollerud, E. 2018, ApJ, 859, L5 [NASA ADS] [CrossRef] [Google Scholar]
 MartínezDelgado, D., Läsker, R., Sharina, M., et al. 2016, AJ, 151, 96 [CrossRef] [Google Scholar]
 McGaugh, S. S. 2005, ApJ, 632, 859 [Google Scholar]
 McGaugh, S. S. 2008, ApJ, 683, 137 [NASA ADS] [CrossRef] [Google Scholar]
 McGaugh, S. S. 2012, AJ, 143, 40 [Google Scholar]
 McGaugh, S. S. 2016, ApJ, 832, L8 [NASA ADS] [CrossRef] [Google Scholar]
 McGaugh, S., & Milgrom, M. 2013a, ApJ, 766, 22 [NASA ADS] [CrossRef] [Google Scholar]
 McGaugh, S., & Milgrom, M. 2013b, ApJ, 775, 139 [NASA ADS] [CrossRef] [Google Scholar]
 McGaugh, S. S., & de Blok, W. J. G. 1998, ApJ, 499, 66 [NASA ADS] [CrossRef] [Google Scholar]
 McGaugh, S. S., Schombert, J. M., Bothun, G. D., & de Blok, W. J. G. 2000, ApJ, 533, L99 [Google Scholar]
 McGaugh, S. S., Lelli, F., & Schombert, J. M. 2016, Phys. Rev. Lett., 117, 201101 [NASA ADS] [CrossRef] [Google Scholar]
 McGaugh, S., Lelli, F., Schombert, J., et al. 2021, AJ, 162, 202 [NASA ADS] [CrossRef] [Google Scholar]
 Merritt, A., van Dokkum, P., Danieli, S., et al. 2016, ApJ, 833, 168 [NASA ADS] [CrossRef] [Google Scholar]
 Mihos, J. C., Durrell, P. R., Ferrarese, L., et al. 2015, ApJ, 809, L21 [Google Scholar]
 Mihos, J. C., Harding, P., Feldmeier, J. J., et al. 2017, ApJ, 834, 16 [Google Scholar]
 Milgrom, M. 1983a, ApJ, 270, 384 [Google Scholar]
 Milgrom, M. 1983b, ApJ, 270, 365 [Google Scholar]
 Milgrom, M. 1983c, ApJ, 270, 371 [Google Scholar]
 Milgrom, M. 1984, ApJ, 287, 571 [CrossRef] [Google Scholar]
 Milgrom, M. 1986, ApJ, 306, 9 [NASA ADS] [CrossRef] [Google Scholar]
 Milgrom, M. 1994, Ann. Phys., 229, 384 [NASA ADS] [CrossRef] [Google Scholar]
 Milgrom, M. 2002, New Astron. Rev., 46, 741 [CrossRef] [Google Scholar]
 Milgrom, M. 2008, New Astron. Rev., 51, 906 [CrossRef] [Google Scholar]
 Milgrom, M. 2010, MNRAS, 403, 886 [NASA ADS] [CrossRef] [Google Scholar]
 Milgrom, M. 2011, Acta Phys. Pol. B, 42, 2175 [Google Scholar]
 Milgrom, M. 2012, Phys. Rev. Lett., 109, 131101 [NASA ADS] [CrossRef] [Google Scholar]
 Milgrom, M. 2014, MNRAS, 437, 2531 [NASA ADS] [CrossRef] [Google Scholar]
 Milgrom, M. 2015, MNRAS, 454, 3810 [NASA ADS] [CrossRef] [Google Scholar]
 Milgrom, M. 2016, ArXiv eprints [arXiv:1609.06642] [Google Scholar]
 Montes, M., InfanteSainz, R., MadrigalAguado, A., et al. 2020, ApJ, 904, 114 [Google Scholar]
 Morrissey, P., Matuszewski, M., Martin, C., et al. 2012, in Groundbased and Airborne Instrumentation for Astronomy IV, eds. I. S. McLean, S. K. Ramsay, & H. Takami, SPIE Conf. Ser., 8446, 844613 [NASA ADS] [Google Scholar]
 Morrissey, P., Matuszewski, M., Martin, D. C., et al. 2018, ApJ, 864, 93 [Google Scholar]
 Muñoz, R. P., Eigenthaler, P., Puzia, T. H., et al. 2015, ApJ, 813, L15 [Google Scholar]
 Müller, O., Jerjen, H., & Binggeli, B. 2018, A&A, 615, A105 [Google Scholar]
 Müller, O., Famaey, B., & Zhao, H. 2019, A&A, 623, A36 [Google Scholar]
 Müller, O., Marleau, F. R., Duc, P.A., et al. 2020, A&A, 640, A106 [Google Scholar]
 Müller, O., Durrell, P. R., Marleau, F. R., et al. 2021, ApJ, 923, 9 [CrossRef] [Google Scholar]
 Natarajan, P., & Zhao, H. 2008, MNRAS, 389, 250 [NASA ADS] [CrossRef] [Google Scholar]
 Nusser, A. 2018, ApJ, 863, L17 [CrossRef] [Google Scholar]
 Nusser, A. 2019, MNRAS, 484, 510 [CrossRef] [Google Scholar]
 Oria, P. A., Famaey, B., Thomas, G. F., et al. 2021, ApJ, 923, 68 [NASA ADS] [CrossRef] [Google Scholar]
 Pandya, V., Romanowsky, A. J., Laine, S., et al. 2018, ApJ, 858, 29 [Google Scholar]
 Papastergis, E., Adams, E. A. K., & Romanowsky, A. J. 2017, A&A, 601, L10 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Pawlowski, M. S., & McGaugh, S. S. 2014, MNRAS, 440, 908 [NASA ADS] [CrossRef] [Google Scholar]
 Pawlowski, M. S., McGaugh, S. S., & Jerjen, H. 2015, MNRAS, 453, 1047 [NASA ADS] [CrossRef] [Google Scholar]
 Peng, E. W., & Lim, S. 2016, ApJ, 822, L31 [NASA ADS] [CrossRef] [Google Scholar]
 Pittordis, C., & Sutherland, W. 2019, MNRAS, 488, 4740 [NASA ADS] [CrossRef] [Google Scholar]
 Pointecouteau, E., & Silk, J. 2005, MNRAS, 364, 654 [NASA ADS] [CrossRef] [Google Scholar]
 Prole, D. J., van der Burg, R. F. J., Hilker, M., & Davies, J. I. 2019, MNRAS, 488, 2143 [NASA ADS] [Google Scholar]
 Reiprich, T. H. 2001, PhD Thesis, MaxPlanckInstitut für Extraterrestrische Physik, Germany [Google Scholar]
 Román, J., & Trujillo, I. 2017, MNRAS, 468, 703 [CrossRef] [Google Scholar]
 Rong, Y., Guo, Q., Gao, L., et al. 2017, MNRAS, 470, 4231 [Google Scholar]
 Roshan, M., Banik, I., Ghafourian, N., et al. 2021, MNRAS, 503, 2833 [NASA ADS] [CrossRef] [Google Scholar]
 Rubin, V. C., Ford, W. K., Jr., & Thonnard, N. 1978, ApJ, 225, L107 [NASA ADS] [CrossRef] [Google Scholar]
 RuizLara, T., Beasley, M. A., FalcónBarroso, J., et al. 2018, MNRAS, 478, 2034 [NASA ADS] [CrossRef] [Google Scholar]
 Sales, L. V., Navarro, J. F., Peñafiel, L., et al. 2020, MNRAS, 494, 1848 [Google Scholar]
 Sandage, A., & Binggeli, B. 1984, AJ, 89, 919 [Google Scholar]
 Sanders, R. H. 1996, ApJ, 473, 117 [NASA ADS] [CrossRef] [Google Scholar]
 Sanders, R. H. 1999, ApJ, 512, L23 [NASA ADS] [CrossRef] [Google Scholar]
 Sanders, R. H. 2003, MNRAS, 342, 901 [NASA ADS] [CrossRef] [Google Scholar]
 Sanders, R. H. 2019, MNRAS, 485, 513 [NASA ADS] [CrossRef] [Google Scholar]
 Sanders, R. H., & McGaugh, S. S. 2002, ARA&A, 40, 263 [NASA ADS] [CrossRef] [Google Scholar]
 Sanders, R. H., & Noordermeer, E. 2007, MNRAS, 379, 702 [NASA ADS] [CrossRef] [Google Scholar]
 Sanders, R. H., & Verheijen, M. A. W. 1998, ApJ, 503, 97 [NASA ADS] [CrossRef] [Google Scholar]
 Scott, T. C., Sengupta, C., Lagos, P., Chung, A., & Wong, O. I. 2021, MNRAS, 503, 3953 [NASA ADS] [CrossRef] [Google Scholar]
 Shi, D. D., Zheng, X. Z., Zhao, H. B., et al. 2017, ApJ, 846, 26 [Google Scholar]
 Skordis, C., & Złośnik, T. 2019, Phys. Rev. D, 100, 104013 [NASA ADS] [CrossRef] [Google Scholar]
 Skordis, C., & Złośnik, T. 2021, Phys. Rev. Lett., 127, 161302 [NASA ADS] [CrossRef] [Google Scholar]
 Swaters, R. A., Sanders, R. H., & McGaugh, S. S. 2010, ApJ, 718, 380 [NASA ADS] [CrossRef] [Google Scholar]
 Thomas, G. F., Famaey, B., Ibata, R., et al. 2018, A&A, 609, A44 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Tiret, O., & Combes, F. 2007, A&A, 464, 517 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Tiret, O., & Combes, F. 2008, A&A, 483, 719 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Toloba, E., Sand, D. J., Spekkens, K., et al. 2016, ApJ, 816, L5 [Google Scholar]
 Toloba, E., Lim, S., Peng, E., et al. 2018, ApJ, 856, L31 [CrossRef] [Google Scholar]
 Trujillo, I., Beasley, M. A., Borlaff, A., et al. 2019, MNRAS, 486, 1192 [Google Scholar]
 van der Burg, R. F. J., Muzzin, A., & Hoekstra, H. 2016, A&A, 590, A20 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 van der Marel, R. P., & Franx, M. 1993, ApJ, 407, 525 [Google Scholar]
 van Dokkum, P. G., Abraham, R., Merritt, A., et al. 2015a, ApJ, 798, L45 [NASA ADS] [CrossRef] [Google Scholar]
 van Dokkum, P. G., Romanowsky, A. J., Abraham, R., et al. 2015b, ApJ, 804, L26 [NASA ADS] [CrossRef] [Google Scholar]
 van Dokkum, P., Abraham, R., Brodie, J., et al. 2016, ApJ, 828, L6 [Google Scholar]
 van Dokkum, P., Abraham, R., Romanowsky, A. J., et al. 2017, ApJ, 844, L11 [NASA ADS] [CrossRef] [Google Scholar]
 van Dokkum, P., Danieli, S., Cohen, Y., et al. 2018, Nature, 555, 629 [Google Scholar]
 van Dokkum, P., Danieli, S., Abraham, R., Conroy, C., & Romanowsky, A. J. 2019a, ApJ, 874, L5 [Google Scholar]
 van Dokkum, P., Wasserman, A., Danieli, S., et al. 2019b, ApJ, 880, 91 [Google Scholar]
 Venhola, A., Peletier, R., Laurikainen, E., et al. 2017, A&A, 608, A142 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Villaume, A., Romanowsky, A. J., Brodie, J., et al. 2021, ApJ, submitted [arXiv:2101.02220] [Google Scholar]
 Wasserman, A., Romanowsky, A. J., Brodie, J., et al. 2018, ApJ, 863, L15 [NASA ADS] [CrossRef] [Google Scholar]
 Wasserman, A., van Dokkum, P., Romanowsky, A. J., et al. 2019, ApJ, 885, 155 [NASA ADS] [CrossRef] [Google Scholar]
 Wittmann, C., Lisker, T., Ambachew Tilahun, L., et al. 2017, MNRAS, 470, 1512 [Google Scholar]
 Wolf, J., Martinez, G. D., Bullock, J. S., et al. 2010, MNRAS, 406, 1220 [NASA ADS] [Google Scholar]
 Wu, X., & Kroupa, P. 2015, MNRAS, 446, 330 [NASA ADS] [CrossRef] [Google Scholar]
 Wu, X., Zhao, H., Famaey, B., et al. 2007, ApJ, 665, L101 [NASA ADS] [CrossRef] [Google Scholar]
 Yagi, M., Koda, J., Komiyama, Y., & Yamanoi, H. 2016, ApJS, 225, 11 [CrossRef] [Google Scholar]
 Yozin, C., & Bekki, K. 2015, MNRAS, 452, 937 [NASA ADS] [CrossRef] [Google Scholar]
 Zaritsky, D., Gonzalez, A. H., & Zabludoff, A. I. 2006a, ApJ, 642, L37 [NASA ADS] [CrossRef] [Google Scholar]
 Zaritsky, D., Gonzalez, A. H., & Zabludoff, A. I. 2006b, ApJ, 638, 725 [Google Scholar]
 Zaritsky, D., Zabludoff, A. I., & Gonzalez, A. H. 2008, in Formation and Evolution of Galaxy Disks, eds. J. G. Funes, & E. M. Corsini, ASP Conf. Ser., 396, 381 [NASA ADS] [Google Scholar]
 Zaritsky, D., Donnerstein, R., Dey, A., et al. 2019, ApJS, 240, 1 [Google Scholar]
 Zhao, H. S. 2005, A&A, 444, L25 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Zhao, H. S., & Famaey, B. 2006, ApJ, 638, L9 [CrossRef] [Google Scholar]
 Zhao, H., & Famaey, B. 2012, Phys. Rev. D, 86, 067301 [NASA ADS] [CrossRef] [Google Scholar]
 Zhao, H., & Tian, L. 2006, A&A, 450, 1005 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
Appendix A: Numerical test of the EFE formula
In order to test the accuracy of Eq. (29), we compared the acceleration it yields to the QUMOND gravitational acceleration obtained by numerical computation on a given test case. For this purpose, we considered a UDG modelled as a Plummer sphere of baryonic mass M_{0} = 3.9 × 10^{8} M_{⊙} and characteristic radius R_{P} = 3.6 kpc (i.e. a halfmass radius of 4.7 kpc) located at distances d = 0.5, 1, 2, 5, 10, and 20 Mpc from a point mass M_{C}(d) representing the centre of the Coma cluster. The mass and radius of the UDG are comparable to those of DF44 (cf. Table 1). The MOND potential Φ was obtained by numerical integration of the QUMOND Poisson equation:
where ρ is the baryonic density and ρ_{PDM} the phantom dark matter (PDM) density (e.g. Milgrom 1986, 2010). The PDM is the theoretical matter that would be needed to obtain MOND effects with Newtonian gravity; its density can be expressed by combining the Newtonian Poisson equation,
with the QUMOND generalised Poisson equation,
yielding
where Φ_{N} is the Newtonian potential. This formula was discretised on a grid pattern using finite differences and computed numerically, as, for example, in Lüghausen et al. (2015) and Oria et al. (2021).
Figure A.1 compares the resulting average QUMOND acceleration profile ⟨g_{r}⟩ with the analytic expression from Eq. (29) for different distances from the point mass modelling the Coma cluster. For comparison, we also considered g_{∥} corresponding to the case where g and g_{e} are aligned (Eq. (60) of Famaey & McGaugh 2012, or Eq. (19) above), and g_{⊥} corresponding to the case where the two fields are perpendicular (Eq. (30) above). The points where g_{N} = a_{0}/100 in each profile (namely, at radii r = 1.1 and 5.2 kpc, cf. the dotted curve in the upper left panel of Fig. A.1) are reported in Fig. 8, which highlights how QUMOND and the different analytic expressions transition from the deepMOND regime (g = ν(g_{N}/a_{0})g_{N}) to the EFEdominated Newtonian regime (g = g_{N}). The different plots confirm that ⟨g_{r}⟩≈g_{⊥} and show that ⟨g_{r}⟩ fits the numerically computed acceleration much better than g_{∥}, which systematically overestimates the EFE and thus underestimates the acceleration field.
Fig. A.1. Comparison between the MOND acceleration with EFE resulting from a QUMOND numerical integration and the analytical expressions ⟨g_{r}⟩ from Eq. (29), g_{∥} from Eq. (60) of Famaey & McGaugh (2012), and g_{⊥} from Eq. (30), considering a Plummer sphere (M_{0} = 3.9 × 10^{8} M_{⊙}, R_{P} = 3.6 kpc) at different distances d from a point mass equal to that of the Coma cluster at that distance given the mass model presented in Sect. 4.1. The plots show that ⟨g_{r}⟩≈g_{⊥} provides a good approximation of the MOND acceleration with EFE. 
Appendix B: Tidal susceptibility
We estimated the tidal susceptibility of a UDG as the ratio η = r_{1/2}/r_{2} between its deprojected halfmass radius r_{1/2} ≈ (4/3)R_{e} and the radius r_{2} of the Roche lobe perpendicularly to the axis linking the UDG and the centre of the cluster. This latter radius r_{2} is the minimum radius associated with the teardropshaped Roche lobe bounded by the gravitational equipotential where the external tidal field from the Coma cluster balances the internal field g of the UDG in the frame of reference of the UDG (e.g. Fig. 1 of Zhao 2005). The apex of the Roche lobe pointing towards the cluster centre is the inner Lagrange point L_{1}, situated at radius r_{1} > r_{2} from the UDG centre.
We determined the radius of the inner Lagrange point r_{1} for a UDG at a distance R from the cluster centre by numerically solving the force equilibrium equation:
where the lefthand term corresponds to the internal gravitational force and the righthand term to the tidal force, involving the spatial derivative of the external field at R. We note that r_{1} is not defined when this derivative is negative, since that would correspond to a compressive rather than disruptive tidal force. The internal acceleration g assumes a spherical mass profile for the UDG (cf. Section 3.1.2) and takes into account the EFE (cf. Section 4.2). The radial derivative of the spherically symmetric external Newtonian acceleration g_{Ne}(R) = − GM_{C}(R)/R^{2} is
and that of the corresponding MOND acceleration g_{e} = ν(g_{Ne}/a_{0})g_{Ne},
For the mass model of the Coma cluster derived by Sanders (2003) and presented in Section 4.1, Eq. (16) for ρ_{C} yields
and
with M_{0} = kT_{C}r_{C}/Gμm_{p} and X = R/r_{C}. The derivatives
and
enable us to calculate
This expression stemming from Eq. (13) can be plugged into Eqs. (B.2) and (B.3) to retrieve dg_{e}/dR, and subsequently to obtain the radius r_{1} of the inner Lagrange point by solving Eq. (B.1) numerically.
Approximating the satellite as a point mass, Zhao (2005) and Zhao & Tian (2006) determined
with
from which we obtained r_{2}. Eq. (B.10) yields Δ_{1} = 2 in the weak acceleration regime of MOND (g_{Ne} ≪ a_{0}), Δ_{1} = 1 in the strong acceleration regime (g_{Ne} ≫ a_{0}), and we obtained, on average, an intermediate Δ ≈ 1.6 for the Coma cluster UDGs at their average distance d_{mean}. The ratio r_{2}/r_{1} is thus equal to 0.471 in the weak acceleration regime, to 0.667 in the strong acceleration regime, and somewhere in between for the Coma cluster UDGs. Although Eq. (B.9) may not be exact for a given satellite mass distribution m(r), it does give the right order of magnitude for the tidal susceptibility and converges to the exact value when r_{1/2} ≪ r_{2} (i.e. η ≪ 1).
We note that Zhao (2005) and Zhao & Tian (2006) also provided an expression for the radius of the inner Lagrange point for a point mass m situated at a distance R from the centre of a spherical mass distribution M(R),
where is the angular frequency for a circular orbit at radius R. Given Eq. (B.9), this yields r_{2} = 0.374R(m/M)^{1/3} in the weak acceleration regime (ζ_{1} = 2, Δ_{1} = 2) and r_{2} = 0.462(m/M)^{1/3} in the strong acceleration regime (ζ_{1} = 3, Δ_{1} = 1). In the present work, we instead used Eq. (B.1) to determine r_{1}, taking into account both the spatial extent of the UDGs and the EFE.
Appendix C: Higher stellar M/L ratio
In Fig. C.1, we show the MOND predicted velocity dispersion with EFE at the average distance d_{mean} inferred from the Einasto UDG distribution of van der Burg et al. (2016) (cf. Section 4.1) for different M/L ratio values, in the case of a radially biased anisotropy β = +0.5. The bestfit values recovering the luminosityweighted effective stellar velocity dispersion σ_{eff} range from 1.5 to 16.5, with an average of 7.0.
Fig. C.1. Comparison between the measured lineofsight velocity dispersion of the sample UDGs (grey error bars and/or shaded area) and the MOND prediction with EFE at the average distance d_{mean} for different values of the stellar masstolight ratio M/L, assuming a radially biased anisotropy β = +0.5. The fiducial M/L ratio from Table 1, the bestfit M/L ratio recovering the effective velocity dispersion σ_{eff}, and the large bust still acceptable M/L = 4 are compared. As in Fig. 9, the plain and dotted light red curves recall the predicted radially biased velocity dispersion profiles in the isolated MOND and Newtonian cases, respectively (cf. Fig. 4). 
Fig. C.1. Continued. 
Appendix D: EMOND profiles
In Fig. D.1, we show the velocity dispersion profiles predicted when emulating EMOND with a critical acceleration A_{0} = 10 × a_{0} at the average distance d_{mean} inferred from the Einasto UDG distribution of van der Burg et al. (2016) (cf. Section 4.1). The external field was determined from the dominant hot gas mass M_{gas} of Eq. (33) and its effect using Eq. (29). EMOND with A_{0} = 10 × a_{0} provides a good match to the velocity dispersion measurements. Galaxy J130005.40+285333.0, for which the agreement is less convincing, is the closest to the cluster centre; hence, it is exposed to the strongest tidal forces and resides the deepest in the potential well of the cluster.
Fig. D.1. Comparison between the measured lineofsight velocity dispersion of the sample UDGs, as in Fig. 9, to the predicted velocity dispersion at the average distance d_{mean} obtained by emulating EMOND with A_{0} = 10 × a_{0}, shown as dashed black (β = 0) and magenta (β = +0.5) lines. 
Fig. D.1. Continued. 
Appendix E: Baryonic EFE profiles
In Fig. E.1, we show the MOND predicted velocity dispersion when the EFE is calculated solely from the intracluster hot gas distribution M_{gas} (Eq. (33)) instead of M_{C} from hydrostatic equilibrium (Eq. (13)). Although M_{gas} is 1 dex below M_{C} at 1 Mpc, the difference is not sufficient to reach the high measured velocity dispersions.
Fig. E.1. Similar to Fig. 9, with the EFE calculated solely from the intracluster hot gas distribution M_{gas} (magenta dashed lines) instead of the mass distribution M_{C} inferred from hydrostatic equilibrium (black dashed lines, similar to Fig. 9). Only the case of a uniform anisotropy parameter β = 0 is shown. 
Fig. E.1. Continued. 
All Tables
Observed and baryonic accelerations at the halflight radius for the radial acceleration relation.
All Figures
Fig. 1. Spatial distribution of the Coma cluster UDGs modelled in this article projected on a colour image of the cluster created from the Dragonfly g and r images (Abraham & van Dokkum 2014; van Dokkum et al. 2015a). The adopted centre of the Coma cluster is indicated by the white cross. DF44 and DFX1, which are the furthest away from the centre, are indicated in red, the Chilingarian et al. (2019) UDGs are displayed in orange. 

In the text 
Fig. 2. Radial acceleration relation for the Coma cluster UDGs (magenta error bars) together with that obtained for SPARC local spiral galaxies by McGaugh et al. (2016, blue colour scale, red squares, and dashed lines for the binned mean and standard deviation) and for dSphs of the Local Group (Lelli et al. 2017, light grey error bars for their least reliable data and dark grey error bars for their better quality sample). The plain black line shows the fit using Eq. (6) obtained by McGaugh et al. (2016) on the SPARC sample of spiral galaxies. The Coma cluster UDGs fall within the scatter obtained for the dSphs and in reasonable agreement with the empirical function of McGaugh et al. (2016), which reproduces MOND phenomenology. 

In the text 
Fig. 3. Comparison between the three MOND interpolating functions considered in this work, namely ν from Eq. (10), the simple interpolating function from Eq. (12), and ν_{2} from Eq. (11) with α = 2. The fiducial ν is very similar to the simple interpolating function. 

In the text 
Fig. 4. Comparison between the measured lineofsight velocity dispersion of the sample UDGs (grey error bars and/or shaded area), the MOND model in isolation (plain lines), and the Newtonian model (dotted lines). Black corresponds to a uniform anisotropy parameter β = 0 (isotropic), magenta to β = +0.5 (radial), and blue to β = −0.5 (tangential). Data points stem from van Dokkum et al. (2019b) for DF44, van Dokkum et al. (2017) for DFX1, and from Chilingarian et al. (2019) for the other UDGs. The shaded area corresponds to the effective velocity dispersion or, when not available (in the case of DFX1), central stellar velocity dispersion. 

In the text 
Fig. 4. continued. 

In the text 
Fig. 5. Coma cluster cumulative dynamical mass derived by Sanders (2003) from the hydrostatic equilibrium of the Xrayemitting hot gas in Newtonian dynamics (dotted black line) and in MOND (plain black line). While the Newtonian dynamical mass (M_{NC}) continues to increase above 1 Mpc, the MOND mass (M_{C}) converges towards 5.6 × 10^{14} M_{⊙}. 

In the text 
Fig. 6. Distribution of the projected distance of the sample UDGs (plain red) and of the parent Yagi et al. (2016) Coma cluster UDG catalogue (plain black) compared to the average Einasto profile derived by van der Burg et al. (2016, dashed black) and a distribution with no UDGs within 10 Mpc (dotted blue). A constant UDG density within 10 Mpc would closely follow the latter, given the xaxis range. The van der Burg et al. (2016) Einasto distribution provides a reasonable fit to the Coma cluster UDGs; a distribution without UDGs within 10 Mpc or with a constant UDG density are disfavoured. 

In the text 
Fig. 7. Spherical selfgravitating system embedded in a constant external field. 

In the text 
Fig. 8. MOND gravitational acceleration g stemming from the different analytic formulae as a function of the Newtonian external field g_{Ne}, for a given internal Newtonian acceleration g_{N} = a_{0}/100. The magenta data points correspond to the radial acceleration field where g_{N} = a_{0}/100 in the QUMOND numerical integrations plotted in Fig. A.1, which consider a Plummer sphere at distances d = 0.5, 1, 2, 5, 10, and 20 Mpc (from right to left) from a point mass representing the Coma cluster. We refer the reader to Fig. A.1 for the values of g_{Ne} corresponding to each distance d within the Coma cluster. The grey line corresponds to a numerical integration of the average radial acceleration of Eq. (24). 

In the text 
Fig. 9. Comparison between the measured lineofsight velocity dispersion of the sample UDGs (grey error bars and/or shaded area) and the MOND prediction at different distances from the cluster centre. Black lines correspond to a uniform anisotropy parameter β = 0 (isotropic), decreasing magenta lines to β = +0.5 (radially biased). The plain and dotted light grey (red) lines recall the predicted isotropic (radiallybiased) velocity dispersion profiles in the isolated MOND and Newtonian cases, respectively (cf. Fig. 4). Data points stem from van Dokkum et al. (2019b) for DF44, van Dokkum et al. (2017) for DFX1, and from Chilingarian et al. (2019) for the other UDGs. The shaded area corresponds to the effective velocity dispersion or, when not available (in the case of DFX1), the central stellar velocity dispersion. 

In the text 
Fig. 9. continued. 

In the text 
Fig. 10. Tidal susceptibility η as a function of the distance R from the centre of the Coma cluster for the sample UDGs. The distance is assumed to be larger than the projected distance d, i.e., R > d. Points correspond to the average distance d_{mean} assuming the van der Burg et al. (2016) Einasto UDG distribution (cf. Sect. 4.1): plain magenta for MOND with EFE, open black for MOND in isolation. 

In the text 
Fig. 11. Bestfit M/L ratio matching the effective velocity dispersion σ_{eff}, assuming MOND with EFE at the average distance d_{mean}. Crosses correspond to the fiducial M/L ratio values from Table 1, plain black circles to the bestfit values with an anisotropy parameter β = 0, and open magenta circles to the bestfit values with β = 0.5. Plain lines correspond to linear fits, dashed and dotted lines to the averages. The average bestfit M/L ratio is 11.5 with β = 0 and 7.0 with β = +0.5. The fits with β = +0.5 are displayed in Fig. C.1. 

In the text 
Fig. 12. Bestfit acceleration A_{0} matching the effective velocity dispersion when emulating EMOND at the average distance d_{mean}, with the external field derived from the dominant hot gas distribution. Plain black circles correspond to an anisotropy parameter β = 0 and the open magenta circles to β = 0.5. Plain lines correspond to linear fits and dashed lines to the averages. A cluster acceleration A_{0} = 10 × a_{0} provides reasonably good fits to the velocity dispersion measurements, especially when some radial anisotropy is allowed (cf. also Fig. D.1). 

In the text 
Fig. A.1. Comparison between the MOND acceleration with EFE resulting from a QUMOND numerical integration and the analytical expressions ⟨g_{r}⟩ from Eq. (29), g_{∥} from Eq. (60) of Famaey & McGaugh (2012), and g_{⊥} from Eq. (30), considering a Plummer sphere (M_{0} = 3.9 × 10^{8} M_{⊙}, R_{P} = 3.6 kpc) at different distances d from a point mass equal to that of the Coma cluster at that distance given the mass model presented in Sect. 4.1. The plots show that ⟨g_{r}⟩≈g_{⊥} provides a good approximation of the MOND acceleration with EFE. 

In the text 
Fig. C.1. Comparison between the measured lineofsight velocity dispersion of the sample UDGs (grey error bars and/or shaded area) and the MOND prediction with EFE at the average distance d_{mean} for different values of the stellar masstolight ratio M/L, assuming a radially biased anisotropy β = +0.5. The fiducial M/L ratio from Table 1, the bestfit M/L ratio recovering the effective velocity dispersion σ_{eff}, and the large bust still acceptable M/L = 4 are compared. As in Fig. 9, the plain and dotted light red curves recall the predicted radially biased velocity dispersion profiles in the isolated MOND and Newtonian cases, respectively (cf. Fig. 4). 

In the text 
Fig. C.1. Continued. 

In the text 
Fig. D.1. Comparison between the measured lineofsight velocity dispersion of the sample UDGs, as in Fig. 9, to the predicted velocity dispersion at the average distance d_{mean} obtained by emulating EMOND with A_{0} = 10 × a_{0}, shown as dashed black (β = 0) and magenta (β = +0.5) lines. 

In the text 
Fig. D.1. Continued. 

In the text 
Fig. E.1. Similar to Fig. 9, with the EFE calculated solely from the intracluster hot gas distribution M_{gas} (magenta dashed lines) instead of the mass distribution M_{C} inferred from hydrostatic equilibrium (black dashed lines, similar to Fig. 9). Only the case of a uniform anisotropy parameter β = 0 is shown. 

In the text 
Fig. E.1. Continued. 

In the text 
Current usage metrics show cumulative count of Article Views (fulltext 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 4896 hours after online publication and is updated daily on week days.
Initial download of the metrics may take a while.