On the absence of backsplash analogues to NGC 3109 in the $\Lambda$CDM framework

The dwarf galaxy NGC 3109 is receding 105 km/s faster than expected in a $\Lambda$CDM timing argument analysis of the Local Group and external galaxy groups within 8 Mpc (Banik \&Zhao 2018). If this few-body model accurately represents long-range interactions in $\Lambda$CDM, this high velocity suggests that NGC 3109 is a backsplash galaxy that was once within the virial radius of the Milky Way and was slingshot out of it. Here, we use the Illustris TNG300 cosmological hydrodynamical simulation and its merger tree to identify backsplash galaxies. We find that backsplashers as massive ($\geq 4.0 \times 10^{10} M_\odot$) and distant ($\geq 1.2$ Mpc) as NGC 3109 are extremely rare, with none having also gained energy during the interaction with their previous host. This is likely due to dynamical friction. Since we identified 13225 host galaxies similar to the Milky Way or M31, we conclude that postulating NGC 3109 is a backsplash galaxy causes $>3.96\sigma$ tension with the expected distribution of backsplashers in $\Lambda$CDM. We show that the dark matter only version of TNG300 yields much the same result, demonstrating its robustness to how the baryonic physics is modelled. If instead NGC 3109 is not a backsplasher, consistency with $\Lambda$CDM would require the 3D timing argument analysis to be off by 105 km/s for this rather isolated dwarf, which we argue is unlikely. We discuss a possible alternative scenario for NGC 3109 and the Local Group satellite planes in the context of MOND, where the Milky Way and M31 had a past close flyby $7-10$ Gyr ago.


INTRODUCTION
The Universe was expanding rather homogeneously at early times (Planck Collaboration XXVII 2014), yet the present velocities of galaxies in the Local Group (LG) deviate significantly from a pure Hubble flow. This is due to the gravity they exert on each other. However, the large distance between the Milky Way (MW) and Andromeda (M31) galaxies implies only a rather weak gravitational attraction if we consider the Newtonian gravity of their baryons alone. This is insufficient to turn around their initial expansion and cause them to approach each other at the observed rate of ≈ 110 km/s (van der Marel et al. 2012). Clearly, there must be an extra source of gravity between the MW and M31. This 'timing argument' was one of the oldest arguments for missing gravity on galactic scales (Kahn & Woltjer 1959).
The now standard Lambda-Cold Dark Matter (ΛCDM) cosmological paradigm (Efstathiou et al. 1990;Ostriker & Steinhardt 1995) proposes that galaxies like the MW and M31 formed and evolved within haloes (White & Rees 1978) of non-baryonic cold dark matter (CDM) particles, accounting for the missing gravity. In this framework, the CDM haloes of the MW and M31 must be massive enough to turn around their initial expansion to the observed extent within the available 13.8 Gyr. This constrains their total mass (e.g. Carlesi et al. 2017). However, this is not a strong test of the ΛCDM model because there are as many data points as model parameters − the total mass and initial separation are varied to match the present separation and radial velocity (RV). This degeneracy can be broken by including data on more distant galaxies in the LG. Such an analysis was attempted by Sandage (1986), who found it difficult to simultaneously explain all the data then available. Using the catalogue of LG dwarfs in McConnachie (2012), a similar study was attempted by Peñarrubia et al. (2014), who found that the observations require an additional source of uncertainty with magnitude 35 ± 5 km/s. Following on from these spherically symmetric dynamical models, Banik & Zhao (2016) constructed an axisymmetric model of the LG consistent with the almost radial MW-M31 orbit (van der Marel et al. 2012) and the close alignment of Centaurus A with this line (Ma et al. 1998). The nearly radial nature of the MW-M31 orbit was later confirmed by van der Marel et al. (2019) and Salomon et al. (2021), which will be important to our discussion in Section 7.2. Treating LG dwarfs as test particles in the gravitational field of these three massive moving objects, Banik & Zhao (2016) investigated a wide range of model parameters using a full grid search. None of the models produced a good fit, even when reasonable allowance was made for inaccuracies in the model as a representation of ΛCDM based on the scatter about the Hubble flow in detailed N -body simulations (Aragon-Calvo et al. 2011). This is because several LG dwarfs have receding RVs much higher than expected in the best-fitting model, though the opposite is rarely the case (figure 9 of Banik & Zhao 2016). Banik & Zhao (2017) used the algorithm described in Shaya & Olling (2011) to test whether this remains the case when using a 3D model of the LG. The typical mismatch between observed and predicted RVs in the best-fitting model is even higher than in the 2D case, with a clear tendency persisting for faster outward motion than expected (figures 7 and 9 of Banik & Zhao 2017). These results are comparable to those obtained by Peebles (2017) using a similar algorithm.  borrowed this algorithm from Peebles and made some significant improvements in order to maximize the chance of finding trajectories consistent with the timing argument (see their section 4.1). Nevertheless, the results remained almost unchanged, with the only major difference being that Tucana became consistent with the model. An important clue is that nearly all the high-velocity galaxies (HVGs) are part of the NGC 3109 association, which was previously identified as having properties that are difficult to understand in ΛCDM (Pawlowski & McGaugh 2014a). The heliocentric RV of NGC 3109 is 403 km/s, which translates to 170 km/s in the Galactocentric frame, slightly below the expected value for a pure Hubble flow (without gravity) centred on the LG barycentre. However, taking into account the effect of Newtonian gravity, this is still 105 km/s too high in the best-fitting model .
Such a high RV reduces the LG timing argument mass inferred from the kinematics of non-satellite dwarf galaxies outside the MW and M31 virial volumes. This might well explain the unusually low LG mass of (1.6 ± 0.2) × 10 12 M found in this manner by Kashibadze & Karachentsev (2018), with their table 4 indicating that their analysis included the NGC 3109 association. Zhai et al. (2020) obtained a much higher timing argument mass of 4.4 +2.4 −1.5 × 10 12 M by searching cosmological simulations for analogues to the LG based on properties of the MW and M31 alone, especially with regards to their relative separation and velocity. This mass is in line with earlier results and simple analytic estimates neglecting information on LG galaxies other than the MW and M31 (Li & White 2008). The mass of M31 alone has been estimated at 1.9 +0.5 −0.4 × 10 12 M based on its giant southern tidal stream (Fardal et al. 2013). The total LG mass is certainly higher as it also includes the MW and material outside the major LG galaxies. Thus, several timing argument analyses of the whole LG found it difficult to explain the high RVs of some dwarf galaxies, with the tension phrased in some works as an anomalously low LG timing argument mass.
The timing argument calculations in Peebles (2017) and  are however not perfect representations of ΛCDM. They should handle long-range interactions between galaxies rather well, but can potentially miss important details due to the lack of dynamical friction between DM haloes. They consequently lack simulated mergers, during which galaxies can temporarily have a high relative velocity that could slingshot a nearby dwarf outwards at high speed in a three-body interaction. This leads to the existence of so-called 'backsplash galaxies' (backsplashers), defined as objects on rather extreme orbits that were once within the virial radius of their host but were subsequently carried outside of it. This backsplash process was studied in detail by Sales et al. (2007), who found it very difficult to get backsplashers at the 1.30 ± 0.02 Mpc distance of NGC 3109 (Soszyński et al. 2006). This is also evident in figure 3 of Teyssier et al. (2012), which additionally showed that backsplashers are very rarely more massive than 10 10 M regardless of their present position. Though NGC 3109 is more massive (Section 3.2), they argued that it is most likely a backsplasher given its position and RV. In a ΛCDM context, it might be very difficult to obtain such a massive and distant backsplasher due to the expected dynamical friction during any encounter with the MW or M31 (e.g. section 4.2 of Kroupa 2015). This would entail ejecting a galaxy as massive as NGC 3109 out of the inner regions of the MW halo against the inevitable dynamical friction. However, a more distant interaction with less dynamical friction would be rather weak, thus having little effect on the trajectory of NGC 3109 beyond that included in a few-body model.
In this contribution, we revisit the ΛCDM-predicted distribution of dwarfs around analogues to the MW or M31 using the much larger volume of the Illustris TNG300 cosmological hydrodynamical simulation (Pillepich et al. 2018). We also use the corresponding dark matter-only simulation to check how our results are affected by modelling of the baryonic physics. Our main goal is to find simulated backsplashers with NGC 3109-like properties today, but whose trajectories are likely to be seriously mis-modelled in the few-body analyses of Peebles (2017) and . If no such trajectories exist, this would improve our confidence in how well those models represent the underlying ΛCDM paradigm, thereby confirming the challenge posed by NGC 3109.
In Section 2, we describe the essential characteristics of the simulation for the present work. We then review the observed properties of NGC 3109 and our selection criteria in Section 3. In Section 4, we search for analogues of NGC 3109 in the simulation without requiring it to be a backsplasher, and review the timing argument analysis for the observed LG and its expected reliability. We then search for analogue backsplashers in Section 5, and present the results in Section 6. We discuss our results and an alternative scenario in Section 7, and conclude in Section 8.

COSMOLOGICAL SIMULATION
To explore whether the few-body models of Peebles (2017) and  might miss trajectories that explain the anomalous kinematics of NGC 3109 within a ΛCDM framework, we use the Illustris TNG300-1 hydrodynamical cosmological simulation (Pillepich et al. 2018;Nelson et al. 2019). This investigates ΛCDM in a cubic region with side length 205 h −1 co-moving Mpc (cMpc). The Hubble constant H0 in units of 100 km/s/Mpc is h = 0.6774, so the simulation box has a side length of 302.6 cMpc (hence the name TNG300). The suffix '-1' indicates that we use the highest available resolution setting for this box size within the Illustris suite. These simulations assume a standard flat cosmology in which the present fraction of the cosmic critical density in matter is 0.3089 (Planck Collaboration XIII 2016). The low H0 in this cosmology is also required by the early Universe observations of Aiola et al. (2020).
We use TNG300 because the larger simulation box compared to TNG100 or TNG50 allows for better statistics. All these simulations can adequately resolve objects much less massive than NGC 3109, as will become apparent in Section 6. However, we expect that we must search through many host galaxies analogous to the MW or M31 to find any backsplashers with properties similar to NGC 3109, or to set a stringent upper limit on their occurrence rate. In what follows, we will refer to host galaxies simply as 'MW analogues' even though the allowed range of properties are extended to allow M31-like galaxies, leaving open the possibility that NGC 3109 is backsplash from M31 (Section 3.3). However, we will see that this is much less plausible than backsplash from the MW.
The Illustris catalogues contain 100 snapshots going back from the present epoch to when the cosmic scale factor was a = 0.0475. The catalogues distinguish between groups and subhaloes. We use the redshift z = 0 group catalogue to identify isolated or LG-like host galaxies, though with some additional checks based on the z = 0 subhalo catalogue (Section 3.3). The position of the MW analogue at any epoch is found using the subhalo catalogue for that epoch, while its virial radius is found using the group catalogue as this is the only one to list virial radii. We use the subhalo catalogues to obtain properties of candidate backsplashers at various epochs. The Illustris sublink merger tree (Gómez et al. 2015) allows us to trace back MW analogues and to trace forward subhaloes within their virial volume, and finally to trace back any candidate backsplasher to better understand its trajectory.
We start by compiling in Section 3 the observed properties of NGC 3109, and our implemented criteria when selecting analogues in the Illustris TNG300 simulation. In Section 4, we search for analogues of NGC 3109 in this simulation without requiring it to be a backsplasher, fol-lowed by a review of the timing argument analysis in the detailed context of the LG. The rest of this paper concerns the backsplash analysis of the Illustris TNG300 simulation.

OBSERVED PROPERTIES OF NGC 3109
To select simulated galaxies analogous to NGC 3109, we first review its observed properties. Due to its small distance and fairly high mass, the uncertainties are rather small.

Distance
The Galactocentric distance of NGC 3109 is one of the most important observational inputs to our analysis. It was measured to be 1.30 ± 0.02 Mpc (Soszyński et al. 2006). Similar results were obtained by Dalcanton et al. (2009) and several other studies. To be conservative, we adopt a distance at the 5σ lower limit of the observationally allowed range. Thus, we require that analogues to NGC 3109 be 1.2 Mpc from their host.
It is possible that NGC 3109 is a backsplasher from M31, whose merger history appears to have been more active than that of the MW (e.g. Hammer et al. 2010;D'Souza & Bell 2018). According to table 2 of McConnachie (2012), the separation between NGC 3109 and M31 is currently 1.99 Mpc. The larger distance arises because the whole NGC 3109 association is in the opposite hemisphere on our sky compared to M31 (e.g. see figure 16 of ). We will see later that this makes a backsplash event in M31 a much less plausible scenario than backsplash from the MW. Thus, we focus almost exclusively on the hypothesis that NGC 3109 was once within the virial radius of the MW in a ΛCDM context. For Illustris host galaxies in a paired LG-like configuration, we require that the backsplasher lies 1.2 Mpc from both hosts.

Mass
In addition to the large distance of NGC 3109, we also need to explain its large mass for a backsplasher in the ΛCDM framework. Its virial mass can be estimated using rotation curve fits that add various halo profiles to the observed baryonic components. Several such fits were conducted by Li et al. (2020) based on the Spitzer Photometry and Accurate Rotation Curves dataset (SPARC; Lelli et al. 2016). The results are summarized in Figure 1, which shows only those halo profiles with a reduced χ 2 below 9. Based on these results, we conservatively assume that the mass of NGC 3109 is at least 10 10.6 M = 4.0 × 10 10 M , since no acceptable fits were obtained with a lower mass. Several studies give larger values, with Valenzuela et al. (2007) estimating a virial mass of 8.1×10 10 M (see their table 4). Moreover, the total mass (as recorded in the Illustris catalogue) is expected to exceed the virial mass (Peñarrubia & Fattahi 2017). For our purposes, using a lower mass is more conservative as we expect less massive subhaloes to be flung out further via the backsplash process, making it easier to match the large distance of NGC 3109. These mass estimates do not consider the rest of the NGC 3109 association, which consists of galaxies at a similar position with a similar anomalously high RV (Pawlowski &  . The estimated virial mass of NGC 3109 according to Newtonian rotation curve fits with different halo profiles (text labels), shown against the reduced χ 2 of the model if this is < 9 (Li et al. 2020). Based on these results, we conservatively assume that NGC 3109 is more massive than 10 10.6 M = 4.0 × 10 10 M .
McGaugh 2014a). As discussed in Section 7 and in , it is very unlikely that multiple dwarf galaxies were flung out in the same direction at a similar time despite initially being independent of each other. A more plausible explanation is that they formed a gravitationally bound association, though this is likely not bound any more (Kourkchi & Tully 2017). If so, the mass of NGC 3109 must be much higher, with Bellazzini et al. (2013) estimating a mass of 3.2 × 10 11 M . The relatively low stellar fraction (see below) could be due to tidal or ram pressure effects during a past interaction with the MW. Note that a past interaction could have resulted in loss of dark matter, so this estimate should be compared with the pre-interaction mass of each backsplasher. A conservative approach would be to compare with the maximum mass of each backsplasher at any snapshot in the Illustris simulation, which we discuss in Section 7.
Since Illustris is a hydrodynamical simulation, we may instead compare the baryonic mass of each backsplasher to that of NGC 3109. This was estimated at 2.1 × 10 9 M by applying rotation curve fitting techniques to high-resolution N -body models (table 4 of Valenzuela et al. 2007). Their estimated neutral hydrogen mass of (6 − 8)×10 8 M is similar to the (4.6 ± 0.5) × 10 8 M reported in table 8 of Carignan et al. (2013). Using stellar population synthesis modelling, Valenzuela et al. (2007) estimated that ≈ 5 × 10 8 M is contributed by stars, with the rest coming from gas.
The more recent SPARC database gives a 3.6 µm luminosity for NGC 3109 of (1.94 ± 0.02) × 10 8 L , which suggests a stellar mass of only 1.0 × 10 8 M for a mass to light ratio of 0.5 (Schombert & McGaugh 2014). Combining this with 1.33× their estimated neutral hydrogen mass of 4.77 × 10 8 M to account for primordial helium, we get a minimum possible baryonic mass of 7.3 × 10 8 M . Table 1 summarizes our adopted mass estimates for NGC 3109, where we have erred on the low side to be conservative. We focus on comparing the virial mass estimate with the total subhalo mass in Illustris backsplashers, Component of NGC 3109 Mass (M ) Stars 1.0 × 10 8 Baryons (in disc) 7.3 × 10 8 Virial (minimum) 4.0 × 10 10 Virial (maximum) 3.2 × 10 11 Table 1. Parameters of NGC 3109 used in this study. The different virial mass estimates refer to whether we consider the kinematics of NGC 3109 alone, or require it and its associated galaxies to be gravitationally bound (Bellazzini et al. 2013). We assume a Galactocentric distance of 1.2 Mpc to be conservative (Section 3.1).
as this should be least affected by uncertainties regarding subgrid baryonic feedback processes. When using the total mass, we still require that analogues to NGC 3109 have a non-zero stellar mass and thus a non-zero baryonic mass (for safety, we require both). We show later that our results remain much the same if we select backsplash analogues to NGC 3109 based on its stellar or baryonic mass.

Isolation conditions and host properties
We consider two kinds of host galaxy − LG-like and isolated.
In both cases, we identify appropriate hosts by considering the z = 0 group catalogue based on the Friends of Friends (FoF) approach. We require each FoF group to have at least 1 subhalo. If it has 2 subhaloes, we require the second most massive subhalo to comprise < 20% of the group virial mass (section 4.1 of Pawlowski & Kroupa 2020).
LG-like hosts consist of two FoF groups with a separation of (0.75 − 1.5) Mpc and total virial mass of (2 − 5) × 10 12 M . This mass range covers the LG mass estimated in various ways, e.g. it is similar to the range reported by González et al. (2014) and Zhai et al. (2020) based on LG analogues in cosmological simulations with a similar separation and relative velocity to the MW-M31 system. The 1D timing argument analyses of Peñarrubia et al. (2014) and Peñarrubia et al. (2016) give values near the lower end of this range. To have a reasonable mass ratio between the galaxies, we require that where Mmax (Mmin) is the virial mass of the heavier (lighter) member of the candidate pair. The lower limit on their separation is based on the 783 kpc distance to M31 (McConnachie 2012). Pairs with such a large separation are unlikely to have turned around and undergone an interaction within a Hubble time. Thus, requiring a present separation > 750 kpc implicitly imposes the condition that the MW and M31 did not undergo a past close interaction, which is correct in a ΛCDM context given their nearly radial orbit (van der Marel et al. 2012(van der Marel et al. , 2019Salomon et al. 2021) and the consequent very strong dynamical friction in any close encounter (Privon et al. 2013). Even without this consideration, a past flyby in Newtonian gravity would entail a very high timing argument mass (Benisty et al. 2019). Including any LG analogues in Illustris which had such an interaction could seriously compromise our analysis as there could be backsplashers from the interaction, which as argued above would not be a viable scenario in ΛCDM. The upper limit on the separation prevents interference from neighbouring groups beyond 3 Mpc (see below).
We ensure a sufficient level of isolation by requiring there to be no other group within 3 Mpc that is more massive than Mmin/3. We also remove pairs where there is another group more massive than 5 (Mmax + Mmin) within 5 Mpc, with the latter condition based on table 3 of Banik & Zhao (2017). This avoids massive nearby groups interfering with the dynamics of the LG, e.g. by pulling a backsplasher out to a much greater distance than it would otherwise reach.
For consistency with the above criteria, we require isolated hosts to have a virial mass M in the range (0.5 − 3.75) × 10 12 M . Their isolated nature is assured by requiring there to be no other group more massive than M/3 within 3 Mpc or more massive than 5M within 5 Mpc.
These selection criteria yield 13225 host galaxies, of which 640 are found in 320 LG-like paired configurations.

NGC 3109 ANALOGUES IN TNG300 AND THE
TIMING ARGUMENT

Frequency of NGC 3109 analogues ignoring the detailed environment of the LG
We now start our analysis by determining the RVs of nearby galaxies with respect to our selected hosts after imposing all the conditions compiled in the previous section (in Section 5.1 we will add to this the backsplash condition). This allows us to focus on the observed properties of galaxies without assumptions about their dynamical history or status as a backsplasher, thereby ignoring for now the detailed observed environment of the LG.
Since the RV will depend on distance, we restrict to a narrow distance range of (1.1 − 1.5) Mpc, which is wide enough to allow good statistics but narrow enough that there is little trend in RV with distance. The selected range is centred on the observed 1.30 Mpc distance to NGC 3109 (Soszyński et al. 2006), allowing also a ±10σ uncertainty. For LG-like hosts, we only consider the above-mentioned distance range relative to the less massive galaxy, and require a distance > 1.5 Mpc from the more massive galaxy. This resembles the LG somewhat more closely, though the statistics are dominated by isolated hosts.
The resulting distribution of RVs is shown in Figure 2, truncated to the range ±300 km/s for clarity. The observed Galactocentric RV of NGC 3109 is 170 km/s (figure 11 of Banik & Zhao 2016), which we show using a vertical grey line. This lies above the vast majority of the distribution, which is consistent with figure 6 of Teyssier et al. (2012) − though their results were based on a much smaller sample size. Although we need to allow a wide enough distance range to build up the statistics, it is clear that the dispersion in RV is much larger than can be explained by variation of the Hubble flow velocity over the narrow distance range considered, thus demonstrating the power of a large cosmological simulation. Importantly, our results show that it is quite possible to have dwarf galaxies receding from the MW as fast as NGC 3109 at its observed distance and with a comparable mass. 1.09% of our tracer galaxies have a higher RV, so the tension is mild. This decreases to 0.72% if imposing isolation conditions on the NGC 3109 analogue as described in Section 5.2.  Figure 2, but after subtracting the mean RV of all tracer galaxies around the same host. We only show results for hosts with 5 tracers. The vertical grey line at 105 km/s shows the RV excess of NGC 3109 relative to the best-fitting 3D ΛCDM timing argument analysis of . 1.28% of our tracers have a higher RV excess defined in the above sense. This is almost unchanged if restricting to only LG-like hosts (1.34%; not shown). Note that the mean RV calculation for each host is allowed to take advantage of tracer galaxies with any mass, but only galaxies more massive than NGC 3109 are shown here (see text).
Our results are based on stacking all host galaxies (with the above restriction on LG-like hosts). Since the hosts are not all equally massive, the expected RV at fixed distance will differ between hosts. To alleviate this, we next consider only hosts with 5 tracer galaxies of any mass, allowing us to calculate their mass-weighted mean RV, and thus the distribution of tracer galaxy RVs around that mean. To improve the accuracy with which this mean is calculated, we relax the condition on the tracer galaxy's mass to compute the mean.
We then subtract the mean RV from the tracer RV, and thereby determine the RV excess. Since not all hosts have 5 tracers with the appropriate position and at the same time at least one tracer more massive than NGC 3109, the statistics are somewhat noisier (Figure 3). For NGC 3109, we show a vertical grey line at 105 km/s, its RV excess compared to the  timing argument analysis. In this case, 1.28% of the distribution lies beyond 105 km/s. This decreases to 0.66% if imposing isolation conditions on the NGC 3109 analogue as described in Section 5.2. We conclude that ignoring the detailed observed environment of the LG, the 105 km/s RV excess of NGC 3109 is rare, but seemingly allowed at the percent level.

Including the LG environment: timing argument calculations and their reliability
At a frequency of ≈ 1%, the RV of NGC 3109 is already uncommon, but not necessarily severely problematic if one neglects the environment around the LG. Nonetheless, considering the 3D positions of perturbers outside the LG in more detail should in principle account for much of the RV variation of the above-discussed analysis. This is precisely what was done in the 3D timing argument calculations of Banik & Zhao (2017), Peebles (2017), and . A list of external perturbers taken into account can be found in table 3 of Banik & Zhao (2017). Despite including all these perturbers and letting their masses vary, Peebles (2017) and  were nevertheless unable to account for the large observed RV of NGC 3109. It is therefore worth discussing whether their models can be trusted to accurately represent ΛCDM expectations for its RV. Figure 9 of Banik & Zhao (2017) shows that if we suppose the model has a 25 km/s uncertainty, then it provides a good fit to the RVs of galaxies when the observed RV lies below the model prediction. The only exception is NGC 4163, but Peebles (2017) argued in his section 6.6 that it is part of the M94 group, and so excluded it from the timing argument analysis. As argued in section 5 of Banik & Zhao (2017), NGC 4163 may well be a backsplasher flung towards us, as suggested by its RV being ≈ 100 km/s lower than that of surrounding galaxies − roughly the amount by which its RV falls below the model prediction. Excluding this problematic galaxy 3.0 Mpc away, a 25 km/s model uncertainty would nicely explain discrepancies between the model and observations in cases where the latter give a lower RV. It is also reasonable on theoretical grounds − the dispersion in RV with respect to the LG barycentre at fixed distance from there should be ≈ 30 km/s (Aragon-Calvo et al. 2011). This suggests that even a 1D model of the LG should be accurate to about this much, so a 3D model can be expected to have an accuracy of ≈ 25 km/s if not better (see also section 4 of Banik & Zhao 2016).
Moreover, the timing argument is mostly sensitive to forces at late times, and thus mainly depends on the matter distribution today (figure 4 of Banik & Zhao 2016). The model includes the Large Magellanic Cloud at a mass of 2.03×10 11 M . 1 This automatically accounts for the induced reflex motion of the MW, which affects how we perceive the velocity field of the LG (Peñarrubia et al. 2016).
Ideally, one should apply the same method to all our selected hosts and their environment in Illustris. While beyond the scope of the present paper, we note that when applied to 32 LG analogues in the Millennium simulation (Springel et al. 2005), the method of Peebles et al. (2011) gave reliable results for the total mass, i.e. the deviations between true and inferred masses were consistent with the inferred uncertainty (Phelps et al. 2013).
Therefore, with an expected accuracy of 25 km/s, it is not at all clear why the overall best-fitting model should underpredict the RV of a fairly massive isolated galaxy like NGC 3109 by 105 ± 5 km/s (table 3 of . The discrepancy could thus be much more severe than the 1% frequency we found in Section 4.1. One explanation could be that some of the NGC 3109 analogues with a high RV are actually backsplashers, objects on rather extreme orbits that were once within the virial radius of their host but were subsequently carried outside of it. Indeed, one important aspect missing from the timing argument analyses is that they do not allow for dynamical friction on the extended dark matter haloes of galaxies, and the resulting mergers. They also do not account for the possibility of significant energy gain by a third galaxy near the spacetime location of the merger. To assess whether such backsplashers might resemble NGC 3109 today, we will in the rest of this paper investigate the distribution of backsplashers around our identified hosts in the Illustris TNG300 simulation.

NGC 3109 AS A BACKSPLASH GALAXY
We now require that the analogues of NGC 3109 in the Illustris TNG300 simulation are backsplash galaxies.

The backsplash condition
A backsplasher must have been within the virial radius of an MW analogue at some time in the past, but is by definition beyond its virial volume at the present time. To avoid making assumptions about the hypothesis being tested, we must allow for the possibility that the backsplasher is currently very far from the MW analogue that it interacted with. We keep the computational cost low by focusing on the virial volume of the MW analogue in all past snapshots for which z 5.22, thus going slightly further back than Teyssier et al. (2012).
We use the Illustris merger tree (Gómez et al. 2015) to trace back the MW analogue subhalo identified at the present epoch. We find the virial radius rvir of the group to which it belongs at each timestep with z 5.22, and then search through all the subhaloes within rvir of the subhalo corresponding to the MW analogue. A backsplasher candidate is defined as a subhalo in this list whose present-day root descendant lies beyond 2 rvir from the MW analogue, with rvir being time-dependent. Some subhaloes within the virial volume at a past epoch are absent from the merger tree or have a root descendant at a previous epoch, i.e. they do not survive up to the present. We reject such cases from further consideration. To improve the efficiency of our algorithm, and since subhaloes within the virial volume of an MW analogue several Gyr ago may well have merged with it by now, we first check if the root descendant is the present MW analogue itself, and reject such cases.
It is quite possible that the same backsplasher is identified within the virial volume of an MW analogue at several past epochs. To avoid double-counting, we keep track of all subhaloes in the z = 0 subhalo catalogue. Once some subhalo S in this catalogue has been recorded as a backsplash candidate, we ignore any subhalo at a previous epoch with root descendant S.
Since galaxies and subhaloes can merge in the hierarchical ΛCDM paradigm, we require that the main progenitor of a backsplasher was once within the virial radius of the host galaxy. Thus, each backsplash candidate is traced back along its main progenitor line to ensure that it was once within the virial radius of its host. In other words, we require that the bulk of the z = 0 subhalo has experienced a past backsplash encounter with the MW analogue. Without this restriction, we could have 'backsplashers' which mostly consist of material that never passed within the virial radius of a massive galaxy. We thus avoid situations where a low-mass backsplasher subsequently merges with a nearby massive galaxy A, causing that the root descendant of the backsplasher is A. Such scenarios are not a viable explanation for NGC 3109 because it is highly unlikely for an even lower mass dwarf to be flung out at a very high speed, only to subsequently catch up and merge with NGC 3109. While the latter's RV would be somewhat affected, the scenario would not explain the anomalous kinematics of other galaxies in the NGC 3109 association. These are most likely not currently bound to NGC 3109 (Kourkchi & Tully 2017) apart from Antlia, which is likely a satellite of NGC 3109 (section 6.4 of , and references therein). Thus, raising the RV of NGC 3109 alone would not be sufficient to raise the RVs of other association members. It also seems unlikely that a merger with NGC 3109 could raise its RV by 100 km/s without seriously disrupting its disc.

Isolation of the backsplasher
A particularly problematic aspect of NGC 3109 is its isolation, which implies that it was not substantially pulled away from the MW or M31 by a nearby massive group. Thus, NGC 3109 should have reached a Galactocentric distance of 1.2 Mpc without a significant 'helping hand' from large scale structure. To avoid selecting dwarf galaxies in Illustris which did receive such a helping hand, we impose various isolation criteria on both the host galaxy and the NGC 3109 analogue. This is motivated by the observed distribution of matter in and around the LG. To get a similarly isolated object as NGC 3109, we require there to be no subhalo with M > 5 × 10 11 M within 3 Mpc, with the mass threshold equal to the lowest allowed mass for an isolated host. An exception is made for the present-day descendant of the 'host' subhalo that the backsplasher once interacted with. 1 In cases where this host represents one member of an LG-like pair, we allow both members to be within 3 Mpc of the backsplasher. This imposes the condition that a 3 Mpc sphere centred on NGC 3109 contains no MW-mass galaxies other than the MW and M31, which is correct observationally ( Table 2).

Requiring energy gain
Our main purpose is to find trajectories with a similar final position to NGC 3109, but which would not be correctly modelled by the 3D timing argument analyses of Peebles (2017) and . Merely passing through the virial radius of an MW-like host galaxy is not sufficient to invalidate especially the latter analysis, since it should have enough time resolution to correctly model the interaction − partly because a softened force law was used to avoid singularities (Banik & Zhao 2017), leading to a smooth trajectory. Regardless of whether the MW mass distribution is modelled perfectly, a dwarf galaxy passing through would typically leave with the same energy as it came in because the model lacks dynamical friction. If such an energy-conserving encounter also happens in Illustris, then we can be fairly confident that it would be appropriately modelled in the Peebles (2017) analysis and in that of , which was very similar but had 10× better temporal resolution.
During galaxy-galaxy encounters, dynamical friction plays an important role (e.g. Privon et al. 2013;Kroupa 2015). This would cause the backsplasher to lose energy, reducing its final RV and making it even more difficult to explain the anomalously high RV of NGC 3109. Therefore, only trajectories with significant energy gain might explain the anomalous kinematics of the HVGs. Such trajectories would very likely be mis-modelled in a few-body timing argument analysis, so they could represent a viable ΛCDMbased explanation.
To focus on such trajectories, we extract the hostbacksplasher separation d (t) using the merger tree (Section 2). We define vin and vout as the backsplasher-host relative velocity at the times tin and tout, respectively, when the backsplasher enters and leaves the region within 2 r vir,mid of the host galaxy. We take r vir,mid to be the average of the virial radii at times tin and tout, with a similar definition used for M vir,mid and the mid-point cosmic scale factor a mid . We look backwards through the trajectory d (t) until d/rvir < 2, with tout being the snapshot when this first happens (looking backwards in time). tout is thus when the backsplasher crossed out of the twice-virial volume of its host, with dout being the separation at that time. To find a suitable choice for tin, we look back even further to find the snapshot when d/rvir is lowest, which is approximately when the subhalo has its closest approach to the host. We stop looking back if d/rvir > 2 again, ensuring that only the most recent encounter is considered in situations with multiple close encounters. We then consider all the snapshots 2 Gyr prior to the point of closest approach (minimum d/rvir), or back to the epoch when a = 0.1. Subject to these limits, we go backwards in time through the trajectory until d > dout again, choosing this to be our tin. If this condition is not satisfied within our allowed time window, then we pick tin to be the snapshot in this time interval with the greatest d, thereby minimizing the gap with dout and allowing as fair a comparison as possible between vin and vout. This scenario can arise if a backsplasher was originally a satellite orbiting well within the virial radius of its host since early times. Once we have found tin, the separation at that time is defined as din.
The above procedure minimizes the difference between din and dout. Nonetheless, some difference remains, making it inaccurate to directly compare the relative velocities at those snapshots. We therefore define ∆v ef f as our measure of the energy gain, where where Φ is the estimated specific potential energy of the backsplasher, with the subscripts indicating if the value is at time tin or tout. To find Φ, we assume that the host galaxy has a Navarro-Frenk-White density profile (Navarro et al. 1997) with the mass-concentration relation given in equation 4 of Duffy et al. (2008). Also adding an allowance for dark energy and dropping the in/out subscripts, we get that where the pivot mass M0 = 2×10 12 h −1 M , and the present fraction of the cosmic critical density in dark energy is ΩΛ,0 = 0.6911. In cases where energy has been lost such that Equation 2 does not yield a real square root, we set To minimize random fluctuations in our estimated potential, we use M vir,mid when calculating both Φin and Φout, neglecting the possibility of a change in host mass over the period in which a backsplasher is inside the twice-virial volume.
Although one can envisage more sophisticated schemes like considering the potential energy of each particle, this would involve handling a very large amount of particle-level data rather than the halo-level data used in our analyses, greatly increasing the computational cost. However, this would not much affect the results for a typical backsplash trajectory due to the good temporal resolution of the Illustris snapshots. Indeed, our results in Section 6 show that the potential adjustment term in Equation 2 is not very important. Moreover, our very conservative choice of threshold on ∆v ef f leaves a significant allowance for uncertainty in how it is calculated (see below).

Toy model
We now estimate the minimum ∆v ef f required to explain the anomalous RV of NGC 3109. For this purpose, we construct an idealized simulation in which a test particle moves under the influence of a galaxy with mass M = 2 × 10 12 M , which we refer to as the MW. As derived in section 2.1 of Banik & Zhao (2016) from General Relativity, the equation of motion for the particle position r relative to the galaxy contains a cosmological acceleration term in addition to the galaxy's gravity.
where r ≡ |r|, and an overdot denotes a time derivative. A force law of this form yields an extended region with a flat rotation curve of amplitude v f , which fixes rs = GM/v 2 f . We set rc = 0.01 rs to prevent a singularity at the centre. To obtain an MW-like galaxy, we use v f = 180 km/s (Kafle et al. 2012). With these values, rs ≈ rvir at z = 0.
We start our fourth-order Runge-Kutta integration when a = 0.1, with the test particle having a peculiar velocity towards the galaxy in addition to some tangential velocity vtan. Both parameters are varied to explore the parameter space. In each case, we must also choose the initial distance di of the test particle, which sets its Hubble flow velocity. Our goal is to find trajectories which turn around and undergo a close encounter with the MW. At that time, the particle's speed v is instantaneously increased as follows: This causes the particle to reach a larger distance than at first turnaround, which is typically at 500 kpc. As di is varied, perigalacticon occurs at different times. Increasing di causes the particle to have more energy, which in turn causes it to encounter the MW later, and to leave with higher v. When di is very small, increasing it significantly raises the post-encounter velocity while not much affecting the amount of time between the encounter and the present epoch, when our simulations end. This raises the present distance d f . However, raising di eventually causes the encounter to occur so late that d f starts decreasing again. We use a gradient ascent method (Fletcher & Powell 1963) to maximize d f as a function of di.
We repeat this procedure for a grid of initial radial and tangential peculiar velocities and ∆v ef f . The range of vtan is limited above by the requirement to have a close encounter with the MW. At the lower limit, we expect results to depend very little on vtan as the orbit is essentially radial. We estimate that the radial component of the peculiar velocity has a ±5σ uncertainty of ±250 km/s (equation 16 of Banik & Zhao 2017). Within this range, we map out the minimum ∆v ef f that allows the particle to reach a post-encounter separation of 1.2 Mpc from the MW within a Hubble time. Our results are shown in Figure 4. It is apparent that under conservative assumptions, we need ∆v ef f 150 km/s.
Assuming that such trajectories can be found in the Illustris simulation, we use our idealized setup to estimate the impact on the final RV. Without the impulse at pericentre, it is completely impossible for the particle to undergo a close approach to the MW and then reach a distance of 1.2 Mpc within the available timeframe. To facilitate a comparison, we construct a control non-impulsed (∆v ef f = 0) trajectory which never turns around and closely approaches the MW. This requires the use of a larger di.
An object ending up at larger d f generally has a larger RV, so this can be fairly compared between the trajectories only if they reach the same d f . Thus, we vary di for both the impulsed and the non-impulsed trajectories to ensure that d f = 1.2 Mpc. With the impulsed trajectory, this implies the maximum possible d f > 1.2 Mpc, so there are two possible choices for di. We choose the larger di since this causes the perigalacticon to occur later, implying the particle must have a larger final RV to reach the same d f . This lets us find how much the final RV could differ between the impulsed and control trajectories. In both cases, we use a Newton-Raphson algorithm to vary di in order to precisely achieve d f = 1.2 Mpc. Figure 5 shows an example where the impulsed trajectory has ∆v ef f = 180 km/s, vtan = 50 km/s, and an initial radial peculiar velocity towards the MW of 100 km/s. For the non-impulsed trajectory, we use the same initial peculiar velocity but larger di. The final RV is 9.73 km/s for the non-impulsed trajectory and 119.64 km/s for the impulsed trajectory. It is clear that the impulse has provided an alternative high-velocity route to reaching the presently observed distance of NGC 3109. The RV excess of this route compared to the 'traditional' (non-impulsed) route is 110 km/s. This would nicely explain why the RV of NGC 3109 exceeds the prediction of the Banik & Zhao (2018) model by 105 km/s. with ∆v ef f = 180 km/s, vtan = 50 km/s, and a radial peculiar velocity towards the MW of 100 km/s when a = 0.1. For comparison, we show a similar trajectory with larger initial separation and no close encounter with the MW such that ∆v ef f = 0 (red). Both trajectories reach the same present distance of 1.2 Mpc, but the present RV of the impulsed trajectory is 109.9 km/s higher. Another impulsed trajectory can also be constructed with the same initial peculiar velocity and final distance but with lower d i , an earlier encounter, and lower final RV (not shown).
We next consider whether backsplashers in the Illustris TNG300 cosmological simulation with the mass and present distance of NGC 3109 ever have trajectories with a similarly large ∆v ef f . This is possible if a dwarf galaxy closely encounters the MW while it is undergoing a minor merger − a significant amount of energy could be gained in a threebody interaction. But such scenarios could prove too rare, or dynamical friction on the dwarf could slow it down such that there is a net loss of orbital energy (∆v ef f < 0). Our simplified analysis in this section neglected the role of dynamical friction, which could be important at the mass of NGC 3109 (Section 7).

RESULTS
We begin by showing the distribution of backsplasher total mass and present distance from the host (Figure 6). In case of an LG-like host, we show the minimum distance from either of the host galaxies. For now, we do not impose any restriction on ∆v ef f . Without this restriction, we identify 1438 backsplashers. The vast majority of these lie at distances 1 Mpc and mass 10 10.5 M . Both the distance and mass of NGC 3109 are individually highly unlikely if it is drawn from the distribution of backsplashers in ΛCDM. Only a handful of backsplashers match NGC 3109 in both respects, and even then only if we use the lower possible NGC 3109 mass of 10 10.6 M . This is only just allowed in some Newtonian rotation curve fits (Figure 1).
Our results in Figure 6 suggest that it is sometimes possible to get backsplashers with the mass and distance of NGC 3109. However, we have not yet considered whether the trajectories of these three backsplashers truly represent behaviour that would not be captured by the modelling of Figure 6. Distribution of total mass and distance from host for the backsplashers we identify in the Illustris TNG300 simulation. No restriction is imposed on the energy gain ∆v ef f during the encounter. Outside the high-density region (contours), we show individual backsplashers, with the colour indicating the encounter time when d/r vir was lowest (Section 5.3). The distance and virial mass of NGC 3109 are shown as blue stars for two possible masses ( Table 1). The total number of backsplashers is indicated at the top. The trajectory of the circled backsplasher slightly right of centre is shown in solid black in Figure 7. . To investigate this further, we use Figure 7 to show the distance and relative velocity between these backsplashers and their hosts as a function of time. It is evident that in all three cases (curves except that in solid black), the trajectory appears symmetric around the time of closest approach to the host. This is borne out by the values of ∆v ef f in km/s with (without) the potential energy adjustment term in Equation 2, which in increasing order of final distance are −160 (−151), −90 (−90), and −121 (−104), indicating energy loss in all cases. For illustrative purposes, we also show a fourth backsplasher (solid black) with a lower mass of 2.26 × 10 10 M . In this case, the backsplasher has clearly gained energy during the encounter, as also evident in that its ∆v ef f = 225 km/s. The good time resolution of the Illustris snapshots allows us to measure the relative velocity at essentially the same separation from the host before and after the encounter, minimizing the potential energy adjustment (if we neglect this, we get ∆v ef f = 247 km/s). Thus, the Illustris TNG300 simulation contains genuine backsplash trajectories with ∆v ef f > 0, but not at the high mass and distance of NGC 3109.
Our simplified model in Section 5.3 suggests that trajectories with ∆v ef f 150 km/s are unable to reach a present distance of 1.2 Mpc. In a cosmological simulation, the effects of large scale structure allow a dwarf to reach this distance despite having a close encounter with an MW analogue which yields zero or even negative ∆v ef f . However, as discussed in Section 5.3, such trajectories should be modelled correctly in the 3D timing argument analyses of Peebles (2017) and . Thus, the failure of their model to correctly represent NGC 3109 cannot be understood using Illustris trajectories with ∆v ef f < 0. If anything, loss of energy during a past close interaction with the MW (e.g. due to dynamical friction) would make Notice that the red and blue trajectories appear symmetric around pericentre, and thus show little sign of energy gain while interacting with their host. This is borne out by their somewhat negative ∆v ef f (see text). For comparison, we also show the trajectory of the most distant backsplasher in our sample with ∆v ef f 0 (solid black), even though its mass is less than that of NGC 3109 (circled object in Figure 6). In this case, energy gain is apparent (∆v ef f = 225 km/s).
it even more difficult to explain the anomalously high RV of NGC 3109.
To be very conservative, we first consider the distribution of backsplashers if we merely require that ∆v ef f 0. Our results show no backsplash analogues to NGC 3109 ( Figure 8). We can find one analogue if we reduce the required mass to 2.26×10 10 M . However, the rotation curve of NGC 3109 reaches an amplitude of ≈ 80 km/s at a distance of 12 kpc, and is likely still rising there (figure 13 of Carignan et al. 2013). This implies a Newtonian dynamical mass of 1.8 × 10 10 M within 12 kpc, making it highly unlikely that the virial mass of NGC 3109 is only 2.3 × 10 10 M . Plausible rotation curve fits in a ΛCDM context yield significantly larger values, with none suggesting a mass below 10 10.6 M = 4.0 × 10 10 M (Figure 1).  Our results in Figure 4 suggest that ∆v ef f 150 km/s for a backsplasher to reach the present distance of NGC 3109 along a substantially different route (and hence different final RV) to a non-backsplash galaxy at the same present position. A smaller impulse would mean the object had more of a helping hand from large scale structure, weakening the case that its trajectory would not be correctly modelled by . Thus, we use Figure 9 to show the effect of requiring ∆v ef f 150 km/s. The overall distribution of backsplashers remains similar, but their number is reduced more than four-fold. NGC 3109 is now much further from the backsplashers' distance and mass distribution. This is especially true if we assign NGC 3109 a mass of 10 11.5 M , as would be required to once have bound the whole NGC 3109 association (section 3 of Bellazzini et al. 2013). We discuss this issue further in Section 7.  (Table 1).
Since Illustris is a hydrodynamical simulation, we can also consider the baryonic and stellar masses of the backsplashers we have identified. The results confirm that the distance of NGC 3109 and its baryonic or stellar mass are indeed significantly higher than expected for backsplashers in ΛCDM (Figure 10). In particular, the bottom panel shows that the well-constrained stellar mass of NGC 3109 is quite unusual for a backsplasher at any distance. This could be related to ram pressure stripping of the backsplasher's gas while it closely encounters a more massive galaxy (see also Teyssier et al. 2012).
So far, we have not distinguished between whether the host galaxy of a backsplasher is isolated or in an LG-like pair. This allows us to build up much better statistics, since we only have 640 host galaxies in 320 LG-like paired configurations. However, this is sufficient to get a good idea if the mass-distance distribution of backsplashers is similar around LG-like host galaxies. We therefore conduct a similar analysis to Figure 8 but only for LG-like hosts, with the result shown in Figure 11. It is clear that the overall distributions are very similar, though the smaller number of objects in the latter case causes the tail to be sampled less well. As a result, there are now no backsplashers as distant as NGC 3109 for any mass, even with our very conservative distance estimate of 1.2 Mpc. The similarity of results between LG-like hosts and the full sample (with mostly isolated hosts) indicates that the presence of M31 does not make it easier to explain how NGC 3109 could be a backsplasher in a conventional gravity context.
Our results allow us to consider whether NGC 3109 could be a backsplasher from M31 rather than the MW. This would require a present distance from the host of 2.0 Mpc (Section 3.1). However, none of the ∆v ef f 0 backsplashers associated with LG-like hosts reach a separation of even 1.2 Mpc, and generally also have a much lower mass than NGC 3109. Clearly, a 2 Mpc separation with the host would make NGC 3109 significantly more of an outlier from the expected backsplasher distribution of distance and mass. Therefore, the (very small) probability of NGC 3109 being a ΛCDM backsplasher arises mostly from the chance that a suitable backsplash event occurred near the MW.

Comparison with dark matter-only simulation
Our analysis thus far has focused exclusively on the Illustris TNG300 simulation. This can easily resolve haloes with the mass of NGC 3109 (Figure 6), while the larger simulation volume than e.g. TNG100 should allow for better statistics. The backsplash process mainly revolves around the motions of fairly massive galaxies, so baryonic physics should play only a small role.
To check this, we compare our results with the dark matter-only version of Illustris TNG300. The analogous results to Figure 8 are shown in Figure 12. The overall massdistance distribution of backsplashers is quite similar in the dark matter-only run, but there are many more backsplashers in this case. This could be related to the much stronger tides upon closely approaching the host galaxy. In a hydrodynamical simulation, this would typically contain a baryonic component that is much more centrally concentrated than the dark matter. Particularly strong tides would arise  Figure 6). The lone backsplasher with higher distance and mass than NGC 3109 (for the lower mass estimate) has ∆v ef f = 105 km/s, which we argued in Section 5.3 is insufficient to explain its anomalously high RV. For clarity, we have omitted a handful of low-mass backsplashers at large distances − these are much less massive than NGC 3109.  Figure 12, but now requiring ∆v ef f 150 km/s. There are no longer enough backsplashers to reliably draw contours of their number density, so they are shown individually. Notice the rather similar result to the corresponding hydrodynamical simulation (Figure 9).
if the host galaxy develops a thin disc, which can efficiently disrupt haloes passing close to it (Pawlowski et al. 2019). Our results in Figure 12 indicate that there is one backsplasher with a marginally larger distance and mass than NGC 3109 for our very conservative choices of these parameters (Section 3). However, this backsplasher has ∆v ef f = 105 km/s, which suggests that it did not gain enough energy during the interaction to explain the anomalously high RV of NGC 3109 (Section 5.3). As discussed there, a more realistic picture can be obtained by requiring ∆v ef f 150 km/s, which yields the results shown in Figure 13. This demonstrates that NGC 3109 is a significant outlier also in the dark matter-only version of TNG300, with the results being rather similar to the standard hydrodynamical version used elsewhere in this contribution (Figure 9). Therefore, it is clear that baryonic physics has only a small effect on our conclusion that NGC 3109 is too distant and massive to be a backsplasher from the MW or M31 in a ΛCDM context.

DISCUSSION
We showed that no backsplashers in the Illustris TNG300 simulation have the right mass and distance to be considered analogues of NGC 3109 even under conservative assumptions. This is consistent with the analytic estimate that backsplashers should not be found at d 2.5 rvir from their host (Mamon et al. 2004). Due to the large number of MW analogues, we are able to get some backsplashers at even larger distances. However, this is quite rare − we found only 781 backsplashers with ∆v ef f 0 from 13225 host galaxies (Figure 8). Some of these probably have d < 2.5 rvir as we only require d > 2 rvir. Thus, our results broadly support the analytic estimate of Mamon et al. (2004).
For a paired host configuration, numerical simulations show that backsplashers can reach up to ≈ 5 rvir from their host (Ludlow et al. 2009;Wang et al. 2009), broadly consistent with our results in Figure 11 for rvir ≈ 200 kpc. This was also demonstrated in figure 1 of Teyssier et al. (2012). However, this figure demonstrates that the distribution of backsplashers is significantly elongated along the axis connecting the two main galaxies. In the orthogonal direction, the extent is similar to the analytic estimate of 2.5 rvir = 660 kpc for Mvir = 2 × 10 12 M . In this regard, it is worth mentioning that NGC 3109 lies 706 kpc from the MW-M31 axis (765 kpc for a distance of 1.3 Mpc), and appears on our sky in the opposite hemisphere to M31 (e.g. see figure 10 of Banik & Zhao 2016). Most of the backsplashers in Teyssier et al. (2012) are located quite close to the two main host galaxies, possibly because their combined gravity makes it difficult to reach a large distance from their barycentre. Their figure 1 shows that it is very difficult to find a backsplasher whose minimum distance from either host is 1.2 Mpc and which lies 700 kpc from the axis between the hosts.
In addition, figure 4 of Teyssier et al. (2012) indicates that regardless of the position, backsplashers more massive than 10 10.2 M are very rare. The mass of NGC 3109 is thus unusually high for a backsplasher even if we assume that its mass is 10 10.6 M , which is the minimum required in Newtonian rotation curve fits with plausible dark matter haloes (Figure 1). Our results agree that this is unusually massive for a ΛCDM backsplasher at any distance ( Figure  8).
To better understand the properties of backsplashers in ΛCDM, we use Figure 14 to show the distribution of ∆v ef f and present distance from the host. As expected, backsplashers from less massive hosts are generally still quite close to their host and did not gain much energy when interacting with it. Since most of our host galaxies are more massive than 10 12 M , our results should not be much affected by slight adjustments to the lower limit on the allowed host mass. Changing the upper limit would have a more significant impact on the statistics, but would preferentially remove those backsplashers which get closest to reproducing the observed properties of NGC 3109.
In a ΛCDM context, an important reason for our lack of backsplash analogues to NGC 3109 is that dynamical friction would be quite significant during any past close interaction with the MW or M31 (Privon et al. 2013). To understand the effect of dynamical friction, Teyssier et al. (2012) used their figure 4 to show how the mass distribution of backsplashers changes with distance from the host. The statistics were limited as only one LG-like host was considered. We revisit this issue in a slightly different way using our sample of 13225 host galaxies. The energy gain ∆v ef f should be smaller at higher mass due to the effect of dynamical friction. Thus, we use Figure 15 to show the relation between ∆v ef f and backsplasher mass for the cases where ∆v ef f 0, with the host mass indicated using the colour. As expected, the upper limit to ∆v ef f declines with mass. There are very few backsplashers more massive than NGC 3109 with ∆v ef f 150 km/s, which we argued is required to reach 1.2 Mpc from the host (Section 5.3). The handful of such massive backsplashers all have d < 1.2 Mpc (Figure 9), even if we relax the energy gain requirement to a very conservative ∆v ef f 0 ( Figure  8).
Another reason for this lack of backsplash analogues to NGC 3109 could be that backsplashers lose mass during their encounter with a more massive galaxy (e.g. Smith et al. 2016). This is apparent by comparing the top and bottom panels of our Figure 15, which show the maximum and present mass, respectively, of each backsplasher. The loss of mass during the encounter is also apparent upon closer examination of the only trajectory shown in Figure 7 where the backsplasher significantly gained energy during the interaction. The high ∆v ef f of 225 km/s comes at the price of the backsplasher mass decreasing from 3.14 × 10 10 M at tin Figure 15. ∆v ef f as a function of maximum (top) or present (bottom) backsplasher total mass, showing only those objects with ∆v ef f 0. The colour indicates the present mass of the host galaxy. Notice the upper limit to ∆v ef f declines with mass. The most massive backsplasher in both panels is the same object, which is currently at a distance of 1.16 Mpc and thus is not analogous to NGC 3109. to 1.48 × 10 10 M at tout. Neither mass would be enough to explain the observed rotation curve of NGC 3109, but the situation is substantially worse post-encounter.
As discussed in Section 3.2, our adopted mass for NGC 3109 should also consider the rest of the NGC 3109 association. While it may well be unbound today (Kourkchi & Tully 2017), it should have been bound in the past in order to explain the filamentary nature of the NGC 3109 association (Bellazzini et al. 2013). Moreover, the other galaxies in this association are likely also backsplashers since their RVs are also too high to be accounted for by the timing argument analyses of Peebles (2017) and . Although the Illustris simulation may well struggle to resolve galaxies like Leo P, our results strongly suggest that it would be very difficult to find a backsplasher of any mass at its observed distance of 1.62 ± 0.15 Mpc (McQuinn et al. 2015). Even if this were possible, we would have to explain why several unassociated dwarf galaxies should be flung away from the LG in much the same direction and end up at a similar distance, suggesting a similar encounter time with the host galaxy. If these dwarfs were falling along a filament onto the MW while it was undergoing a minor merger, then the short dynamical timescale of the merger implies that even slight differences in the infall time of the dwarf could substantially alter the direction in which it is ultimately flung out. Moreover, galaxies with different mass should have experienced different amounts of dynamical friction. This means that even if several galaxies fell into the MW along a filament at much the same time and were ejected outwards in a similar direction, they would end up at rather different distances, e.g. Sextans A would still end up much further ahead than NGC 3109. In reality, both have a similar distance, with Sextans A only ≈ 12% further away (McQuinn et al. 2017).
It is also worth mentioning that essentially all the HVGs identified by  lie in the NGC 3109 association, even though they considered 33 galaxies in addition to the MW and M31. Their figure 10 shows that there are at best three other HVGs in addition to those in the NGC 3109 association. Of these, the distance to HIZSS 3 is seriously questionable due to observational difficulties caused by its extremely low Galactic latitude of 0.09 • (section 6.3 of . Meanwhile, the RVs of KKH 98 and DDO 190 are marginally consistent with the dynamical model if we allow a model uncertainty of 25 km/s, which is also suggested by focusing on only those galaxies whose RVs lie below the model prediction. Thus, postulating that the NGC 3109 association was never bound would still leave us with the task of explaining the anomalously high RVs of NGC 3109, Sextans A, Sextans B, and very probably Leo P. Whatever process is responsible for these HVGs, it is not very common. In the relatively small fraction of cases where the unknown process operates, the resulting HVGs should not end up in the same direction at a similar distance if the HVGs were flung out in individual events. This was discussed in great detail by , who suggested that the HVGs must have been correlated in some way based purely on how they define a thin plane. A correlation becomes almost inevitable when we consider that most if not all of the HVGs are actually located quite close to a line (Pawlowski & McGaugh 2014a).
The most plausible solution is that the NGC 3109 association was once a gravitationally bound group. The mass required to bind it would be rather large, with Bellazzini et al. (2013) estimating that the required mass was 10 11.5 M . Such a high mass could well alleviate the above-mentioned issues regarding the NGC 3109 association, but would also increase the amount of dynamical friction during any close encounter with a massive galaxy. Our results in Figure 9 show that a backsplasher of this mass is highly implausible in a model where galaxies have dark matter haloes that would inevitably create significant dynamical friction during interactions (Privon et al. 2013;Kroupa 2015). Since there are no analogues to NGC 3109 for an assumed mass of just 10 10.5 M , it is clear that the galaxy and the rest of its association pose severe problems for ΛCDM if their high RVs indicate that they are backsplash from the MW or M31, as argued here and in previous works (Teyssier et al. 2012;Pawlowski & McGaugh 2014a;.
So far, we have mostly focused on comparing backsplashers to the present mass of NGC 3109. However, our preceding discussion suggests that it should have been much more massive in the past to bind the NGC 3109 associ-ation. Since mass could be lost during a past encounter with the MW, a conservative approach would be to consider the maximum mass of each backsplasher at any time in its past, as traced by the Illustris merger tree. This is shown in the bottom panel of Figure 15 against the backsplasher's ∆v ef f . It is evident that only one backsplasher with ∆v ef f 150 km/s ever had a mass 10 11.5 M , but it is too close to its host to resemble NGC 3109.

Broader context: the satellite planes challenge
The HVGs in the NGC 3109 association should be considered together with the LG satellite planes because these could all have a common origin, as suggested in Section 7.2. Indeed, Pawlowski & McGaugh (2014a) showed that the high RVs of galaxies in the NGC 3109 association strongly suggest a past close interaction with the MW, even though the association currently lies outside the zero-velocity surface of the LG. For this reason, we do not combine the probabilities of these challenges, but do consider the level of tension with ΛCDM.
Flattened distributions of very likely co-orbiting satellite galaxies around the MW (Lynden-Bell 1976, 1982Kroupa et al. 2005;Pawlowski et al. 2012) and M31 (Metz et al. 2007;Ibata et al. 2013) have long posed a challenge to our understanding of galaxy formation in the ΛCDM context. Recent proper motion data confirm that most of the classical MW satellites do indeed have a common orbital plane (Pawlowski et al. 2013;Pawlowski & Kroupa 2020) aligned with the plane normal defined by the satellite positions alone (Santos-Santos et al. 2020). Their velocities show a very significant bias towards the tangential direction, as occurs for a rotating disc (Cautun & Frenk 2017). Proper motions of two M31 satellite plane members indicate that this structure is likely also coherently rotating (Sohn et al. 2020), as suggested by the RVs of satellites in this nearly edge-on structure (Ibata et al. 2013). After careful consideration of several proposed scenarios for how primordial CDMrich satellites might end up in a thin plane, Pawlowski et al. (2014) concluded that none of them agree with observations for either the MW or M31. Structures as extreme as those observed are exceedingly rare in cosmological simulations (Ibata et al. 2014;Pawlowski & McGaugh 2014b), including hydrodynamical simulations (Ahmed et al. 2017;Shao et al. 2019;Pawlowski & Kroupa 2020) and simulations which model the effects of a central disc galaxy (Pawlowski et al. 2019). The arguments raised by Metz et al. (2009) andPawlowski et al. (2014) against the group infall and filamentary accretion scenarios were later independently confirmed by Shao et al. (2018) using the EAGLE hydrodynamical cosmological simulation (Crain et al. 2015;Schaye et al. 2015). For a recent review of the satellite plane problem, we refer the reader to Pawlowski (2018), who considered both LG satellite planes and the recently discovered one around Centaurus A (Müller et al. 2018(Müller et al. , 2021. To help our discussion, we quantify the level of tension that each challenge represents for ΛCDM, and compare to the one found here. Since we found no NGC 3109 analogues around 13225 hosts, we conservatively assign a frequency of 1/13225 to the HVG challenge. The equivalent number of standard deviations χ corresponding to this frequency can  Table 3. The level of tension between ΛCDM and various challenges it faces within the LG. We have chosen these challenges because they have a common explanation in an alternative framework (Section 7.2). The frequencies for the MW and M31 satellite planes come from section 4.2 of Pawlowski & Kroupa (2020) and figure 2 of Ibata et al. (2014), respectively. Values in the final column are obtained by applying Equation 9 to the frequencies.
be found using We solve this using the Newton-Raphson algorithm with initial guess (3 − log 10 P ) for events with P < 0.001 (the tension is not very significant otherwise, making the initial guess less important for numerical convergence). In this way, we estimate that the HVG challenge corresponds to a 3.96σ event. This is based on allowing the LG mass to lie in the range (2 − 5) × 10 12 M , with consequent implications for the allowed mass range of isolated hosts (Section 3.3) and the isolation condition on backsplashers (Section 5.2). If the LG mass is restricted to the range (2 − 4) × 10 12 M by reducing Mmax in Equation 1, the slightly reduced number of hosts raises the frequency to 1/12187 (3.94σ). If instead we restrict the LG mass to (3 − 5) × 10 12 M by raising Mmin while keeping the nominal Mmax, the number of hosts decreases to 10089, reducing the significance to 3.89σ. Using the nominal LG mass range of (2 − 5) × 10 12 M but focusing on our dark matter-only simulation, we get 12960 hosts, yielding a significance of 3.95σ. In all these cases, the estimated statistical significances should be treated as lower limits because we did not identify any backsplash analogues to NGC 3109. Table 3 summarizes the statistical significance of the HVG challenge for ΛCDM and that of the LG satellite planes. These probabilities are a frequentist interpretation of the "number of trials" for each individual test (e.g. Bayer & Seljak 2020). Since the timing argument analysis of  considered the kinematics of 33 LG non-satellite dwarf galaxies of which NGC 3109 was not the only HVG, the challenge to ΛCDM posed by NGC 3109 is difficult to understand merely via the look-elsewhere effect. Moreover, only around the MW, M31, and Centaurus A do we have information on the 3D distribution of satellite galaxies and at least one component of their velocity. A satellite plane is also evident around Centaurus A, with properties that are likely to arise in ΛCDM only 0.2% of the time (Müller et al. 2021). As a result, it would be difficult to repeat the above-mentioned analyses further afield. In particular, distance uncertainties would be larger for more distant objects, creating significant uncertainty on the peculiar velocity and making it very tricky to do a timing argument analysis. Distance uncertainties also make it difficult to determine the 3D distribution of satellite galaxies, which in addition are very faint and not easy to observe even at the distance of Centaurus A. Thus, the above-mentioned challenges for ΛCDM arise in the only cases where the paradigm can be tested in detail based on the timing argument and the phase space distribution of satellites. One can consider other tests of ΛCDM beyond the LG, some of which we briefly discuss next.

An alternative scenario
Anisotropically distributed satellite galaxies are known to form out of the tidal debris expelled during the interaction between galaxies, as observed in the Antennae (Mirabel et al. 1992). Therefore, the MW and M31 co-orbiting planes of satellite galaxies may have formed as tidal dwarf galaxies (TDGs; Pawlowski et al. 2011), implying a past major interaction. The HVGs in the NGC 3109 association may then be backsplash from this event, even if backsplash events with the required properties do not occur in ΛCDM.
Any second-generation origin for the MW and M31 satellite planes runs into the issue that such satellites would be free of dark matter (Barnes & Hernquist 1992;Wetzstein et al. 2007;Haslbauer et al. 2019) − this was discussed in more detail by Kroupa (2012). But since CDM haloes have never been detected independently of their presumed gravitational effects (e.g. Hoof et al. 2020) and require particles beyond the well-tested standard model of particle physics, it is prudent to consider alternative paradigms without such haloes (e.g. Kroupa 2015). The most promising such paradigm is Milgromian dynamics (MOND; Milgrom 1983). In MOND, the gravitational field strength g at distance r from an isolated point mass M transitions from the Newtonian GM/r 2 law at short range to g = GM a 0 r for g a 0 .
MOND introduces a 0 as a fundamental acceleration scale of nature below which the deviation from Newtonian dynamics becomes significant. Empirically, a 0 ≈ 1.2 × 10 −10 m/s 2 to match galaxy rotation curves (Begeman et al. 1991;Gentile et al. 2011). With this value of a 0 , MOND continues to fit galaxy rotation curves very well using only their directly observed baryonic matter (e.g. Kroupa et al. 2018;Li et al. 2018;Sanders 2019). In particular, observations confirm the prior MOND prediction of very large departures from Newtonian dynamics in low surface brightness galaxies (e.g. de Blok & McGaugh 1997;McGaugh & de Blok 1998). More generally, there is a very tight empirical 'radial acceleration relation' (RAR) between the gravity inferred from rotation curves and that expected from the baryons alone in Newtonian dynamics , with the relation also extending to ellipticals (Lelli et al. 2017;Chae et al. 2020;Shelest & Lelli 2020). This confirms the central prediction of Milgrom (1983). The evidence for MOND on galaxy scales goes beyond the observed tightness of the RAR. For instance, the dynamical friction experienced by galactic bars rotating through a CDM halo is problematic because it would cause the bar to slow down (Debattista & Sellwood 2000), conflicting with observations (Algorry et al. 2017;Peschken & Lokas 2019) − the tension is at the 8σ level (Roshan et al. 2021). In addition, bar-halo angular momentum exchange would cause a resonant effect leading to a quite strong bar after only a few Gyr (Athanassoula 2002), making it difficult to ex-plain rather isolated galaxies like M33 with only a weak bar (Sellwood et al. 2019). This is naturally accounted for in a hydrodynamical MOND simulation of M33, which bears good overall resemblance to observations . The lack of massive CDM haloes and the resulting dynamical friction in close interactions causes a reduced major merger rate (Nipoti et al. 2007;Renaud et al. 2016), which might better explain the high prevalence of thin disc galaxies in the local Universe with little or no bulge (Kormendy et al. 2010;Peebles & Nusser 2010). This continues to challenge the latest ΛCDM cosmological simulations (Peebles 2020). For a review of MOND including its strengths and weaknesses, we refer the reader to Famaey & McGaugh (2012), while Milgrom (2015) provides a more theoretical review.
In the LG, Equation 10 implies a much stronger MW-M31 mutual attraction than the Newtonian inverse square law. As a result, applying MOND to the almost radial MW-M31 orbit (van der Marel et al. 2012(van der Marel et al. , 2019Salomon et al. 2021) implies that they underwent a close encounter 9 ± 2 Gyr ago, as first put forward by Zhao et al. (2013). This is approximately when the MW bar formed and its disc underwent the buckling instability (Grady et al. 2020), which could be due to the interaction if it was 8 − 9 Gyr ago. Due to the high MW-M31 relative velocity around the time of their flyby, they would likely have gravitationally slingshot several LG dwarfs out at high speed. This could well explain the unusually high RV of NGC 3109 − it might have been near the spacetime location of the flyby, thereby gaining a significant amount of energy from the time-dependent LG potential ). Their figure 5 shows that backsplashers from such a highly energetic flyby can easily reach the 1.3 Mpc distance of NGC 3109, and even the 1.6 Mpc distance of Leo P. This is because in addition to the fast MW-M31 relative velocity of ≈ 700 km/s at pericentre, dynamical friction would be greatly reduced as galaxies would not have dark matter haloes (Bílek et al. 2018).
In this scenario, NGC 3109 must have closely approached the MW and/or M31. However, figure 6 of  shows that it is quite possible for the MW-M31 interaction to efficiently slingshot a tracer particle out to > 1.6 Mpc even if it never approached within 40 disc scale lengths of either galaxy. Thus, it is easy to envisage the NGC 3109 association being flung out in this way with negligible dynamical friction on its constituents. It is likely that the association as a whole would be tidally disrupted, such that it is likely unbound today (Kourkchi & Tully 2017). This would make the association analogous to a tidal stream traced by dwarf galaxies rather than stars. However, tidal effects on individual galaxies in the association may have been rather small due to the large pericentric distance and the short duration of any such interaction. Even if there were tidal signatures imprinted at pericentre, the long time since then would make it nearly impossible to identify them today.
During the MW-M31 flyby, tidal tails would likely have formed and might later have condensed into TDGs (Zhao et al. 2013). This phenomenon occurs in some observed galactic interactions like the Antennae (Mirabel et al. 1992) and in MOND simulations of them (Tiret & Combes 2008;Renaud et al. 2016). Due to the way in which such TDGs form out of a thin tidal tail, they would end up lying close to a plane and co-rotating within that plane (Wetzstein et al. 2007;Haslbauer et al. 2019), though a small fraction might well end up counter-rotating depending on the exact details (Pawlowski et al. 2011). The possibility of explaining the LG satellite planes in this way was investigated with MOND Nbody simulations of the MW-M31 encounter (Bílek et al. 2018). Those authors demonstrated the formation of tidal tails connecting the galaxies.  investigated a much wider range of orbital geometries using a restricted N -body approach where the MW and M31 were treated as point masses surrounded by test particle discs. The tidal debris around each galaxy were generally distributed in a thin plane, as evidenced by a sharp concentration of orbital poles. In some models, the preferred direction aligned with the corresponding observed satellite plane for both the MW and M31 (see their figure 5). One reason for this success is that the MW and M31 satellite planes rotate in the same sense, with their orbital poles separated by only ≈ 50 • (Pawlowski et al. 2014). While some mismatch is expected due to the orientations of the MW and M31 discs differing by ≈ 65 • (table 1 of Banik & Zhao 2018), a much larger angle would be difficult to accommodate if both satellite planes condensed out of a common tidal tail.
Since the encounter would have been very long ago, the metallicities and other internal properties of the M31 satellite plane members might be rather similar to those of primordial dwarfs (Recchi et al. 2015), especially in a model where both TDGs and primordial dwarfs are purely baryonic and thus lack any fundamental difference. This might explain the similarity in internal properties between on-and off-plane satellites of M31 (Collins et al. 2015). While those authors interpreted their results as evidence against the TDG hypothesis, field dwarfs (which are presumably mostly primordial in a ΛCDM context) follow a similar mass-radius relation to confirmed TDGs (Dabringhausen & Kroupa 2013). However, a clear splitting is expected in cosmological ΛCDM simulations (figure 12 of Haslbauer et al. 2019). The similarity between primordial and tidal dwarfs is expected in MOND as both would be purely baryonic.
Therefore, the MOND scenario of a past MW-M31 flyby could well explain the LG satellite planes and the high internal velocity dispersions of their members while also accounting for the unusually high RV of NGC 3109 for its position. This should be seriously considered as an alternative to the standard approach of treating the HVG and satellite plane problems as separate statistical flukes in the ΛCDM paradigm (Table 3). It should be the topic of further detailed simulations in a MOND context.

CONCLUSIONS
A detailed 3D Newtonian timing argument calculation of the LG and its surroundings underpredicts the RV of NGC 3109 by 105 ± 5 km/s (table 3 of  . This is despite the significantly more exhaustive search through parameter space described in their section 4.1 compared to the similar analysis of Peebles (2017), who reached similar conclusions. No simple trajectory can be found for these galaxies that respects the Newtonian timing argument and matches available observations at z = 0.
However, the analyses of Peebles (2017) and  are not cosmological simulations. In such a sim-ulation, there could be processes which are not correctly handled by the above-mentioned analyses. In particular, mergers between galaxies can temporarily lead to high relative velocities. A nearby dwarf could then be flung outwards at high speed, possibly explaining the anomalously high RV of NGC 3109. This would make it a backsplasher, as previously suggested by Teyssier et al. (2012) and Pawlowski & Mc-Gaugh (2014a). Using a simplified calculation, we found that this scenario requires an energy gain of ∆v ef f 150 km/s during a past interaction with the MW (Section 5.3). Such trajectories can increase the RV by ≈ 110 km/s compared to a ∆v ef f = 0 trajectory that reaches the same present Galactocentric distance of 1.2 Mpc. Thus, backsplash can in principle explain the anomalously high RV of NGC 3109.
To find out if such trajectories are expected in ΛCDM, we investigated the Illustris TNG300 hydrodynamical cosmological simulation. We identified 13225 host galaxies similar to the MW or M31, and used the merger tree to trace them back in time. At each snapshot, we identified all subhaloes within their virial volume, and traced them forwards as far as possible. Backsplashers are those with a recognizable root descendant at the present epoch that lies beyond 2 rvir from the associated host (Section 5.1).
We found that backsplashers with a larger distance and mass than NGC 3109 are very rare. In the handful of cases where they do occur, ∆v ef f < 0, probably due to dynamical friction. These backsplashers must have received a significant helping hand from large scale structure to reach their present distance, since during the encounter they actually lost energy. However, the timing argument analyses of Peebles (2017) and  include the major galaxy groups outside the LG up to a distance of almost 8 Mpc (table 3 of Banik & Zhao 2017). Therefore, the Illustris cosmological simulation does not reveal trajectories with NGC 3109-like final states that might be significantly mis-modelled by the above-mentioned 3D timing argument analyses. As these neglect dynamical friction, including this process would if anything make it even more difficult to explain the high RV of NGC 3109.
Since the backsplash process concerns the motions of fairly massive galaxies with significant CDM haloes, it should not be affected much by baryonic physics in galaxies. We tested this by comparing our results to the dark matter-only version of TNG300 (Section 6.1). There were many more backsplashers in this case, perhaps due to the lack of strong disruptive tides from e.g. a baryonic disc in the host (Pawlowski et al. 2019). Nonetheless, we found no backsplash analogues to NGC 3109 in the dark matter-only simulation, which yielded a similar overall distance-mass distribution for backsplashers compared to the hydrodynamical TNG300. We therefore conclude that this distribution does not extend to the observed properties of NGC 3109 regardless of precisely how the baryonic physics is treated.
To explain the anomalous kinematics of the NGC 3109 association via the backsplash process, we would need several backsplashers to be flung out in nearly the same direction at a similar time. This strongly suggests that the whole association was once a bound group which closely approached the MW or M31 and was subsequently flung out (Bellazzini et al. 2013). The high mass required for the NGC 3109 group in this scenario renders it infeasible in the ΛCDM context because of the inevitable very strong dynamical friction during the encounter. This is apparent in the lack of sufficiently massive and distant backsplashers in the Illustris TNG300 simulation (Figure 8). The situation remains the same if we trace each backsplasher back in time and consider its maximum mass (bottom panel of Figure  15).
Our null detection of backsplash analogues to NGC 3109 allows us to place an upper limit on their frequency of 1/13225, implying ΛCDM is in > 3.96σ tension with the observed properties of NGC 3109 if it is a backsplasher (Section 7.1). We argue that this is more probable than a 105 km/s error in the timing argument analysis of  for an isolated dwarf galaxy 1.3 Mpc away (Section 4.2) that is also quite far from any major galaxy outside the LG (Table 2). This problem may be related to the phase space correlated distribution of satellite galaxies around the MW and M31, each of whose satellite planes are in 3.55σ tension with ΛCDM (Table 3). These should also be combined with the severe tensions that ΛCDM faces on cosmological scales with regards to the locally measured expansion rate (Riess 2020;Di Valentino 2021), the unusually low matter density within 300 Mpc (Keenan et al. 2013;Haslbauer et al. 2020), and the too-rapid formation of observed galaxy clusters like El Gordo .
We therefore propose an alternative scenario in which the unusual kinematics of the NGC 3109 association might bear witness to a past close MW-M31 flyby in the MOND context (Section 7.2). Tidal debris from the flyby could have formed into the LG satellite planes Bílek et al. 2018). Fitting this picture into a broader cosmological context (as suggested by Haslbauer et al. 2020) would require a relativistic MOND theory such as that of Skordis & Z lośnik (2019), which may well enhance the growth of structure sufficiently to address the above-mentioned issues.

DATA AVAILABILITY
The data underlying this article are available in the article. The algorithms used to download and analyze the Illustris simulations have been made publicly available. 1