The Relevance of Dynamical Friction for the MW/LMC/SMC Triple System

Simulations of structure formation in the standard cold dark matter cosmological model quantify the dark matter halos of galaxies. Taking into account dynamical friction between the dark matter halos, we investigate the past orbital dynamical evolution of the Magellanic Clouds in the presence of the Galaxy. Our calculations are based on a three-body model of rigid Navarro-Frenk-White profiles for the dark matter halos, but were verified in a previous publication by comparison to high-resolution N-body simulations of live self-consistent systems. Under the requirement that the LMC and SMC had an encounter within 20 kpc between 1 and 4 Gyr ago, in order to allow the development of the Magellanic Stream, and using the latest astrometric data, the dynamical evolution of the MW/LMC/SMC system is calculated backwards in time. With the employment of the genetic algorithm and a Markov-Chain Monte-Carlo method, the present state of this system is unlikely with a probability of<10^{-9} (6 sigma complement), because solutions found do not fit into the error bars for the observed plane-of-sky velocity components of the Magellanic Clouds. This implies that orbital solutions that assume dark matter halos according to cosmological structure formation theory to exist around the Magellanic Clouds and the Milky Way are not possible with a confidence of more than 6 sigma


Introduction
The standard ΛCDM cosmological model [1][2][3] requires about 81 per cent of the matter-content of the Universe to be comprised of pressure-less cold dark matter (CDM) particles which are not described by the standard model of physics, and about 68 per cent of the total mass-energy density of the universe to be composed of dark energy represented by the cosmological constant Λ [4,5].The related ΛWDM model is based on dark matter (DM) being made of warm dark matter particles that have a smaller mass than CDM particles but lead to very similar DM halos within which the observed galaxies reside (e.g.[6][7][8][9]).Detecting DM particles has become a major world-wide effort which has so far yielded null results (e.g.[10][11][12][13][14][15] and references therein).Finding evidence for DM particles assuming these have a finite interaction cross section with standard-model particles may lead to null detection if the interaction cross section is too small to be measurable.However, arguments based on the observed strong correlations between standard particles and dark matter suggest a significant cross section [16].This indicates an impasse in the experimental verification of the existence of DM particles.
An alternative and robust method to establishing the existence of DM particles is to develop techniques which rely solely on their gravitational interaction with ordinary matter.A suitable approach is given by Chandrasekhar dynamical friction (e.g.[5]): when a satellite galaxy with its own DM halo enters the DM halo of a host galaxy, it's orbit decays as a result of dynamical friction and the satellite merges with the host.The decay of the orbit does not depend on the mass of the DM particle but only on the mass density of DM particles which is fixed by the cosmological parameters.This is the primary reason for the existence of the merger tree in standard cosmology (e.g.[17,18]), and for interacting galaxies to be thought of as merging galaxies.Effectively, a DM halo, being 10-20 times more extended and 50-100 times as massive as the observable part of a galaxy, works like a spider's web.In contrast, if there were to be no DM halos around galaxies, then the galaxygalaxy encounters would be significantly less dissipative, galaxies would encounter each other multiple times and mergers would be rare.In the eventuality that DM halos were not to exist, however, non-relativistic gravitational theory would need to be non-Newtonian and, given the correlations that galaxies are observed to obey, this theory would need to be effectively Milgromian [19][20][21][22].Galaxy-galaxy encounters and mutual orbits in Milgromian dynamics and without DM halos have been studied (for example the interacting Antennae galaxy pair [23], and the Milky Way-Andromeda binary explaining the mutually correlated planes of satellites around both hosts [24][25][26]).The law of gravitation fundamentally defines the formation and evolution of galaxies.For example, by mergers being much rarer in Milgromian gravitation (i.e.without DM), galaxies evolve largely in isolation [27][28][29], and the formation of cosmological structures proceeds differently than in the DM-based models [5,30,31].It is thus of paramount importance to test for the existence of DM halos around galaxies.
The test based on Chandrasekhar dynamical friction for the existence of DM halos has been introduced by Angus, Diaferio & Kroupa [32] by addressing the question if the present-day Galactocentric distances and motion vectors of observed satellite galaxies of the Milky Way (MW) conform to their putative infall many Gyr ago.The solutions for those satellite galaxies for which proper motions were available imply tension with the existence of DM halos since no in-fall solutions were found.In contrast, without the DM component, the satellite galaxies would be orbiting about the MW having been most likely born as tidal-dwarf galaxies during the MW-Andromeda encounter about 10 Gyr ago [25,26,33].A further test for the presence of the DM component applying dynamical friction was achieved by Roshan et al. [34] who found the observed bars of disk galaxies to be too long and rotating too fast in comparison with the theoretically expected bars in the presence of DM halos that absorb the bar's angular momentum.The reported discrepancy amounts to significantly higher than the 5 sigma threshold such that the observations are incompatible with the existence of DM halos.Another independent application of the dynamical friction test is available on the basis of the observed distribution of matter within the M81 group of galaxies.The extended and connected tidal material implies multiple past close encounters of the group members.But the dynamical evolution of the M81 group of galaxies is difficult to understand theoretically if the galaxies are contained in the DM halos that are expected in standard cosmology ( [35] and references therein).The problem is that in the presence of DM halos the group merges too rapidly to allow the tidal material to be dispersed as observed.
The interesting aspect of these results is that the independent three analyses of the orbital dynamics of the MW satellite galaxies, galactic bars and of the M81 group agree in the conclusion that the data are difficult if not impossible to understand in the presence of DM halos.Because the implications of the non-existence of DM bears major implications for theoretical cosmology and galaxy formation and evolution, further tests are important to ascertain the above conclusions, or indeed to question them.
Given the high-accuracy and high-precision position and velocity data available today through the astrometric Gaia mission, the problem of how the existence of the massive and extended DM halos in which galaxies reside can be tested for is revisited here using the process of Chandrasekhar dynamical friction on the triple-galaxy system comprised of the MW, the Large (LMC) and Small (SMC) Magellanic Clouds including the Magellanic Stream.This system is comparable to the M81 system as it too has a large gaseous structure stemming from the Magellanic Clouds, the Magellanic Stream, which constrains the past orbits of the components of the system.In the following we assume that the standard cosmological model is correct and that each, the MW, the LMC and the SMC are contained in DM halos that are consistent with those obtained from the theory of structure formation in the ΛCMD model of cosmology.That is, we associate with the baryonic mass of each galaxy a DM halo profile as predicted by the theory in order to test the theory.
Given that the LMC and SMC are about 50 and 60 kpc distant from the centre of the MW and about 20 kpc distant from each other and that the radii of the DM halos are such that all three galaxies are immersed in the dark matter halo of the MW and that the SMC is immersed in the dark matter halo of the LMC and vice versa, it is apparent that Chandrasekhar dynamical friction is likely to play a very significant role in establishing the orbits of the three galaxies relative to each other.Recent work ( [36][37][38][39][40][41]) suggests the Magellanic Clouds to be on their first pericentre passage such that the bulk properties of their putative DM halos will not be significantly stripped.Detailed dynamical modelling of the Magellanic Clouds problem has, until now, not come to the conclusion that there is a problem.
For example, Besla et al. ( [37]) performed simulations of the LMC and SMC in the first-infall scenario by assuming that dynamical friction from the MW DM halo can be neglected.This is not correct and significantly helps the LMC/SMC pair to exist longer before merging.Indeed, Lucchini et al. ( [42]) study the formation of the Magellanic Stream using hydrodynamical simulations in the first-infall scenario for the LMC and SMC.Their results, obtained in fully live DM halos of all components such that dynamical friction is correctly computed self-consistently, show that the LMC/SMC binary orbit shrinks significantly more rapidly (their fig. 1) than calculated by Besla et al. ([37]).But Lucchini et al. adopted DM halo masses of the LMC (1.8 × 10 11 M ⊙ ), of the SMC (1.9 × 10 10 M ⊙ ) and of the MW (1.0 × 10 12 M ⊙ ) that are significantly less massive than predicted by the ΛCDM model (see Table 1 below).Their calculations also do not consider if the in-falling LMC/SMC binary would have survived for sufficiently long before infall because it is implausible for the present-day LMC/SMC binary to have formed during or shortly before infall, and if the infall velocity vector is consistent with the Hubble flow 3.46 Gyr ago at the start of their simulation.Vasiliev ([41]) studies a scenario in which the LMC is on its second passage past the MW but ignores the SMC.The calculated orbits of the LMC decay due to dynamical friction consistently to the results presented here.Thus, this past work either neglects important contributions to dynamical friction, members in the problem or uses light-weight DM halos that are not consistent with the ΛCDM cosmological model.Essentially, the past work has demonstrated that the MW/LMC/SMC triple system can broadly be understood in the context of orbital dynamics by effectively suppressing Chandrasekhar dynamical friction, but it does not demonstrate that the calculated orbits are consistent with the standard DM-based cosmological theory.
In order to assess if the standard DM-based models can account for the existence of the observed MW/LMC/SMC triple system, the present work studies their orbital history by strictly constraining the DM halo of each of the three involved galaxies to be consistent with the ΛCDM-allowed DM halo masses.The dark-matter component is thus assumed to be pressureless dust.The calculations furthermore impose the condition that the LMC and SMC had an encounter within 20 kpc of each other between 1 and 4 Gyr ago in order to allow the Magellanic Stream to form, and they take into account cosmological expansion to constrain the relative positions of the galaxies 5 and more Gyr ago.To check for consistency of the solutions, two independent statistical methods are applied to search for all allowed initial conditions within the error bars for the transverse velocity components of the LMC and SMC.The two tests yield indistinguishable results.The calculations based on rigid Navarro-Frenk-White (NFW) DM halo profiles are conservative because equivalent simulations with live DM halos lead to faster merging [35].These calculations show there to be no orbital solutions and that the observed system comprised of the MW/LMC/SMC plus Magellanic Stream cannot be understood in the presence of DM halos in the context of the ΛCDM model.The problem which gravitational theory leads to the observed system is left for a future contribution.

The Model 2.1. NFW Profiles
When testing a theory for consistency with data it is of paramount importance to not mix the two: purely theoretically calculated properties need to be compared with empirical data that have not been modulated by the very theory to be tested, and vice versa.It is therefore not permissible to choose sub-massive DM halos that allow a solution in order to argue that ΛCDM or ΛWDM theory account for the MW/LMC/SMC system.Thus, here we are adamant at requiring to use the theoretically predicted NFW dark matter halo profiles ( [43]).These arise from self-consistent cosmological structure formation simulations such as documented in [44] for the range of baryonic galaxy masses considered here.More generally, the NFW DM halo profile is a standardised outcome of structure formation simulations in ΛCDM theory (e.g.[45]), and cold or warm DM halos have been shown to have similar overall profiles and densities [6][7][8][9].While the inner and outer DM halo profiles can differ from the NFW form, the latter well-represents the distribution of theoretically-predicted DM particles around the baryonic component of a galaxy especially for galaxies that are on a first-infall orbit (e.g.[46][47][48][49]).In seeking orbital solutions the stellar masses are allowed to vary within ±30 per cent of their nominal values with the correspondingly different DM halo masses (Table 1 below).The range of theoretically predicted DM halo masses (−41 to +68 per cent for the MW, −42 to +14 per cent for the LMC and −17 to +16 per cent for the SMC) is thus accounted for.
The DM halo of either galaxy is treated as a rigid halo with a density profile according to [43] (NFW-profile), truncated at the radius R 200 : with R s = R 200 /c, R 200 denoting the radius yielding an average density of the halo of 200 times the cosmological critical density, and the concentration parameter c [50], The DM halo masses are derived from the stellar masses of the galaxies by means of fig.7 of [44].Note that using rigid DM halo profiles in such orbital computations is admissible because these have been verified to give conservative solutions (slower merging times) than self-consistent simulations with live DM halos [35] and in particular because the MW/LMC/SMC system has only a recent encounter history.

Dynamical Friction and Equations of Motion
Exploring the dynamics of bodies orbiting in the interior of DM halos implies that dynamical friction 1 has to be taken into account in an appropriate manner [51].Here the formulation of Chandrasekhar's dynamical friction as used in [35], their eq. 1 and 2, is applied.
Chandrasekhar's formula gives a quick-to-compute estimate for orbital decay which is needed for the sake of establishing statistical statements about merger rates between galaxies.High-resolution simulations of live self-consistent systems confirm our approach of employing this semi-analytical formula in our three-body calculations (see sec. 7 as well as figures 13 and 14 in [35]).Using the computationally significantly more time intensive live simulations in fact leads to more rapid orbital decay such that our here-employed semi-analytical approach is conservative by allowing a larger range of solutions than would be the case if live DM halos were to be used.
The equations of motion for the individual galaxies are as given in appendix C of [35] 2 , and are here augmented with an additional term taking into account the Hubble flow as follows: Based on the assumption of a flat universe curvature parameter k = 0) we extended the equations of motion by the cosmic acceleration term s s ⃗ r i caused by the Hubble flow of the expanding universe.Here s(t) is the scale-factor of the universe and ⃗ r i is the position of a galaxy in the centre-of-mass frame of the group.In Cartesian coordinates the equations of motion are then given by (k = 1, 2, 3 for the MW, LMC and SMC, respectively). with Here ⃗ F i is the total force acting on galaxy i, V ij is the potential energy between the galax- ies i and j, and ⃗ F DF ij is the dynamical friction force acting on galaxy i caused by the overlap of the DM halos of galaxy i and j (according to actio est reactio ⃗ F DF ji is taken into account, too.) Assuming a flat dark energy dominated cosmology (Ω m,0 + Ω Λ,0 = 1), the second Friedmann equation becomes where Ω m,0 and Ω Λ,0 are the matter and dark energy densities at the present time, respectively, scaled by ρ crit , and H 0 is the Hubble constant.Setting the conditions s(0) = 0 and ṡ = H 0 at s(t) = 1 (present time) yields an analytical expression for the cosmic-scale factor where t is the age of the Universe.For our calculations we used the following values: Ω m,0 = 0.315, Ω Λ,0 = 0.685 and H 0 = 67.3km/s/Mpc.However, the cosmic acceleration term does not play a significant role for the time frame considered here.The order of magnitude contribution, compared to the forces between the DM halos in Eq. 4, is in the range from 10 −3 to 10 −2 within the time range [−5 Gyr, today].All numerical computations were performed using SAP's ABAP development workbench.

Observational Data
For the sake of easy accessibility, we list the observational constraints in Tables 1-7, and the references regarding the basic observational data in Table 2.
As mentioned in Sec.2.1 we derived the DM halo masses of MW, LMC and SMC according to [44], based on the stellar masses extracted from the references given in Table 2. Our approach, therefore, is in full correspondence with the predictions of ΛCDM regarding the DM halo masses.However, in order to overcome possible uncertainties in that regard, and to achieve more independent results, we decided to apply our statistical evaluations of Sec. 4 individually to 27 mass combinations by varying the stellar mass of each galaxy by plus minus 30%, as displayed in Table 1.The 27 mass combinations are then denoted by the indicators "o" for the original masses, "m" for −30%, and "p" for +30%, for instance "o-m-p" for the MW/LMC/SMC mass triple.The relevant coordinates and velocities are presented in Tables 3, 4 and 5 and the transformation to the Cartesian heliocentric equatorial system is shown in Table 6 and 7.
The Magellanic Stream consists of H I gas which trails behind both the LMC and the SMC across a large fraction of the sky.It is thought to be the result of a combination of tidal forces and ram pressure stripping through the orbit of the LMC and SMC within the hot gaseous halo of the MW.Studies of the origin of the Magellanic Stream have shown that it was most likely created when the LMC and the SMC had a close encounter during which some of the gas of the LMC and SMC became less bound to be subsequently removed from the pair through ram pressure stripping [39,40,62,63].

Statistical Methods
In order to cross-check the solutions and to enhance the confidence in these we utilize two independent statistical methods, a Markov-Chain Monte-Carlo method (MCMC, see Sec. 4.1) and the genetic algorithm (GA, see Sec. 4.2).The reason, why we do so, is given in the introductory statement of Sec.4.2.This is the initial situation for each method MCMC and GA: Due to the error bars of the transverse velocity components in RA and DEC for the Magellanic Clouds (see Table 4) we have to consider four open parameters P i regarding the calculation of the three-body orbits of the dark matter halos, as displayed in Table 8.This means that employing the MCMC method and the GA method separately for each of the 27 mass combinations (see Sec. 3), we search for solutions of the three-galaxy orbits with best fits to the transverse velocity components of the Magellanic Clouds.It is important to note that our goal is to achieve results independent of the particular masses of the galaxies, and not to find a best fit mass combination.
Further, based on the radio-astronomical observational data concerning the Magellanic H I -Stream explained in Sec. 3, we specify the following broad condition we here from refer to as condition COND regarding admissible past orbits of the Magellanic Clouds: LMC and SMC encountered each other within the past time interval of [−4 Gyr, −1 Gyr] at a pericentre distance of less than 20 kpc.Incorporating this condition, the algorithms searched for solutions by integrating Eq. 4 backwards in time up to −5 Gyr.

Markov-Chain Monte-Carlo (MCMC)
We follow a methodology proposed by [64] employing an affine-invariant ensemble sampler for the Markov-Chain Monte-Carlo method.A detailed guideline for an imple-  mentation can be found in [65].The basics of applying this formalism to our situation are outlined in appendix D of [35].

Definition of the Posterior Probability Density
First of all, concerning the open parameters, we need to account for the error bars of the transverse velocity components v i of the Magellanic Clouds.This is ensured by an appropriate definition of the prior distribution p(X) where X symbolizes the parameter vector (P 1 ...P 4 ).With the values from Table 8, the four contributions to the prior distribution read (i = 1, ..., 4) Exploiting the minimal distance d 23 between LMC and SMC within the time period [−4 Gyr,−1 Gyr], the condition COND implies the likelihood function with d per = 20 kpc and d 0 = 5 kpc.The posterior probability density is then given by the product of Eqs. 8 and 9: Here the normalizing constant for an overall probability of 1 is neglected because it is clear from the comparative search algorithm that the absolute values of π(X) are irrelevant.
with varying n, to cater for extended error bars as explained in Sec.4.3.

Genetic Algorithm (GA)
General aspects of the genetic algorithm are explained in detail in [66].As pointed out there, the major advantage of this method is perceived to be its capability of avoiding local maxima in the process of searching the global maximum of a given function (here the fitness function).This feature makes it worthwhile to employ the GA method as a second independent statistical approach besides the MCMC method.
A precise description of the algorithm, especially instructions for the implementation, can also be found in [67].To put it in a nutshell: The open parameters are related to genes, which are concatenated to genotypes.The set of a number of genotypes is considered to be a population where the members pairwise produce a follow-up generation with new features due to randomly performed cross-over and mutation.Winners per generation are determined by means of the fitness function, and by comparison between the generations, the overall winner is found.
Comparing to the more familiar MCMC method we have the following parallels: walker <-> genotype, ensemble <-> population, set of ensembles <-> generations, posterior probability density <-> fitness function.

The GA Generations
Each open parameter from Table 8 is mapped to a 4-digit string ("gene") [abcd] i (i ∈ {1, ..., 4}), again with varying n like in Sec.4.1.2,in order to cater for extended error bars as explained in Sec.4.3.All genes together define the genotypes as 4×4-digit strings, [abcd] 1 ...[abcd] 4 , which generate a population of N pop genotypes.The first generation is established by randomly creating N pop 4×4-digit strings.

Definition of the Fitness Function
Our goal is to establish GA evaluations that are compatible to the evaluations by means of the Markov-Chain Monte-Carlo method.Therefore we choose the definition of the fitness function to be identical to the definition of the posterior probability density (Eq.10), where the individual components are taken from Eq. 8 and 9.

Results
Both statistical methods, MCMC and GA, deliver mutually consistent, and in fact indistinguishable results.

Step I
As a first step, we tried to find solutions for which all four open parameters P i lie within the 1σ error bars of the transverse velocity components of the Magellanic Clouds (see Table 8).Employing the methods MCMC and GA as search engines, using 1000 ensembles with 100 walkers (MCMC) and accordingly 1000 generations with a population of 100 genotypes (GA), we repeated the search 100 times for both statistical methods for each of the 27 mass combinations according to Table 1.This attempt, to find a solution, failed.That is, within the here probed 1σ uncertainty range of the transverse velocity components, neither algorithm found orbits that fulfil the condition COND.In other words, the MW/LMC/SMC plus Magellanic Stream system cannot exist in its observed configuration in the presence of dark matter halos.
To proceed further, we gradually extended the allowed intervals for the parameters P i to be v i ± nσ i with increasing integers n, see Eqs. 11 and 12. Neither the MCMC nor the GA statistical method found a solution for σ = 1, 2, 3.For σ = 4 both methods delivered solutions for a single mass combination only, namely p-m-p.The details regarding all 27 mass combinations are given in App. A.

Step II
How to go about quantifying the first results obtained in Sec.4.3.1?First of all, for a mass combination with first solutions based on nσ error intervals we choose the (n + 1)σ intervals to be allowed for the open parameters.This is motivated by the thought that relaxing the overall nσ condition for all open parameters may result in solutions with individual smaller-than 1 σ deviations of the parameters.
Furthermore, we constructed a "probability-sigma-grid" in the following sense: If a parameter P i lies within the error bar (< 1σ i ) then its probability weighting factor is 1, a conservative cautious setting.If a parameter P i lies outside the error bar then its deviation from v i is weighted with the remainder of the probability function, based on 0.1σ intervals.For instance, if we have P i = v i + 3.74σ i then the corresponding probability weighting factor for this parameter is calculated according to the remainder of the probability function for 3.7σ.
Based on the same search method as in Sec.4.3.1,we searched for best-fit solutions by calculating an overall probability for the combination of the open parameters.The results are displayed in Tables 9 and 10 for the methods MCMC and GA, respectively.Regarding the MCMC method, it is interesting to note that in some cases individual parameters are not confined to the mentioned (n + 1)σ intervals because the stretch moves can place walkers outside of the intervals specified by Eq. 11, thus supporting the overall process of maximizing the posterior probability density consisting of all open parameters.This is not possible for the GA method as the genotypes are apriori confined to the preset intervals.
The result of Sec.4.3.1 that the mass combination p-m-p provides the best-fit solution is confirmed by either statistical method, MCMC and GA.Furthermore, Fig. 1 shows the trajectories, displayed as pairwise distances, of the individual best-fit solutions for the mass combinations o-o-o and p-m-p.Even in that detail, both methods deliver identical results.Table 9. Results using the MCMC method: Probabilities regarding the evaluation of the error intervals of the plane-of-sky velocity components of LMC and SMC for the individual best-fit solutions of each mass combination considered, as explained in Sec.4.3.2.

Combination Deviation Deviation Deviation Deviation Probability of masses
of Table 10.Same as for Table 9, but for the GA method.
Combination Deviation Deviation Deviation Deviation Probability of masses of  2), and creating sets of 1000 solutions in each case 3 .For the GA method we undertook broader search runs establishing 1000 solutions for the mentioned mass combinations.
The autocorrelation functions are presented in Figs. 2 and 3, demonstrating that convergence is achieved quickly for the MCMC method, and the improved probabilities are displayed in Table 11 and 12 for the MCMC and GA methods, respectively.the uncertainty range was therefore increased to ±(n + 1) σ with the result that values of n are needed to obtain viable orbital solutions that lie outside the 5 σ confidence range of the quantities such that the orbital solutions are likely with probabilities smaller than 10 −9 (Tables 9 and 10).
Step III (Sec.4.3.3)demonstrates that the MCMC and GA methods yield indistinguishable results, thus confirming the conclusion that orbital solutions do not exist for the observed MW/LMC/SMC plus Magellanic Stream system in the presence of the theoretically expected DM halos.The dynamical behaviour of the MW/LMC/SMC system demonstrates the importance of dynamical friction in the context of interacting galaxies, as dynamical friction significantly influences the orbits.Fig. 4 shows that the forces due to dynamical friction between the LMC and the SMC are comparable to the pure gravitational force between the overlapping DM halos.
As a thought experiment, we calculated the orbits back in time to −7 Gyr with and without dynamical friction for the overall best-fit solution, derived for the mass combination p-m-p, and for the best-fit solution for the original mass combination o-o-o.Dynamical friction can be turned of by omitting the terms ⃗ F DF in Eq. 5.This retains the gravitational pull of the DM halo but avoids the deceleration through Chandrasekhar dynamical friction.This Gedanken experiment is thus a rough approximation of the situation in Milgromian dynamics according to which a galaxy generates a Milgromian gravitational potential that can be viewed as stemming from a Newtonian plus phantom DM halo that does not generate dynamical friction.With this approximation to MOND, the GA method as the search engine immediately delivers 17 mass combinations already within the 1σ intervals (see also App.A).In other words, solutions matching the error intervals for the transverse velocities of the LMC and SMC are found readily without Chandrasekhar's dynamical friction but with the potential generated by the DM halo.The result is displayed in Fig. 5.With the absence of Chandrasekhar dynamical friction on DM halos the Magellanic Clouds would have a long orbital lifetime as satellites of the MW, possibly being massive tidal dwarf galaxies formed during the MW-Andromeda encounter about 10 Gyr ago.
These results thus suggest that the MW/LMC/SMC plus Magellanic Stream system may have a straight-forward orbital solution in Milgromian dynamics, and it is thus of much interest to simulate this system and its possible origin in this non-Newtonian framework.

Conclusions
The conclusions reached by this analysis are very robust, since both the MCMC and GA algorithms lead to indistinguishable results.Taking the observed configuration in six-dimensional phase space of the MW/LMC/SMC plus Magellanic Stream system as the necessary boundary condition, it is impossible to find orbital solutions backwards in time that fulfil the very liberal condition COND (Sec.4).The orbits accelerate too rapidly (backwards in time) such that the LMC and SMC could not have remained bound long enough to have the required encounter that is needed to have occurred to produce the Magellanic Stream.In other words, the system merges too rapidly forward in time to allow a close encounter as defined through COND and to still be visible today as two distinctly separated galaxies next to the MW.Solutions do not even appear if the DM halo masses of the three galaxies are allowed to be larger or smaller by up to 42 per cent than those of the standard assumptions (see Table 1).The possibility that the LMC and SMC fell into the MW DM halo independently of each other but at a similar time in order to allow them to pair up to the observed binary is arbitrarily unlikely because the LMC and SMC would have to have had relative velocities to each other and to the MW that oppose the Hubble flow.Such an unlikely solution is in any case not possible because the condition COND requires the LMC and SMC to have had a close encounter at a time in the past such that the two would have merged today due to the mutual dynamical friction on their respective DM halos.In any case, neither the MCMC nor the GA methods found such solutions.This work thus shows that the observed configuration of the MW/LMC/SMC plus Magellanic Stream system is not possible in the presence of the theoretically expected DM halos.
The results based on the Chandrasekhar dynamical friction test applied to the MW/-LMC/SMC triple system arrived at here corroborate the previous evidence based on the same test but other systems, noted in the Introduction, that question the existence of dark matter particles.The independently documented problems [5,[68][69][70][71][72][73][74] of fitting the standard model of cosmology to the observed Universe on most probed scales is consistent with these results.
The thought experiment in Sec.4.3.4 in which the potentials of the DM halos are kept but the Chandrasekhar dynamical friction term is set to zero naturally leads to solutions.This experiment is an approximation to the situation in Milgromian dynamics and demonstrates that the origin and evolution of the MW/LMC/SMC plus Magellanic Stream system needs to be studied in this non-Newtonian framework.
There is a printing mistake in the equation for Φ i (s i ), the third equation of appendix C in [35]: the − sign in the first row is incorrect.
3 As explained in [35] the autocorrelation functions deliver an estimate for the convergence of the search algorithm for sufficiently large numbers of ensembles.Therefore, in order to establish full convergence for the approximate calculation of the autocorrelation functions, we created 50 000 ensembles for both mass combinations considered, o-o-o and p-m-p.The mentioned 1000 solutions fulfilling the conditions (n + 1)σ intervals and COND were extracted from follow-up runs based on the set of walkers of ensemble 50 000 in either mass combination case.

Figure 1 .
Figure1.Orbits, displayed as pairwise distances, calculated backwards in time up to −7 Gyr for the best fit solutions for the mass combinations o-o-o and p-m-p for either statistical method MCMC and GA.Left panel: Orbits obtained using the MCMC method.Right panel: Orbits obtained using the GA method.4.3.3.Step IIIEspecially due to the MCMC method not being utilized as a pure search engine only, we established its validity by checking the results obtained in the previous Section 4.3.2 in terms of the autocorrelation function for the mass combinations o-o-o (original masses) and p-m-p (mass combination with the overall best-fit solutions according to Sec. 4.3.1 and 4.3.2),and creating sets of 1000 solutions in each case3 .For the GA method we undertook broader search runs establishing 1000 solutions for the mentioned mass combinations.The autocorrelation functions are presented in Figs.2 and 3, demonstrating that convergence is achieved quickly for the MCMC method, and the improved probabilities are displayed in Table11and 12 for the MCMC and GA methods, respectively.

4. 3 Figure 2 .Figure 3 .
Figure 2. The MCMC method: Autocorrelation functions related to the open parameters for the primordial mass combination o-o-o.The calculation is based on a set of 50 000 ensembles.

Figure 4 .Figure 5 .
Figure 4. Forces between the Magellanic Clouds, based on today's observational data with a distance of 22.5 kpc between the CM's of the dark matter halos, and the present-day relative velocity of 89.8 pc/Myr (87.8 km/s), shown as the vertical dotted line.

Table 2 .
List of references for the basic observational data.

Table 3 .
Observational data for LMC and SMC, and in parts for the Galactic centre.

Table 5 .
Circular velocity and proper motion of the Sun (Galactic coordinates).

Table 6 .
Cartesian equatorial coordinates and heliocentric equatorial velocities for the Magellanic Clouds and the centre of the Galaxy.

Table 7 .
Transformation of the error bars in RA and DEC for the transverse velocity components of the Magellanic Clouds to Cartesian equatorial coordinates.

Table 8 .
Open parameters for the statistical methods and the corresponding 1σ uncertainty values from the observational data.∈ {0, ..., 1000} separately for each walker, and for each open parameter of a given walker, the first ensembles are created according to (i ∈ {1, ..., 4}) RA ) (LMC: v DEC ) (SMC: v RA ) (SMC: v DEC )

Table 11 .
Results from the MCMC method: Probabilities regarding the evaluation of the error intervals of the plane-of-sky velocity components of the LMC and SMC for the individual best-fit solutions of the mass combinations o-o-o and p-m-p, based on sets of 1000 solutions as explained in Sec.4.3.3.

Table 12 .
Same as for Table11, but for the GA method.