https://orcid.org/0000-0002-3951-8718
ABSTRACT
Due to their low surface brightness, dwarf galaxies are particularly susceptible to tidal forces. The expected degree of disturbance depends on the assumed gravity law and whether they have a dominant dark halo. This makes dwarf galaxies useful for testing different gravity models. In this project, we use the Fornax Deep Survey (FDS) dwarf galaxy catalogue to compare the properties of dwarf galaxies in the Fornax Cluster with those predicted by the Lambda cold dark matter (ΛCDM) standard model of cosmology and Milgromian dynamics (MOND). We construct a test particle simulation of the Fornax system. We then use the Markov Chain Monte Carlo (MCMC) method to fit this to the FDS distribution of tidal susceptibility η (half-mass radius divided by theoretical tidal radius), the fraction of dwarfs that visually appear disturbed as a function of η, and the distribution of projected separation from the cluster centre. This allows us to constrain the η value at which dwarfs should get destroyed by tides. Accounting for an r′-band surface brightness limit of 27.8 magnitudes per square arcsec, the required stability threshold is |$\eta _{\textrm {destr}} = 0.25^{+0.07}_{-0.03}$| in ΛCDM and |$1.88^{+0.85}_{-0.53}$| in MOND. The ΛCDM value is in tension with previous N-body dwarf galaxy simulations, which indicate that ηdestr ≈ 1. Our MOND N-body simulations indicate that ηdestr = 1.70 ± 0.30, which agrees well with our MCMC analysis of the FDS. We therefore conclude that the observed deformations of dwarf galaxies in the Fornax Cluster and the lack of low surface brightness dwarfs towards its centre are incompatible with ΛCDM expectations but well consistent with MOND.
1 INTRODUCTION
Dwarf galaxies are the smallest and most common type of galaxy. They are characterized by their low mass (|$M \lt\ 10^9 \, \mathrm{M}_{\odot }$|) and low metallicity. Most dwarfs are found in galaxy clusters or near a larger galaxy, making them potentially susceptible to the gravitational effect of these larger structures. The currently standard Lambda cold dark matter (ΛCDM) cosmological model (Efstathiou, Sutherland & Maddox 1990; Ostriker & Steinhardt 1995) provides two different scenarios by which dwarf galaxies can be formed (the dual dwarf galaxy theorem; Kroupa 2012):
From the collapse of dark matter particles into haloes, which then accrete baryonic matter into their potential wells (White & Rees 1978). Such dwarfs are known as ‘primordial dwarf galaxies’ and are expected to be dark matter-dominated; and
From the collapse of overdense regions in tidal tails generated by an interaction between larger, gas-rich galaxies. These so-called ‘tidal dwarf galaxies’ (TDGs) must be free of dark matter as the velocity dispersion of the dark matter particles surrounding the host galaxy is too high to allow for their efficient capture by the shallow potential wells of substructures in the tidal tail (Barnes & Hernquist 1992; Wetzstein, Naab & Burkert 2007). In recent years, cosmological ΛCDM simulations have advanced to the point where they can resolve TDGs (Ploeckinger et al. 2018; Haslbauer et al. 2019b).
Dwarf galaxies can also be classified according to their morphology into early and late types depending on whether they have star-forming regions, which are present only for late-type dwarfs. This category includes blue compact dwarfs and dwarf irregular galaxies like the Magellanic Clouds, while early-type dwarfs include dwarf elliptical (dE) and dwarf spheroidal (dSph) galaxies, with dSphs generally having a lower stellar mass (M⋆). The lowest M⋆ dwarfs tend to have velocity dispersions (σ) which are too high if one assumes virial equilibrium, with σ sometimes even exceeding the escape velocity (Aaronson 1983; Grebel 2001).
This discrepancy relies on the validity of General Relativity and our ability to detect nearly all the matter. ΛCDM is a cosmological model based on General Relativity in which the addition of the dark matter component was motivated by the mismatch between the observed baryonic mass and the mass calculated dynamically from the observed σ assuming the virial theorem (Zwicky 1933). Such acceleration discrepancies are also apparent in the gravity between the Milky Way (MW) and Andromeda (M31; Kahn & Woltjer 1959) and in the outer rotation curves of galaxies (e.g. Babcock 1939; Rubin & Ford 1970; Rogstad & Shostak 1972; Roberts & Whitehurst 1975; Bosma 1978, 1981), as reviewed in Faber & Gallagher (1979). Therefore, the natural ΛCDM explanation for dSphs having such high σ is to assume that most of their mass is in the form of dark matter, in which case they must be primordial dwarfs.
ΛCDM predicts that primordial dwarfs should be distributed nearly isotropically around galaxies (Moore et al. 1999; Gao et al. 2004). However, the dwarf satellite galaxies of the MW, M31, and Centaurus A preferentially align in flattened planes (Lynden-Bell 1976; Ibata et al. 2013; Tully et al. 2015; Müller et al. 2018). This is in significant tension with the ΛCDM model (Kroupa, Theis & Boily 2005). While it was later shown that the distribution of dark matter subhaloes is not supposed to be exactly isotropic due to the preferential accretion of subhaloes along cosmic filaments and the intrinsic triaxiality of dark matter haloes (Libeskind et al. 2005; Zentner et al. 2005), the mild expected flattening is not sufficient to explain the strong correlation in position and velocity space observed in nearby satellite systems (Ibata et al. 2014; Pawlowski et al. 2014; Pawlowski & Kroupa 2020; Müller et al. 2021; Pawlowski & Tony Sohn 2021). The satellite plane problem is reviewed in Pawlowski (2021b), which also considers tentative evidence for more satellite planes beyond the three mentioned above. The Local Group (LG) satellite planes are each in 3.55σ tension with ΛCDM (table 3 of Banik et al. 2021, and references therein), while the satellite plane around Centaurus A is only 0.2 per cent (3.09σ) likely to arise in this paradigm (Müller et al. 2021). These are the only three host galaxies near enough for us to reliably know the phase-space distribution of their satellites. We can approximately combine their low likelihoods in ΛCDM using Gaussian statistics. Since we effectively have χ2 = 3.552 + 3.552 + 3.092 = 34.75, the combined tension can be estimated as the likelihood of the χ2 statistic exceeding this value for three degrees of freedom. This suggests that the LG and Centaurus A satellite planes combined cause a tension of 1.40 × 10−7 (5.27σ). A new interpretation is thus needed to explain the origin of the observed satellite galaxy planes.
Another less widely known problem is the distorted morphologies of MW satellites, which strongly imply that they have been affected by tidal forces (Kleyna et al. 1998; Walcher et al. 2003; Sand et al. 2012). Because the inner region of a satellite galaxy can hardly be affected by tides if it is protected by a dominant dark matter halo (Kazantzidis et al. 2004), |$\lesssim\!{10}{{\ \rm per\ cent}}$| of the MW satellites are expected to be distorted in this paradigm (Kroupa 2012). However, McGaugh & Wolf (2010) found that the majority of the MW satellites present signs of being disturbed, both in their elevated σ and in their observed ellipticity. More recently, Hammer et al. (2020) pointed out that the high σ of dSphs surrounding the MW and their proximity to perigalacticon makes it extremely unlikely for them to be dark matter dominated.
An alternative explanation for the planar distribution of the satellite galaxies is that they are of tidal origin. This is because TDGs are expected to be phase-space correlated (Pawlowski, Kroupa & de Boer 2011; Kroupa 2012; Pawlowski 2018; Haslbauer et al. 2019b). But if the observed satellites are of tidal origin, they would be dark matter free, in which case their high σ for their low M⋆ should be explained in a different way. Kroupa (1997) proposed that due to close encounters of the TDGs with their parent galaxy, the TDGs are highly perturbed. As a result, they should be significantly anisotropic both in terms of their internal structure and their velocity dispersion tensor. More generally, they should not be in dynamical equilibrium, making it incorrect to directly apply the virial theorem to infer the mass from σ as this could cause a significant overestimate. However, purely baryonic dwarfs would be very fragile and easily destroyed, making it unlikely that so many of them exist in the LG right now (Haslbauer et al. 2019a,b). Even if this scenario can explain the high σ of all observed dSphs, ΛCDM would still struggle to explain why almost all observed dwarf satellites of the MW, M31, and Centaurus A are of tidal origin − the quenching mechanisms invoked to solve the missing substructure problem are not expected to be so destructive as to get rid of all observable primordial dwarfs (Kim, Peter & Hargis 2018; Read & Erkal 2019; Webb & Bovy 2020).
Given these difficulties, it is important to note that the properties of both primordial and tidal dSphs can be explained without resorting to the assumption of a surrounding dark matter halo. This entails discarding the ΛCDM cosmological model and using instead an alternative framework, the currently leading contender being Milgromian dynamics (MOND; Milgrom 1983). MOND proposes that the deviations from Newtonian behaviour in the rotation curves of galaxies should be attributed to a departure from Newtonian gravity in the regime of weak gravitational fields (|$g \lesssim a_{_0} = 1.2 \times 10^{-10}$| m s−2 = 3.9 pc Myr−2; Begeman, Broeils & Sanders 1991; Gentile, Famaey & de Blok 2011; McGaugh, Lelli & Schombert 2016). The gravity boost that dwarf galaxies experience in this regime would explain their high σ (McGaugh & Wolf 2010; McGaugh & Milgrom 2013a,b; McGaugh et al. 2021). It would also make the dwarfs less vulnerable to tides and stellar feedback than Newtonian TDGs, which are expected to be extremely fragile. Moreover, MOND offers an elegant scenario for the origin of the LG satellite planes by means of a past flyby encounter between M31 and the MW 9 ± 2 Gyr ago, which is required in MOND (Zhao et al. 2013) and seems to reproduce important aspects of their satellite planes (Banik, O’Ryan & Zhao 2018; Bílek et al. 2018, 2021; Banik et al. 2022a). Therefore, we will focus mainly on ΛCDM and MOND in this contribution.
The planes of satellites problem is one of the most well-known challenges to ΛCDM on galaxy scales (Kroupa et al. 2005; Pawlowski 2018, 2021a,b). It provides a compelling motivation to further investigate dwarf galaxies and question their very nature. Fortunately, the properties of dwarf galaxies make them very suitable for testing different gravity theories. Due to their low mass and especially their low surface brightness, dwarf galaxies can be very susceptible to the effects of gravitational tides. Depending on whether we assume the ΛCDM or MOND model to be valid significantly affects the expected influence of tides on dwarfs. These expectations can then be compared with observations to try and distinguish the models.
Since MOND is a non-linear theory of gravity, the internal dynamics of an object can be affected by the presence of an external field (Bekenstein & Milgrom 1984). This is because the enhancement to the self-gravity depends on the total strength of g, including any external sources. In a dwarf galaxy that experiences a strong gravitational field (usually from a nearby massive galaxy), the MOND boost to the self-gravity will be limited by the dominant external field from the larger central galaxy. This effect becomes stronger as the dwarf gets closer to the central galaxy, to the point that the dwarf can become almost fully Newtonian. Because of this, dwarfs are expected to be more vulnerable to tides in MOND than in ΛCDM, where they would be shielded by their dark matter halo throughout their whole trajectory (Brada & Milgrom 2000).1
In this project, we use the Fornax Deep Survey (FDS) dwarf galaxy catalogue (Venhola et al. 2018, 2019) to compare the observed morphological properties of Fornax Cluster dwarf galaxies with the properties predicted by ΛCDM and MOND. Our aim is to find out if the observed level of disturbance in the Fornax dwarfs is similar to that expected in ΛCDM or MOND, or if neither model works well. ΛCDM could provide too much protection against tides such that it underpredicts the observed level of disturbance in the Fornax dwarfs population. Meanwhile, the lack of protective dark matter haloes around all dwarf galaxies and their reduced self-gravity due to the background cluster gravity could mean that in the MOND scenario, dwarfs are too fragile to survive in the harsh Fornax Cluster environment. Determining which of these scenarios is more likely would help to clarify the physics governing the formation and dynamics of galaxies, whose dominant source of gravity remains unknown.
The layout of this paper is as follows: In Section 2, we describe the FDS dwarf galaxy catalogue and the selection criteria that we apply to it (Section 2.1). In Section 3, we explain the relevant types of gravitational interactions that dwarfs might experience in this cluster: disruption from cluster tides (Section 3.1) and galaxy–galaxy harassment (Section 3.2). These sections consider only Newtonian gravity − the generalization to MOND is presented in Section 3.3. In Section 4, we provide the equations describing the susceptibility of dwarfs to tidal forces in the ΛCDM and MOND models, obtain the tidal susceptibility of the dwarfs in the FDS catalogue for each model (Section 4.1), and show how this theoretical quantity is related to the distribution of the dwarfs (Section 4.2) and whether their observed morphology appears disturbed or undisturbed (Section 4.3). In Section 5, we construct a test particle simulation of the orbits of Fornax dwarfs and, using the Markov Chain Monte Carlo (MCMC) method, fit it to the real Fornax system using the FDS catalogue. In Section 6, we present the results obtained from our MCMC analysis and how they compare to the results of N-body simulations, which we complement with our own N-body simulations of a typical Fornax dwarf in MOND (Section 7). We then discuss our results in Section 8 before concluding in Section 9.
2 THE FDS
The Fornax Cluster is one of the nearest galaxy clusters (dFornax = 20.0 ± 0.3 Mpc; Blakeslee et al. 2009). It is named after its sky position in the Southern hemisphere constellation of Fornax. The cluster is structured into two main components: the main Fornax Cluster centred on NGC 1399, and an infalling subcluster (Fornax A) centred 3° to the south-west in which NGC 1316 is the central galaxy (Drinkwater, Gregg & Colless 2001). The Fornax Cluster contains a significant number of dwarf galaxies with different luminosities, colours, shapes, sizes, and distances to the cluster centre, making it very valuable for studying the properties of dwarf galaxies.
The FDS is the most recent survey of the Fornax Cluster. It includes the main Fornax Cluster and part of the Fornax A subcluster, with a total sky coverage of 26 deg2 (Venhola et al. 2018). The FDS represents a significant improvement in resolution and image depth with respect to the previous spatially complete Fornax Cluster Catalog (FCC; Ferguson 1989). This has allowed the FDS to identify a large number of previously unknown faint galaxies, which can be useful to test the effects of the cluster environment on smaller, more vulnerable galaxies. The FDS reaches the 50 per cent completeness limit at an apparent (absolute) magnitude in the red band of |$M_{r^{\prime }} = -10.5$| (|$m_{r^{\prime }} = 21$|), while the corresponding surface brightness limit is |$\mu _{\mathrm{ e},r^{\prime }} = 26 ~\textrm {mag} ~\textrm {arcsec}^{-2}$|. However, the FDS can still clearly detect some dwarf galaxies down to |$M_{r^{\prime }} = -9$| and |$\mu _{\mathrm{ e},r^{\prime }} = 27.8 ~\textrm {mag} ~\textrm {arcsec}^{-2}$| (Venhola et al. 2018).
The FDS catalogue of dwarf galaxies (Venhola et al. 2017, 2018, 2019) includes 564 dwarf galaxies with 2 × 105 < M⋆/M⊙ < 2 × 109, some in the main Fornax Cluster and others in the infalling subcluster. As in other galaxy clusters, dEs and dSphs are the most common types of dwarf galaxy that can be found in the Fornax Cluster. These are estimated to have an age of tFornax = 10 ± 1 Gyr (Rakos et al. 2001), where tFornax is the age of the elliptical galaxies in Fornax, which we assume to have a similar age to that of the dwarf galaxies. Because of the similarities in some of their morphological properties, the FDS classifies dE and dSph galaxies as one single type, dE. The FDS catalogue also provides information about other properties of the dwarfs. The ones which are relevant for this project are: M⋆, the effective radius, the right ascension (RA) and declination (Dec), the apparent surface brightness in the r′ band, the Sérsic index of the surface brightness profile (Sérsic 1963), the morphological type, the nucleated flag indicating if the dwarf is nucleated or non-nucleated, and the tidal morphology (undisturbed, possibly/mildly disturbed, very disturbed, or unclear; Venhola et al. 2022). The effective radius, the Sérsic index, and the apparent brightness in the r′-band are obtained by fitting the data to a two-dimensional (2D) Sérsic profile (Venhola et al. 2018) using the galfit software (Peng et al. 2002). M⋆ is obtained from the empirical relation between the g′ − i′ colour and mass-to-light (M/L) ratio (Taylor et al. 2011; for further details, see Venhola et al. 2019). The morphological classifications such as the nucleated flags, the Hubble type (Venhola et al. 2018, 2019), and the tidal morphologies are done visually. The tidal morphology is classified in Venhola et al. (2022) based on the following criteria:
Undisturbed: Dwarf galaxies that do not present irregularities, distortions to their shape, or tidal tails;
Possibly/mildly disturbed: Hints of irregularities are present in the outskirts of the dwarf galaxy;
Very disturbed: Dwarf galaxies with tidal tails and/or very clear distortion in the shape; and
Unclear: Nearby bright objects or data artefacts make the classification difficult.
Fig. 1 shows some illustrative examples of dwarfs in these categories.
Images of three FDS dwarfs presenting different levels of disturbance in different colour bands and filters. Each row shows the same dwarf as a red-green-blue colour image (left-hand column) and in the r′ band with a filter enhancing the dwarf’s low surface brightness features (right-hand column). The dwarf in the first, second, and third row is classified as ‘undisturbed’, ‘mildly disturbed’, and ‘very disturbed’, respectively. The horizontal red lines show an angular scale of 10 arcsec, which corresponds to 970 pc at the 20 Mpc distance to the Fornax Cluster.
2.1 Data selection
We remove dwarfs with Rsky > 800 kpc as dwarfs further out mostly belong to the subcluster Fornax A, so including these would contaminate our sample of dwarfs belonging to the main Fornax Cluster (see fig. 4 of Venhola et al. 2019). This leaves us with 353 dwarf galaxies.
3 EFFECTS OF GRAVITATIONAL INTERACTIONS ON DWARFS
Before discussing the gravitational perturbations experienced by Fornax Cluster dwarf galaxies, we first discuss why non-gravitational forces are not expected to perturb Fornax Cluster dwarfs today. Old dwarf galaxies in a cluster environment are expected to be gas-poor. Most dwarfs in the FDS catalogue are classified as early-type galaxies, implying that they are dominated by old stellar populations and are not currently forming new stars. The scarcity of star-forming dwarfs in the Fornax Cluster is consistent with the fact that they are likely to be gas-poor. One important reason for this is ram pressure stripping (Gunn & Gott 1972). This takes place when a galaxy containing a large amount of cold gas moves through a galaxy cluster full of hot gas. The temperature difference and motion between the two gas components generate a pressure gradient that strips the cold gas from the galaxy. Venhola et al. (2020) estimated in the left-hand panel of their fig. 21 that ram pressure stripping of Fornax Cluster dwarfs at the low masses relevant to our analysis should have been quite efficient − the vast majority of the dwarfs in our sample have |$M_{\star } \lt\ 10^8 \, \mathrm{M}_{\odot }$| (Section 2.1). The fact that the Fornax dwarfs are gas-poor has been observationally confirmed by Zabel et al. (2019), who studied the molecular gas in the Fornax Cluster and showed that its dwarfs are gas deficient. Loni et al. (2021) showed the same for neutral hydrogen in FDS dwarfs with M⋆ down to a few times |$10^7 \, {\mathrm{M}_\odot}$|, below which theoretical arguments indicate that the gas reservoir should have been ram pressure stripped by now (see section 7.3.1 of Venhola et al. 2019). Moreover, the colours of the FDS dwarfs also suggest a lack of recent star formation (see their fig. 18). Ongoing gas loss is thus very unlikely to explain the observed disturbances to the structures of some Fornax Cluster dwarfs. We therefore conclude that their internal structure is to a good approximation only affected by gravity from surrounding structures.
The main types of gravitational interaction that can disturb and transform the structure of a dwarf galaxy in the Fornax Cluster are tidal disruption from the cluster’s tidal field and galaxy–galaxy harassment due to encounters with the cluster’s massive elliptical galaxies (see section 7 of Venhola et al. 2019). In the following, we discuss these processes in the context of Newtonian gravity before deriving their generalization to MOND (Section 3.3).
3.1 Disruption from cluster tides
Priors for the free parameters in our model of the Fornax Cluster dwarf galaxy population.
| Parameter | Minimum | Maximum |
|---|---|---|
| |$\textrm {Slope}_{P_\mathrm{r}}$| | −9 | −3 |
| |$\textrm {Slope}_{P_\mathrm{e}}$| | −2 | 2 |
| rcore (Mpc) | 0.01 | 3 |
| Pdist, floor | 0 | 1 |
| Pdist, ceiling | 0 | 1 |
| ηmin, dist | 0 | 5 |
| ηdestr | ηmin, dist | 5 |
| Parameter | Minimum | Maximum |
|---|---|---|
| |$\textrm {Slope}_{P_\mathrm{r}}$| | −9 | −3 |
| |$\textrm {Slope}_{P_\mathrm{e}}$| | −2 | 2 |
| rcore (Mpc) | 0.01 | 3 |
| Pdist, floor | 0 | 1 |
| Pdist, ceiling | 0 | 1 |
| ηmin, dist | 0 | 5 |
| ηdestr | ηmin, dist | 5 |
Priors for the free parameters in our model of the Fornax Cluster dwarf galaxy population.
| Parameter | Minimum | Maximum |
|---|---|---|
| |$\textrm {Slope}_{P_\mathrm{r}}$| | −9 | −3 |
| |$\textrm {Slope}_{P_\mathrm{e}}$| | −2 | 2 |
| rcore (Mpc) | 0.01 | 3 |
| Pdist, floor | 0 | 1 |
| Pdist, ceiling | 0 | 1 |
| ηmin, dist | 0 | 5 |
| ηdestr | ηmin, dist | 5 |
| Parameter | Minimum | Maximum |
|---|---|---|
| |$\textrm {Slope}_{P_\mathrm{r}}$| | −9 | −3 |
| |$\textrm {Slope}_{P_\mathrm{e}}$| | −2 | 2 |
| rcore (Mpc) | 0.01 | 3 |
| Pdist, floor | 0 | 1 |
| Pdist, ceiling | 0 | 1 |
| ηmin, dist | 0 | 5 |
| ηdestr | ηmin, dist | 5 |
3.2 Galaxy–galaxy harassment
To obtain rh, dwarf from the projected effective radius re containing half of the dwarf’s total stellar mass, we use equation (B3) of Wolf et al. (2010), though a good approximation is that rh, dwarf ≈ (4/3)re. Our adopted Mp, * = 1010 M⊙ is the median stellar mass of the large galaxies catalogued in table C1 of Iodice et al. (2019) and in the FCC. In the ΛCDM case, the contribution of the dark halo should be added to this mass. Unlike with the dwarf galaxies, the full extent of the dark halo is considered for the large galaxies because these are expected to be quite robust to cluster tides, so the full halo mass should be considered when estimating the perturbation to a passing dwarf. Venhola et al. (2019) found Mp, ΛCDM = 1011.6 M⊙ following this procedure, which we also verified.
Using a single Mp value for all perturbers gives only an approximate estimate of the heating rate. A more accurate calculation should use the power-law distribution of all the galaxies and make predictions based on that, but this would be extremely difficult. Moreover, the other simplifications assumed throughout the whole calculation of td have a larger impact on the result than taking into account the right distribution of perturbing galaxy masses. Fortunately, we will see that td greatly exceeds a Hubble time, a conclusion which should remain valid even with small adjustments to the calculation. In particular, we will show that considering the mass spectrum of perturbers should affect the estimated heating rate by only a small factor such that td remains very long (Section 4.1).
The rh, p value of the large galaxies is also obtained from the median of all the documented large galaxies (perturbers) in the cluster, yielding rh, p = 4 kpc based on the luminous matter. This is applicable to MOND, but in the ΛCDM case, the rh, p of the large galaxies should account for half of the perturber’s total mass, not only the stellar mass given in the catalogues. This is because the gravitational effect of the dark matter halo also contributes to perturb the stellar content of a passing dwarf. To find out the relation between rh, p and rh, p, ΛCDM, Venhola et al. (2019) looked into the Illustris cosmological simulations (Pillepich et al. 2018) to infer the relation between these two quantities in simulated large galaxies in a galaxy group with a similar mass to the Fornax Cluster, yielding rh, p, ΛCDM/rh, p = 3.6. Therefore, the half-mass radius of the perturbers in ΛCDM is taken to be rh, p, ΛCDM = 14.4 kpc.
To summarize, the disruption time-scale in ΛCDM can be found by directly applying equation (13) once we include the contribution of the dark matter halo to Mdwarf, Mp, and rh, p. In Section 3.3.2, we describe how to obtain the corresponding disruption time-scale expression in MOND.
3.3 Generalization to MOND
The MOND model proposes that Newtonian gravity breaks down in the limit of low accelerations such that the actual gravitational field g is related to the Newtonian field |$g_{_\mathrm{N}}$| according to |$g = \sqrt{a_{_0} g_{_\mathrm{N}}}$|. Milgrom’s constant |$a_{_0} = 1.2 \times 10^{-10}$| m s−2 is a new fundamental acceleration scale added by MOND. Its value has been empirically determined by matching observed galaxy rotation curves (Begeman et al. 1991; Gentile et al. 2011; McGaugh et al. 2016), which MOND does extremely well (Famaey & McGaugh 2012; Lelli et al. 2017; Li et al. 2018). Due to the very small numerical value of |$a_{_0}$| (which may be related to the quantum vacuum, see Milgrom 1999; Senay, Mohammadi Sabet & Moradpour 2021), the behaviour of gravity has never been directly tested in the deep-MOND regime (|$g \ll a_{_0}$|). Indeed, Solar system tests are typically only sensitive to the behaviour of gravity in the regime where g exceeds |$a_{_0}$| by many orders of magnitude (though for a proposed Solar system test in the MOND regime, see Penner 2020).
It is well known that although MOND is capable of fitting the rotation curves of galaxies without dark matter (see the review by Famaey & McGaugh 2012), it cannot fit the temperature and density profiles of galaxy clusters using only their visible mass − MOND still needs an additional contribution to the gravitational field (Sanders 1999; Aguirre, Schaye & Quataert 2001). The central galaxy of the Fornax Cluster (NGC 1399) is no exception (Samurović 2016). To solve this discrepancy and to account for other observations hinting at the presence of collisionless matter in galaxy clusters (most famously in the Bullet cluster; Clowe et al. 2006), it has been proposed that MOND should be supplemented by sterile neutrinos with a rest energy of 11 eV, a paradigm known as the neutrino hot dark matter (νHDM) cosmological model (Angus 2009). νHDM can fit observations of virialized galaxy clusters using the MOND gravity of their directly detected baryons plus the sterile neutrinos (Angus, Famaey & Diaferio 2010). It can also fit the power spectrum of anisotropies in the cosmic microwave background (CMB) because the typical gravitational field at the epoch of recombination was |$\approx\!{20} \, a_{_0}$| and the cosmic expansion history would be standard. Neutrino free streaming reduces the power on small scales compared to ΛCDM, but this is consistent with CMB observations provided the rest energy of the neutrinos exceeds 10 eV (see section 6.4.3 of Planck Collaboration XIII 2016). The gravitational fields from density perturbations would enter the MOND regime only when the redshift |$\lesssim\!{50}$|, before which the MOND corrections to General Relativity should be small (for a more detailed explanation of this model, see Haslbauer, Banik & Kroupa 2020). νHDM relies on the existence of eV-scale sterile neutrinos, but these are also hinted at by several terrestrial experiments (for a recent review, see Berryman et al. 2022).
Equation (16) shows that unlike Newtonian gravity, MOND is a non-linear theory of gravity. A physical consequence of this non-linearity is the so-called external field effect (EFE; Milgrom 1986). This implies that the internal gravity of a system can be weakened by a constant gravitational field from its external environment even if this is completely uniform, violating the strong equivalence principle. The reason is that the MOND boost to the Newtonian gravity is approximately given by ν, which is damped due to the external field. In MOND, the EFE explains why some galaxies like NGC 1052-DF2 have a very low observed velocity dispersion (Famaey, McGaugh & Milgrom 2018; Kroupa et al. 2018; van Dokkum et al. 2018; Haghi et al. 2019a), even though other galaxies like DF44 with similar properties but in a more isolated environment have a much higher velocity dispersion (Bílek, Müller & Famaey 2019; Haghi et al. 2019b; van Dokkum et al. 2019).2 Strong evidence for the EFE has recently been obtained based on the outer rotation curves of galaxies showing a declining trend if the galaxy experiences a significant EFE, while galaxies in more isolated environments have flat outer rotation curves (Haghi et al. 2016; Chae et al. 2020, 2021). For a discussion of observational evidence relating to the EFE, we refer the reader to section 3.3 of Banik & Zhao (2022).
3.3.1 Tidal radius
3.3.2 Galaxy–galaxy harassment
When a dwarf interacts with a massive galaxy in the Fornax Cluster environment, we need to consider both the gravity from the elliptical and the background EFE due to the cluster potential. As in Section 3.2, we estimate the perturbation to the dwarf by assuming it is a collection of test particles that receive some impulse |$\boldsymbol {u}$| from the elliptical, with the heating rate of the dwarf proportional to the square of |$\vert \Delta \boldsymbol {u} \vert$|, the spread in |$\boldsymbol {u}$| across the dwarf. Once it has moved away from the elliptical, the binding energy of the dwarf is given by equation (12) but with G → Geff as discussed above. The main difficulty lies in estimating the energy gained by the dwarf due to interactions with impact parameter b, which for a high-velocity encounter is approximately the same as the closest approach distance between the dwarf and the elliptical.
We need to consider encounters in two different regimes:
The QN regime in which |$g_\mathrm{c} \ll a_{_0}$| dominates over gravity from the elliptical; and
The isolated deep-MOND (IDM) regime in which the gravity from the elliptical dominates over gc but is still much weaker than |$a_{_0}$|.
We do not need to consider the Newtonian regime because the perturbers have a radius that is numerically similar to their MOND radius for the parameters given in Section 3.2. This is not unique to the Fornax Cluster: Elliptical galaxies generally have a size similar to their MOND radius (Sanders 2000). This is because if the initial radius was much smaller and the system is nearly isothermal, then a significant proportion of the mass in the outskirts would be moving faster than the Newtonian escape velocity, causing the system to expand to its MOND radius (Milgrom 1984, 2021).
4 TIDAL SUSCEPTIBILITY
If a dwarf has strong self-gravity (e.g. due to being surrounded by a dark matter halo or being in the deep-MOND regime), then the point at which the tidal force of the cluster will start to dominate over the self-gravity of the dwarf will be far from the centre of the dwarf. Therefore, the dwarf’s rtid will be large and its tidal susceptibility will be low. Such a dwarf should be little disturbed by the cluster tides. If instead the dwarf has only weak self-gravity (e.g. because it is a TDG with little dark matter or because it is a MONDian dwarf but the EFE from the cluster is very significant), then the point at which the tidal force of the cluster will start to dominate over the self-gravity of the dwarf will be close to the dwarf’s centre. Its rtid will then be small and its tidal susceptibility high. Such a dwarf would be significantly disturbed by tides. In the extreme case that rtid ≪ rh (ηrtid ≫ 1), the dwarf will be destroyed within a few dynamical times. As a result, we need to consider the maximum value of ηrtid attained throughout the trajectory, i.e. we need to evaluate ηrtid at pericentre.
Although our definitions for ηrtid and ηhar differ somewhat because the former is a ratio of radii while the latter is a ratio of time-scales, both definitions share the feature that low values of η indicate that a dwarf should be little affected by the process under consideration. In principle, there should not be any dwarf galaxies for which either η ≫ 1. It is possible to have η slightly above 1 due to projection effects and other subtleties like the time required to achieve destruction, which can be significant for ηrtid as multiple pericentre passages may be required and the orbital period can be long (Section 3.1). However, we should very seriously doubt the validity of any theory which tells us that a significant fraction of the dwarf galaxies in a galaxy cluster have ηrtid ≫ 1 or ηhar ≫ 1. It is harder to falsify a theory in the opposite limit where it yields very low values for both measures of η for all the dwarfs in a galaxy cluster. In this case, we could gain evidence against the theory if there is strong evidence that the dwarf galaxy population has been significantly affected by tides. In this project, we apply these considerations to the dwarf galaxy population in the Fornax Cluster.
4.1 Tidal susceptibility of the Fornax dwarfs
Our first quantitative result is the susceptibility of dwarfs in the FDS catalogue to cluster tides, which we calculate in ΛCDM and MOND using equations (11) and (20), respectively. We show the results as histograms in the top row of Fig. 2, with ΛCDM shown in the left-hand panel and MOND in the right-hand panel. The ηrtid values are ≈5× higher in MOND than in ΛCDM. Since an isolated dwarf has a similar amount of self-gravity in both frameworks by construction, the difference in ηrtid values is primarily caused by the EFE weakening the self-gravity of a MONDian dwarf as it approaches the cluster centre (Section 1). This effect does not exist for a ΛCDM dwarf, which would retain the same dark matter fraction within its baryonic extent throughout its trajectory.
Histogram of the tidal susceptibility values of Fornax Cluster dwarfs in ΛCDM (left-hand column) and MOND (right-hand column) to tides from the overall cluster potential (top row) and harassment by interactions with individual massive galaxies (bottom row). The bin widths are: 0.01 (top-left panel), 0.05 (top-right panel), 0.01 (bottom-left panel), and 0.005 (bottom-right panel). Notice the different ηrtid scales for ΛCDM and MOND. In both theories, we typically have ηhar ≪ ηrtid.
The bottom row of Fig. 2 shows the susceptibility of FDS dwarfs to galaxy–galaxy harassment according to ΛCDM (equation 13) and MOND (equation 28). In both theories, the histogram of ηhar peaks at very low values such that ηhar ≪ ηrtid and ηhar ≪ 1. Therefore, both frameworks predict that the FDS dwarfs should be little affected by interactions with massive elliptical galaxies in the Fornax Cluster.
This implies at face value that in ΛCDM, the observed signs of tidal disturbance (section 7.4 of Venhola et al. 2022) cannot be assigned to either cluster tides or to harassment. Since we explore the impact of cluster tides more carefully later in this contribution, we briefly reconsider our calculation of ηhar. As explained in Section 3.2, one simplifying assumption we made is that there are 48 equal mass and equal size perturbers within the 0.77 Mpc virial radius of the Fornax Cluster. However, the heating rate due to any individual perturber scales as |$\dot{E} \propto \left(M_\mathrm{p}/r_{\mathrm{h},\mathrm{ p}} \right)^2$| (equation 23). We can use this to find the ratio |$\dot{E}/\dot{E}_{\textrm {fid}}$| between the heating rate due to individual perturbers and the assumed heating rate |$\dot{E}_{\textrm {fid}}$| for an ‘average’ perturber with |$M_{\star } = 10^{10} \, \mathrm{M}_{\odot }$| and rh, p = 4 kpc, taking into account that the actual mass and size are larger in ΛCDM and assuming a de Vaucouleurs profile for the stars (de Vaucouleurs 1948). We obtain that in descending order of M⋆, the ratio |$\dot{E}/\dot{E}_{\textrm {fid}}$| for the perturbers listed in table C1 of Iodice et al. (2019) is 14.7 (FCC 219), 42.7 (FCC 167), 10.3 (FCC 184), 4.76 (FCC 161), 5.41 (FCC 147), 11.8 (FCC 170), 1.06 (FCC 276), 1.93 (FCC 179), and 0.13 (FCC 312). Other perturbers have |$M_{\star } \lt\ 10^{10} \, \mathrm{M}_{\odot }$|, so we assume their contribution to the heating rate is small. Adding up the above ratios and averaging over 48 perturbers (many of which are too low in mass to appreciably harass Fornax dwarfs), we get that |$\dot{E}/\dot{E}_{\textrm {fid}}$| is on average 1.9. Therefore, using a more accurate treatment of the heating rate would not change our conclusion that the FDS dwarfs are not really susceptible to galaxy–galaxy harassment: Doubling all the ηhar values would still lead to its distribution having a mode <0.1 and all the dwarfs having ηhar < 1.
Moreover, using tFornax as the time-scale for interactions is an optimistic assumption − dwarfs in ΛCDM may have been accreted by the cluster long after they formed, while in MOND they could be TDGs that formed more recently (Renaud, Famaey & Kroupa 2016). This implies that the dwarfs would not have experienced that many encounters with elliptical galaxies, which themselves might only have been accreted ≪10 Gyr ago. As an example, we may consider the case of FCC 219 ≡ NGC 1404, the most massive perturber listed in table C1 of Iodice et al. (2019) in terms of M⋆. Its radial velocity exceeds that of the brightest cluster galaxy NGC 1399 by 522 km s−1, but modelling indicates that the relative velocity could be higher still as most of it should lie within the sky plane (Machacek et al. 2005). Moreover, NGC 1404 appears to lie in front of the Fornax Cluster: Its heliocentric distance is only 18.7 ± 0.3 Mpc (Hoyt et al. 2021), whereas the distance to NGC 1399 is 20.0 ± 0.3 Mpc (Blakeslee et al. 2009). Detailed modelling in a ΛCDM context indicates that although NGC 1404 is not on a first infall, it has likely spent |$\lesssim\!{3}$| Gyr within the cluster (Sheardown et al. 2018). During this time, the high relative velocity would have reduced the heating rate on any dwarf galaxy that it came near (equation 23). It is therefore clear that ηhar is overestimated by assuming that both all the dwarfs and all 48 massive ellipticals were in the virial volume of the Fornax Cluster over the last 10 Gyr.
Based on this, we will neglect the role of harassment in what follows and focus on cluster tides.4 Thus, η will be used to mean ηrtid unless stated otherwise. An important example of this is our discussion of Newtonian TDGs that are purely baryonic, where ηhar plays an important role (Appendix D).
4.2 Testing the effect of cluster tides on Fornax dwarfs
A significant fraction of the FDS dwarfs appear disturbed in a manual visual classification (Fig. 1, see also Venhola et al. 2022). To check if cluster tides are truly the main mechanism responsible for the apparent disturbance of the Fornax dwarfs − as our results in Fig. 2 seem to suggest − in Fig. 3 we plot the projected distance of the selected Fornax dwarfs against the ratio between their effective radius re and rmax, where rmax is the maximum re that the dwarf could have to remain detectable given its M⋆ and the FDS detection limit of 27.8 mag arcsec−2. Dwarfs with larger size at fixed stellar mass − i.e. lower surface brightness dwarfs − are more susceptible to tides and will be more easily destroyed, especially near the cluster centre where the tides are stronger. In Fig. 3, we can see a deficit of low surface brightness dwarfs near the cluster centre. The absence of dwarfs in this region of the parameter space cannot be explained by the survey detection limit as we find an increasing number of dwarfs with the same or lower surface brightness at larger Rsky, e.g. if we consider a horizontal line at re/rmax = 0.4. This tendency is highlighted in Fig. 3 using a sloped dotted line that appears to be a tidal edge. Further from the cluster, its tides become weaker, so it is quite possible that dwarfs in this region are not much affected by tides.
Distribution of the projected distances of Fornax Cluster dwarfs against re/rmax, where re is the projected half-light radius and rmax is the maximum re at fixed M⋆ that the dwarf can have to remain detectable given the surface brightness limit of the survey of 27.8 mag arcsec−2 in the r′ band. Dwarfs visually classified as ‘undisturbed’ are shown in blue, while those classified as ‘disturbed’ are shown in red. Notice the lack of low surface brightness dwarfs near the cluster centre. We have emphasized this by drawing a dashed grey line for illustrative purposes, which we interpret as a tidal edge. This interpretation is bolstered by the lack of dwarfs above this line and the high proportion of disturbed dwarfs just below this line.
Additional evidence for the importance of tides towards the cluster centre comes from the colours of the dots in Fig. 3, which indicate whether the dwarf visually appears disturbed (red) or undisturbed (blue). Just below the claimed tidal edge, we would expect that the dwarfs are much more likely to appear disturbed as they should be close to the threshold of being destroyed altogether. This is indeed apparent: The proportion of disturbed galaxies is much higher in this part of the parameter space.5
To emphasize this trend further, we use Fig. 4 to show the observed fraction of disturbed dwarfs (fd) in different Rsky bins. This is found as fd = S/T, with the uncertainties calculated using binomial statistics as |$\sqrt{S \left(T - S \right)/T^3}$|, where T is the number of galaxies in each Rsky bin and S ≤ T is the number of these galaxies which appear disturbed.6 As expected from our previous results, fd is very high in the central 200 kpc of the Fornax Cluster. Although fd is very low further out, it is still non-zero and remains so out to the largest distances covered by our data set. We attribute this to the complexities of visually assessing whether a dwarf is tidally disturbed: If a dwarf appears asymmetric due to observational difficulties or due to a dense star cluster on one side, this could lead to a false positive. It is also possible that the dwarf is genuinely disturbed due to a recent close encounter with a massive galaxy in the cluster, which could happen even in the cluster outskirts. When we construct a detailed model of the Fornax Cluster dwarf galaxy population in Section 5.2.4, we will need to allow a non-zero likelihood that a dwarf appears disturbed even if it is unaffected by cluster tides.
The proportion of Fornax dwarfs that appear disturbed in different projected separation bins of width 200 kpc. The error bars show the binomial uncertainty assuming the likelihood of appearing disturbed is the same as the proportion of disturbed dwarfs.6
4.3 Correlating tidal susceptibility with the observed level of disturbance
We use Fig. 5 to plot the mean and standard deviation obtained in this way against the central η value for the bin under consideration. In both ΛCDM and MOND, a clear trend is apparent whereby dwarfs with higher η are more likely to appear disturbed. We quantify this by dividing the FDS sample into two subsamples where η is below or above some threshold ηt, thereby assuming only a monotonic relation between fd and η that is not necessarily linear. Appendix C explains how we obtain the likelihood that the same fd can explain the number of disturbed dwarfs and the total number of dwarfs in both subsamples given binomial uncertainties. Using this method, we find that the ‘signal’ is maximized in ΛCDM if we use ηt = 0.36, in which case the null hypothesis of fd being the same in both subsamples can be rejected at a significance of P = 4.1 × 10−3 (2.87σ). If instead we use MOND, the optimal ηt = 0.85 and the significance rises to P = 4.4 × 10−4 (3.52σ).7 Though both theories imply that fd is higher in the high η subsample, fd starts rising at a much lower value of η in ΛCDM than in MOND, as clearly shown by the optimal ηt values. We may expect that dwarfs start to look disturbed when their half-mass radius is about the same as their tidal radius, so fd should start rising only when |$\eta \gtrsim 0.5$|. This is not the case in ΛCDM, which implies that dwarfs are more likely to be classified as disturbed once their |$\eta \gtrsim 0.1{-}0.2$|. A dwarf with such a low η should be little affected by tides, indicating a problem for this framework. In the MOND case, we see that dwarfs start being classified as disturbed more often once their |$\eta \gtrsim 1{-}1.5$|, which is much more plausible physically.
The likelihood that a dwarf appears disturbed in each tidal susceptibility bin in ΛCDM (orange) and MOND (blue). The value and uncertainty are calculated using binomial statistics (equation 32) and plotted at the centre of each bin. The bin width is 0.5 for MOND and 0.1 for ΛCDM. In both cases, the last bin also includes all dwarfs with higher η. Notice that the likelihood of a dwarf appearing disturbed rises with η. The higher uncertainties at high η are due to a small sample size (see the η distribution in Fig. 2).
Another important aspect is the overall distribution of η, whose decline towards the highest bin is responsible for a larger uncertainty in the probability of appearing disturbed. The distribution of η is shown explicitly in the top row of Fig. 2. There are no ΛCDM dwarfs with η > 0.7, even though a dwarf with η = 0.7 should still be tidally stable. In MOND, the maximum η ≈ 3, though there are very few dwarfs with η > 1.5. The high calculated η for these dwarfs could indicate that they lie very close to the cluster centre in projection but not in reality. To handle such projection effects and other uncertainties like the unknown orbital eccentricity distribution of the dwarfs, we next construct a test mass simulation of the dwarf galaxy population in the Fornax Cluster.
5 TEST MASS SIMULATION OF THE FORNAX CLUSTER
In order to quantify the aforementioned trends and thereby obtain the range of values that the minimum η required for disturbance and the η required for destruction can have to be consistent with observations − both in ΛCDM and in MOND − we need to construct a forward time-evolution model of the Fornax Cluster. With this forward model, we can also account for projection effects that can make dwarfs appear closer to the cluster centre than they actually are. In this section, we describe the set-up of the simulated Fornax system with test masses, as well as the methods that we use to quantify the properties of the Fornax dwarfs and their orbits. Here, we focus only on those dwarfs classified as ‘non-nucleated’ as this type of galaxy is more numerous than the ‘nucleated’ type. Moreover, having the same deprojection method (Appendix A) for all dwarfs will simplify the analysis. Removing the nucleated dwarfs from the sample leaves us with 279 dwarfs.
5.1 Orbit integration
The first step in building a simulation of test masses orbiting in the observed cluster potential is to generate a grid of orbits for a wide range of semimajor axis (Ri) and eccentricity (e) values, with the integrations started at R = Ri. The initial radii have a range of values from 15 to 2015 kpc, while the eccentricities cover the full range of values for an ellipse (0 < e < 1). The grid is divided into 100 × 100 cells. Initially, we assign the test mass a mass and half-mass radius which are typical for a Fornax dwarf (Mdwarf = 3.16 × 107 M⊙ and rh = 0.84 kpc), but these values are not relevant as the results will be rescaled later according to the distribution of dwarf densities in the system (Section 5.2.3).
5.2 Assigning probabilities to the orbits
The orbital and internal properties of the Fornax dwarfs (e.g. the radial profile of the orbits, the distribution of their eccentricities, and the likelihood of appearing perturbed) follow certain probability distributions. Because of this, we assign probabilities to each of our simulated orbits by fitting them to a few crucial observed properties (next subsection) in order to make our simulated system as similar as possible to the observed Fornax dwarf galaxy system. The parameters governing these probability distributions are described below.
5.2.1 Number density of dwarfs
5.2.2 The eccentricity distribution
5.2.3 Distribution of dwarf densities
The density ρ of each Fornax dwarf within its rh can be inferred from the data in the FDS catalogue using |$\rho = 3 M_{\star }/ \left(8 \mathrm{\pi } r_\mathrm{h}^3\right)$|. Fig. 6 shows a histogram of the so-obtained densities of these dwarfs, from which it can be seen that the FDS distribution of |$\log _{10} \, \rho$| follows a Gaussian distribution with mean −2.74 in units of M⊙ pc−3. Therefore, when we assign a density to each of the simulated dwarfs obtained in Section 5.1, we associate a probability to this density according to a lognormal distribution. This is assumed to be independent of Ri since the central region of a cluster should be able to accrete dwarfs that formed further out, leading to mixing of dwarfs that formed in different positions within the cluster.
The distribution of each dwarf galaxy’s mean baryonic density ρ within its half-mass radius. The orange vertical line shows the sample mean, while the magenta lines offset by ±0.57 dex show the standard deviation around it. The grey line shows the density of a dwarf corresponding to the observational surface density detection limit of Σmin = 0.26 M⊙ pc−2 assuming the mean |$M/L_{r^{\prime }} = 1.10$| and mean |$\rho / \Sigma = 0.59 \, \textrm {kpc}^{-1}$|. If instead we add (subtract) the standard deviation in the |$M/L_{r^{\prime }}$| ratios, we have that |$\Sigma _{\textrm {min}} = 0.35 \left(0.17 \right) \, \mathrm{M}_{\odot }\,\textrm {pc}^{-2}$|. From these and by adding (subtracting) the standard deviation in the ρ/Σ ratios, we obtain the ρ value given by the dashed (solid) black line.
The relation between mass and luminosity in the r′ band for the non-nucleated sample of dwarfs in the FDS catalogue, with both shown on a log10 scale. The assumed M/L ratios are based on empirical relations with the colour (Section 2). The dotted grey line shows the linear regression, while the solid grey line shows our adopted fit assuming a slope of 1. The dashed (solid) black line shows one standard deviation above (below) the mean |$M_{\star }/L_{r^{\prime }}$|. The horizontal red line at |$10^{7.2} \, \mathrm{M}_{\odot }$| shows the stellar mass below which core formation is inefficient in ΛCDM (Section 8).
The relation between the average 3D mass density of baryons within their half-mass radius and their surface mass density within their projected half-light radius for our galaxy sample. The dotted grey line shows the linear regression, while the solid grey line shows our adopted fit assuming a slope of 1. The dashed (solid) black line shows one standard deviation above (below) the mean ρ/Σ. The red star shows the values for the dwarf galaxy used in our N-body simulations (Section 7).
In the ΛCDM case, we need to include the halo mass within the baryonic extent of each dwarf (equation 10), leading to higher volume densities. This causes a steeper slope and a larger amount of scatter in the mass–luminosity relation, making it difficult to follow the above-mentioned method. To keep the procedure similar, we set ρmin to a value 0.09 dex below ρmin, FDS as this is the gap assumed for MOND. The steps involved with this model are shown in Appendix E.
5.2.4 Disturbance to the dwarf structure
5.3 Comparison with observations
The observed parameters of the Fornax dwarfs that we aim to reproduce in our simulation are:
The distribution of sky-projected distances (Rsky) to the cluster centre;
The distribution of apparent η values at pericentre (ηobs); and
The disturbed fraction of dwarfs as a function of ηobs.
To do the comparison, we start by dividing the range of Rsky and ηobs into several bins. We then classify the observed dwarfs into these bins according to their values of projected distance or estimated η at pericentre. To obtain the probability for a dwarf to have a projected distance or ηobs which falls in the range of values delimited by each of these bins, we count the number of dwarfs in each bin and compare it to the total number of dwarfs. To obtain the probability of disturbance, we count the number of dwarfs classified as disturbed in each ηobs bin and compare it to the total number of dwarfs in that bin.
For the simulated sample (i.e. the dwarfs generated for all possible combinations of Ri, e, ρ, and θ), we consider the same bins as for the observed sample. For each bin, we add the probability that each simulated dwarf has Rsky or ηobs values that fall in the range given by the bin. We then normalize this by the sum of all the probabilities in all bins. For the probability of disturbance, we apply an additional factor of Pdist to the likelihood of each (Ri, e, ρ, θ) combination and add this to the appropriate ηobs bin. We then divide this sum by the probability of ηobs falling in that bin (i.e. without considering Pdist).
In order to maximize this Ptotal, we leave as free parameters: rcore, |$\textrm {Slope}_{P_\mathrm{r}}$|, |$\textrm {Slope}_{P_\mathrm{e}}$|, ηmin, dist, ηdestr, Pdist, floor, and Pdist, ceiling. We explore this set of parameter values using the MCMC method discussed below.
5.3.1 MCMC analysis
The MCMC method generates a sequence of parameter values in such a way that their frequency distribution matches the posterior inference on the model parameters. The basic idea is to start with some initial guess for the parameters with likelihood Ptotal and generate a proposal by adding Gaussian random perturbations to the parameters, leading to a likelihood of Pnext with the revised parameters. The proposal is accepted if Pnext > Ptotal or if a random number drawn uniformly from the range (0–1) is <Pnext/Ptotal. If the proposal is rejected, the parameter perturbations are not applied but the previous parameters must be recorded once more.
We run a total of 105 trials in each chain and check that the acceptance fraction is close to 0.234, the optimal acceptance rate for an efficient MCMC algorithm (Gelman, Gilks & Roberts 1997). This is achieved by rerunning the chain a few times to determine the optimal step sizes for the parameter perturbations. To ensure that the algorithm chooses physically reasonable parameter values, we impose the priors listed in Table 1. If the algorithm chooses a value for any of these parameters outside the specified range, it is asked to draw another proposal, but this does not count as a new MCMC trial. We let the algorithm consider a sufficiently large number of proposals at each stage in the chain that we are sure to obtain a physically plausible proposal for the parameter combination to try next, even if this is rejected because it fits the observations poorly.
To prevent the MCMC algorithm from starting with a set of values which is too far away from the optimal set, we first fit the simulation’s free parameters to the observations using a gradient ascent algorithm (Fletcher & Powell 1963). This maximizes Ptotal by increasing or decreasing the step size according to how much Ptotal increased or decreased with respect to the previous set of parameter values that it tested. This is done until the step size becomes very small, indicating that the algorithm cannot increase Ptotal any more. Then the algorithm converges and returns the optimal set of parameter values.
6 RESULTS OF THE STATISTICAL ANALYSIS
We present our best-fitting model in each theory (Section 6.1) before discussing the parameter uncertainties obtained with the MCMC method (Section 6.2).
6.1 The best-fitting model
The optimal set of parameters found by the gradient ascent algorithm are given in columns 2 and 3 of Table 2 for ΛCDM and MOND, respectively. These are the initial values at which we start the MCMC chains. Due to the use of 105 trials, the MCMC method provides a set of parameter values (a model) that fits the observations slightly better (higher Ptotal in equation 49) than we achieved with gradient ascent. The best-fitting parameter values in the MCMC chain are also given in Table 2 (columns 4–5) along with the goodness of fit to the observations (last row). In this regard, there is little difference between the theories, though the optimal parameters are rather different. We will return to this later (Section 8).
| Gradient ascent | MCMC | |||
|---|---|---|---|---|
| Parameter | ΛCDM | MOND | ΛCDM | MOND |
| |$\textrm {Slope}_{P_\mathrm{r}}$| | −3.77 | −3.67 | −5.85 | −4.55 |
| |$\textrm {Slope}_{P_\mathrm{e}}$| | −1.55 | 0.34 | −1.98 | −1.70 |
| rcore | 0.62 | 0.65 | 1.35 | 0.90 |
| Pdist, floor | 0.09 | 0.04 | 0.10 | 0.02 |
| Pdist, ceiling | 0.65 | 0.76 | 0.54 | 0.53 |
| ηmin, dist | 0.11 | 0.24 | 0.12 | 0.10 |
| ηdestr | 0.24 | 1.88 | 0.23 | 1.24 |
| log10Ptotal | −30.69 | −32.46 | −30.53 | −32.25 |
| Gradient ascent | MCMC | |||
|---|---|---|---|---|
| Parameter | ΛCDM | MOND | ΛCDM | MOND |
| |$\textrm {Slope}_{P_\mathrm{r}}$| | −3.77 | −3.67 | −5.85 | −4.55 |
| |$\textrm {Slope}_{P_\mathrm{e}}$| | −1.55 | 0.34 | −1.98 | −1.70 |
| rcore | 0.62 | 0.65 | 1.35 | 0.90 |
| Pdist, floor | 0.09 | 0.04 | 0.10 | 0.02 |
| Pdist, ceiling | 0.65 | 0.76 | 0.54 | 0.53 |
| ηmin, dist | 0.11 | 0.24 | 0.12 | 0.10 |
| ηdestr | 0.24 | 1.88 | 0.23 | 1.24 |
| log10Ptotal | −30.69 | −32.46 | −30.53 | −32.25 |
| Gradient ascent | MCMC | |||
|---|---|---|---|---|
| Parameter | ΛCDM | MOND | ΛCDM | MOND |
| |$\textrm {Slope}_{P_\mathrm{r}}$| | −3.77 | −3.67 | −5.85 | −4.55 |
| |$\textrm {Slope}_{P_\mathrm{e}}$| | −1.55 | 0.34 | −1.98 | −1.70 |
| rcore | 0.62 | 0.65 | 1.35 | 0.90 |
| Pdist, floor | 0.09 | 0.04 | 0.10 | 0.02 |
| Pdist, ceiling | 0.65 | 0.76 | 0.54 | 0.53 |
| ηmin, dist | 0.11 | 0.24 | 0.12 | 0.10 |
| ηdestr | 0.24 | 1.88 | 0.23 | 1.24 |
| log10Ptotal | −30.69 | −32.46 | −30.53 | −32.25 |
| Gradient ascent | MCMC | |||
|---|---|---|---|---|
| Parameter | ΛCDM | MOND | ΛCDM | MOND |
| |$\textrm {Slope}_{P_\mathrm{r}}$| | −3.77 | −3.67 | −5.85 | −4.55 |
| |$\textrm {Slope}_{P_\mathrm{e}}$| | −1.55 | 0.34 | −1.98 | −1.70 |
| rcore | 0.62 | 0.65 | 1.35 | 0.90 |
| Pdist, floor | 0.09 | 0.04 | 0.10 | 0.02 |
| Pdist, ceiling | 0.65 | 0.76 | 0.54 | 0.53 |
| ηmin, dist | 0.11 | 0.24 | 0.12 | 0.10 |
| ηdestr | 0.24 | 1.88 | 0.23 | 1.24 |
| log10Ptotal | −30.69 | −32.46 | −30.53 | −32.25 |
Using these parameters, Fig. 9 shows the simulated and observed probability distributions of Rsky, ηobs, and disturbed fraction versus ηobs, revealing a good overall fit to the observations in both theories.8 In particular, the rising likelihood of a dwarf appearing disturbed as a function of ηobs is nicely reproduced by the best-fitting models.
Comparison between observations and the best-fitting simulation in ΛCDM (left-hand column) and MOND (right-hand column) in terms of the distribution of projected separation Rsky from the cluster (first row), tidal susceptibility ηobs (second row), and likelihood that a dwarf appears disturbed as a function of ηobs (third row). The observations (blue points with error bars) and the best-fitting simulation in each theory (red points) are plotted at the centre of each bin, but dithered slightly along the x-axis for clarity in case the model works well. The bin width in Rsky is 100 kpc in both theories. For ηobs, the bin width is 0.15 in ΛCDM and 0.65 in MOND.
6.2 Parameter uncertainties
To fit the test mass simulation of the Fornax dwarf galaxy system to its observed properties, we require several free parameters in the model (Section 5). Having discussed the values of these parameters in the most likely model (Table 2), we now find the most likely value of each parameter and its uncertainty. This is somewhat different because instead of considering the most likely model, we use the MCMC chain to obtain the posterior inference on each model parameter, which we then characterize using its mode and 1σ confidence interval. The results are shown in Table 3.
The most likely value and 1σ confidence interval of each model parameter in our test mass simulation of the Fornax Cluster dwarf galaxy population, based on 105 MCMC trials.
| Parameter | ΛCDM | MOND |
|---|---|---|
| |$\textrm {Slope}_{P_\mathrm{r}}$| | |$-7.43_{-0.99}^{+2.24}$| | |$-7.58_{-0.88}^{+2.18}$| |
| |$\textrm {Slope}_{P_\mathrm{e}}$| | |$-1.65_{-0.30}^{+1.80}$| | |$0.75_{-1.22}^{+1.20}$| |
| rcore | |$2.00_{-0.98}^{+0.34}$| | |$2.02_{-0.88}^{+0.52}$| |
| Pdist, floor | |$0.10_{-0.03}^{+0.03}$| | |$0.07_{-0.03}^{+0.04}$| |
| Pdist, ceiling | |$0.49_{-0.15}^{+0.30}$| | |$0.79_{-0.20}^{+0.15}$| |
| ηmin, dist | |$0.11_{-0.06}^{+0.05}$| | |$0.24_{-0.19}^{+0.24}$| |
| ηdestr | |$0.25_{-0.03}^{+0.07}$| | |$1.88_{-0.53}^{+0.85}$| |
| Parameter | ΛCDM | MOND |
|---|---|---|
| |$\textrm {Slope}_{P_\mathrm{r}}$| | |$-7.43_{-0.99}^{+2.24}$| | |$-7.58_{-0.88}^{+2.18}$| |
| |$\textrm {Slope}_{P_\mathrm{e}}$| | |$-1.65_{-0.30}^{+1.80}$| | |$0.75_{-1.22}^{+1.20}$| |
| rcore | |$2.00_{-0.98}^{+0.34}$| | |$2.02_{-0.88}^{+0.52}$| |
| Pdist, floor | |$0.10_{-0.03}^{+0.03}$| | |$0.07_{-0.03}^{+0.04}$| |
| Pdist, ceiling | |$0.49_{-0.15}^{+0.30}$| | |$0.79_{-0.20}^{+0.15}$| |
| ηmin, dist | |$0.11_{-0.06}^{+0.05}$| | |$0.24_{-0.19}^{+0.24}$| |
| ηdestr | |$0.25_{-0.03}^{+0.07}$| | |$1.88_{-0.53}^{+0.85}$| |
The most likely value and 1σ confidence interval of each model parameter in our test mass simulation of the Fornax Cluster dwarf galaxy population, based on 105 MCMC trials.
| Parameter | ΛCDM | MOND |
|---|---|---|
| |$\textrm {Slope}_{P_\mathrm{r}}$| | |$-7.43_{-0.99}^{+2.24}$| | |$-7.58_{-0.88}^{+2.18}$| |
| |$\textrm {Slope}_{P_\mathrm{e}}$| | |$-1.65_{-0.30}^{+1.80}$| | |$0.75_{-1.22}^{+1.20}$| |
| rcore | |$2.00_{-0.98}^{+0.34}$| | |$2.02_{-0.88}^{+0.52}$| |
| Pdist, floor | |$0.10_{-0.03}^{+0.03}$| | |$0.07_{-0.03}^{+0.04}$| |
| Pdist, ceiling | |$0.49_{-0.15}^{+0.30}$| | |$0.79_{-0.20}^{+0.15}$| |
| ηmin, dist | |$0.11_{-0.06}^{+0.05}$| | |$0.24_{-0.19}^{+0.24}$| |
| ηdestr | |$0.25_{-0.03}^{+0.07}$| | |$1.88_{-0.53}^{+0.85}$| |
| Parameter | ΛCDM | MOND |
|---|---|---|
| |$\textrm {Slope}_{P_\mathrm{r}}$| | |$-7.43_{-0.99}^{+2.24}$| | |$-7.58_{-0.88}^{+2.18}$| |
| |$\textrm {Slope}_{P_\mathrm{e}}$| | |$-1.65_{-0.30}^{+1.80}$| | |$0.75_{-1.22}^{+1.20}$| |
| rcore | |$2.00_{-0.98}^{+0.34}$| | |$2.02_{-0.88}^{+0.52}$| |
| Pdist, floor | |$0.10_{-0.03}^{+0.03}$| | |$0.07_{-0.03}^{+0.04}$| |
| Pdist, ceiling | |$0.49_{-0.15}^{+0.30}$| | |$0.79_{-0.20}^{+0.15}$| |
| ηmin, dist | |$0.11_{-0.06}^{+0.05}$| | |$0.24_{-0.19}^{+0.24}$| |
| ηdestr | |$0.25_{-0.03}^{+0.07}$| | |$1.88_{-0.53}^{+0.85}$| |
We also use Fig. 10 to show the results of the MCMC analysis by plotting the probability distribution of each parameter and showing contour plots for all possible parameter pairs. The parameters |$\textrm {Slope}_{P_\mathrm{r}}$|, rcore, Pdist, floor, and Pdist, ceiling cover a similar range of values in both theories. This is to be expected because the distribution of dwarfs in the Fornax Cluster is known observationally such that |$\textrm {Slope}_{P_\mathrm{r}}$| and rcore are not strong tests of the gravity law, while Pdist, floor and Pdist, ceiling are set by the proportion of dwarfs in different ηobs bins that appear disturbed (Fig. 3). Unlike these four parameters, |$\textrm {Slope}_{P_\mathrm{e}}$|, ηmin, dist, and ηdestr cover very different ranges in these two models. As discussed below, these are the parameters which can help us discern between ΛCDM and MOND, allowing us to assess which model performs better when compared to observations.
The 1σ confidence regions for the parameters in our model of the Fornax Cluster dwarf galaxy population using ΛCDM (orange) and MOND (blue), based on the priors listed in Table 1. The top panel in each column shows the inference on a single parameter, while the other panels show the 1σ confidence region for a pair of parameters. The results shown in this ‘triangle plot’ are based on 105 MCMC trials (Section 5.3.1). All the triangle plots shown in this contribution were generated using the pygtc package (Bocquet & Carter 2016).
The inference on |$\textrm {Slope}_{P_\mathrm{e}}$| (shown in the top panel of column 2 of Fig. 10) peaks close to the minimum allowed value of −2 in ΛCDM. The opposite happens in MOND, where the peak is close to 1. Negative slopes in equation (41) assign higher probabilities to nearly circular orbits. However, according to Ambartsumian (1937), we expect the eccentricity distribution to be thermal and thus have |$\textrm {Slope}_{P_\mathrm{e}} \approx 2$| (for a derivation, see section 4.2 of Kroupa 2008). In this regard, MOND performs better than ΛCDM.
The major differences between ΛCDM and MOND are in the parameters ηmin, dist and ηdestr, whose posterior inferences are shown in detail in Figs 11 and 12 due to their importance to our argument. The low values in ΛCDM arise because dwarfs have quite strong self-gravity by virtue of being embedded in a dominant dark matter halo throughout their trajectory. This makes them less susceptible to the effect of tides (stronger self-gravity raises rtid and thus reduces η, see equation 29). As a result, the algorithm needs to set ηmin, dist and ηdestr to very low values in order to match the observed fact that many dwarfs are morphologically disturbed and we do not observe dwarfs beyond a certain limiting η. MOND also boosts the baryonic self-gravity of a dwarf, but this boost is damped due to the EFE of the cluster’s gravitational field. This effect gets stronger as dwarfs approach the pericentre of their orbits, to the point that dwarfs which are sufficiently close to the cluster centre can become almost Newtonian despite a very low internal acceleration. Because of this, MONDian dwarfs are significantly more susceptible to tides than their ΛCDM counterparts. This causes the algorithm to choose significantly higher ηmin, dist and ηdestr values in the MOND case.
The probability distribution of the tidal susceptibility above which a dwarf is more likely to appear disturbed (Section 5.2.4). Notice that the MCMC analysis prefers significantly higher values for MOND (blue) than for ΛCDM (orange).
The probability distribution of the tidal susceptibility at which a dwarf is destroyed (Section 5.2.4) according to ΛCDM (orange) and MOND (blue).
N-body simulations of dwarf galaxies show that ηdestr should be ≈1 in ΛCDM (Peñarrubia et al. 2009; van den Bosch et al. 2018). However, fitting the observations with our MCMC method gives a much lower value of |$\eta _{\textrm {destr}} = 0.25^{+0.07}_{-0.03}$|. This implies an important discrepancy between model expectations in ΛCDM and actual observations of dwarf galaxies in the Fornax Cluster.
Turning to MOND, comparing the ηdestr value inferred from observations with that obtained using simulations is not so straightforward given that the best available N-body simulations studying the resilience of Milgromian dwarf galaxies to tides is by now very old and poorly suited to the present study (Brada & Milgrom 2000). Because of this, we perform our own N-body simulations of a typical Fornax Cluster dwarf galaxy, as described next.
7 N-BODY SIMULATIONS OF A FORNAX DWARF
For the dwarf, we use a half-mass radius of rh = 0.84 kpc and a total mass of Mdwarf = 3.16 × 107 M⊙ represented by 105 particles, making the mass resolution |$316 \, \mathrm{M}_{\odot }$|. These are typical parameters for a dwarf in the Fornax Cluster (see the red star in Fig. 8). Setting the velocity dispersion σ is non-trivial because we need to account for the cluster EFE when we initiate the simulation. We do this by using the Fornax dwarf templates kindly provided by Prof. Xufen Wu, who used a similar method to that described in section 3.3 of Haghi et al. (2019a) to generate these templates. The idea is to take a Newtonian template and then enhance the velocities by the factor needed to ensure virial equilibrium given the enhanced gravity (Wu & Kroupa 2013).
To set up the dwarf, we apply a Galilean transformation to the template whereby the Cartesian positions of all particles are boosted by (x0 = Ri, y0 = 0, z0 = 0) and the velocities are boosted depending on the circular velocity at Ri and the orbital eccentricity e, as described in Section 5.1. We start the simulation with the dwarf at the semimajor axis of its orbit and receding from the cluster. We then evolve the system until shortly after the dwarf reaches apocentre for the second time so that there is ample time to assess the impact of the pericentre passage. The code generates an output of the mass, position, and velocity of every particle every 20 Myr, allowing us to analyse the structure of the dwarf and find out if it has been destroyed.
Our main objective is to find the threshold value of η at pericentre beyond which the dwarf gets destroyed in the simulation. This requires us to perform multiple simulations with different eccentricities in order to obtain different η values at pericentre. To guide our choice of parameters, we use a simple MOND Runge–Kutta orbit integrator of a point mass orbited by a test particle in 2D. This is also very helpful when deciding the appropriate duration for each simulation, which we keep fixed for models with the same Ri.
7.1 Analysis
We extract the particle positions |$\boldsymbol {r}_i$|, velocities |$\boldsymbol {v}_i$|, and masses mi using extract_por (Nagesh et al. 2021), with the index i used in what follows to distinguish the particles. To assess if a dwarf has been destroyed, we infer three properties of the dwarf from the output at each snapshot: its half-mass radius, velocity dispersion, and aspect ratio. Unlike in Newtonian gravity, the time-varying EFE implies that these quantities are expected to vary around the orbit even if the dwarf is completely tidally stable (η ≪ 1), perhaps most famously for the velocity dispersion (Kroupa et al. 2018). To assess tidal stability, we check whether the dwarf responds adiabatically to the time-varying EFE. Tidal stability requires the dwarf to recover the initial values for these parameters after the pericentre passage, at least by the time of the next apocentre. If this is not the case, then the dwarf is either destroyed or unstable, in which case several pericentre passages may be required to destroy the dwarf. However, it is beyond the scope of this project to simulate multiple pericentre passages.
7.1.1 Finding the barycentre
7.1.2 Velocity dispersion
7.1.3 Half-mass radius
7.1.4 Aspect ratio
7.2 Results
The results of our por simulations are shown in Fig. 13. Unlike in the Newtonian case, even dwarfs with a very low tidal susceptibility exhibit significant variations in their properties due to the time-varying EFE. We can see that in the cases with low e, the dwarf manages to recover the properties it had before pericentre. However, in the cases with higher e, these properties do not regain their initial values, indicating that the dwarf is tidally unstable.9 This was expected because dwarfs with more eccentric orbits have closer pericentre passages and thus higher η values at pericentre.
Evolution of the half-mass radius (first column), 3D velocity dispersion (second column), and aspect ratio (third column) of the simulated dwarfs over time starting from an initial distance of Ri = 150 kpc (first row), Ri = 300 kpc (second row), and Ri = 450 kpc (third row). In each panel, the different curves show different orbital eccentricities as indicated in the legend, which gives their corresponding η values at pericentre (solid grey line) based on the EFE and tidal stress there but with the mass and half-mass radius found at the next apocentre (see the text). The mass at that time is similar to the initial value. The vertical dashed grey lines represent the first and second apocentre of the orbit. The solid (dashed) coloured lines represent those dwarfs which do (do not) recover their initial properties. The dotted lines that repeat one of the eccentricities in each panel correspond to a higher resolution simulation (8 × 105 particles), indicating that resolution hardly affects our results. The horizontal black line in the lower middle panel represents the expected velocity dispersion of the dwarf in the IDM limit (equation 58).
To assess whether a dwarf is destroyed in the simulation, the criterion that we apply is to consider destroyed those dwarfs which have a higher rh at the second apocentre than at pericentre. Since the dwarf is likely to expand even further as it heads towards its next pericentre, this implies that the dwarf has been too destabilized by tides to contract back to its size at its first pericentre passage. As a result, the dwarf would have an even higher tidal susceptibility at subsequent pericentres. This makes it very likely that the dwarf would not be able to survive multiple pericentre passages. On the other hand, if a dwarf that experiences a pericentre passage has a smaller rh at the subsequent apocentre and is contracting further, then it may well get back to its size at first pericentre by the time it reaches its second pericentre. This should allow it to survive multiple pericentre passages, which in the Fornax Cluster case should allow survival over a Hubble time.
To fairly compare our N-body results with our MCMC analysis, we should consider how observers calculate ηobs. The rh entering into equation (29) is the observed size, so ideally we would calculate η at pericentre using the EFE and tidal stress there but using the presently observed size. As a proxy for this, we use the size at apocentre since this is the orbital phase at which we are most likely to observe the dwarf. Physically, the tidal stability of a dwarf depends on the ratio between its size and tidal radius at pericentre. Using the ratio between the tidal radius at pericentre and the half-mass radius at apocentre may seem somewhat counter-intuitive. However, the ηdestr values obtained in this way are much more comparable to those obtained from our MCMC analysis of the Fornax Cluster for the reasons discussed above. In what follows, we will use η to mean the value calculated in this way, though Table 4 also shows results based on the size at pericentre.
Summary of our MOND N-body simulation results for a Fornax dwarf with an initial distance of Ri = 150 kpc and different orbital eccentricities (first column). The tidal susceptibility is calculated assuming the EFE and tidal stress at pericentre but using the half-mass radius of the dwarf at pericentre (second column) or at the subsequent apocentre (third column), which we argue in the text is more comparable to our MCMC results. The fourth column gives our assessment of the simulation based on the top-left panel of Fig. 13.
| η using rh at … | |||
|---|---|---|---|
| e | Pericentre | Apocentre | Outcome |
| 0.03 | 0.6 | 0.5 | Stable |
| 0.29 | 0.9 | 0.6 | Stable |
| 0.45 | 1.2 | 1.0 | Stable |
| 0.48 | 1.3 | 1.4 | Marginal |
| 0.51 | 1.4 | 2.0 | Unstable |
| 0.52 | 1.4 | 2.3 | Unstable |
| 0.53 | 1.4 | 2.7 | Destroyed |
| η using rh at … | |||
|---|---|---|---|
| e | Pericentre | Apocentre | Outcome |
| 0.03 | 0.6 | 0.5 | Stable |
| 0.29 | 0.9 | 0.6 | Stable |
| 0.45 | 1.2 | 1.0 | Stable |
| 0.48 | 1.3 | 1.4 | Marginal |
| 0.51 | 1.4 | 2.0 | Unstable |
| 0.52 | 1.4 | 2.3 | Unstable |
| 0.53 | 1.4 | 2.7 | Destroyed |
Summary of our MOND N-body simulation results for a Fornax dwarf with an initial distance of Ri = 150 kpc and different orbital eccentricities (first column). The tidal susceptibility is calculated assuming the EFE and tidal stress at pericentre but using the half-mass radius of the dwarf at pericentre (second column) or at the subsequent apocentre (third column), which we argue in the text is more comparable to our MCMC results. The fourth column gives our assessment of the simulation based on the top-left panel of Fig. 13.
| η using rh at … | |||
|---|---|---|---|
| e | Pericentre | Apocentre | Outcome |
| 0.03 | 0.6 | 0.5 | Stable |
| 0.29 | 0.9 | 0.6 | Stable |
| 0.45 | 1.2 | 1.0 | Stable |
| 0.48 | 1.3 | 1.4 | Marginal |
| 0.51 | 1.4 | 2.0 | Unstable |
| 0.52 | 1.4 | 2.3 | Unstable |
| 0.53 | 1.4 | 2.7 | Destroyed |
| η using rh at … | |||
|---|---|---|---|
| e | Pericentre | Apocentre | Outcome |
| 0.03 | 0.6 | 0.5 | Stable |
| 0.29 | 0.9 | 0.6 | Stable |
| 0.45 | 1.2 | 1.0 | Stable |
| 0.48 | 1.3 | 1.4 | Marginal |
| 0.51 | 1.4 | 2.0 | Unstable |
| 0.52 | 1.4 | 2.3 | Unstable |
| 0.53 | 1.4 | 2.7 | Destroyed |
To constrain ηdestr, we focus mainly on models with Ri = 150 kpc as dwarfs with a larger semimajor axis would typically be observed much further out than the region contributing to the apparent tidal edge in Fig. 3, especially if the eccentricity is significant. The results of these models are summarized in Table 4. The models with η ≤ 1.0 respond adiabatically. We choose ηdestr = 1.4 as the lowest value at which a dwarf can get destroyed in MOND since dwarfs with this η still seem to be marginally capable of contracting their rh back to their pericentre value by the time they reach apocentre.10 For the upper limit to ηdestr at pericentre, we choose a value of 2.0 because for this η, the dwarfs in our simulations are clearly larger at apocentre than at pericentre and are still expanding at the end of the simulation, indicating irreversible behaviour. We therefore infer that ηdestr = 1.70 ± 0.30 if rh is measured at the second apocentre. If instead we obtain rh at pericentre, then ηdestr has a slightly lower value of 1.35 ± 0.05.
As expected, ηdestr is of order unity because the main physics should be captured by analytic arguments (Zhao 2005; Zhao & Tian 2006). Our numerical results suggest that it would be more accurate to drop the factor of |$\frac{2}{3}$| in equation (20), which would also reconcile the numerical pre-factor with that in the Newtonian tidal radius formula (equation 11) for the case α = 1 and |$g \gg a_{_0}$|. This seems to indicate that we should identify the tidal radius with the distance to the L1 Lagrange point in the derivation of Zhao & Tian (2006) − their equation (36) introduces a factor of |$\frac{2}{3}$| in the Newtonian limit because the Roche Lobe extends to a shorter distance in the two non-radial directions than in the radial direction by about this factor. However, it could be that for somewhat eccentric orbits, the Roche Lobe’s extent along the radial direction is the limiting factor to the dwarf’s size.11
Our simulations also show that the higher the initial distance to the cluster, the more resilient the dwarf is to the effect of cluster tides. This is because a more eccentric orbit implies a shorter amount of time spent near pericentre, so the dwarf is exposed to a high η value for only a very brief period, allowing it to recover. Therefore, we would probably still be able to observe dwarfs which have η = 2.4 (or higher) at pericentre if these have sufficiently large apocentric distances. Given that in our analysis we considered dwarfs up to 800 kpc from the cluster centre, it is likely that there are several dwarfs in our sample which experienced a somewhat higher η at some point in their past − but for a sufficiently brief period that the dwarf remained intact. This is fairly consistent with the results of our MCMC analysis, which found that |$\eta _{\textrm {destr}} = 1.88_{-0.53}^{+0.85}$|.
The observed shape of a dwarf is one of the indicators for whether it has been perturbed. Therefore, to estimate the η at which simulated dwarfs should start appearing morphologically disturbed, we look at the evolution of their aspect ratio (equation 62). We need to bear in mind that even a uniform external field can cause a MONDian dwarf to become deformed because the potential of a point mass is not spherical once the EFE is considered (Banik & Zhao 2018a). N-body simulations of dwarfs experiencing the EFE but not tides explicitly show that this process can yield axis ratios of ≈0.7 (Wu et al. 2017). This is very much in line with our lowest eccentricity orbit with Ri = 150 kpc, so the mild degree of flattening evident here is not necessarily indicative of tidal effects. We find that models with Ri = 150 kpc start to acquire significantly elongated morphologies throughout most of their trajectories only when |$\eta \gtrsim 0.6$| (see column 3 in Fig. 13). Therefore, we take ηmin, dist ≈ 0.6. This is slightly higher than what our MCMC analysis requires (|$\eta _{\textrm {min,~dist}} = 0.24_{-0.19}^{+0.24}$|). One possible explanation is that dwarfs with higher Ri start acquiring elongated morphologies at lower η.
To check if increasing the resolution would affect our results, we perform a high-resolution rerun of one of our models for each Ri. This is shown using the dotted line in each panel of Fig. 13. The only resolution-related effect which we can observe is that the half-mass radius of a distant dwarf expands less than at lower resolution. Because of this, we obtain slightly lower pericentric η values for the same orbit with higher resolution. However, the evolution of the dwarf properties as a function of η at pericentre remains almost the same as for the low-resolution model. Therefore, our conclusions should barely be affected by the resolution of the simulation.
8 DISCUSSION
Observations of Fornax Cluster dwarf galaxies show that some of them present a detectable level of disturbance in their morphology. Among the environmental effects inside a galaxy cluster that could be causing this disturbance, we found that gravitational tides from the cluster are the most likely cause (Section 3). The condition for a dwarf galaxy in a galaxy cluster to be tidally stable is approximately the same as the requirement that the dwarf’s density exceed the average density of the cluster interior to the dwarf’s orbit (equation 7).12 This should be the case for a ΛCDM dwarf in a cluster because we expect the dwarf to be dominated by dark matter and to have formed much earlier than the cluster, at which time the cosmic mean density was higher. Therefore, in this paradigm, the dwarf galaxies in the Fornax Cluster should be little affected by the tides it raises. This is indeed what our calculations show (Fig. 2).
In MOND, the enhancement to the Newtonian gravity of an isolated dwarf is similar to that provided by the dark matter halo in ΛCDM. However, MONDian dwarfs in a galaxy cluster are also affected by the resulting EFE, which weakens their self-gravity. As a result, they are more susceptible to tides than dwarfs in ΛCDM, which has no EFE due to the strong equivalence principle. Therefore, observations of Fornax dwarfs can be used to compare which of the two models performs better.
To check if tides might be important in the Fornax Cluster, we plotted the projected separation (Rsky) of each FDS dwarf against a measure of its surface brightness (Fig. 3). This revealed a lack of low surface brightness dwarfs in the central ≈200 kpc even though such dwarfs are evident further out, indicating that selection effects are not responsible for the tentative tidal edge marked on this figure as a grey line. Just below this, the proportion of apparently disturbed dwarfs is also much higher than elsewhere in the cluster (see Fig. 4). We quantified this trend by plotting the disturbed fraction as a function of the tidal susceptibility η of each dwarf (equation 29), revealing a clear rising trend detected at 2.9σ significance in ΛCDM and 3.5σ in MOND (Fig. 5). These arguments suggest that the dwarf galaxy population in the FDS catalogue has been significantly shaped by tides, as previously argued by Venhola et al. (2022).
However, the overall distribution of η only goes up to ≈0.5 in ΛCDM (Fig. 2). We expect a dwarf to be destroyed or severely disturbed only if η ≈ 1, as indicated by ΛCDM N-body simulations (Peñarrubia et al. 2009; van den Bosch et al. 2018). We quantified this discrepancy using our MCMC analysis, which shows that the tidal stability limit of the Fornax dwarfs should be |$\eta _{\textrm {destr}} = 0.25^{+0.07}_{-0.03}$| to match observations. Therefore, ΛCDM dwarfs should be destroyed when the tidal force that they experience is ≈0.253 = 1.56 × 10−2 times smaller than their internal gravity (tidal force/internal gravity ≈ η3). Not only is this unrealistic, but also such a low ηdestr is in >5σ tension with the ηdestr value of 1 inferred from ΛCDM N-body simulations (Fig. 14). The highest ηdestr value achieved with our MCMC analysis for ΛCDM is only 0.60. This corresponds to the 4.42σ upper limit because we ran 105 MCMC trials. Since the uncertainty on ηdestr towards higher values from the mode is only 0.07, it is clear that ηdestr = 1 is strongly excluded by the observations if the tidal susceptibilities are calculated within the ΛCDM framework.
Joint inference on ηmin, dist and ηdestr (Section 5.2.4). We show the 1σ (inner solid line), 3σ (dashed line), and 5σ (outer dotted line) confidence region for MOND (blue) and ΛCDM (orange). The thick orange line shows the ΛCDM expectation that ηdestr ≈ 1. For MOND, the corresponding expectation from our N-body simulations (Section 7) is that ηdestr = 1.70 ± 0.30 (horizontal blue stripe). The grey shaded region below the line of equality is not allowed by our choice of prior because it is unphysical.
In MOND, we obtained a tidal stability limit with the MCMC analysis of |$\eta _{\textrm {destr}} = 1.88^{+0.85}_{-0.53}$|, which is closer to the expected value of ≈1 based on analytic arguments (equation 20). To check if this limit is accurate, we performed several N-body simulations of a dwarf orbiting a central potential similar to the Fornax Cluster (Section 7). These simulations suggest that cluster tides would make Fornax dwarfs appear disturbed when |$\eta _{\textrm {min,dist}} \gtrsim 0.6$| and destroy them at ηdestr = 1.70 ± 0.30, which is in good agreement with our MCMC results (see Fig. 14).
We considered several possible explanations for the discrepancy between the low tidal susceptibility values of ΛCDM dwarfs and the fact that some of the observed Fornax dwarfs appear disturbed. This could be due to the fact that cluster tides are not the main effect responsible for the observed morphological disturbances. However, there are several trends in the FDS that suggest exactly this. These trends are as follows:
There are fewer low surface brightness dwarfs towards the centre of the cluster, where they are most susceptible to tides (Fig. 3). Since such dwarfs are detectable further out, this feature cannot be ascribed to selection effects. A related finding is that FDS dwarfs are typically larger towards the cluster centre, which could be related to tidal heating (for a more detailed discussion, see section 7.4 of Venhola et al. 2022); and
The algorithm in charge of fitting the simulated Fornax system to the observations clearly noticed a rising trend between η and the probability of disturbance (Pdist). This is shown by the fact that the algorithm chose Pdist, ceiling > Pdist, floor with ≈3σ confidence in both ΛCDM and MOND (see Fig. 15), even though we did not impose this condition a priori.
We have seen that these trends cannot be understood in ΛCDM as a direct consequence of cluster tides given the very low η values. Moreover, the other major environmental effect that could be causing the observed disturbance (galaxy–galaxy harassment) also presents very low η values (see Section 4.1).
Joint inference on Pdist, floor and Pdist, ceiling (Section 5.2.4). We show the 1σ (inner solid line), 2σ (dashed line), and 3σ (outer dotted line) confidence region for MOND (blue) and ΛCDM (orange). Physically, we expect to get values above the solid red line of equality (Pdist, ceiling ≥ Pdist, floor), though this is not imposed as a prior. Even so, this is favoured by the MCMC analysis, which gives a likelihood that Pdist, ceiling ≤ Pdist, floor of only 3.14 × 10−3 (2.95σ) in MOND and 6.43 × 10−3 (2.73σ) in ΛCDM. Both theories prefer a non-zero false positive rate of Pdist, floor ≈ 0.1, which is related to the similar fraction of dwarfs classified as disturbed in the outskirts of the Fornax Cluster where tides should be unimportant (Fig. 4).
Another possibility is that our results could be affected by some of the assumptions or choices that we made during the analysis. To check if this is the case, we repeat the procedures described in Section 5 but change some of the assumed conditions and/or parameters in the following ways:
Considering that the FDS dwarfs could have a lower dark matter fraction within their optical radii: We consider the possibility that the dark matter fraction of the FDS dwarfs is lower than assumed in our nominal case (this is motivated in Section 8.1.1). Assuming that ΛCDM explains the properties of isolated dwarfs, we use the velocity dispersions of nearby isolated dwarfs to estimate their typical dark matter fraction, which returns a somewhat lower value than assumed in our nominal analysis. Substituting this fit (equation 65) into our MCMC chain raises ηdestr slightly, but it is still only |$0.33^{+0.04}_{-0.05}$|. We then consider a very conservative scenario in which there is only 10× as much dark matter as stars within the optical extent of each dwarf, which requires altering equation (10) to Mdwarf, ΛCDM = 11 M⋆. For this very low dark matter fraction, we obtain that |$\eta _{\textrm {destr}} = 0.54^{+0.19}_{-0.09}$|, which reduces the tension between observations and ηdestr = 1 (as expected from N-body simulations) to 2.29σ (the triangle plot for this analysis is shown in Fig. F1). While this is a significant improvement with respect to the >5σ tension in the nominal case, we see that even when considering one of the most conservative assumptions for the amount of dark matter contained within the optical radius of a dwarf, ηdestr ≥ 1 is still excluded at 97.8 per cent confidence. Moreover, we show in Section 8.1.1 that in a recent high-resolution cosmological ΛCDM simulation, the dark matter fraction within the stellar rh of a dwarf is far higher than this at the relevant M⋆, and is actually quite close to our nominal assumption;
Changing the lower limit to the distribution of dwarf densities in the test mass simulation: To check if the adopted detection limit to the density of the Fornax dwarfs significantly affects the results, we repeat the analysis using a density threshold ρt that is 5σ below the mean logarithmic density. We also consider a density limit of ρmean (grey line in Fig. 6). For reference, the nominal ρt in MOND is 2.88σ below the mean logarithmic density, while ρmean is 1.91σ below. The corresponding values in ΛCDM are 3.58σ and 2.56σ, respectively (see Appendix E). Fig. F2 shows the triangle plots comparing the results obtained using these two density limits with the nominal one for ΛCDM and MOND. From these plots (described further in Appendix F), we can see that choosing a lower ρt worsens the tension for ΛCDM while maintaining consistency in MOND. Using a higher ρt helps to increase the estimated values for ηdestr in ΛCDM. However, even if we use ρt = ρmean, the inferred ηdestr is still significantly below the threshold of ≈1 required in N-body simulations, while the inference on ηmin, dist hardly changes. Thus, choosing even higher ρt could perhaps help ΛCDM to reach a reasonable ηdestr. However, taking such high values for ρt would be in disagreement with observations as the whole point of ρt is that dwarfs are not detectable if they have a lower density, but dwarfs with a lower density are clearly observed if we adopt such a high ρt;
Changing the values of the deprojection parameters (see Appendix A): The deprojection parameters in our nominal analysis were offset = 0.4° and nnucfloor = 1.2° based on fig. 6 of Venhola et al. (2019). We repeat our analysis using deprojection values at the upper limit of the envelope in this figure: offset = 0.5° and nnucfloor = 1.5°. Fig. F3 shows the triangle plots comparing the results for these two different deprojections in ΛCDM and MOND. From these plots, we can see that these two deprojections give almost the same results in either theory;
Changing the ratio between present and pericentre distances (see Appendix B): A related change we could make is to consider altering the assumed ratio of 0.29 between the average R and the pericentre distance. This is valid for a thermal eccentricity distribution with |$\textrm {Slope}_{P_\mathrm{e}} = 2$|, which is expected theoretically but is the highest possible value (equation 41). With a lower |$\textrm {Slope}_{P_\mathrm{e}}$|, the ratio would rise as orbits would typically be more circular, reducing the calculated tidal susceptibility at pericentre. This would worsen the problem for ΛCDM; and
Increasing the resolution: In Section 5.1, we created a grid of 100 × 100 cells for different values of the orbital eccentricity (e) and initial distance to the cluster centre (Ri). We increase the resolution to 200 × 200 and repeat the analysis to check if this has any effect on the results. The triangle plots showing the results in ΛCDM and MOND for these two resolutions are shown in Fig. F4. From these plots, we can see that the results are nearly identical for the high- and low-resolution cases.
From these tests, we infer that our results are not significantly affected by modelling assumptions.
8.1 The dark matter content of dwarf galaxies in ΛCDM
Our conclusion that ΛCDM is inconsistent with the FDS dwarfs relies heavily on their low values of η in this paradigm, which in turn relies on the assumption that they should be dominated by dark matter. We therefore explore whether consistency could be gained by partially relaxing this assumption in a manner consistent with other constraints.
To try and raise η while continuing to use Newtonian gravity, we consider the possibility that the FDS dwarfs are TDGs. Our results are presented in Appendix D. We see that this scenario is also not viable because the elliptical galaxies in the cluster must still contain substantial dark matter haloes, leading to highly efficient disruption of dwarfs through galaxy–galaxy harassment.
It thus seems clear that the FDS dwarfs should be primordial. In this case, we may consider whether the dark matter density in their central regions could be substantially less than assumed here, raising their tidal susceptibility within the ΛCDM framework. The transformation of central cusps in the dark matter density profile into cores is expected to be rather inefficient for dwarfs with |$M_{\star } \lesssim 10^{7.2} \, \mathrm{M}_{\odot }$| (Di Cintio et al. 2014; Dutton et al. 2016; Tollet et al. 2016). Most FDS dwarfs have a lower M⋆ (i.e. they lie below the red line in Fig. 7). This makes it unlikely that baryonic feedback has substantially reduced the central dark matter density of most FDS dwarfs, especially at the low-mass and low surface brightness end important to our argument about tidal stability. Adiabatic contraction could actually raise the central dark matter density (Forouhar Moreno et al. 2022; Li et al. 2022), as could tidal stripping of the dark matter halo (Peñarrubia et al. 2008). The colours of the FDS dwarfs also indicate that star formation stopped early, most likely due to ram pressure stripping of the gas (Section 3). Thus, it would only be possible for strong feedback to substantially reduce the baryonic potential depth once. This is insufficient to substantially affect the central dark matter density even in the extreme case that the entire gas disc is instantaneously removed (Gnedin & Zhao 2002). Multiple bursts of star formation would be required to substantially affect the dominant dark matter halo (Pontzen & Governato 2012), but it is very unlikely that this occurred in most FDS dwarfs. Consequently, they should still have a significant amount of dark matter in their central regions, as is the case with Galactic satellites whose star formation ended early (Read, Walker & Steger 2019). Moreover, the low surface brightness nature of the FDS dwarfs considered here implies an atypically large size at fixed M⋆, causing the baryonic portion of the dwarf to enclose a larger amount of dark matter than for the more typical Illustris galaxies considered by Díaz-García et al. (2016).
Another way in which FDS dwarfs could lose dark matter is through interactions with a massive elliptical galaxy. This scenario has been shown to lead to a dwarf like DF2 with an unusually low dark matter content (Shin et al. 2020). However, such examples are rare in cosmological simulations (Haslbauer et al. 2019a; Moreno et al. 2022). In addition, the possibility that most FDS dwarfs lack dark matter altogether runs into severe difficulties based on simple analytic arguments: Newtonian TDGs would be very fragile and easily disrupted by interactions with massive cluster ellipticals, which must have substantial dark matter haloes in a ΛCDM context (Appendix D). MOND seems to offer the right level of tidal stability: neither too much such that all the dwarfs are completely shielded from tides and the observed signs of tidal disturbance remain unexplained, nor too little such that the dwarfs would have been destroyed by now in the harsh cluster environment studied here. The FDS dwarfs behave just as they ought to in MOND.
This conclusion is in agreement with the recent work of Keim et al. (2022), which used the observed tidal disturbance of the dwarf galaxies NGC 1052-DF2 and NGC 1052-DF4 to argue that they must be ‘dark matter free’, since otherwise their dark matter halo would have shielded them from tides. Phrased in a less model-dependent way, these observations indicate much weaker self-gravity than for a typical isolated dwarf, which is a clear prediction of MOND due to the EFE (Famaey et al. 2018; Kroupa et al. 2018; Haghi et al. 2019a). In the more isolated galaxy DF44, the self-gravity is stronger despite a similar baryonic content (van Dokkum et al. 2019), but this too is in line with MOND expectations (Bílek et al. 2019; Haghi et al. 2019b). Strong evidence for the EFE has also been reported from the outer rotation curves of spiral galaxies, which tend to be flat for isolated galaxies but have a declining trend for galaxies in a more crowded environment (Haghi et al. 2016; Chae et al. 2020, 2021).
Our results with the FDS are similar to those of Chilingarian et al. (2019) and Freundlich et al. (2022), who also report signs of tidal disturbance in some of the dwarf galaxies in the Coma cluster. Another case in point is the recent study of the dwarf galaxy population in the Hydra I cluster, where the proximity to the cluster centre seems to be affecting the morphology of the dwarfs in a manner suggestive of tidal effects (e.g. larger half-mass radii for dwarfs closer to the cluster centre; La Marca et al. 2022). Closer to home, the MW satellites also show signs of tidal disturbance like elliptical isophotes (McGaugh & Wolf 2010). There is a good correlation between these features and the value of η in MOND, which moreover has a maximum value very close to 1 (see their fig. 6). However, the maximum η in ΛCDM is |$\lesssim\!{0.2}$|, making it difficult to understand the observations in this framework.
8.1.1 Revised dark matter fraction in ΛCDM dwarfs
Throughout our analysis, we followed the Díaz-García et al. (2016) prescription that |$4{{\ \rm per\ cent}}$| of the total dark matter halo of each dwarf lies within its optical radius, with the total halo mass Mhalo following from M⋆ through the Moster et al. (2010) abundance matching relation. The factor of |$4{{\ \rm per\ cent}}$| was obtained by fitting to the dynamically inferred dark matter masses MDM within the optical radii of S4G galaxies, as shown in fig. 6 of Díaz-García et al. (2016). In this figure, we can see that for low-mass galaxies (|$M_{\star } \lesssim 10^9 \mathrm{M}_{\odot }$|), the MDM/M⋆ versus M⋆ relation seems to flatten at |$M_{\textrm {DM}} \approx 10 \, M_{\star }$|. However, this is unclear because S4G has very few well-observed galaxies with such a low mass.
To check the consistency between the assumed dark matter fraction and observations of isolated dwarfs, we use Fig. 16 to plot Mdyn/M⋆ of the galaxies in four different galaxy surveys (semi-transparent coloured dots), assuming the Díaz-García et al. (2016) result for the dark matter fraction as used in our nominal analysis (black line; equation 10), and assuming conservatively that |$M_{\textrm {DM}} = 10 \, M_{\star }$| (blue line). We can see that it is rather unlikely that the FDS dwarfs generally have much less dark matter in their baryonic region than we assumed, since the linear regression to the survey data over the M⋆ range of the FDS dwarfs (dashed black line) is quite close to our nominal dark matter fraction. We can also use the Illustris TNG50 cosmological simulation (Pillepich et al. 2018, 2019; Nelson et al. 2019a,b) to check the dark matter fraction that we expect dwarfs to have in the ΛCDM paradigm. We do this in Fig. 16, where we show the mean and standard deviation of MDM/M⋆ + 1 within the stellar rh in M⋆ bins of width 0.25 dex (cyan dots with error bars). The trend followed by these simulated dwarfs is even steeper than that given by the observed dwarfs, though both give a similar dark matter fraction at the low-mass end crucial to our analysis (the median M⋆ of the FDS dwarfs is shown by the vertical dashed green line at log10(M⋆/M⊙) = 6.96). This further supports our nominal choice for the dark matter fraction of FDS dwarfs. One reason for their high expected dark matter fraction is that the vast majority of them have too little stellar mass for efficient core formation, the threshold for which is shown by the red vertical line at log10(M⋆/M⊙) = 7.2 for the reasons discussed above. All these arguments highlight that the |$M_{\textrm {DM}} = 10 \, M_{\star }$| case is clearly very conservative given the steep relation followed by low-mass galaxies that we expect from abundance matching arguments, Illustris TNG50 results, and the velocity dispersions of nearby dwarfs.
The relation between stellar mass and Newtonian dynamical mass (equation 64). The semi-transparent coloured dots represent galaxies from four different galaxy surveys as indicated in the legend (Falcón-Barroso et al. 2011; McConnachie 2012; Ryś, van de Ven & Falcón-Barroso 2014; Toloba et al. 2014). The cyan dots with error bars represent the logarithmic mean and dispersion of the total/stellar mass ratio within the stellar rh of the dwarfs in each M⋆ bin in the ΛCDM cosmological simulation Illustris TNG50 (Pillepich et al. 2018, 2019; Nelson et al. 2019a,b). These bins have a width of 0.25 dex and cover the mass range log10(M⋆/M⊙) = 4.5–12. The solid black line represents the expected trend in ΛCDM with the nominal dark matter fraction from abundance matching (equation 10). The dashed black line represents the fit to the dwarfs from the aforementioned galaxy surveys in the mass range covered by the FDS dwarfs. The horizontal blue line at log1011 shows our conservative assumption that |$M_{\textrm {DM}} = 10 \, M_{\star }$|. In these three cases, the relations are only plotted over the M⋆ range of the FDS dwarfs. Their median mass is shown with the dashed green vertical line. The solid red vertical line corresponds to the stellar mass below which core formation is inefficient in ΛCDM (see the text).
For completeness, we repeat our analysis with the very conservative assumption that |$M_{\textrm {DM}} = 10 \, M_{\star }$|. In this case, the distribution of dwarf densities is similar to that in MOND (Fig. 6) but scaled up 11×. Thus, the logarithmic dispersion remains σ = 0.57 dex and the density threshold ρt = 4.66 × 10−4 M⊙ pc−3 is still 2.88σ below the mean |$\log _{10} \, \rho$|, which is now −1.69 in these units. As expected, the ρt value is 11× higher than in the MOND model − and thus much less than in our nominal ΛCDM analysis. We found that in this reduced density case, |$\eta _{\textrm {destr}} = 0.54^{+0.19}_{-0.09}$| and the probability that ηdestr ≥ 1 is 2.23 × 10−2 (2.29σ).
Appendix F shows the complete triangle plot with the distributions of the model parameters and parameter pairs for the nominal ΛCDM analysis and the two revised cases described above. There is little impact to the inferences on parameters other than ηdestr, ηmin, dist, and |$\textrm {Slope}_{P_\mathrm{e}}$|.
Therefore, it is clear that assuming a lower dark matter fraction for the ΛCDM dwarfs helps to alleviate the tension between observations and N-body simulations only if this fraction is reduced significantly. However, having a dark matter fraction of MDM/M⋆ = 10 within the optical radius is a very conservative assumption at odds with many other lines of evidence, including cosmological simulations. Even with this assumption, ηdestr ≥ 1 is still excluded by our MCMC analysis of the FDS at 97.8 per cent confidence.
9 CONCLUSIONS
We studied the tidal susceptibility of dwarf galaxies in the Fornax Cluster to gravitational effects of the cluster environment in both ΛCDM and MOND. In both theories, we found cluster tides to be the main effect. Thus, cluster tides should be able to explain the observed morphological disturbance of some Fornax dwarfs and the lack of low surface brightness dwarfs towards the cluster centre (Fig. 3). By constructing a test mass simulation of the Fornax system and performing a statistical analysis using the MCMC method, we constrained the tidal susceptibility (η ≡ rh/rtid) value at which a Fornax dwarf should get destroyed in order to match the observations, which we call ηdestr. We found that |$\eta _{\textrm {destr}} = 0.25_{-0.03}^{+0.07}$| in ΛCDM and |$1.88_{-0.53}^{+0.85}$| in MOND.
The ηdestr value in ΛCDM falls significantly below analytic expectations (equation 11) and is in >5σ tension with N-body simulation results, which indicate that ηdestr ≈ 1 (Peñarrubia et al. 2009; van den Bosch et al. 2018). In other words, the very low η values of FDS dwarfs imply that they should be unaffected by cluster tides, contradicting the observed signs of tidal disturbance. We also found that the other major environmental influence of interactions with individual massive galaxies in the cluster should not be a significant process in ΛCDM (see also section 7.3.3 of Venhola et al. 2019). We discarded the possibility that the above-mentioned discrepancy is due to the minimum allowed density of the simulated sample of dwarfs being too low, the deprojection parameters being different from our nominal ones, the resolution of the test mass simulation not being high enough to get reliable results, and the dwarfs having less dark matter than we assumed (Section 8). In particular, the velocity dispersions of nearby isolated dwarfs suggest a slightly lower dark matter fraction (dashed line in Fig. 16). Using this only slightly raises ηdestr to |$0.33^{+0.04}_{-0.05}$|. Even if we conservatively assume that the FDS dwarfs have only 10× as much dark matter as stars within their optical radius, we still get a 2.29σ tension with expectations (equation 11). Therefore, our results reliably show that the ΛCDM paradigm is in serious tension with observations of perturbed dwarf galaxies in the Fornax Cluster (observations which are strongly suggestive of tidal effects, see also section 7.4 of Venhola et al. 2022).
An alternative model that assumes different properties for the dark matter particles could perhaps reconcile the basics of the ΛCDM cosmology with the observed morphological disturbances of some Fornax dwarfs. One of the most popular alternatives is the ‘superfluid dark matter’ model (Berezhiani & Khoury 2015; Hossenfelder & Mistele 2020). Like most hybrid models, it attempts to reconcile the successes of MOND on galaxy scales with the advantages of dark matter on larger scales, especially with regards to the CMB anisotropies and galaxy cluster dynamics. However, this model also presents its own problems, including orbital decay of stars in the Galactic disc from Cherenkov radiation (Mistele 2021) and that the LG satellite planes extend beyond the estimated superfluid core radii of the MW and M31, making it difficult to explain the high observed internal velocity dispersions of the satellites in these planes (see section 5.6 of Roshan et al. 2021a). There are also difficulties explaining the observed regularities in rotation curves consistently with gravitational lensing results in a theory where baryons feel extra non-gravitational forces that do not affect photons (Mistele, McGaugh & Hossenfelder 2022). Another possibility is that the dark matter particles are fuzzy with a low mass and thus a long de Broglie wavelength, reducing their density in the central region of a dwarf galaxy. However, ultralight bosons (Hu, Barkana & Gruzinov 2000; Hui et al. 2017) are in significant tension with observations of the Lyman-α forest (Rogers & Peiris 2021). More generally, reducing the ability of dark matter to cluster on small scales would make it difficult to form dwarf galaxies at high redshift and to explain their high Newtonian dynamical M/L ratios.
This brings us to the MOND case, in which the inferred ηdestr is much more consistent with analytic expectations (equation 20). In order to compare ηdestr with the results of N-body simulations as we did for ΛCDM, we had to perform numerical MOND simulations ourselves (though one pioneering study exists, see Brada & Milgrom 2000). From our simulations tailored to the properties of a typical dwarf galaxy in the Fornax Cluster, we obtained that ηdestr = 1.70 ± 0.30, in excellent agreement with the value required to fit the observational data according to the MCMC method. We therefore conclude that MOND performs significantly better than ΛCDM and is clearly the preferred model in all the tests that we conducted throughout this work, even though it was not designed with the FDS in mind. Nevertheless, MOND still needs an additional ingredient to explain some of the observations on larger scales, especially the temperature and pressure profiles in galaxy clusters and the CMB power spectrum (Famaey & McGaugh 2012). For this, several models have been proposed that complement MOND. Some of the most promising ones are the relativistic MOND theory which can fit the speed of gravitational waves and the CMB anisotropies but likely cannot explain the dynamics of virialized galaxy clusters (Skordis & Złośnik 2021); and the νHDM model that assumes MOND gravity and 11 eV sterile neutrinos (Angus 2009). These proposed particles would play the role of a collisionless component that only aggregates at the scale of galaxy clusters, helping to explain the Bullet Cluster (Angus et al. 2007) and other virialized galaxy clusters (Angus et al. 2010), where the MOND corrections to Newtonian gravity are generally small. MOND has also proved capable of explaining several physical phenomena that ΛCDM has been failing to describe, including the planes of satellite galaxies in the LG and beyond (Pawlowski 2021a,b), the weakly barred morphology of M33 (Sellwood, Shen & Li 2019; Banik et al. 2020), and the pattern speeds of galaxy bars (Roshan et al. 2021a,b). Using the νHDM extension, MOND can also explain the CMB (Angus & Diaferio 2011), the KBC void and Hubble tension (Haslbauer et al. 2020), and the early formation of the interacting galaxy cluster El Gordo (Katz et al. 2013; Asencio, Banik & Kroupa 2021). Therefore, this later model is capable of explaining both the CMB and the dynamics of galaxy clusters while preserving the successes of MOND at galaxy scales (Banik & Zhao 2022, and references therein). In this study, we have shown that it should also be capable of resolving the problem faced by ΛCDM with regards to the observed signs of tidal disturbance in Fornax Cluster dwarf galaxies.
ACKNOWLEDGEMENTS
EA is supported by a stipend from the Stellar Populations and Dynamics Research Group at the University of Bonn. IB is supported by Science and Technology Facilities Council grant ST/V000861/1, which also partially supports HZ. IB acknowledges support from a ‘Pathways to Research’ fellowship from the University of Bonn. PK acknowledges support through the Deutscher Akademischer Austauschdienst-Eastern European Exchange Programme. EA would like to thank Prof. Xufen Wu for providing the initial conditions templates of the dwarf galaxy used in the MOND N-body simulations. The authors are grateful to Sara Eftekhari for providing the table of literature data shown in Fig. 16. They are also grateful to the referee for comments which substantially improved this paper.
DATA AVAILABILITY
The results presented can be reproduced by using the data available in the Vizier catalogue13 and following the methods described in this paper. For a user guide describing how to install por and providing links from which it can be downloaded, we refer the reader to Nagesh et al. (2021).
Footnotes
For an isolated dwarf, the dark matter halo in ΛCDM and the correction to Newtonian gravity in MOND both provide a similar enhancement to the self-gravity.
This is consistent with the previous ΛCDM result that harassment is not very significant for dwarfs in a Virgo-like cluster (Smith et al. 2015).
This is not expected if the disturbances are due to harassment because dwarfs subject to this would be well mixed throughout the cluster (Smith et al. 2015).
This is based on the binomial uncertainty in S assuming that the probability of a galaxy appearing disturbed in each Rsky bin is fd = S/T. In reality, fd is not precisely constrained by the observations − we handle this complexity later (equation 32).
Section 8 provides a more rigorous quantification of how confident we can be that fd rises with η.
The low values of Ptotal arise due to the large sample size.
This seems to be the case for the MW satellite Crater II (Torrealba et al. 2016), whose low surface brightness, small pericentre (Li et al. 2021), and low velocity dispersion for ΛCDM (Caldwell et al. 2017) suggest that it is the remnant of an originally smaller object that got severely disrupted by tides during its perigalacticon passage (Borukhovetskaya et al. 2022; Errani et al. 2022). It is also expected to be tidally unstable in MOND (see section 3.3 of Banik & Zhao 2022).
This certainly appears to be the case for the η = 1.5 model with Ri = 300 kpc.
Without the |$\frac{2}{3}$| factor in equation (20), the tidal susceptibility threshold is ηdestr = 1.13 ± 0.20 when using rh at apocentre and ηdestr = 0.90 ± 0.03 when using rh at pericentre. Note that the MOND tidal susceptibilities of FDS dwarfs would also be reduced by a factor of |$\frac{2}{3}$| in this case, which would affect the inferred ηdestr posterior.
The tidal stress Δgc/Δr is related to the cluster mass profile Mc(<R) by Δgc/Δr = GMc(2 − α)/R3, from which it follows that |$r_{\textrm {tid}}^3/M_{\textrm {dwarf}} \approx R^3/M_\mathrm{ c}$|. Thus, a dwarf with rh ≈ rtid has |$M_{\textrm {dwarf}}/r_\mathrm{ h}^3 \approx M_\mathrm{ c}/R^3$|.
Results are also shown for dwarf irregulars, but we removed these from our sample.
REFERENCES
APPENDIX A: DEPROJECTING DISTANCES IN THE SKY PLANE TO 3D DISTANCES
APPENDIX B: OBTAINING Rper FROM A 3D DISTANCE
APPENDIX C: DO TWO EXPERIMENTS HAVE THE SAME PROPORTION OF SUCCESSES?
In Section 4.3, we encountered the problem that one experiment gives Sobs, 1 ‘successes’ out of T1 trials while another experiment gives Sobs, 2 successes out of T2 trials, with a success defined as a dwarf galaxy that appears disturbed. The problem is to test the null hypothesis that the proportion of successes (x) is the same in both experiments assuming that T1 and T2 are set in advance independently of the actual number of successes. We consider this problem in following two stages as follows:
Keeping x fixed, we evaluate the likelihood Px of obtaining data as bad as or worse than the observed combination (Sobs, 1, Sobs, 2) for the null hypothesis; and
We then obtain a weighted mean value for Px by considering all plausible x, each time weighting by the likelihood that the observed (Sobs, 1, Sobs, 2) arises with that x.
If we know x, we can use binomial statistics (equation 48) to find the likelihood of obtaining any combination (S1, S2). We obtain Px by adding the probabilities of all (S1, S2) combinations which are as likely as or less likely than the observed combination (Sobs, 1, Sobs, 2). This follows the usual principle that if the data seems unlikely given the null hypothesis, we should consider all the ways in which it could look as bad or even worse.
In the particular case of Section 4.3, calculating the significance P in this way only tells us how plausible it is that fd is the same in the low η and high η subsamples, which is the null hypothesis. Our alternative hypothesis specifies that fd should be higher in the high η subsample on physical grounds, not merely that fd should have some sort of correlation with η. Since the inferred fd indeed rises with η, we should bear in mind that the low likelihood of the null hypothesis is caused by a deviation in just the sense expected on physical grounds under the alternative hypothesis where tides are relevant. On the other hand, we tried all possible choices of ηt to maximize the significance of the signal, leading to a look-elsewhere effect.
APPENDIX D: TIDAL SUSCEPTIBILITY OF NEWTONIAN TDGs
As discussed in Section 8, our results indicate a higher level of tidal susceptibility than is expected in ΛCDM. This could be a sign that the Fornax dwarfs lack dark matter altogether, which is possible in this framework if the FDS dwarfs are mostly TDGs. These are expected to be rather rare in ΛCDM, so the scenario is not very plausible (Haslbauer et al. 2019b). We none the less consider it for completeness.
If the dwarfs are of tidal origin, they would be free of dark matter (Barnes & Hernquist 1992; Wetzstein et al. 2007). However, the massive cluster galaxies would still be surrounded by a dark matter halo. In this scenario, the mass ratio between the dwarfs and the massive galaxies would be rather extreme, suggesting a serious problem with the stability of the dwarfs.
To quantify this, we obtain the tidal radius of a dwarf by applying equation (11) considering only its baryonic mass. Similarly, we can obtain the disruption time-scale by applying equation (13) and accounting for the fact that the terms referring to the dwarf (those labelled with a subindex ‘dwarf’) should be purely baryonic while the terms referring to the large galaxies (labelled with a subindex ‘p’) should still account for the dark matter contribution to the mass and half-mass radius. We can then substitute in these results to obtain the susceptibility to cluster tides (equation 29) and galaxy–galaxy harassment (equation 31). The results are shown in Fig. D1.
The distribution of tidal susceptibility values of Fornax Cluster dwarfs in a Newtonian TDG scenario to cluster tides (top panel) and harassment (bottom panel), with a bin width of 0.05 and 0.5, respectively.
As expected, the dwarfs are now much more susceptible to cluster tides (higher ηrtid than in Fig. 2). The distribution of ηrtid becomes very similar to MOND, suggesting that maybe the Newtonian TDG scenario is plausible. However, the tidal susceptibility to harassment (ηhar) is very large in this scenario and greatly exceeds 1 for the vast majority of the dwarfs. The high ηhar values arise because the dwarfs are completely unprotected: They do not have a boost to their self-gravity either from MOND or from a dark matter halo. Given their low surface brightness, this leads to very weak self-gravity. However, in a ΛCDM universe, the large galaxies must still have dark matter haloes. As a result, purely baryonic dwarfs governed by Newtonian gravity should have already been destroyed by encounters with the massive cluster galaxies. Therefore, we can consider that the TDG scenario in ΛCDM is extremely unlikely. Note that in MOND, our analysis is not sensitive to whether the dwarfs are TDGs or formed primordially − they are purely baryonic in either case.
APPENDIX E: DISTRIBUTION OF DWARF DENSITIES IN ΛCDM
Our MCMC analysis relies on an assumed distribution for the dwarf densities, which are crucial to their tidal stability. We therefore need to repeat the steps discussed in Section 5.2.3 for the case of ΛCDM. For this model, we show the mass–luminosity relation (Fig. E1), the surface density–volume density relation (Fig. E2), and the histogram of volume densities of the dwarfs in the FDS catalogue (Fig. E3). The main difference is that the mass of the dwarfs is higher since it includes the contribution of the dark matter component within the optical radius (equation 10). This raises their surface and volume density. We found that |$M/L_{r^{\prime }} = 74.92 \pm 52.38~\mathrm{ M}_{\odot }/\mathrm{ L}_{\odot , r^{\prime }}$|, indicating a rather high dispersion. Moreover, we can no longer approximate that the slope of the relation is 1 on logarithmic axes, indicating non-linearity.
Similar to Fig. 7, but for ΛCDM instead of MOND and showing only the linear regression.
Similar to Fig. 8, but for ΛCDM instead of MOND and showing only the linear regression.
Due to these difficulties, we found that it would be unsuitable to repeat the steps described in Section 5.2.3. To enable a fair comparison with MOND, we none the less used as similar a procedure as possible. For this, we fixed the logarithmic offset between the density of the least dense dwarf in our sample (ρmin, FDS) and the adopted density threshold of the survey (ρt). As a result, the minimum observational limit (black line in Fig. E3) is 0.09 dex below ρmin, FDS, the mean observational limit (grey line in this figure) is 0.46 dex above ρmin, FDS, and the maximum observational limit (dashed black line in this figure) is 0.79 dex above ρmin, FDS. As in the MOND case, we choose the minimum observational limit (black line in Fig. E3) as our nominal density limit for the distribution since it is the only one of these three choices that implies ρt < ρmin, FDS, which is required of a realistic detection threshold. Assuming instead the mean observational limit would make us lose two observed dwarf galaxies from the low-density tail of the distribution. Note also that these dwarfs have a clear tidal morphology because we removed any dwarfs where this is unclear (Section 2.1).
The distribution of each dwarf galaxy’s mean density within its half-mass radius, accounting for both baryonic and dark matter. The orange vertical line at −0.99 shows the sample mean, while the magenta lines offset by ±0.54 dex show the standard deviation around it. Other lines have a similar meaning to Fig. 6 and have been obtained similarly to the MOND case to allow a fair comparison (see the text).
To summarize, our nominal ρt in ΛCDM is 3.58σ below the mean logarithmic density, while ρmean is 2.56σ below.
APPENDIX F: TRIANGLE PLOTS WITH ALTERNATIVE MODELLING CHOICES
In this appendix, we rerun our MCMC analysis with different modelling assumptions and show their impact using triangle plots similar to Fig. 10. Instead of showing ΛCDM and MOND results on the same graph as done there, our approach will be that each graph shows results for different modelling assumptions but within the context of the same theory. We will use different panels for the different theories. As before, we show only the 1σ contour for each pair of parameters, though the full probability distribution is shown when considering the posterior on one parameter marginalized over all others. The results presented here are discussed in more detail in Section 8.
In Fig. F1, we check how decreasing the dark matter fraction within the optical radius of the FDS dwarfs affects the results. In particular, we consider the revised dark matter fraction given in equation (65) based on the observed velocity dispersions of nearby dwarfs (Section 8.1.1). As discussed there, we also consider the very conservative case |$M_{\textrm {DM}} = 10 \, M_{\star }$|. The main impact is on the parameters ηdestr and ηmin, dist. The inference on the slope of the eccentricity distribution is rather different for the case |$M_{\textrm {DM}} = 10 \, M_{\star }$|, but otherwise the posteriors are not much different to the nominal ΛCDM case in both revised analyses shown here.
Triangle plot showing the inferred parameter values in the ΛCDM model constrained by our statistical analysis assuming the nominal dark matter fraction (equation 10; blue), our fit to the empirically determined dark matter fractions of nearby isolated dwarfs (equation 65; orange), and a very conservative scenario in which the mass of dark matter within the optical radius of each dwarf is only 10× that of the baryons (green).
In Fig. F2, we compare the parameter inferences resulting from the MCMC analysis assuming three different lower limits (ρt values) to the density distribution of the dwarfs:
The lowest considered ρt is set at 5σ below the mean logarithmic density;
The second-lowest considered ρt is the nominal value used in the main analysis; and
The highest considered ρt is the mean observational limit (ρmean), which we obtained in Section 5.2.3 and Appendix E for MOND and ΛCDM, respectively.
In Fig. F3, we compare ΛCDM and MOND while assuming two different values for the deprojection parameters ‘offset’ and ‘non-nucleated floor’ (see Appendix A). In addition to the nominal values used in the main analysis, we also consider a higher set of values corresponding to the highest plausible 3D distance given the sky-projected distance (see fig. 6 of Venhola et al. 2019). This entails setting nnucfloor = 1.5° and offset = 0.5° instead of the nominal nnucfloor = 1.2° and offset = 0.4°.
Triangle plot similar to Fig. 10, but this time showing ΛCDM (left-hand panel) and MOND (right-hand panel) in separate panels. Each panel shows the difference in the final results when using three different lower limits to the density distribution of the simulated Fornax dwarfs. The results for the 5σ lower limit are shown in green, those with the nominal density limit are shown in blue, and results for the mean observational limit (ρmean) are shown in orange.
Similar to Fig. F2, but showing the difference in the final results for ΛCDM (left-hand panel) and MOND (right-hand panel) when two different sets of values for the parameters ‘offset’ and ‘non-nucleated floor’ are used to deproject distances (Appendix A). Results for the nominal deprojection (offset = 0.4°, non-nucleated floor = 1.2°) are shown in blue, while results with the revised deprojection parameters (offset = 0.5°, non-nucleated floor = 1.5°) are shown in orange.
In Fig. F4, we check if increasing the resolution of the orbital elements in the test mass simulation affects the results for ΛCDM and MOND. The nominal resolution used is a grid of size 100 × 100 for the eccentricity e and initial distance from the cluster centre (Ri), which also corresponds to the semimajor axis. In the higher resolution case, this is raised to 200 × 200.
Similar to Fig. F2, but showing the difference in the final results for ΛCDM (left-hand panel) and MOND (right-hand panel) with two different resolutions in orbital elements. The nominal resolution case (blue) uses a grid of 100 × 100 bins to generate orbits with different eccentricities and initial positions/semimajor axes. The high-resolution case (orange) uses a grid of 200 × 200 bins.
