pdftexdestination with the same hyperrefOption hyperrefToken pdftex(dest)

# The distribution and morphologies of Fornax Cluster dwarf galaxies suggest
they lack dark matter

Elena Asencio1, Indranil Banik2, Steffen Mieske3, Aku Venhola4, Pavel
Kroupa1,5 and Hongsheng Zhao2
1Helmholtz-Institut für Strahlen und Kernphysik (HISKP), University of Bonn,
Nussallee 14$-$16, D-53115 Bonn, Germany
2Scottish Universities Physics Alliance, University of Saint Andrews, North
Haugh, Saint Andrews, Fife, KY16 9SS, UK
3European Southern Observatory, Alonso de Cordova 3107, Vitacura, Santiago,
Chile
4Space Physics and Astronomy Research Unit, University of Oulu, Pentti
Kaiteran katu 1, 90014 Oulu, Finland
5Astronomical Institute, Faculty of Mathematics and Physics, Charles
University, V Holešovičkách 2, CZ-180 00 Praha 8, Czech Republic Email:
<EMAIL_ADDRESS>(Elena Asencio),
<EMAIL_ADDRESS>(Indranil Banik)

###### 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 ($\Lambda$CDM) standard model of cosmology and
Milgromian dynamics (MOND). We construct a test particle simulation of the
Fornax system. We then use the MCMC method to fit this to the FDS distribution
of tidal susceptibility $\eta$ (half-mass radius divided by theoretical tidal
radius), the fraction of dwarfs that visually appear disturbed as a function
of $\eta$, and the distribution of projected separation from the cluster
centre. This allows us to constrain the $\eta$ value at which dwarfs should
get destroyed by tides. Accounting for an $r^{\prime}$-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 $\Lambda$CDM and
$1.88^{+0.85}_{-0.53}$ in MOND. The $\Lambda$CDM value is in tension with
previous N-body dwarf galaxy simulations, which indicate that
$\eta_{\textrm{destr}}\approx 1$. Our MOND N-body simulations indicate that
$\eta_{\textrm{destr}}=1.70\pm 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 $\Lambda$CDM expectations but well
consistent with MOND.

###### keywords:

gravitation – dark matter – galaxies: clusters: individual: Fornax – galaxies:
dwarf – galaxies: interactions – galaxies: statistics

††pubyear: 2022††pagerange: The distribution and morphologies of Fornax
Cluster dwarf galaxies suggest they lack dark matter–24

##  1 Introduction

Dwarf galaxies are the smallest and most common type of galaxy. They are
characterized by their low mass ($M<10^{9}\,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 ($\Lambda$CDM)
cosmological model (Efstathiou et al., 1990; Ostriker & Steinhardt, 1995)
provides two different scenarios by which dwarf galaxies can form (the Dual
Dwarf Galaxy Theorem; Kroupa, 2012):

  1. 1.

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

  2. 2.

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 et al., 2007). In recent
years, cosmological $\Lambda$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_{\star}$). The lowest
$M_{\star}$ dwarfs tend to have velocity dispersions ($\sigma$) which are too
high if one assumes virial equilibrium, with $\sigma$ 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. $\Lambda$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 $\sigma$ 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
$\Lambda$CDM explanation for dSphs having such high $\sigma$ is to assume that
most of their mass is in the form of dark matter, in which case they must be
primordial dwarfs.

$\Lambda$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
$\Lambda$CDM model (Kroupa et al., 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;
Pawlowski & Tony Sohn, 2021; Müller et al., 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\sigma$ tension with $\Lambda$CDM (table 3
of Banik et al., 2021, and references therein), while the satellite plane
around Centaurus A is only 0.2% ($3.09\sigma$) 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 $\Lambda$CDM
using Gaussian statistics. Since we effectively have
$\chi^{2}=3.55^{2}+3.55^{2}+3.09^{2}=34.75$, the combined tension can be
estimated as the likelihood of the $\chi^{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\times 10^{-7}$
($5.27\sigma$). 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), $\la
10\%$ 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
$\sigma$ and in their observed ellipticity. More recently, Hammer et al.
(2020) pointed out that the high $\sigma$ 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 et al., 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 $\sigma$ for
their low $M_{\star}$ 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 $\sigma$ 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 $\sigma$ of all
observed dSphs, $\Lambda$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 et al., 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
$\Lambda$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\la a_{{}_{0}}=1.2\times
10^{-10}$ m/s${}^{2}=3.9$ pc/Myr2; Begeman et al., 1991; Gentile et al., 2011;
McGaugh et al., 2016). The gravity boost that dwarf galaxies experience in
this regime would explain their high $\sigma$ (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\pm 2$ Gyr ago, which is required in MOND (Zhao et
al., 2013) and seems to reproduce important aspects of their satellite planes
(Banik et al., 2018; Bílek et al., 2018, 2021; Banik et al., 2022a).
Therefore, we will focus mainly on $\Lambda$CDM and MOND in this contribution.

The planes of satellites problem is one of the most well-known challenges to
$\Lambda$CDM on galaxy scales (Kroupa et al., 2005; Pawlowski, 2018, 2021a,
2021b). 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 $\Lambda$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 $\Lambda$CDM, where they would be shielded
by their dark matter halo throughout their whole trajectory (Brada & Milgrom,
2000).111For an isolated dwarf, the dark matter halo in $\Lambda$CDM and the
correction to Newtonian gravity in MOND both provide a similar enhancement to
the self-gravity.

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 $\Lambda$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 $\Lambda$CDM or MOND, or if
neither model works well. $\Lambda$CDM could provide too much protection
against tides such that it under-predicts 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 $\Lambda$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 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 Fornax Deep Survey (FDS)

The Fornax Cluster is one of the nearest galaxy clusters
($d_{\textrm{Fornax}}=20.0\pm 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^{\circ}$ 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\,\textrm{deg}^{2}$ (Venhola et al., 2018). The FDS represents a
significant improvement in resolution and image depth with respect to the
previous spatially complete Fornax Cluster Catalogue (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%
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_{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_{e,r^{\prime}}=27.8\,\textrm{mag}\,\textrm{arcsec}^{-2}$ (Venhola et al.,
2018).

Figure 1: 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 column) and in the $r^{\prime}$
band with a filter enhancing the dwarf’s low surface brightness features
(right column). The dwarf in the first, second, and third row is classified as
“very disturbed”, “mildly disturbed”, and “undisturbed”, 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.

The FDS catalogue of dwarf galaxies (Venhola et al., 2017, 2018, 2019)
includes 564 dwarf galaxies with $2\times 10^{5}<M_{\star}/M_{\odot}<2\times
10^{9}$, 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 $t_{\textrm{Fornax}}=10\pm 1$ Gyr (Rakos et al.,
2001), where $t_{\textrm{Fornax}}$ 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_{\star}$, the effective radius, the right
ascension and declination, the apparent surface brightness in the $r^{\prime}$
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^{\prime}$-band
are obtained by fitting the data to a 2D Sérsic profile (Venhola et al., 2018)
using the galfit software (Peng et al., 2002). $M_{\star}$ is obtained from
the empirical relation between the $g^{\prime}-i^{\prime}$ 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:

  1. 1.

Undisturbed: Dwarf galaxies that do not present irregularities, distortions to
their shape, or tidal tails;

  2. 2.

Possibly/mildly disturbed: Hints of irregularities are present in the
outskirts of the dwarf galaxy;

  3. 3.

Very disturbed: Dwarf galaxies with tidal tails and/or very clear distortion
in the shape; and

  4. 4.

Unclear: Nearby bright objects or data artefacts make the classification
difficult.

Fig. 1 shows some illustrative examples of dwarfs in these categories.

###  2.1 Data selection

From the 564 FDS dwarfs, we remove those which are classified as late-type as
there is a high chance that these are not physically in the cluster but
instead represent line of sight contamination (Venhola et al., 2019). We also
remove dwarfs which have an ‘unclear’ tidal morphology because they are not
useful for the analysis. This leaves us with 456 dwarfs. We then obtain the
angular distance between each dwarf and the centre of the Fornax Cluster based
on the right ascension (RA) and declination (Dec) of the dwarf and that of the
Fornax Cluster, whose sky coordinates are
$\textrm{RA}_{\textrm{centre}}=54.6^{\circ}$,
$\textrm{Dec}_{\textrm{centre}}=-35.5^{\circ}$ (table D1 of Watson et al.,
2009).

$\displaystyle\\!\\!\\!\\!\\!\\!\Delta\textrm{RA}$ $\displaystyle\equiv$
$\displaystyle\textrm{RA}-\textrm{RA}_{\textrm{centre}}\,,$ (1)
$\displaystyle\\!\\!\\!\\!\\!\\!\Delta\textrm{Dec}$ $\displaystyle\equiv$
$\displaystyle\textrm{Dec}-\textrm{Dec}_{\textrm{centre}}\,,$ (2)
$\displaystyle\\!\\!\\!\\!\\!\\!\Delta^{\prime}\textrm{RA}$ $\displaystyle=$
$\displaystyle\Delta\textrm{RA}\cdot\cos\left(\frac{\textrm{Dec}+\textrm{Dec}_{\textrm{centre}}}{2}\right),$
(3) Angular distance $\displaystyle=$
$\displaystyle\sqrt{\left(\Delta^{\prime}\textrm{RA}\right)^{2}+\left(\Delta\textrm{Dec}\right)^{2}}\,.$
(4)

Expressing this angular distance in radians and multiplying it by the 20 Mpc
distance to Fornax (Blakeslee et al., 2009) then gives the dwarf’s sky-
projected distance $R_{\textrm{sky}}$ from the centre of the Fornax Cluster.

$\displaystyle
R_{\textrm{sky}}~{}=~{}d_{\textrm{Fornax}}\times\left(\textrm{Angular
distance}\right)\,.$ (5)

We remove dwarfs with $R_{\textrm{sky}}>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 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}<10^{8}\,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_{\star}$ down to a few times
$10^{7}\,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

In this type of interaction, the structure of a dwarf with mass
$M_{\textrm{dwarf}}$ is affected by gravitational tides coming from the
overall cluster potential, i.e., from the difference in the cluster gravity
across the finite size of the dwarf. We quantify the influence of cluster
tides on a dwarf using the concept of its tidal radius $r_{\textrm{tid}}$.
This is defined such that if $r_{\textrm{tid}}$ were the dwarf’s actual size,
then the tidal force of the cluster and the self-gravity of the dwarf would
have the same strength. We can intuitively see that

$\displaystyle\frac{GM_{\textrm{dwarf}}}{r^{2}_{\textrm{tid}}}~{}$
$\displaystyle\approx$
$\displaystyle~{}r_{\textrm{tid}}\overbrace{\left(\frac{\Delta g_{c}}{\Delta
R}\right)}^{\textrm{Tidal stress}}\,,$ (6) $\displaystyle\Rightarrow
r_{\textrm{tid}}~{}$ $\displaystyle\approx$
$\displaystyle~{}\left(\frac{GM_{\textrm{dwarf}}}{\Delta g_{c}/\Delta
R}\right)^{1/3}\,,$ (7)

where $G$ is the Newtonian constant of gravitation and $\Delta g_{c}/\Delta R$
is the tidal stress from the cluster potential, with $g_{c}$ and $R$ being the
cluster gravity and the 3D distance to the cluster centre, respectively. Since
we want to find out the maximum degree of disturbance that a dwarf can
experience due to the cluster potential, we obtain $g_{c}$ and its gradient
when the dwarf is at pericentre ($R=R_{\textrm{per}}$). In order to obtain
$R_{\textrm{per}}$ for each dwarf from its projected distance in the FDS, we
use $R_{\textrm{per}}=0.29\,R$ (see Appendix B), with $R$ obtained by
deprojecting $R_{\textrm{sky}}$ using the method described in Appendix A.

As in Venhola et al. (2019), we assume that the galaxy number density and
cluster potential have remained constant over time. This approximation is
reasonable because the orbital periods of galaxies in the Fornax Cluster are
typically much shorter than a Hubble time: The estimated 1D velocity
dispersion of 370 km/s (Drinkwater et al., 2001) combined with a maximum size
of 800 kpc (Section 2.1) implies a crossing time of only 1.2 Gyr. We assign
the cluster a Newtonian dynamical mass profile given by

$\displaystyle
M_{c}\left(<\theta_{\textrm{3D}}\right)~{}=~{}M_{\textrm{norm}}\left(\frac{\theta_{\textrm{3D}}}{\theta_{\textrm{norm}}}\right)^{\alpha}\,,$
(8)

where $\theta_{\textrm{3D}}\equiv R/d_{\textrm{Fornax}}$ is the 3D angular
distance to the Fornax Cluster centre. The parameters are:
$M_{\textrm{norm}}=3\times 10^{10}~{}M_{\odot}$,
$\theta_{\textrm{norm}}=10\arcsec$, and $\alpha=1.1$. This radial mass
dependency is obtained from fitting the above power-law to the mass profile
derived in fig. 17b of Paolillo et al. (2002), which uses the X-ray surface
brightness distribution of the central Fornax galaxy and its gas temperature
profile to find the gas density distribution. The mass profile is then derived
assuming hydrostatic equilibrium by applying the spherical Jeans equation.
Note that the mass derived here is a Newtonian dynamical mass. A more model-
independent way to describe the observations is in terms of the cluster
gravity $g_{c}\equiv GM_{c}/R^{2}$. This method of obtaining $g_{c}$ relies on
the well-understood physical process of thermal X-ray emission from hot gas.
Its temperature and density profile require a particular radial run of $g_{c}$
regardless of the gravity law. Therefore, it is not relevant whether $g_{c}$
has been enhanced by a dark matter halo or by MOND (or indeed by some elements
of both, as argued in Section 3.3). Consequently, $g_{c}$ will be the same in
the $\Lambda$CDM and MOND scenarios, as will the resulting tidal stress on
each dwarf.

This is not the case for $M_{\textrm{dwarf}}$. The FDS catalogue gives only
$M_{\star}$ for each dwarf. This can be equated with $M_{\textrm{dwarf}}$ in
MOND, but not in $\Lambda$CDM where each dwarf is expected to have a
substantial dark halo of mass $M_{\textrm{halo}}$. We find this using the same
abundance matching procedure as Venhola et al. (2019). We first find
$M_{\textrm{halo}}$ from the relation between $M_{\star}$ and
$M_{\textrm{halo}}$ given in equation 2 of Moster et al. (2010):

$\displaystyle\frac{M_{\star}}{M_{\textrm{halo}}}$ (9)
$\displaystyle=2\left(\frac{M_{\star}}{M_{\textrm{halo}}}\right)_{0}\left[\left(\frac{M_{\textrm{halo}}}{M_{1}}\right)^{-\beta}+\left(\frac{M_{\textrm{halo}}}{M_{1}}\right)^{-\gamma}\right]^{-1}\,.$

Their table 1 clarifies that the parameters in this equation are:
$\left(\frac{M_{\star}}{M_{\textrm{halo}}}\right)_{0}=0.0282$,
$M_{1}=10^{11.884}~{}M_{\odot}$, $\beta=1.057$, and $\gamma=0.556$. As the
dark halo of each dwarf is not observable and remains hypothetical, we are
only interested in whether tides are perturbing the dwarf’s stellar component
(which they might not be even if its dark matter halo is being stripped; see
Smith et al., 2016). For this, the Shell Theorem indicates that we only need
to consider the dark matter within the dwarf’s optical radius. Following
Venhola et al. (2019), we assume that this is only 4% of the total halo mass
$-$ Díaz-García et al. (2016) found this fraction to be consistent with the
dark matter masses within the optical radii of S4G galaxies (Sheth et al.,
2010). Adding the halo contribution to $M_{\star}$, the total mass of the
dwarf in $\Lambda$CDM for the purposes of our analysis is therefore:

$\displaystyle M_{\textrm{dwarf,
}\Lambda\textrm{CDM}}~{}=~{}M_{\star}+0.04\,M_{\textrm{halo}}\,.$ (10)

In Section 8.1.1, we consider other possible choices for the fraction of the
halo mass within the optical radius of a dwarf.

Equation 7 is only a very crude estimate for the tidal radius of a dwarf.
While it should capture the essential physics, we expect a more careful
treatment to yield an additional factor of order unity. Numerical simulations
are required to capture the details of mass loss from a dwarf undergoing tidal
disruption, which is expected to substantially distort its shape. To account
for this, the $\Lambda$CDM expression for $r_{\textrm{tid}}$ in equation 1 of
Baumgardt et al. (2010) includes an extra factor of $2^{-1/3}$. Taking this
into consideration, we adopt the following expression for $r_{\textrm{tid}}$
in $\Lambda$CDM:

$\displaystyle r_{\textrm{tid,
}\Lambda\textrm{CDM}}~{}=~{}\left(\frac{G\,M_{\textrm{dwarf,
}\Lambda\textrm{CDM}}}{2\Delta g_{c}/\Delta R}\right)^{1/3}\,.$ (11)

This is based on using their study to obtain the numerical pre-factor in
Equation 7 for circular orbits in a central potential with a flat rotation
curve ($\alpha=1$) $-$ other approaches are discussed below Equation 63.
Notice that $g_{c}$ itself does not directly affect the tidal radius: The
cluster gravity only affects the dwarf through the tidal stress it creates on
the dwarf. This is not so in the corresponding expression for MOND (Equation
20), which we derive in Section 3.3.

###  3.2 Galaxy-galaxy harassment

The morphology of the Fornax Cluster dwarf galaxies can also be disrupted by
gravitational interactions with individual large galaxies in the cluster. This
effect is called harassment (Venhola et al., 2019). Assuming a high relative
velocity between the dwarf galaxy and the larger galaxy, we can use the
impulse approximation to estimate the impact of each encounter on the internal
structure of the dwarf. We then need to combine the effects of many such
interactions, each time adding the squares of the velocity perturbations as
these would generally be in random directions, leading to a process resembling
a diffusive random walk. Equivalently, we should add the energy gained by the
dwarf from each encounter, leading to the concept of a heating rate $\dot{E}$
(equation 8.52 of Binney & Tremaine, 2008). The disruption time-scale
$t_{d,\Lambda\textrm{CDM}}$ is the time-scale over which putting energy into
the dwarf at the presently calculated $\dot{E}$ would cause it to become
unbound given its present gravitational binding energy per unit dwarf mass of

$\displaystyle\lvert E\rvert~{}=~{}\frac{GM_{\textrm{dwarf,
}\Lambda\textrm{CDM}}}{2r_{h,\textrm{dwarf}}}\,,$ (12)

where $r_{h,\textrm{dwarf}}$ is the half-mass radius of the dwarf. Since only
the baryons are visible, we again restrict our attention to the baryonic
component of each dwarf, so $r_{h,\textrm{dwarf}}$ refers to only its visible
component and $M_{\textrm{dwarf, }\Lambda\textrm{CDM}}$ is again found using
Equation 10. Dividing the magnitude of the binding energy by the heating rate
gives the disruption time-scale (equation 8.54 of Binney & Tremaine, 2008):

$\displaystyle t_{d,\Lambda\textrm{CDM}}\equiv\frac{\lvert
E\rvert}{\dot{E}}=\frac{0.043}{W_{p}}\frac{\sqrt{2}\sigma
M_{\textrm{dwarf,}\Lambda\textrm{CDM}}r_{h,p,{\Lambda\textrm{CDM}}}^{2}}{GM_{p,\Lambda\textrm{CDM}}^{2}n_{p}r_{h,\textrm{dwarf}}^{3}}\,.$
(13)

The ‘p’ subscript denotes the massive galaxy (perturber), while ‘dwarf’ refers
to the dwarf galaxy that is being perturbed. $W_{p}$ is a factor accounting
for the shape of the perturber galaxy’s mass distribution. We choose $W_{p}=1$
as an intermediate value between that of the Plummer and Hernquist models
(chapter 8.2 of Binney & Tremaine, 2008). $n_{p}$ is the number density of
perturbers, which Venhola et al. (2019) estimated to be
$25\,\textrm{Mpc}^{-3}$ by counting 48 large galaxies inside the virial volume
of the Fornax Cluster ($R_{\textrm{vir}}=0.77$ Mpc). Its 1D velocity
dispersion is $\sigma=370$ km/s (Drinkwater et al., 2001), with the extra
factor of $\sqrt{2}$ accounting for the fact that we need to consider the
dwarf-perturber relative velocity. $M_{p,\Lambda\textrm{CDM}}$ and
$r_{h,p,{\Lambda\textrm{CDM}}}$ are the perturber galaxy’s mass and half-mass
radius, respectively. Note that we use $r_{h}$ for the deprojected half-mass
radius of the baryonic component. $r_{h}$ does not include the dark matter
halo unless we explicitly say so and label it accordingly as
$r_{h,\Lambda\textrm{CDM}}$. Venhola et al. (2019) use $r_{h}$ for the radius
containing half of the total mass including dark matter, so our notation is
different in this respect.

To obtain $r_{h,\textrm{dwarf}}$ from the projected effective radius $r_{e}$
containing half of the dwarf’s total stellar mass, we use equation B3 of Wolf
et al. (2010), though a good approximation is that
$r_{h,\textrm{dwarf}}\approx\left(4/3\right)r_{e}$. Our adopted
$M_{p,*}=10^{10}~{}M_{\odot}$ is the median stellar mass of the large galaxies
catalogued in table C1 of Iodice et al. (2019) and in the FCC. In the
$\Lambda$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
$M_{p,\Lambda\textrm{CDM}}=10^{11.6}~{}M_{\odot}$ following this procedure,
which we also verified.

Using a single $M_{p}$ 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 $t_{d}$ have a larger impact on
the result than taking into account the right distribution of perturbing
galaxy masses. Fortunately, we will see that $t_{d}$ 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 $t_{d}$ remains very long (Section 4.1).

The $r_{h,p}$ value of the large galaxies is also obtained from the median of
all the documented large galaxies (perturbers) in the cluster, yielding
$r_{h,p}=4$ kpc based on the luminous matter. This is applicable to MOND, but
in the $\Lambda$CDM case, the $r_{h,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 $r_{h,p}$ and $r_{h,p,\Lambda\textrm{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 $r_{h,p,\Lambda\textrm{CDM}}/r_{h,p}=3.6$. Therefore, the half-mass
radius of the perturbers in $\Lambda$CDM is taken to be
$r_{h,p,\Lambda\textrm{CDM}}=14.4$ kpc.

To summarize, the disruption time-scale in $\Lambda$CDM can be found by
directly applying Equation 13 once we include the contribution of the dark
matter halo to $M_{\textrm{dwarf}}$, $M_{p}$, and $r_{h,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_{{}_{N}}$ according to $g=\sqrt{a_{{}_{0}}g_{{}_{N}}}$.
Milgrom’s constant $a_{{}_{0}}=1.2\times 10^{-10}$ m/s2 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 et al., 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).

For an isolated spherically symmetric problem, the expression for the MOND
gravitational field $g$ as a function of the Newtonian field $g_{{}_{N}}$ can
be written as

$\displaystyle g~{}=~{}g_{{}_{N}}\nu\left(g_{{}_{N}}\right)\,,$ (14)

where $\nu$ is the interpolating function with argument $g_{{}_{N}}$. To
satisfy Solar system constraints and the observed flat rotation curves in the
outskirts of galaxies, this function must have the following asymptotic
limits:

$\displaystyle\nu\to\begin{cases}1\,,&\textrm{if}~{}g_{{}_{N}}\gg
a_{{}_{0}}\,,\\\
\sqrt{\frac{a_{{}_{0}}}{g_{{}_{N}}}}\,,&\textrm{if}~{}g_{{}_{N}}\ll
a_{{}_{0}}\,.\end{cases}$ (15)

The first case is the Newtonian regime in which $\nu=1$ and $g=g_{{}_{N}}$ to
a very good approximation. In the MOND regime,
$g=\sqrt{a_{{}_{0}}g_{{}_{N}}}$. This causes the gravity from an isolated
point mass $M$ to decline as $1/r$ beyond its MOND radius
$r_{{}_{\textrm{MOND}}}\equiv\sqrt{GM/a_{{}_{0}}}$, which is necessary to
explain the rotation curve data using only the luminous matter. Several forms
of the MOND interpolating function have been proposed (Kent, 1987; Hees et
al., 2014, 2016; McGaugh et al., 2016). Among these, the simple interpolating
function (Famaey & Binney, 2005) seems to work better with recent observations
(Iocco et al., 2015; Banik & Zhao, 2018c; Chae et al., 2018). Therefore, we
will use the simple interpolating function:

$\displaystyle\nu\left(g_{{}_{N}}\right)~{}=~{}\frac{1}{2}+\sqrt{\frac{1}{4}+\frac{a_{{}_{0}}}{g_{{}_{N}}}}\,.$
(16)

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 et al., 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 ($\nu$HDM) cosmological model (Angus,
2009). $\nu$HDM can fit observations of virialized galaxy clusters using the
MOND gravity of their directly detected baryons plus the sterile neutrinos
(Angus et al., 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 $\Lambda$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
$\la 50$, before which the MOND corrections to General Relativity should be
small (for a more detailed explanation of this model, see Haslbauer et al.,
2020). $\nu$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 $\nu$, 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 (van Dokkum et al., 2018; Famaey et al.,
2018; Kroupa 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 (van Dokkum et al., 2019; Bílek et al., 2019;
Haghi et al., 2019b).222In a conventional gravity context, the very low
observed velocity dispersion of NGC 1052-DF2 implies a lack of dark matter,
which however is not easily explained in $\Lambda$CDM (Haslbauer et al.,
2019a; Moreno et al., 2022). 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).

The EFE is also important to Fornax Cluster dwarfs because their low surface
brightness implies rather little self-gravity, allowing the gravitational
field of the cluster to dominate over that of the dwarf. As a result, the
dwarf is in the quasi-Newtonian (QN) regime where its internal dynamics are
similar to a Newtonian dwarf but with a renormalized gravitational constant
$G_{\textrm{eff}}>G$. We need to determine $G_{\textrm{eff}}$ from the cluster
gravitational field $g_{c}$. We do this by writing Equation 16 in the inverse
form:

$\displaystyle g_{{}_{N}}~{}$ $\displaystyle=$
$\displaystyle~{}g\mu\left(g\right)\,,\text{ where}$ (17)
$\displaystyle\mu\left(g\right)~{}$ $\displaystyle=$
$\displaystyle~{}\frac{g}{g+a_{{}_{0}}}\,.$ (18)

As the cluster gravity is dominant over the self-gravity of the dwarf, we can
set $g=g_{c}$, with $g_{c}$ obtained from observations as described in Section
3.1. Since the Newtonian gravity of the cluster is directly proportional to
the Newtonian gravitational constant ($g_{c,N}\propto G$), the effective
gravity of a dwarf in the cluster will be directly proportional to an
analogous constant parameter $G_{\textrm{eff}}$ defined such that
$g_{c}=\left(G_{\textrm{eff}}/G\right)g_{c,N}$. From Equation 18, we infer
$G_{\textrm{eff}}$ to be:

$\displaystyle
G_{\textrm{eff}}~{}=~{}\left(\frac{a_{{}_{0}}+g_{c}}{g_{c}}\right)G\,.$ (19)

Note that replacing $G\to G_{\textrm{eff}}$ can only be applied if the dwarf’s
self-gravity is dominated by the external field of the cluster, so that the
combined gravitational field of the dwarf and cluster will remain
approximately constant with increasing distance with respect to the dwarf’s
centre.

####  3.3.1 Tidal radius

At the tidal radius of a dwarf, the difference in cluster gravity across the
dwarf is comparable to its self-gravity. Therefore, the total cluster gravity
$g_{c}$ dominates over the dwarf’s self-gravity. Thus, the MOND tidal radius
of any dwarf is necessarily in the EFE-dominated/QN regime where its dynamics
are approximately Newtonian but with $G\to G_{\textrm{eff}}$. Substituting
this into Equation 7 gives an approximate expression for the MOND tidal
radius. Accounting for additional details like the non-spherical nature of the
point mass potential in the QN regime (discussed further in section 2.4 of
Banik & Zhao, 2022), the MOND tidal radius can be expressed as (equations 26
and 36 of Zhao & Tian, 2006):

$\displaystyle r_{\textrm{tid, MOND}}$ (20)
$\displaystyle=\frac{2}{3}\left.\sqrt{\frac{\partial\ln g}{\partial\ln
g_{{}_{N}}}}\right|_{g=g_{c}}\left[\left(\frac{2-\alpha}{3-\alpha}\right)\frac{G_{\textrm{eff}}\,M_{\textrm{dwarf}}}{\Delta
g_{c}/\Delta R}\right]^{1/3},$

where the factor of order unity is the MOND Roche lobe scaling factor
accounting for such subtleties. Note that we have generalized their equation
26 to write the result in terms of $G_{\textrm{eff}}$ and the tidal stress.
The parameter $\alpha\equiv 2+\frac{d\ln g_{c}}{d\ln r}$ has the same meaning
as in Equation 8, so its value remains 1.1. For the case of a dwarf orbiting a
point mass in the deep-MOND limit ($\alpha=1$), the numerical factors combine
to give $2^{1/6}/3$, matching equation 44 of Zhao (2005).333Equation 20 is the
extent of the Roche Lobe in the tangential direction within the orbital plane.
The extent along the orbital pole is similar, and in both cases is smaller
than the extent along the radial direction (see section 4.2 of Zhao & Tian,
2006).

####  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 $\bm{u}$ from the elliptical, with the heating rate
of the dwarf proportional to the square of $\lvert\Delta\bm{u}\rvert$, the
spread in $\bm{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\to G_{\textrm{eff}}$ 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:

  1. 1.

The QN regime in which $g_{c}\ll a_{{}_{0}}$ dominates over gravity from the
elliptical; and

  2. 2.

The isolated deep-MOND (IDM) regime in which the gravity from the elliptical
dominates over $g_{c}$ 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).

The QN and IDM regimes are separated by encounters with
$b=r_{{}_{\textrm{EFE}}}$, the distance from the elliptical beyond which the
cluster gravity dominates.

$\displaystyle r_{{}_{\textrm{EFE}}}~{}=~{}\sqrt{\frac{GM_{p}}{g_{c,N}}}\,,$
(21)

where $g_{c,N}\equiv g_{c}\mu\left(g_{c}\right)$ is the Newtonian gravity of
the cluster at the location of the dwarf-elliptical encounter. These
encounters would generally not occur when the dwarf is at the pericentre of
its orbit around the Fornax Cluster. However, encounters at this point would
be more damaging because the dwarf’s self-gravity would be weaker. We
therefore assume that the encounters with ellipticals take place at a typical
distance from the cluster of $R_{\textrm{enc}}=0.5\,R$, which is slightly more
than the pericentre distance of $0.29\,R$ (Appendix B) but less than the
present distance.

We will first consider the heating rate $\dot{E}_{\textrm{{QN}}}$ from
encounters in the QN regime before turning to the heating rate
$\dot{E}_{\textrm{{IDM}}}$ from encounters in the IDM regime. The total
heating rate is then

$\displaystyle\dot{E}_{\textrm{MOND}}~{}\equiv~{}\dot{E}_{\textrm{Newt}}\times\textrm{CF}~{}=~{}\dot{E}_{\textrm{{QN}}}+\dot{E}_{\textrm{IDM}}\,,$
(22)

where CF is the correction factor that needs to be applied to the Newtonian
$\dot{E}$ to make it MONDian. Our approach is to assume a sharp transition
between the QN and IDM regimes such that the EFE is completely dominant in the
former and completely negligible in the latter. This approximate approach
should be accurate to within a factor of order unity, which we will argue
later is sufficient for our purposes.

In all regimes, the heating rate due to encounters with an impact parameter in
the range $b\pm db/2$ ($db\ll b$) is
$\dot{E}_{b}=\dot{C}\,\langle\Delta\widetilde{E}\rangle$, where
$\dot{C}\propto b\,db$ is the average rate of such encounters and
$\langle\Delta\widetilde{E}\rangle$ is the average energy gain of the dwarf
per unit mass due to each such encounter. Since accelerating the dwarf as a
whole does not alter its internal structure, we only need to consider the
variation in the impulse $\bm{u}$ across the dwarf, so
$\langle\Delta\widetilde{E}\rangle\propto{\lvert\Delta\bm{u}\rvert}^{2}$. In
Newtonian dynamics, the magnitude of the impulse on a passing test particle is
$u\propto 1/b$, so $\lvert\Delta\bm{u}\rvert\propto 1/b^{2}$ (equation 8.41 of
Binney & Tremaine, 2008) and $\langle\Delta\widetilde{E}\rangle\propto
1/b^{4}$. This explains the $1/b^{3}$ scaling in the integrand in equation
8.53 of Binney & Tremaine (2008), which states that the Newtonian heating rate
per unit dwarf mass is:

$\displaystyle\dot{E}_{\textrm{Newt}}~{}=~{}\underbrace{\frac{14}{3}\sqrt{2\mathrm{\pi}}\frac{G^{2}M^{2}_{p}n_{p}r^{2}_{h,\textrm{dwarf}}}{\sqrt{2}\sigma}}_{A}\int^{\infty}_{r_{h,p}}\frac{db}{b^{3}}~{}=~{}\frac{A}{2r^{2}_{h,p}}\,,$
(23)

where $A$ is a constant.

We are now in a position to MONDify this result for the QN regime. Both the
dwarf’s self-gravity and the elliptical’s gravity on the dwarf are similar to
the Newtonian result but with $G\to G_{\textrm{eff}}$. The heating rate in the
QN regime is thus similar to Equation 23, but using $G_{\textrm{eff}}$ instead
of $G$ in the calculation of the normalization constant. To distinguish this
result from the Newtonian case, we call the QN normalization constant
$A^{\prime}=A\left(G_{\textrm{eff}}/G\right)^{2}$. Since by definition the QN
regime involves only those encounters with $b>r_{{}_{\textrm{EFE}}}$, the
total heating rate from encounters in this regime is

$\displaystyle\dot{E}_{\textrm{QN}}~{}=~{}A^{\prime}\int^{\infty}_{r_{{}_{\textrm{EFE}}}}\frac{db}{b^{3}}~{}=~{}\frac{A^{\prime}}{2r^{2}_{\textrm{EFE}}}\,.$
(24)

In the IDM regime, the scalings are different because the gravity from the
elliptical follows an inverse distance law. Since the interaction time-scale
rises linearly with the closest approach distance, the impulse becomes
independent of this ($u\propto b^{0}$). However, as the direction from the
elliptical to the dwarf is still different for different parts of the dwarf,
the variation in the impulse across it scales as
$\lvert\Delta\bm{u}\rvert\propto 1/b$, implying that the energy gain per
encounter scales as
$\langle\Delta\widetilde{E}\rangle\propto{\lvert\Delta\bm{u}\rvert}^{2}\propto
1/b^{2}$. Since the encounter rate again behaves as $\dot{C}\propto b\,db$ due
to the geometry being the same in both models, we obtain:

$\displaystyle\dot{E}_{\textrm{IDM}}~{}$ $\displaystyle=$
$\displaystyle~{}\frac{A^{\prime}}{r^{2}_{\textrm{EFE}}}\int^{r_{{}_{\textrm{EFE}}}}_{r_{{}_{\textrm{MOND}}}}\frac{db}{b}$
(25) $\displaystyle=$
$\displaystyle~{}\frac{A^{\prime}}{r^{2}_{\textrm{EFE}}}\ln\left(\frac{r_{{}_{\textrm{EFE}}}}{r_{{}_{\textrm{MOND}}}}\right)\,.$
(26)

The normalization of the integrand ensures continuity of the specific heating
rate per unit $b$ between the QN and IDM regimes.

Inserting our results for $\dot{E}_{\textrm{QN}}$ and $\dot{E}_{\textrm{IDM}}$
into Equation 22 and noting that
$A^{\prime}=A\left(G_{\textrm{eff}}/G\right)^{2}$, we obtain that

$\displaystyle\textrm{CF}~{}=~{}\left[1+\ln\left(\frac{a_{{}_{0}}}{g_{c,N}}\right)\right]\left(\frac{G_{\textrm{eff}}\,r_{\textrm{h,p}}}{G\,r_{{}_{\textrm{EFE}}}}\right)^{2}\,.$
(27)

Since $t_{d}\equiv\lvert E\rvert/\dot{E}$ and the MONDian binding energy of
the dwarf exceeds the Newtonian result (Equation 12) by a factor of
($G_{\textrm{eff}}/G$), the effect of the MOND corrections to Newtonian
gravity amount to multiplying the Newtonian $t_{d}$ (Equation 13) by a factor
of $\textrm{CF}^{-1}\left(G_{\textrm{eff}}/G\right)$.

$\displaystyle t_{d,\textrm{MOND}}~{}$ $\displaystyle\equiv$
$\displaystyle~{}\frac{\lvert
E\rvert_{\textrm{MOND}}}{\dot{E}_{\textrm{MOND}}}$ $\displaystyle=$
$\displaystyle\frac{0.043}{W_{p}}\frac{\sqrt{2}\sigma\,M_{\textrm{dwarf}}\,r_{{}_{\textrm{EFE}}}^{2}}{G_{\textrm{eff}}\,M_{p}^{2}\,n_{p}\,r_{h,\textrm{dwarf}}^{3}\left[1+\ln\left(\frac{a_{{}_{0}}}{g_{c,N}}\right)\right]}\,.$

We assume $W_{p}=1$ as in the Newtonian case. Our derivation assumes that
$g_{c,N}\ll a_{{}_{0}}$, which is valid in the Fornax Cluster. In general, we
recommend that the logarithmic term be omitted if $g_{c,N}>a_{{}_{0}}$.

##  4 Tidal susceptibility

Now that we have defined the main effects which can disturb the structure of a
dwarf in a galaxy cluster, we estimate the susceptibility of a dwarf to these
effects in both $\Lambda$CDM and MOND. To quantify the disturbance caused by
tides from the global cluster potential, we define the tidal susceptibility as
the ratio between the half-mass radius $r_{h}$ and the tidal radius
$r_{\textrm{tid}}$ of a dwarf:

$\displaystyle\eta_{\textrm{rtid}}~{}\equiv~{}\frac{r_{h}}{r_{\textrm{tid}}}\,.$
(29)

From the definition of $r_{\textrm{tid}}$ in both $\Lambda$CDM (Equation 11)
and MOND (Equation 20), we have that $r_{\textrm{tid}}\propto M^{1/3}$. This
implies that:

$\displaystyle\eta_{\textrm{rtid}}~{}\propto~{}\frac{r_{h}}{M^{1/3}}~{}\propto~{}\rho^{-1/3}\,.$
(30)

Therefore, only the density $\rho$ of the dwarf is relevant to its tidal
susceptibility in both $\Lambda$CDM and MOND.

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
$r_{\textrm{tid}}$ 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 $r_{\textrm{tid}}$ will then be small and its tidal
susceptibility high. Such a dwarf would be significantly disturbed by tides.
In the extreme case that $r_{\textrm{tid}}\ll r_{h}$ ($\eta_{\textrm{rtid}}\gg
1$), the dwarf will be destroyed within a few dynamical times. As a result, we
need to consider the maximum value of $\eta_{\textrm{rtid}}$ attained
throughout the trajectory, i.e., we need to evaluate $\eta_{\textrm{rtid}}$ at
pericentre.

If the disturbance is caused by interaction with massive galaxies
(harassment), we define the tidal susceptibility as the ratio between the age
of the elliptical galaxies in the Fornax Cluster ($t_{\textrm{Fornax}}\approx
10$ Gyr; Rakos et al., 2001) and the disruption time-scale $t_{d}$ of the
dwarf, which we assume to typically be about as old as the cluster itself.

$\displaystyle\eta_{\textrm{har}}~{}\equiv~{}\frac{t_{\textrm{Fornax}}}{t_{d}}\,.$
(31)

According to this definition, if $t_{\textrm{Fornax}}\ll t_{d}$ for a dwarf,
then it will hardly be susceptible to the effect of galaxy-galaxy harassment.
If instead $t_{\textrm{Fornax}}\gg t_{d}$ for a dwarf, then we expect that it
will be significantly disturbed due to this process.

Although our definitions for $\eta_{\textrm{rtid}}$ and $\eta_{\textrm{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
$\eta$ 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 $\eta\gg 1$. It is possible to have $\eta$ slightly above 1 due to
projection effects and other subtleties like the time required to achieve
destruction, which can be significant for $\eta_{\textrm{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 $\eta_{\textrm{rtid}}\gg 1$ or $\eta_{\textrm{har}}\gg 1$.
It is harder to falsify a theory in the opposite limit where it yields very
low values for both measures of $\eta$ 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

Figure 2: Histogram of the tidal susceptibility values of Fornax Cluster
dwarfs in $\Lambda$CDM (left column) and MOND (right 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),
0.05 (top right), 0.01 (bottom left), and 0.005 (bottom right). Notice the
different $\eta_{\textrm{rtid}}$ scales for $\Lambda$CDM and MOND. In both
theories, we typically have $\eta_{\textrm{har}}\ll\eta_{\textrm{rtid}}$.

Our first quantitative result is the susceptibility of dwarfs in the FDS
catalogue to cluster tides, which we calculate in $\Lambda$CDM and MOND using
Equations 11 and 20, respectively. We show the results as histograms in the
top row of Fig. 2, with $\Lambda$CDM shown on the left and MOND on the right.
The $\eta_{\textrm{rtid}}$ values are $\approx 5\times$ higher in MOND than in
$\Lambda$CDM. Since an isolated dwarf has a similar amount of self-gravity in
both frameworks by construction, the difference in $\eta_{\textrm{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 $\Lambda$CDM dwarf, which would retain the same dark matter
fraction within its baryonic extent throughout its trajectory.

The bottom row of Fig. 2 shows the susceptibility of FDS dwarfs to galaxy-
galaxy harassment according to $\Lambda$CDM (Equation 13) and MOND (Equation
3.3.2). In both theories, the histogram of $\eta_{\textrm{har}}$ peaks at very
low values such that $\eta_{\textrm{har}}\ll\eta_{\textrm{rtid}}$ and
$\eta_{\textrm{har}}\ll 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 $\Lambda$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 $\eta_{\textrm{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_{p}/r_{h,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}\,M_{\odot}$ and $r_{h,p}=4$ kpc, taking into account that
the actual mass and size are larger in $\Lambda$CDM and assuming a de
Vaucouleurs profile for the stars (de Vaucouleurs, 1948). We obtain that in
descending order of $M_{\star}$, 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}<10^{10}\,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 $\eta_{\textrm{har}}$ values
would still lead to its distribution having a mode $<0.1$ and all the dwarfs
having $\eta_{\textrm{har}}<1$.

Moreover, using $t_{\textrm{Fornax}}$ as the time-scale for interactions is an
optimistic assumption $-$ dwarfs in $\Lambda$CDM may have been accreted by the
cluster long after they formed, while in MOND they could be TDGs that formed
more recently (Renaud et al., 2016). This implies that the dwarfs would not
have experienced that many encounters with elliptical galaxies, which
themselves might only have been accreted $\ll 10$ Gyr ago. As an example, we
may consider the case of FCC 219 $\equiv$ NGC 1404, the most massive perturber
listed in table C1 of Iodice et al. (2019) in terms of $M_{\star}$. Its radial
velocity exceeds that of the brightest cluster galaxy NGC 1399 by 522 km/s,
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\pm 0.3$ Mpc (Hoyt et al., 2021), whereas the distance
to NGC 1399 is $20.0\pm 0.3$ Mpc (Blakeslee et al., 2009). Detailed modelling
in a $\Lambda$CDM context indicates that although NGC 1404 is not on a first
infall, it has likely spent $\la 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 $\eta_{\textrm{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.444This is consistent with the previous $\Lambda$CDM
result that harassment is not very significant for dwarfs in a Virgo-like
cluster (Smith et al., 2015). Thus, $\eta$ will be used to mean
$\eta_{\textrm{rtid}}$ unless stated otherwise. An important example of this
is our discussion of Newtonian TDGs that are purely baryonic, where
$\eta_{\textrm{har}}$ plays an important role (Appendix D).

###  4.2 Testing the effect of cluster tides on Fornax dwarfs

Figure 3: Distribution of the projected distances of Fornax Cluster dwarfs
against $r_{e}/r_{\textrm{max}}$, where $r_{e}$ is the projected half-light
radius and $r_{\textrm{max}}$ is the maximum $r_{e}$ at fixed $M_{\star}$ that
the dwarf can have to remain detectable given the surface brightness limit of
the survey of $27.8\,\textrm{mag}\,\textrm{arcsec}^{-2}$ in the $r^{\prime}$
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.

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 $r_{e}$ and $r_{\textrm{max}}$, where
$r_{\textrm{max}}$ is the maximum $r_{e}$ that the dwarf could have to remain
detectable given its $M_{\star}$ and the FDS detection limit of
$27.8\,\textrm{mag}\,\textrm{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 $R_{\textrm{sky}}$, e.g. if we consider a
horizontal line at $r_{e}/r_{\textrm{max}}=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.

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.555This 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).

Figure 4: 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

To emphasize this trend further, we use Fig. 4 to show the observed fraction
of disturbed dwarfs ($f_{d}$) in different $R_{\textrm{sky}}$ bins. This is
found as $f_{d}=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 $R_{\textrm{sky}}$ bin and $S\leq T$ is the number of these
galaxies which appear disturbed.666This is based on the binomial uncertainty
in $S$ assuming that the probability of a galaxy appearing disturbed in each
$R_{\textrm{sky}}$ bin is $f_{d}=S/T$. In reality, $f_{d}$ is not precisely
constrained by the observations $-$ we handle this complexity later (Equation
32). As expected from our previous results, $f_{d}$ is very high in the
central 200 kpc of the Fornax Cluster. Although $f_{d}$ is very low further
out, it is still non-zero and remains so out to the largest distances covered
by our dataset. 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.

###  4.3 Correlating tidal susceptibility with the observed level of
disturbance

Having obtained the tidal susceptibility $\eta$ of each Fornax dwarf in our
sample (Section 4.1), we can compare this to its visual level of disturbance.
We do so using the proportion of dwarfs classified as disturbed in each $\eta$
bin, which is similar to the analysis shown in Fig. 4 but binning in $\eta$
instead of $R_{\textrm{sky}}$. We consider each $\eta$ bin as an experiment
with $T$ trials (dwarfs) out of which $S$ are ‘successes’ (disturbed-looking
dwarfs). We then use binomial statistics to infer the probability distribution
of the disturbed fraction $f_{d}$ assuming a uniform prior over the range
$0-1$ and applying Bayes’ Theorem. The mean and standard deviation of $f_{d}$
are:

mean $\displaystyle=$ $\displaystyle\frac{S+1}{T+2}\,,$ (32) standard
deviation $\displaystyle=$
$\displaystyle\frac{1}{T+2}\sqrt{\frac{\left(S+1\right)\left(T-S+1\right)}{\left(T+3\right)}}\,.$

For the extreme case $S=T=0$, we expect that the probability distribution of
$f_{d}$ is uniform over the range $0-1$ as there is no data. In this case, we
recover the standard result that the mean of this distribution is $1/2$ and
its variance is $1/12$.

Figure 5: The likelihood that a dwarf appears disturbed in each tidal
susceptibility bin in $\Lambda$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
$\Lambda$CDM. In both cases, the last bin also includes all dwarfs with higher
$\eta$. Notice that the likelihood of a dwarf appearing disturbed rises with
$\eta$. The higher uncertainties at high $\eta$ are due to a small sample size
(see the $\eta$ distribution in Fig. 2).

We use Fig. 5 to plot the mean and standard deviation obtained in this way
against the central $\eta$ value for the bin under consideration. In both
$\Lambda$CDM and MOND, a clear trend is apparent whereby dwarfs with higher
$\eta$ are more likely to appear disturbed. We quantify this by dividing the
FDS sample into two subsamples where $\eta$ is below or above some threshold
$\eta_{t}$, thereby assuming only a monotonic relation between $f_{d}$ and
$\eta$ that is not necessarily linear. Appendix C explains how we obtain the
likelihood that the same $f_{d}$ 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
$\Lambda$CDM if we use $\eta_{t}=0.36$, in which case the null hypothesis of
$f_{d}$ being the same in both subsamples can be rejected at a significance of
$P=4.1\times 10^{-3}$ ($2.87\sigma$). If instead we use MOND, the optimal
$\eta_{t}=0.85$ and the significance rises to $P=4.4\times 10^{-4}$
($3.52\sigma$).777Section 8 provides a more rigorous quantification of how
confident we can be that $f_{d}$ rises with $\eta$. Though both theories imply
that $f_{d}$ is higher in the high $\eta$ subsample, $f_{d}$ starts rising at
a much lower value of $\eta$ in $\Lambda$CDM than in MOND, as clearly shown by
the optimal $\eta_{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 $f_{d}$ should start rising only when $\eta\ga 0.5$. This is not the case
in $\Lambda$CDM, which implies that dwarfs are more likely to be classified as
disturbed once their $\eta\ga 0.1-0.2$. A dwarf with such a low $\eta$ 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\ga 1-1.5$, which is much more plausible physically.

Another important aspect is the overall distribution of $\eta$, whose decline
towards the highest bin is responsible for a larger uncertainty in the
probability of appearing disturbed. The distribution of $\eta$ is shown
explicitly in the top row of Fig. 2. There are no $\Lambda$CDM dwarfs with
$\eta>0.7$, even though a dwarf with $\eta=0.7$ should still be tidally
stable. In MOND, the maximum $\eta\approx 3$, though there are very few dwarfs
with $\eta>1.5$. The high calculated $\eta$ 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 $\eta$ required for disturbance and the $\eta$
required for destruction can have to be consistent with observations $-$ both
in $\Lambda$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
semi-major axis ($R_{i}$) and eccentricity ($e$) values, with the integrations
started at $R=R_{i}$. The initial radii have a range of values from 15 kpc to
2015 kpc, while the eccentricities cover the full range of values for an
ellipse ($0<e<1$). The grid is divided into $100\times 100$ cells. Initially,
we assign the test mass a mass and half-mass radius which are typical for a
Fornax dwarf ($M_{\textrm{dwarf}}=3.16\times 10^{7}~{}M_{\odot}$ and
$r_{h}=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).

We initialize the simulated dwarfs for every possible combination of $R_{i}$
and $e$ as described below. We start the simulation at the semi-major axis of
the orbit, where the velocity $v$ satisfies

$\displaystyle v~{}=~{}v_{c}~{}=~{}\sqrt{-\bm{r}\cdot\bm{g}}\,.$ (33)

As discussed in section 2.3.1 of Banik & Zhao (2018c), the eccentricity $e$ is
defined such that

$\displaystyle
e~{}\equiv~{}\lvert\widehat{\bm{r}}\cdot\widehat{\bm{v}}\rvert\,,$ (34)

where $r=R_{i}\widehat{\bm{r}}$ and $\bm{g}=-g_{c}\widehat{\bm{r}}$, with
$\widehat{\bm{v}}$ indicating the unit vector parallel to any vector $\bm{v}$
of length $v$. The modulus is not required in our case because we start with
the dwarf going away from the cluster if $e>0$. Using Cartesian coordinates,
we define the initial positions and velocities of the orbit as:

$\displaystyle x~{}$ $\displaystyle=$ $\displaystyle~{}R_{i}\,,$ (35)
$\displaystyle y~{}$ $\displaystyle=$ $\displaystyle~{}0\,,$ (36)
$\displaystyle\vspace{10pt}v_{x}~{}$ $\displaystyle=$ $\displaystyle~{}ve\,,$
(37) $\displaystyle v_{y}~{}$ $\displaystyle=$
$\displaystyle~{}v\sqrt{1-e^{2}}\,.$ (38)

Equation 33 defines $v$ and Equation 34 sets the component of $\bm{v}$ along
the radial direction. $v_{y}$ is the remaining tangential velocity.

In order to obtain the positions and velocities of the simulated dwarf at each
point of the orbit, we implement a fourth-order Runge-Kutta integrator in 2D.
To ensure that the time-step we use for each iteration is computationally
efficient but also small enough to yield accurate results, we use an adaptive
time-step that depends on the dynamical time-scale at the instantaneous
orbital radius $R$:

$\displaystyle dt~{}=~{}0.01\sqrt{\frac{R}{g_{c}}}\,.$ (39)

We evolve the system for $t_{\textrm{Fornax}}=10$ Gyr, the estimated age of
the system (Rakos et al., 2001). At each time-step, we calculate the tidal
radius of the simulated dwarf at its current position and, by comparing this
with the half-mass radius, we obtain its instantaneous tidal susceptibility
$\eta$. We record the $e$ value of each simulated orbit and its final $R$, the
distance with respect to the cluster centre at which we should be seeing the
dwarf today. We also record two $\eta$ values in each orbit simulation: the
maximum $\eta$ over the whole simulation ($\eta_{\textrm{max}}$), and the
maximum $\eta$ in the last 2 Gyr ($\eta_{\textrm{max, recent}}$). We use
$\eta_{\textrm{max}}$ to decide whether the dwarf is destroyed and should be
removed from our statistical analysis. If not, then $\eta_{\textrm{max,
recent}}$ is used to set the likelihood that the dwarf appears disturbed. This
is because we expect a dwarf to return to a nearly undisturbed appearance if
it experiences only low $\eta$ values along its orbit for over 2 Gyr, provided
$\eta$ is never so high as to destroy the dwarf.

###  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, 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

The number density $n$ of dwarfs is assumed to be a function only of the
distance $R$ from the cluster centre. It is related to the radial probability
distribution $P_{r}$ as: $n\propto P_{r}/R^{2}$. We assume that $P_{r}$ is
described by a double power-law:

$\displaystyle
P_{r}~{}=~{}R^{2}\left(R+r_{\textrm{core}}\right)^{\textrm{Slope}_{P_{r}}}\,,$
(40)

where $r_{\textrm{core}}$ is the radius of the constant density central region
of the Fornax Cluster and $\textrm{Slope}_{P_{r}}$ is the power-law slope of
the radial profile in the cluster outskirts. To obtain a convergent number of
dwarfs, $\textrm{Slope}_{P_{r}}<-3$.

####  5.2.2 The eccentricity distribution

For the probability distribution of the orbital eccentricities, we assume a
linear function as in Banik & Zhao (2018c):

$\displaystyle
P_{e}~{}=~{}1+\textrm{Slope}_{P_{e}}\left(e-\frac{1}{2}\right)\,,$ (41)

where $\textrm{Slope}_{P_{e}}$ is the slope of the eccentricity probability
distribution.

####  5.2.3 Distribution of dwarf densities

The tidal susceptibility of a dwarf depends on both its mass and its radius,
which in general differ from the values assumed in our test mass simulation.
As discussed below Equation 29, the mass and radius of a dwarf affect its
tidal susceptibility only to the extent that they affect its density $\rho$.
Therefore, the $\eta$ values that we recorded in Section 5.1 should be
multiplied by a density-related factor accounting for the difference between
the intended density $\rho$ and the fixed value $\rho_{0}$ assumed in that
section. We therefore set

$\displaystyle\eta_{\textrm{max}}~{}$ $\displaystyle=$
$\displaystyle~{}\eta_{\textrm{max,
0}}\left(\frac{\rho}{\rho_{0}}\right)^{-1/3}\,,$ (42)
$\displaystyle\eta_{\textrm{max, recent}}~{}$ $\displaystyle=$
$\displaystyle~{}\eta_{\textrm{max, recent,
0}}\left(\frac{\rho}{\rho_{0}}\right)^{-1/3}\,,$ (43)

where the ‘0’ subscript denotes values obtained in Section 5.1. The $-1/3$
exponent comes from the fact that $\eta\propto{r_{h}}/M^{1/3}$ in both
theories.

Figure 6: The distribution of each dwarf galaxy’s mean baryonic density $\rho$
within its half-mass radius. The orange vertical line shows the sample mean,
while the magenta lines offset by $\pm 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
$\Sigma_{\textrm{min}}=0.26~{}M_{\odot}/\textrm{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)\,M_{\odot}/\textrm{pc}^{2}$. From
these and by adding (subtracting) the standard deviation in the $\rho/\Sigma$
ratios, we obtain the $\rho$ value given by the dashed (solid) black line.

The density $\rho$ of each Fornax dwarf within its $r_{h}$ can be inferred
from the data in the FDS catalogue using
$\rho=3M_{\star}/\left(8\mathrm{\pi}r_{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_{\odot}/\textrm{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 log-
normal distribution. This is assumed to be independent of $R_{i}$ 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.

Figure 7: The relation between mass and luminosity in the $r^{\prime}$-band
for the non-nucleated sample of dwarfs in the FDS catalogue, with both shown
on a $\log_{10}$ 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}\,M_{\odot}$ shows the stellar mass below which core formation is
inefficient in $\Lambda$CDM (Section 8).

In order to set the lowest density that can be assigned to a dwarf in a way
that is consistent with the observational constraints of the FDS, we check
down to which surface brightness $\mu$ dwarfs can be detected in this survey.
The limiting $\mu$ is given by the $1\sigma$ signal-to-noise threshold per
pixel, which in the FDS is $27.8\,\textrm{mag}\,\textrm{arcsec}^{-2}$ in the
red band (section 4.1 of Venhola et al., 2018). To infer the corresponding
$\rho$, we first convert this $\mu$ value to astronomical units
($L_{\odot}/\textrm{pc}^{2}$):

$\displaystyle\log_{10}\,\mu\left[L_{\odot}/\textrm{pc}^{2}\right]=\frac{\mu\left[\textrm{mag}\,\textrm{arcsec}^{-2}\right]-21.57-\textrm{Mag}_{\odot}}{-2.5}\,,$
(44)

where $\textrm{Mag}_{\odot}=4.65$ is the absolute magnitude of the Sun in the
red band (table 3 of Willmer, 2018). This gives
$\mu_{\textrm{min}}=0.23\,L_{\odot}/\textrm{pc}^{2}$. We then use the mass-
luminosity relation (solid grey line in Fig. 7) to obtain that
$M/L_{r^{\prime}}=1.10\pm 0.38\,M_{\odot}/L_{\odot,r^{\prime}}$. From this we
can convert $\mu_{\textrm{min}}$ to a surface density $\Sigma_{\textrm{min}}$
with some error due to the scatter in $M/L_{r^{\prime}}$, yielding
$\Sigma_{\textrm{min}}=0.26\pm 0.09\,M_{\odot}/\textrm{pc}^{2}$. Finally, we
can convert this $\Sigma_{\textrm{min}}$ to a threshold density $\rho_{t}$ by
plotting the surface density of the Fornax dwarfs against their volume density
and doing a linear regression (Fig. 8). Since the slope is very close to 1, we
fix it to 1 for simplicity, leading to a fixed ratio of $\rho/\Sigma=0.59\pm
0.33\,\textrm{kpc}^{-1}$. The limiting $\rho$ of the Fornax survey that we
obtain with this method is $\rho_{t}=1.51^{+1.67}_{-1.09}\times
10^{-4}\,M_{\odot}/\textrm{pc}^{3}$ considering the $1\sigma$ lower and upper
limits to both $M/L_{r^{\prime}}$ and $\rho/\Sigma$. From Fig. 6, we can see
that the distribution of dwarfs is only included in its entirety if we take
the lower limit and thus adopt a threshold of
$\rho_{t}=\rho_{\textrm{min}}=4.2\times 10^{-5}\,M_{\odot}/\textrm{pc}^{3}$.
Given that the images of the dwarfs have been carefully analysed by observers
and labelled as ‘unclear’ whenever the image was not clear enough, we assume
that all the considered dwarfs were observed without difficulty by the FDS.
Therefore, we consider that a reasonable lower limit for the density
distribution in our statistical analysis should encompass all the dwarfs in
the dataset, so we take $\rho_{t}=\rho_{\textrm{min}}=4.2\times
10^{-5}M_{\odot}/\textrm{pc}^{3}$ (black line in Fig. 6) as our nominal lower
limit to the density distribution. This choice of $\rho_{t}$ is 0.09 dex below
$\rho_{\textrm{min, FDS}}$, the lowest $\rho$ of any considered dwarf in the
FDS. If instead we had assumed that $\rho_{t}=\rho_{\textrm{mean}}=1.51\times
10^{-4}M_{\odot}/\textrm{pc}^{3}$ (grey line in Fig. 6), we would have needed
to discard 7 of the observed dwarfs in the FDS. These and other choices for
$\rho_{t}$ are discussed in Section 8.

Figure 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 $\rho/\Sigma$. The red star shows the values for the dwarf galaxy
used in our N-body simulations (Section 7).

In the $\Lambda$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 $\rho_{\textrm{min}}$ to a value 0.09 dex below
$\rho_{\textrm{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

Assuming that tides are the main cause of the apparent disturbance to the
structure of many Fornax dwarfs, we expect the probability of a dwarf
appearing perturbed to grow with its tidal susceptibility. We assume a linear
relation between $\eta$ and the probability of disturbance
($P_{\textrm{dist}}$) with slope:

$\displaystyle\textrm{Slope}_{P_{\textrm{dist}}}~{}=~{}\frac{P_{\textrm{dist,
ceiling}}-P_{\textrm{dist, floor}}}{\eta_{\textrm{destr}}-\eta_{\textrm{min,
dist}}}\,,$ (45)

where $\eta_{\textrm{min, dist}}$ is the lowest $\eta$ value at which the
dwarf is disturbed by tides, $\eta_{\textrm{destr}}$ is the $\eta$ value at
which the dwarf is destroyed (the algorithm rejects all simulated orbits in
which $\eta_{\textrm{max}}$ surpasses this value), $P_{\textrm{dist,
ceiling}}$ is the probability for a dwarf to appear disturbed right before it
gets destroyed at $\eta=\eta_{\textrm{destr}}$, and $P_{\textrm{dist, floor}}$
is the minimum probability for a dwarf to appear disturbed if
$\eta_{\textrm{max, recent}}<\eta_{\textrm{min, dist}}$. We allow
$P_{\textrm{dist, floor}}>0$ to capture the possibility that a dwarf appears
disturbed for reasons unrelated to cluster tides, e.g. asymmetric star
formation. Similarly, we expect that $P_{\textrm{dist, ceiling}}<1$ because a
significantly perturbed dwarf might be elongated along the line of sight and
thus appear circular. For a dwarf with $\eta_{\textrm{max,
recent}}\geq\eta_{\textrm{min, dist}}$, the probability of disturbance is:

$\displaystyle P_{\textrm{dist}}~{}=~{}P_{\textrm{dist,
floor}}+\textrm{Slope}_{P_{\textrm{dist}}}\left(\eta_{\textrm{max,
recent}}-\eta_{\textrm{min, dist}}\right).$ (46)

###  5.3 Comparison with observations

The observed parameters of the Fornax dwarfs that we aim to reproduce in our
simulation are:

  1. 1.

The distribution of sky-projected distances ($R_{\textrm{sky}}$) to the
cluster centre;

  2. 2.

The distribution of apparent $\eta$ values at pericentre
($\eta_{\textrm{obs}}$); and

  3. 3.

The disturbed fraction of dwarfs as a function of $\eta_{\textrm{obs}}$.

Because these quantities are projected or depend on the deprojection method,
we need to obtain the $R_{\textrm{sky}}$ values of our simulated dwarfs and
then deproject them using the same method that we use for the observed dwarfs.
To obtain the $R_{\textrm{sky}}$ values for each 3D distance $R$ of the
simulated dwarf, we consider the view if it is observed from all possible
angles $0^{\circ}\leq\theta\leq 90^{\circ}$ in steps of $1^{\circ}$, where
$\theta$ is the angle between $\bm{R}$ and the line of sight. The projected
distance is given by

$\displaystyle R_{\textrm{sky}}~{}=~{}R\sin\theta\,.$ (47)

Each value of $\theta$ is statistically weighted by the difference in
$\cos\theta$ across the corresponding bin. We then apply the deprojection
method described in Appendix A and obtain the corresponding distance at
pericentre (Appendix B). With this, we can calculate $R_{\textrm{tid}}$ and
$\eta$ at pericentre in a similar way to that in which we obtain these
parameters for the observed dwarfs. We name the new $\eta$ parameter that we
obtain with this method $\eta_{\textrm{obs}}$. Therefore, the simulated
quantities that we compare to the previously mentioned observables are:
$R_{\textrm{sky}}$, the distribution of $\eta_{\textrm{obs}}$, and the
probability of disturbance at each $\eta_{\textrm{obs}}$.

To do the comparison, we start by dividing the range of $R_{\textrm{sky}}$ and
$\eta_{\textrm{obs}}$ into several bins. We then classify the observed dwarfs
into these bins according to their values of projected distance or estimated
$\eta$ at pericentre. To obtain the probability for a dwarf to have a
projected distance or $\eta_{\textrm{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
$\eta_{\textrm{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 $R_{i}$, $e$, $\rho$, and $\theta$), we consider the same bins
as for the observed sample. For each bin, we add the probability that each
simulated dwarf has $R_{\textrm{sky}}$ or $\eta_{\textrm{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 $P_{\textrm{dist}}$ to the likelihood of each
$\left(R_{i},e,\rho,\theta\right)$ combination and add this to the appropriate
$\eta_{\textrm{obs}}$ bin. We then divide this sum by the probability of
$\eta_{\textrm{obs}}$ falling in that bin (i.e., without considering
$P_{\textrm{dist}}$).

To quantify how closely the properties of the simulated sample of dwarfs
resemble the properties of the observed FDS dwarfs in terms of each of the
above-mentioned observables, we use the binomial probability

$\displaystyle
P_{\textrm{x}}~{}=~{}\prod_{\textrm{Bins}}\frac{T!}{\left(T-S\right)!S!}p^{S}\left(1-p\right)^{T-S}\,,$
(48)

where $T$ is the total number of observed dwarfs, $S$ is the number of
observed dwarfs in a bin, $p$ is the simulated probability that a dwarf is in
that bin, and the ‘x’ subscript refers to the observable under consideration.
If this is the disturbed fraction, $T$ is the total number of observed dwarfs
in a particular $\eta_{\textrm{obs}}$ bin, $S$ is the observed number of
disturbed dwarfs in that bin, and $p$ is the probability given by the
simulation that a dwarf in that bin is disturbed. The total probability is
given by multiplying all probabilities for all the bins and all the
observables:

$\displaystyle
P_{\textrm{total}}~{}=~{}P_{R_{\textrm{sky}}}P_{\eta_{\textrm{obs}}}P_{\left.\textrm{perturbed}\right|\eta_{\textrm{obs}}}\,.$
(49)

In order to maximize this $P_{\textrm{total}}$, we leave as free parameters:
$r_{\textrm{core}}$, $\textrm{Slope}_{P_{r}}$, $\textrm{Slope}_{P_{e}}$,
$\eta_{\textrm{min, dist}}$, $\eta_{\textrm{destr}}$, $P_{\textrm{dist,
floor}}$, and $P_{\textrm{dist, ceiling}}$. We explore this set of parameter
values using the Markov Chain Monte Carlo (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 $P_{\rm{total}}$ and generate a proposal by adding
Gaussian random perturbations to the parameters, leading to a likelihood of
$P_{\rm{next}}$ with the revised parameters. The proposal is accepted if
$P_{\rm{next}}>P_{\rm{total}}$ or if a random number drawn uniformly from the
range $\left(0-1\right)$ is $<P_{\rm{next}}/P_{\rm{total}}$. If the proposal
is rejected, the parameter perturbations are not applied but the previous
parameters must be recorded once more.

Table 1: Priors for the free parameters in our model of the Fornax Cluster dwarf galaxy population. Parameter | Minimum | Maximum
---|---|---
$\textrm{Slope}_{P_{r}}$ | $-9$ | $-3$
$\textrm{Slope}_{P_{e}}$ | $-2$ | 2
$r_{\textrm{core}}$/Mpc | 0.01 | 3
$P_{\textrm{dist, floor}}$ | 0 | 1
$P_{\textrm{dist, ceiling}}$ | 0 | 1
$\eta_{\textrm{min, dist}}$ | 0 | 5
$\eta_{\textrm{destr}}$ | $\eta_{\textrm{min, dist}}$ | 5

We run a total of $10^{5}$ 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 $P_{\rm{total}}$ by increasing or decreasing the step
size according to how much $P_{\textrm{total}}$ 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 $P_{\textrm{total}}$ 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

Table 2: The parameters of our best-fitting model in each theory, obtained with the gradient ascent method (columns $2-3$) and based on $10^{5}$ MCMC trials (columns $4-5$). The last row shows the likelihood of the model (Equations 48 and 49). | Gradient ascent | MCMC
---|---|---
Parameter |  $\Lambda$CDM | MOND |  $\Lambda$CDM | MOND
$\textrm{Slope}_{P_{r}}$ | $-3.77$ | $-3.67$ | $-5.85$ | $-4.55$
$\textrm{Slope}_{P_{e}}$ | $-1.55$ | 0.34 | $-1.98$ | $-1.70$
$r_{\textrm{core}}$ | 0.62 | 0.65 | 1.35 | 0.90
$P_{\textrm{dist, floor}}$ | 0.09 | 0.04 | 0.10 | 0.02
$P_{\textrm{dist, ceiling}}$ | 0.65 | 0.76 | 0.54 | 0.53
$\eta_{\textrm{min, dist}}$ | 0.11 | 0.24 | 0.12 | 0.10
$\eta_{\textrm{destr}}$ | 0.24 | 1.88 | 0.23 | 1.24
$\log_{10}P_{\textrm{total}}$ | $-30.69$ | $-32.46$ | $-30.53$ | $-32.25$

The optimal set of parameters found by the gradient ascent algorithm are given
in columns 2 and 3 of Table 2 for $\Lambda$CDM and MOND, respectively. These
are the initial values at which we start the MCMC chains. Due to the use of
$10^{5}$ trials, the MCMC method provides a set of parameter values (a model)
that fits the observations slightly better (higher $P_{\textrm{total}}$ in
Equation 49) than we achieved with gradient ascent. The best-fit 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).

Figure 9: Comparison between observations and the best-fitting simulation in
$\Lambda$CDM (left column) and MOND (right column) in terms of the
distribution of projected separation $R_{\textrm{sky}}$ from the cluster
(first row), tidal susceptibility $\eta_{\textrm{obs}}$ (second row), and
likelihood that a dwarf appears disturbed as a function of
$\eta_{\textrm{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 $R_{\textrm{sky}}$ is
100 kpc in both theories. For $\eta_{\textrm{obs}}$, the bin width is 0.15 in
$\Lambda$CDM and 0.65 in MOND.

Using these parameters, Fig. 9 shows the simulated and observed probability
distributions of $R_{\textrm{sky}}$, $\eta_{\textrm{obs}}$, and disturbed
fraction vs. $\eta_{\textrm{obs}}$, revealing a good overall fit to the
observations in both theories.888The low values of $P_{\textrm{total}}$ arise
due to the large sample size. In particular, the rising likelihood of a dwarf
appearing disturbed as a function of $\eta_{\textrm{obs}}$ is nicely
reproduced by the best-fitting models.

###  6.2 Parameter uncertainties

Table 3: The most likely value and $1\sigma$ confidence interval of each model parameter in our test mass simulation of the Fornax Cluster dwarf galaxy population, based on $10^{5}$ MCMC trials. Parameter |  $\Lambda$CDM | MOND
---|---|---
$\textrm{Slope}_{P_{r}}$ | $-7.43_{-0.99}^{+2.24}$ | $-7.58_{-0.88}^{+2.18}$
$\textrm{Slope}_{P_{e}}$ | $-1.65_{-0.30}^{+1.80}$ | $0.75_{-1.22}^{+1.20}$
$r_{\textrm{core}}$ | $2.00_{-0.98}^{+0.34}$ | $2.02_{-0.88}^{+0.52}$
$P_{\textrm{dist, floor}}$ | $0.10_{-0.03}^{+0.03}$ | $0.07_{-0.03}^{+0.04}$
$P_{\textrm{dist, ceiling}}$ | $0.49_{-0.15}^{+0.30}$ | $0.79_{-0.20}^{+0.15}$
$\eta_{\textrm{min, dist}}$ | $0.11_{-0.06}^{+0.05}$ | $0.24_{-0.19}^{+0.24}$
$\eta_{\textrm{destr}}$ | $0.25_{-0.03}^{+0.07}$ | $1.88_{-0.53}^{+0.85}$

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\sigma$
confidence interval. The results are shown in Table 3.

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_{r}}$,
$r_{\textrm{core}}$, $P_{\textrm{dist, floor}}$, and $P_{\textrm{dist,
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_{r}}$ and $r_{\textrm{core}}$ are
not strong tests of the gravity law, while $P_{\textrm{dist, floor}}$ and
$P_{\textrm{dist, ceiling}}$ are set by the proportion of dwarfs in different
$\eta_{\textrm{obs}}$ bins that appear disturbed (Fig. 3). Unlike these four
parameters, $\textrm{Slope}_{P_{e}}$, $\eta_{\textrm{min, dist}}$, and
$\eta_{\textrm{destr}}$ cover very different ranges in these two models. As
discussed below, these are the parameters which can help us discern between
$\Lambda$CDM and MOND, allowing us to assess which model performs better when
compared to observations.

Figure 10: The $1\sigma$ confidence regions for the parameters in our model of
the Fornax Cluster dwarf galaxy population using $\Lambda$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\sigma$ confidence region for a pair of parameters. The results shown in
this ‘triangle plot’ are based on $10^{5}$ 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_{e}}$ (shown in the top panel of column 2
of Fig. 10) peaks close to the minimum allowed value of $-2$ in $\Lambda$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_{e}}\approx 2$ (for a derivation,
see section 4.2 of Kroupa, 2008). In this regard, MOND performs better than
$\Lambda$CDM.

Figure 11: 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
$\Lambda$CDM (orange).

The major differences between $\Lambda$CDM and MOND are in the parameters
$\eta_{\textrm{min, dist}}$ and $\eta_{\textrm{destr}}$, whose posterior
inferences are shown in detail in Figs. 11 and 12 due to their importance to
our argument. The low values in $\Lambda$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 $r_{\textrm{tid}}$ and thus reduces
$\eta$; see Equation 29). As a result, the algorithm needs to set
$\eta_{\textrm{min, dist}}$ and $\eta_{\textrm{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 $\eta$. 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 $\Lambda$CDM counterparts.
This causes the algorithm to choose significantly higher $\eta_{\textrm{min,
dist}}$ and $\eta_{\textrm{destr}}$ values in the MOND case.

Figure 12: The probability distribution of the tidal susceptibility at which a
dwarf is destroyed (Section 5.2.4) according to $\Lambda$CDM (orange) and MOND
(blue).

N-body simulations of dwarf galaxies show that $\eta_{\textrm{destr}}$ should
be $\approx 1$ in $\Lambda$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 $\Lambda$CDM and actual
observations of dwarf galaxies in the Fornax Cluster.

Turning to MOND, comparing the $\eta_{\textrm{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

As the last part of this project, we conduct our own N-body simulations of a
typical Fornax dwarf to find out the expected $\eta_{\textrm{destr}}$ in MOND.
The motivation is that while the analytic formula for the tidal radius
(Equation 20) should capture the scalings with the relevant variables like the
tidal stress and the EFE, there could be a constant numerical pre-factor that
arises from a detailed simulation. We investigated this using the Milgromian
N-body code phantom of ramses (por) developed in Bonn by Lüghausen et al.
(2015), who adapted it from the Newtonian N-body code ramses (Teyssier, 2002).
As a result, por inherits many features of ramses, including the adaptive mesh
refinement technique to better resolve denser regions. por can work with both
particle and gas dynamics. It is suited for simulations of isolated galaxies
(Banik et al., 2020; Roshan et al., 2021a; Banik et al., 2022b), interacting
galaxies (Renaud et al., 2016; Thomas et al., 2017, 2018; Bílek et al., 2018;
Banik et al., 2022a), galaxy formation (Wittenburg et al., 2020), and even for
cosmological structure formation (N. Wittenburg et al., in preparation). The
main difference between por and ramses is the fact that por solves the
ordinary Poisson equation twice, with $\bm{g}_{{}_{N}}$ found using standard
techniques in the first stage and the following equation solved in the second
stage to implement the MOND corrections:

$\displaystyle\nabla\cdot\bm{g}~{}=~{}\nabla\cdot\left(\nu\bm{g}_{{}_{N}}\right)\,,$
(50)

where $\nu$ was defined in Equation 16. The boundary condition for the
Milgromian potential $\Phi$ is:

$\displaystyle\Phi~{}=~{}\sqrt{GMa_{{}_{0}}}\ln r\,,$ (51)

where $M$ is the total mass in the simulation volume and $r$ is the distance
from the barycentre in the simulation unit of length, the choice of which has
no bearing on the result.

Since Fornax Cluster dwarfs are expected to contain little gas (Section 3), we
can simplify the set-up greatly by using the ‘particle-only’ version of the
por code. In particular, we use the ‘staticparts’ patch (described in section
4.1 of Nagesh et al., 2021) which allows the use of particles that provide
gravity but do not move if their mass exceeds a user-defined threshold. This
is helpful because we treat the cluster gravity as sourced by a point mass
fixed at the origin, with the dwarf at three possible initial distances
$R_{i}$. To ensure the gravity on the dwarf is the same as in the Fornax
Cluster, we use Equation 8 to obtain $g_{c}$ and then obtain the corresponding
$g_{c,N}$ with the simple interpolating function in the inverse form (Equation
18), from which we get the central mass:

$\displaystyle M_{c}~{}\equiv~{}\frac{g_{c,N}R_{i}^{2}}{G}\,.$ (52)

The different MOND dynamical cluster masses obtained in this way are:
$M_{c}=2.18\times 10^{12}~{}M_{\odot}$ at 150 kpc, $M_{c}=2.89\times
10^{12}~{}M_{\odot}$ at 300 kpc, and $M_{c}=3.31\times 10^{12}~{}M_{\odot}$ at
450 kpc. We use $7-13$ refinement levels and set the box length to $6\,R_{i}$
as the apocentre could be at almost $2\,R_{i}$.

For the dwarf, we use a half-mass radius of $r_{h}=0.84$ kpc and a total mass
of $M_{\textrm{dwarf}}=3.16\times 10^{7}~{}M_{\odot}$ represented by $10^{5}$
particles, making the mass resolution $316\,M_{\odot}$. These are typical
parameters for a dwarf in the Fornax Cluster (see the red star in Fig. 8).
Setting the velocity dispersion $\sigma$ 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
($x_{0}=R_{i}$, $y_{0}=0$, $z_{0}=0$) and the velocities are boosted depending
on the circular velocity at $R_{i}$ and the orbital eccentricity $e$, as
described in Section 5.1. We start the simulation with the dwarf at the semi-
major 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 $\eta$ 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 $\eta$ 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
$R_{i}$.

###  7.1 Analysis

We extract the particle positions $\bm{r}_{i}$, velocities $\bm{v}_{i}$, and
masses $m_{i}$ 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 ($\eta\ll 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

We apply an iterative outlier rejection scheme to accurately obtain the
barycentre position $\overline{\bm{r}}$ and velocity $\overline{\bm{v}}$ based
on the positions and velocities of the particles. In the first iteration, we
consider all the particles and calculate

$\displaystyle\overline{\bm{r}}~{}\equiv~{}\frac{\sum_{i}m_{i}\bm{r}_{i}}{M}\,,\quad
M\equiv\sum_{i}m_{i}\,.$ (53)

We use a similar definition for $\overline{\bm{v}}$. The barycentre position
and velocity are then used to find the root mean square (rms) dispersion in
position and velocity.

$\displaystyle
r^{2}_{\textrm{rms}}~{}\equiv~{}\frac{\sum_{i}m_{i}{\lvert\bm{r}_{i}-\overline{\bm{r}}\rvert}^{2}}{M}\,,$
(54)

with a similar definition used for $v_{\textrm{rms}}$, which we call $\sigma$
for consistency with other workers. This lets us define a $\chi^{2}$ statistic
for each particle based on its position.

$\displaystyle\chi^{2}_{\textrm{pos}}~{}\equiv~{}\left(\frac{\left|\bm{r}_{i}-\overline{\bm{r}}\right|}{r_{\textrm{rms}}}\right)^{2}\,,$
(55)

with a similar definition used for $\chi^{2}_{\textrm{vel}}$ based on the
velocity.

In the second iteration, we repeat the above steps for only those particles
whose $\chi^{2}_{\textrm{pos}}$ and $\chi^{2}_{\textrm{vel}}$ are both below
25, which changes the calculated quantities. In subsequent iterations, we
expect to have pinned down the barycentre more precisely, so we use the
stricter condition that

$\displaystyle\chi^{2}_{\textrm{pos}}+\chi^{2}_{\textrm{vel}}~{}<~{}\chi^{2}_{\textrm{max}}\,,$
(56)

where $\chi^{2}_{\textrm{max}}=11.83$ is set so that the likelihood of the
$\chi^{2}$ statistic for two degrees of freedom exceeding
$\chi^{2}_{\textrm{max}}$ is the same as the likelihood of a Gaussian random
variable deviating from its mean value by $\geq 3\sigma$. Our procedure can
thus be thought of as $3\sigma$ outlier rejection.

We consider the algorithm to have converged once the difference in
$\overline{\bm{r}}$ and $\overline{\bm{v}}$ between successive iterations is
so small that

$\displaystyle\frac{\left|\Delta\overline{\bm{r}}\right|^{2}}{r^{2}_{\textrm{rms}}}+\frac{\left|\Delta\overline{\bm{v}}\right|^{2}}{\sigma^{2}}~{}<~{}10^{-5}\,,$
(57)

with the additional requirement that the number of ‘accepted’ particles
deviates from that in the previous iteration by no more than the Poisson
uncertainty. In the analyses described below, we will only consider those
particles which are accepted on the final iteration.

####  7.1.2 Velocity dispersion

The velocity dispersion $\sigma$ is already available as part of our $3\sigma$
outlier rejection system for finding the barycentre of the dwarf. This 3D
$\sigma$ is found by applying Equation 54 but using velocities rather than
positions. If the dwarf were isolated and unaffected by tides, equation 14 of
Milgrom (1994) tells us to expect that

$\displaystyle\sigma~{}=~{}\left(\frac{4}{9}GMa_{{}_{0}}\right)^{1/4}\,.$ (58)

This assumes dynamical equilibrium and the deep-MOND limit, but does not make
any assumptions concerning whether the orbits are mostly radial or tangential.
If the system is not spherically symmetric, the velocity dispersion would not
be the same along every direction, but the bulk 3D velocity dispersion above
would still hold. Another important caveat is that the system should consist
only of particles with $m_{i}\ll M$.

####  7.1.3 Half-mass radius

To obtain the half-mass radius $r_{h}$, we order the particles in ascending
order of their distance to the above-determined dwarf barycentre
$\overline{\bm{r}}$. We then find the index $p$ such that

$\displaystyle\sum_{i=1}^{p}m_{i}~{}=~{}\frac{M}{2}\,,$ (59)

with the total mass $M$ of all accepted particles in general being slightly
below the initial mass of the dwarf. By definition, $r_{h}$ is the distance of
particle $p$ from the dwarf’s barycentre.

$\displaystyle r_{h}~{}\equiv~{}\lvert\bm{r}_{p}-\overline{\bm{r}}\rvert\,.$
(60)

####  7.1.4 Aspect ratio

To quantify the shape of the simulated dwarf, we obtain its inertia tensor

$\displaystyle\textbf{{I}}_{jk}~{}\equiv~{}\sum_{i}m_{i}\left(\bm{r}-\overline{\bm{r}}\right)_{j}\left(\bm{r}-\overline{\bm{r}}\right)_{k}\,,$
(61)

where the spatial indices $j$ and $k$ take values in the range $1-3$ because
there are three dimensions. We then find the eigenvalues of I. The aspect
ratio of the dwarf is defined as

$\displaystyle\textrm{aspect
ratio}~{}\equiv~{}\sqrt{\frac{\lambda_{\textrm{min}}}{\lambda_{\textrm{max}}}}\,,$
(62)

where $\lambda_{\textrm{min}}$ ($\lambda_{\textrm{max}}$) is the smallest
(largest) eigenvalue.

###  7.2 Results

Figure 13: 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 $R_{i}=150$ kpc (first
row), $R_{i}=300$ kpc (second row), and $R_{i}=450$ kpc (third row). In each
panel, the different curves show different orbital eccentricities as indicated
in the legend, which gives their corresponding $\eta$ 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\times 10^{5}$
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 isolated deep-MOND limit (Equation
58).

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.999This 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 $\Lambda$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 was
expected because dwarfs with more eccentric orbits have closer pericentre
passages and thus higher $\eta$ values at pericentre.

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 $r_{h}$ 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
$r_{h}$ 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 $\eta_{\textrm{obs}}$. The $r_{h}$ entering
into Equation 29 is the observed size, so ideally we would calculate $\eta$ 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 $\eta_{\textrm{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
$\eta$ to mean the value calculated in this way, though Table 4 also shows
results based on the size at pericentre.

Table 4: Summary of our MOND N-body simulation results for a Fornax dwarf with an initial distance of $R_{i}=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. |  $\eta$ using $r_{h}$ at $\ldots$ |
---|---|---
$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 $\eta_{\textrm{destr}}$, we focus mainly on models with
$R_{i}=150$ kpc as dwarfs with a larger semi-major 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 $\eta\leq 1.0$ respond
adiabatically. We choose $\eta_{\textrm{destr}}=1.4$ as the lowest value at
which a dwarf can get destroyed in MOND since dwarfs with this $\eta$ still
seem to be marginally capable of contracting their $r_{h}$ back to their
pericentre value by the time they reach apocentre.101010This certainly appears
to be the case for the $\eta=1.5$ model with $R_{i}=300$ kpc. For the upper
limit to $\eta_{\textrm{destr}}$ at pericentre, we choose a value of 2.0
because for this $\eta$, 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
$\eta_{\textrm{destr}}=1.70\pm 0.30$ if $r_{h}$ is measured at the second
apocentre. If instead we obtain $r_{h}$ at pericentre, then
$\eta_{\textrm{destr}}$ has a slightly lower value of $1.35\pm 0.05$.

As expected, $\eta_{\textrm{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 $\alpha=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.111111Without the $\frac{2}{3}$ factor in
Equation 20, the tidal susceptibility threshold is
$\eta_{\textrm{destr}}=1.13\pm 0.20$ when using $r_{h}$ at apocentre and
$\eta_{\textrm{destr}}=0.90\pm 0.03$ when using $r_{h}$ 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
$\eta_{\textrm{destr}}$ posterior.

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 $\eta$ value for only a very
brief period, allowing it to recover. Therefore, we would probably still be
able to observe dwarfs which have $\eta=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
$\eta$ 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