It is commonly believed that our own Milky Way is on a collision course with the neighbouring Andromeda galaxy. As a result of their merger predicted in around five billion years the two large spiral galaxies that define the present Local Group would form a new elliptical galaxy. Here we consider the latest and most accurate observations by the Gaia and Hubble space telescopes along with recent consensus mass estimates to derive possible future scenarios and identify the major sources of uncertainty in the evolution of the Local Group over the next 10 billion years. We find that the next most massive Local Group member galaxies  namely M33 and the Large Magellanic Cloud  distinctly and radically affect the Milky Way  Andromeda orbit. While including M33 increases the merger probability the orbit of the Large Magellanic Cloud runs perpendicular to the Milky Way  Andromeda orbit and makes their merger less likely. In the full system we find that uncertainties in the present positions motions and masses of all galaxies leave room for drastically different outcomes and a probability of close to 50 that there is no Milky Way  Andromeda merger during the next 10 billion years.
['Astrophysics - Astrophysics of Galaxies']
Apocalypse When No Certainty of a Milky Way  Andromeda Collision
{'Apocalypse When? No Certainty of a Milky Way Andromeda Collision': "Till Sawala 1,2* , Jehanne Delhomelle 1,3 , Alis J. Deason 2 , Carlos S. Frenk 2 , Peter H. Johansson 1 , Atte Keitaanranta 1 , Alexander Rawlings 1 , Ruby Wright 1  \n1 Department of Physics, University of Helsinki, Gustaf Hallstromin katu 2, Helsinki, FI-00014, Finland. 2 Institute for Computational Cosmology, Durham University, South Road, Durham, DH1 3LE, United Kingdom. Universit'e Toulouse III - Paul Sabatier, 118 Route de Narbonne,  \n3 Toulouse, F-31062, France.  \n*Corresponding author(s). E-mail(s): till.sawala@helsinki.fi;", 'Abstract': "It is commonly believed that our own Milky Way is on a collision course with the neighbouring Andromeda galaxy. As a result of their merger, predicted in around five billion years, the two large spiral galaxies that define the present Local Group would form a new elliptical galaxy. Here we consider the latest and most accurate observations by the Gaia and Hubble space telescopes, along with recent consensus mass estimates to derive possible future scenarios and identify the major sources of uncertainty in the evolution of the Local Group over the next 10 billion years. We find that the next most massive Local Group member galaxies - namely, M33 and the Large Magellanic Cloud - distinctly and radically affect the Milky Way - Andromeda orbit. While including M33 increases the merger probability, the orbit of the Large Magellanic Cloud runs perpendicular to the Milky Way - Andromeda orbit and makes their merger less likely. In the full system, we find that uncertainties in the present positions, motions, and masses of all galaxies leave room for drastically different outcomes, and a probability of close to 50% that there is no Milky Way - Andromeda merger during the next 10 billion years.  \nThe Local Group (LG) contains two large spiral galaxies, our own Milky Way (MW) and the Andromeda galaxy (hereafter M31), along with approximately 100 known smaller galaxies [1]. In addition, it likely hosts additional galaxies yet to be discovered [2] and, according to the standard cosmological model, a vast number of completely dark substructures [3].  \nThe negative radial velocity of M31 towards the MW has been known for over a century [4], even before its distance was first accurately measured [5]. However, while indirect methods had since been used to constrain the transverse components of M31's velocity vector [6, 7, 8], direct measurements of the minute proper motions were only achieved much more recently, with the Hubble Space Telescope (HST) [9].  \nThe first numerical studies [10] of a possible MW-M31 merger predate even the early estimates of the transverse velocity. The finding that the MW-M31 motion is close to radial immediately led to the prediction of a likely future collision and merger [11, 12, 13]. This scenario has since become the prevalent narrative [14, 15] and textbook knowledge [16, 17].", 'Predicting the future of the Local Group': 'The MW and M31 both contain remnants of past mergers and interactions with other galaxies [18, 19, 20, 21]. Predicting future mergers requires knowledge about the present coordinates, velocities, and masses of the systems partaking in the interaction. In addition to the gravitational force between galaxies, dynamical friction is the dominant process in the lead-up to galactic mergers. It describes a transfer of orbital kinetic energy to internal energy of the objects involved, and consequently leads to the decay of galactic orbits.  \nIn this study, we parameterise the density profile of each halo as isotropic NavarroFrenk-White profiles [22], use standard assumptions about their concentrations and velocity dispersion profiles [23] and calculate dynamical friction using an analytic formalism (see Methods for a detailed description). While non-gravitational processes, such as gas drag and star formation leading to increased central densities, etc., also shape the final stages of a merger, the phase of the orbit that defines the occurrence of mergers is largely determined by gravity, which in turn, is dominated by the dark matter component in the standard cosmological model [24].  \nThe simplest model for the MW-M31 orbit, as considered by [11], contains only the two main galaxies. Due to the planar symmetry, only five parameters are required: the two masses, along with the initial separation and the two velocity components. In general, an orbit with N > 2 galaxies requires specifying 3 × ( N -1) coordinates and 3 × ( N -1) velocity components, along with the N masses. More recently, [12] and [13] have considered three-body orbits including M33, the third most massive LG galaxy; and [13] also considered the Large Magellanic Cloud (LMC), the fourth most massive LG galaxy. They still conclude that a merger is certain. We will consider two-body, three-body, and four-body systems to study the future evolution of the LG and reveal the distinct effects of the M33 and LMC on the MW-M31 orbit.', 'A fiducial model based on the most accurate values available': 'In the context of the Local Group, it is important to note that apart from the sky positions, all parameters including distance moduli, line-of-sight velocities, proper motions, masses, and concentrations carry non-negligible uncertainties. We will use Monte Carlo (NC) sampling of all values to investigate how these observational errors propagate to uncertainties on the future evolution, and in particular, the probability of a merger between the MW and M31.  \nOur fiducial LG model contains the four most massive LG members: the MW, M31, M33, and the LMC, based on the latest and most accurate available data. The mass of the MW has recently been extensively studied using Gaia data, with a consensus emerging of a total mass close to 10 12 M /circledot . We adopt a value of M 200 = 1 ± 0 . 2 × 10 12 M /circledot (excluding the mass of the LMC, which we treat separately). For all other galaxies, there is significantly more uncertainty. For M31, we adopt M 200 = 1 . 3 ± 0 . 4 × 10 12 M /circledot . For M33, we assume M 200 = 3 ± 1 × 10 11 M /circledot and for the LMC we assume M 200 = 1 . 5 ± 0 . 5 × 10 11 M /circledot . See Methods for a review of mass estimates.  \nThe line-of-sight velocities are well-known and we adopt the values given by [1]. For the distance moduli, µ , we choose the most accurate and precise recent values in the literature: [25] for the LMC and [26] for M33. For M31, we use recent HST Cepheid results [27].  \nFor the proper motions, where we use the notation, µ δ for the proper motion in declination and µ ∗ α = µ α cos δ for the proper motion in right ascension, we use the Hubble Space Telescope (HST) proper motions of [28] for the LMC, and the combined HST and Gaia DR2 proper motions of [13] for M33. Our fiducial M31 proper motions are the most precise values in the literature based on Gaia DR3 astrometry [29]. However, we find very similar results using the combined Gaia DR2 and HST proper motions of [13], as shown in Extended Data Figure 1. All assumed values are listed in Extended Table 1. We assume Galactocentric coordinates of RA GC = 266 . 41 · , Dec GC = -28 . 94 · [30], d GC = 8 . 122 kpc [31], and a velocity of (12 . 9 , 245 . 6 , 7 . 78) kms -1 with respect to the Sun [30, 31, 32].  \nTo account for the fact that the true probability distributions may not be Gaussian and to exclude possible effects caused by unrealistic or even unphysical outliers, we truncate all distributions at ± 2 σ . This truncation increases the probability of a MWM31 merger by only ∼ 10% (see Extended Data Figure 2). For the fiducial model, we use 50000 MC samples while for all other variants, we use 2500 samples, ensuring that all statistical errors on the merger rate are below 1%.', 'The evolution of the MW-M31-M33-LMC system': "The Monte Carlo initial conditions are integrated numerically using a symplectic direct scheme (see Methods for details). Figure 1 shows 100 realisations each of the MW-M31 orbit, in the two-body MW-M31, and the four-body MW-M31-M33-LMC systems. Figure 2 shows the evolution of the distances between the MW and M31 for the same sets of orbits, and additionally, for the MW-M31-M33 and MW-M31-LMC three-body systems.  \nThe MW-M31 two-body orbit evolves in a plane and leads to a merger in slightly less than half of cases, the majority of which occur during the second pericentre. The addition of M33 increases the merger probability to ∼ 2 / 3, with a similar median merger time. However, the addition of the LMC has the opposite effect: the pure MWM31-LMC system experiences a merger in only slightly more than 1 / 3 of cases and the merger probability of the full M31-MW-M33-LMC system is just over 50%.  \nFor each system, we also show the single orbit obtained by adopting the most likely value of each observable, either in the fiducial model or assuming HST+ Gaia DR2 proper motions for M31 [13]. In the case of the MW-M31-M33 system, which is the one considered by [12] and [13], we reproduce similar orbits and times of first pericentre and merger (the differences are likely due to differences in the other model parameters). However, despite using newer and more precise measurements, when we perform a Monte Carlo analysis, we find considerable uncertainty in the outcome that was not previously reported. In particular, in the full MW-M31-M33-LMC system, a merger between the MW and M31 occurs within the next 10 Gyr in only approximately half of the cases. This is in stark contrast to all previous results that had only considered the most likely values without accounting for the numerous and significant uncertainties.  \nIn Figure 3, we show the probability distributions of the merger time and of the minimum distance between the MW and M31, as well as the 'survival' rate of the MW over time, i.e. the probability that no merger with M31 has occurred. In the fiducial model, we consider a merger to occur when the distance between any two galaxies is below 20 kpc, but we find that our results are not sensitive to this particular choice (see Methods for details). Due to the effect of dynamical friction, we find that there are two distinct possibilities for the eventual fate of the MW and M31: orbits that come within less than ∼ 200 kpc eventually merge, which would likely lead to the formation of an intermediate-mass elliptical galaxy [10, 33, 34]. For systems that merge, we find a median time of 7 . 6 Gyr in the fiducial model, or 8 . 0 Gyr adopting a 10 kpc threshold. By contrast, orbits with larger pericentres do not decay due to dynamical friction. In this case, the MW and M31 continue to evolve in isolation. Based on the best current data, both outcomes are almost equally likely.", 'The roles of M33 and the LMC': 'The distinct effects of each satellite on the MW-M31 orbit are illustrated in Figures 4 and 5, where we compare the trajectories of the MW and M31 in two-body systems and in systems that also include either M33 or the LMC. Both satellites provide some additional acceleration in the radial direction of the MW-M31 orbit. However, importantly, both satellites also affect the motion of their respective hosts. Including M33 in the calculation decreases the transverse velocity of M31 with respect to the MW. By contrast, as already pointed out in [35], at its current orbital phase, the recoil due to the LMC results in a lower transverse velocity measured between the MW and M31. During its short orbital period of ∼ 1 . 5 Gyr, the LMC will accelerate the MW to a higher transverse velocity. In addition, the inclusion of M33 largely provides momentum in the original MW-M31 plane, while the inclusion of the LMC also provides significant momentum perpendicular to the MW-M31 plane. In our  \nanalysis, the LMC is certain to merge with the MW, and M33 is highly likely ( ∼ 86%) to merge with M31 before any possible MW-M31 merger. The net effect of adding M33 to the two-body system is to increase the merger probability, while the net effect of adding the LMC is to decrease it.', 'Sources of uncertainty': 'In Figure 6, we show how the merger probability depends on the different observables. In each panel, we show the dependence on two variables in the ± 2 σ ranges, with the remaining variables Monte Carlo sampled. The effects of the concentration parameters are shown in Extended Data Figure 3. We find that all masses, the proper motions of M31 and M33, and the distance moduli of M31 and M33, significantly impact the probability of a merger. The merger probability is positively correlated with the masses of the MW, M31, and M33 and negatively correlated with the mass of the LMC. The impact of the satellite masses is more pronounced for lower combined masses of the MWand M31. The ± 2 σ ranges of the M31 proper motions include values that imply a merger probability above 90%, but also values that imply a merger probability close to zero. The most likely proper motions (assuming no errors) only lead to a merger in ∼ 2 / 3 of cases.  \nFuture more precise proper motion measurements may significantly change the expected outcome, although they could make the merger either more or less likely. However, if they fall within ± 1 σ of the current most likely values, even precise M31 proper motion measurements alone will not suffice to determine the outcome.  \nEven the comparatively high precision of the line-of-sight velocity and distance moduli for M33 and M31 contribute significant uncertainty, with the probability of a merger varying between 40% and over 60% for different combinations in the ± 2 σ ranges around the most likely values.  \nGiven the considerable measurement errors, it is worth noting that cosmological simulations result in a present-day MW-M31 transverse velocity prior of v t = 75 +65 -40 kms -1 [36]. There is thus no reason to assume that the transverse velocity measured using Gaia DR3 (most likely value of v t = 76 kms -1 ) is overestimated or to expect that more precise measurements will result in a lower value. It is also worth noting that a perfectly radial present-day M31-MW motion ( µ ∗ α = 38 . 03 mas yr -1 , µ δ = -21 . 37 mas yr -1 ) is neither compatible with current observations nor does it result in the highest merger probability in the four-body system.', 'Summary': 'Even using the latest and most precise available observational data, the future evolution of the Local Group is uncertain. Intriguingly, we find an almost equal probability for the widely publicised merger scenario (albeit with a later median time to merger) and one where the MW and M31 survive unscathed. We reach this conclusion by including the LMC and importantly, for the first time, considering the relevant uncertainties in the observables.  \nOur results are not sensitive to the necessary choices of gravitational softening (see Extended Data Figure 4), merger threshold (Extended Data Figure 5), or dynamical friction scheme (Extended Data Figure 6)).  \nWhile we have shown that considerable uncertainty results from the proper motion measurements of M31, we also find that a more accurate prediction requires more precise measurements of the positions, motions, and masses of all participating galaxies. The dependence of the evolution of the MW-M31 system on the treatment of other LG galaxies points to further uncertainties. Cosmological simulations suggest that ∼ 25% of the bound mass of the LG is outside the two main haloes [37]. The next most massive individual Local Group galaxies that could impact the MW-M31 orbit are M32 (an M31 satellite and possible merger remnant [38]) and the Small Magellanic Cloud (SMC), a satellite of the MW. Both are at least a factor of five less massive than the LMC, and we find that including the SMC, for which proper motion measurements are available, does not significantly change the MW-M31 merger probability (see Extended Data Figure 7). However, the unaccounted cumulative effects of additional substructures as well as that of the cosmic environment introduce further uncertainty, particularly towards the far future. An accurate prediction of the evolution even from perfectly precise observations may require cosmological constrained simulations [39, 40, 41] that can account for these effects. Meanwhile, the assumptions and simplifications we have made here are likely conservative regarding our central claim, that there is considerable uncertainty about the MW-M31 merger.  \nUpcoming Gaia data releases will improve the proper motion constraints and mass models are continuously refined. However, it is clear that Galactic eschatology is still in its infancy and significant work is required before the eventual fate of the Local Group can be predicted with any certainty. As it stands, proclamations of the impending demise of our Galaxy appear greatly exaggerated.', 'References': "- [1] McConnachie, A. W. The Observed Properties of Dwarf Galaxies in and around the Local Group. AJ 144 , 4 (2012). 1204.1562.\n- [2] Newton, O. et al. The Undiscovered Ultradiffuse Galaxies of the Local Group. ApJ 946 , L37 (2023). 2212.05066.\n- [3] Sawala, T. et al. The APOSTLE simulations: solutions to the Local Group's cosmic puzzles. MNRAS 457 , 1931-1943 (2016). 1511.01098.\n- [4] Slipher, V. M. The radial velocity of the Andromeda Nebula. Lowell Observatory Bulletin 2 , 56-57 (1913).\n- [5] Hubble, E. P. A spiral nebula as a stellar system, Messier 31. ApJ 69 , 103-158 (1929).\n- [6] Peebles, P. J. E., Phelps, S. D., Shaya, E. J. & Tully, R. B. Radial and Transverse Velocities of Nearby Galaxies. ApJ 554 , 104-113 (2001). astro-ph/0010480.\n- [7] Loeb, A., Reid, M. J., Brunthaler, A. & Falcke, H. Constraints on the Proper Motion of the Andromeda Galaxy Based on the Survival of Its Satellite M33. ApJ 633 , 894-898 (2005). astro-ph/0506609.\n- [8] van der Marel, R. P. & Guhathakurta, P. M31 Transverse Velocity and Local Group Mass from Satellite Kinematics. ApJ 678 , 187-199 (2008). 0709.3747.\n- [9] Sohn, S. T., Anderson, J. & van der Marel, R. P. The M31 Velocity Vector. I. Hubble Space Telescope Proper-motion Measurements. ApJ 753 , 7 (2012). 1205.6863.\n- [10] Dubinski, J., Mihos, J. C. & Hernquist, L. Using Tidal Tails to Probe Dark Matter Halos. ApJ 462 , 576 (1996). astro-ph/9509010.\n- [11] Cox, T. J. & Loeb, A. The collision between the Milky Way and Andromeda. MNRAS 386 , 461-474 (2008). 0705.1170.\n- [12] van der Marel, R. P. et al. The M31 Velocity Vector. II. Radial Orbit toward the Milky Way and Implied Local Group Mass. ApJ 753 , 8 (2012). 1205.6864.\n- [13] van der Marel, R. P. et al. First Gaia Dynamics of the Andromeda System: DR2 Proper Motions, Orbits, and Rotation of M31 and M33. ApJ 872 , 24 (2019). 1805.04079.\n- [14] Cowen, R. Andromeda on collision course with the Milky Way. Nature (2012).\n- [15] Harvey-Smith, L. When galaxies collide (Melbourne University Press, 2020).\n- [16] Binney, J. & Tremaine, S. Galactic Dynamics: Second Edition (Princeton University Press, 2008).\n- [17] Eicher, D. J. The New Cosmos: Answering Astronomy's Big Questions (Cambridge University Press, 2015).\n- [18] Helmi, A. et al. The merger that led to the formation of the Milky Way's inner stellar halo and thick disk. Nature 563 , 85-88 (2018). 1806.06038.\n- [19] Belokurov, V., Erkal, D., Evans, N. W., Koposov, S. E. & Deason, A. J. Coformation of the disc and the stellar halo. MNRAS 478 , 611-619 (2018). 1802. 03414.\n- [20] Ruiz-Lara, T., Gallart, C., Bernard, E. J. & Cassisi, S. The recurrent impact of the Sagittarius dwarf on the star formation history of the Milky Way. Nature Astronomy 4 , 965-973 (2020). 2003.12577.\n- [21] Naidu, R. P. et al. Reconstructing the Last Major Merger of the Milky Way with  \nthe H3 Survey. ApJ 923 , 92 (2021). 2103.03251.  \n- [22] Navarro, J. F., Frenk, C. S. & White, S. D. M. A Universal Density Profile from Hierarchical Clustering. ApJ 490 , 493-508 (1997). astro-ph/9611107.\n- [23] Zentner, A. R. & Bullock, J. S. Halo Substructure and the Power Spectrum. ApJ 598 , 49-72 (2003). astro-ph/0304292.\n- [24] Cautun, M., Deason, A. J., Frenk, C. S. & McAlpine, S. The aftermath of the Great Collision between our Galaxy and the Large Magellanic Cloud. MNRAS 483 , 2185-2196 (2019). 1809.09116.\n- [25] Pietrzy'nski, G. et al. A distance to the Large Magellanic Cloud that is precise to one per cent. Nature 567 , 200-203 (2019). 1903.08096.\n- [26] Ou, J.-Y. et al. A Distance Measurement to M33 Using Optical Photometry of Mira Variables. AJ 165 , 137 (2023). 2302.02901.\n- [27] Li, S. et al. A Sub-2% Distance to M31 from Photometrically Homogeneous Near-infrared Cepheid Period-Luminosity Relations Measured with the Hubble Space Telescope. ApJ 920 , 84 (2021). 2107.08029.\n- [28] Kallivayalil, N., van der Marel, R. P., Besla, G., Anderson, J. & Alcock, C. Thirdepoch Magellanic Cloud Proper Motions. I. Hubble Space Telescope/WFC3 Data and Orbit Implications. ApJ 764 , 161 (2013). 1301.0832.\n- [29] Salomon, J. B. et al. The proper motion of Andromeda from Gaia EDR3: confirming a nearly radial orbit. MNRAS 507 , 2592-2601 (2021). 2012.09204.\n- [30] Reid, M. J. & Brunthaler, A. The Proper Motion of Sagittarius A*. II. The Mass of Sagittarius A*. ApJ 616 , 872-884 (2004). astro-ph/0408107.\n- [31] GRAVITY Collaboration et al. Detection of the gravitational redshift in the orbit of the star S2 near the Galactic centre massive black hole. A&A 615 , L15 (2018). 1807.09409.\n- [32] Drimmel, R. & Poggio, E. On the Solar Velocity. Research Notes of the American Astronomical Society 2 , 210 (2018).\n- [33] Forbes, D. A., Masters, K. L., Minniti, D. & Barmby, P. The elliptical galaxy formerly known as the Local Group: merging the globular cluster systems. A&A 358 , 471-480 (2000). astro-ph/0001477.\n- [34] Schiavi, R., Capuzzo-Dolcetta, R., Arca Sedda, M. & Spera, M. The collision between the Milky Way and Andromeda and the fate of their Supermassive Black Holes. In Bragaglia, A., Davies, M., Sills, A. & Vesperini, E. (eds.) Star Clusters: From the Milky Way to the Early Universe, vol. 351, 161164 (2020). 1908.07278.\n- [35] Pe˜narrubia, J., G'omez, F. A., Besla, G., Erkal, D. & Ma, Y.-Z. A timing constraint on the (total) mass of the Large Magellanic Cloud. MNRAS 456 , L54-L58 (2016). 1507.03594.\n- [36] Sawala, T., Teeriaho, M. & Johansson, P. H. The Local Group's mass: probably no more than the sum of its parts. MNRAS 521 , 4863-4877 (2023). 2210.07250.\n- [37] Sawala, T., Pe˜narrubia, J., Liao, S. & Johansson, P. H. The timeless timing argument and the total mass of the Local Group. MNRAS 526 , L77-L82 (2023). 2307.13732.\n- [38] D'Souza, R. & Bell, E. F. The Andromeda galaxy's most important merger about 2 billion years ago as M32's likely progenitor. Nature Astronomy 2 , 737-743  \n(2018). 1807.08819.  \n- [39] Libeskind, N. I. et al. The HESTIA project: simulations of the Local Group. MNRAS 498 , 2968-2983 (2020). 2008.04926.\n- [40] Sawala, T. et al. The SIBELIUS Project: E Pluribus Unum. MNRAS 509 , 1432-1446 (2022). 2103.12073.\n- [41] Wempe, E. et al. Constrained cosmological simulations of the Local Group using Bayesian hierarchical field-level inference. arXiv e-prints arXiv:2406.02228 (2024). 2406.02228.\n- [42] Astropy Collaboration et al. The Astropy Project: Sustaining and Growing a Community-oriented Open-source Project and the Latest Major Release (v5.0) of the Core Package. ApJ 935 , 167 (2022). 2206.14220.\n- [43] Harris, C. R. et al. Array programming with NumPy. Nature 585 , 357-362 (2020). URL https://doi.org/10.1038/s41586-020-2649-2.\n- [44] Virtanen, P. et al. SciPy 1.0: Fundamental Algorithms for Scientific Computing in Python. Nature Methods 17 , 261-272 (2020).\n- [45] Diemer, B. COLOSSUS: A Python Toolkit for Cosmology, Large-scale Structure, and Dark Matter Halos. ApJS 239 , 35 (2018). 1712.04512.  \nFig. 1 Possible future MW-M31 orbits. Coloured lines show probability densities for the positions of the MW and M31 in 100 Monte Carlo samples of the fiducial model, integrated over 10 Gyr or until the merger. On the top and bottom panels, respectively, trajectories are projected in the orbital plane and perpendicular to the orbital plane defined by the initial positions and velocities of the MW and M31. White markers denote MW-M31 mergers. In the left column, we show the MW-M31 two-body system, while in the right, we show the four-body system that includes the MW, M31, M33, and the LMC.  \n<!-- image -->  \n<!-- image -->  \n<!-- image -->  \n<!-- image -->  \nFig. 2 Distance between the MW and M31. On each panel, we show 100 realisations of our fiducial MC model and also state the probability for a MW-M31 merger within 10 Gyr. In the top row, we show the MW-M31 two-body orbits (top-left) and MW-M31-M33-LMC four-body orbits (top-right) shown in Figure 1. In the bottom row, we show the MW-M31-M33 (bottom-left) and MW-M31-LMC (bottom-right) three-body orbits. White markers denote MW-M31 mergers; percentages indicate the fraction of orbits that merge within 10 Gyr. Only slightly more than half of the four-body orbits lead to a merger within 10 Gyr. The inclusion of M33 increases the likelihood of a merger, whereas the inclusion of the LMC decreases it. White lines show the individual orbits using the most likely values of every variable, either assuming the Gaia DR3 proper motions [29] of the fiducial model (solid) or HST + Gaia DR2 proper motions [13] (dashed).  \n<!-- image -->  \n50%  \n20%  \n10%  \n5%  \n1%  \n<5  \ndm  \n=20 kpc  \ndm  \n=10 kpc  \n5  \n6  \n7  \n8  \n9  \n10 >10  \ntm  \n(  \nM  \n31  \nMW  \n) [Gyr]  \n<15 15 25  \n50 100 200 400>400  \nmin(  \ndM  \n31  \nMW  \n)[kpc]  \nFig. 3 Distributions of the merger time, t m , the minimum distance, min( d M 31 -MW ), and the 'survival' rate of the MW in the MW-M31-M33-LMC system. Blue and red colours show results using a distance threshold of 20 kpc (our default) or 10 kpc for a merger, respectively. Mergers happen on average later when the distance threshold is lower, but the fraction of systems that merge within 10 Gyr is similar. The median time to merger is 7.6 Gyr for systems that merge with a 20 kpc threshold and 8.0 Gyr for systems that merge with a 10 kpc threshold. The distributions of minimum distance show a clear bimodality, independently of the threshold: about half of the systems reach the merger threshold within 10 Gyr while the vast majority of the remaining systems do not approach closer than 200 kpc. The 'survival' rate of the MW drops sharply between ∼ 6 -9 Gyr and levels off afterwards.  \n<!-- image -->  \n50%  \n20%  \n10%  \n5%  \n1%  \ndm  \n=20 kpc  \ndm  \n=10 kpc  \nFig. 4 Effects of M33 and the LMC on the trajectory of the MW. As in Figures 1 and 5, panels in the top row are projected in the orbital plane defined by the initial MW and M31 positions and velocities, while those in the bottom row are projected perpendicular to this orbital plane. Arrows point towards the initial position of M31. On the left, we compare MW trajectories in the two-body MW-M31 system to those in the three-body MW-M31-M33 system. On the right, we show the orbit of the LMC, and that of the MW before and after the merger with the LMC in the MW-M31-LMC system. Compared to the two-body orbit, the inclusion of M33 reduces the transverse velocity of the MW relative to M31 and introduces only a small velocity perpendicular to the MW-M31 orbital plane. By contrast, the LMC increases the MW-M31 transverse velocity and causes a larger velocity perpendicular to the MW-M31 orbital plane.  \n<!-- image -->  \nFig. 5 Effects of M33 and the LMC on the trajectory of M31. As in Figures 1 and 4, panels in the top row are projected in the orbital plane defined by the initial MW and M31 positions and velocities, while those in the bottom row are projected perpendicular to this orbital plane. Arrows point towards the initial position of the MW. On the right, we compare M31 trajectories in the two-body MW-M31 system to those in the three-body MW-M31-LMC system. On the left, we show the orbit of M33, and that of M31 before and after a possible merger with M33 in the MW-M31-M33 system. Compared to the 2-body orbit, the inclusion of the LMC increases the transverse velocity of M31 and introduces motion perpendicular to the initial orbital plane. On the other hand, the inclusion of M33 reduces the transverse velocity of M33 with respect to the MW.  \n<!-- image -->  \n<!-- image -->  \n<!-- image -->  \n<!-- image -->  \n<!-- image -->  \n<!-- image -->  \n<!-- image -->  \n<!-- image -->  \n<!-- image -->  \nFig. 6 Dependence of the merger probability on observables of the MW-M31-M33-LMC system. In each panel, we show the probability of a MW-M31 merger within 10 Gyr as a function of two variables, with all other observables sampling the probability distributions of the fiducial model. White lines on the colour bars indicate the minimum and maximum merger probabilities for the range of values shown on each panel, indicating the sensitivity of the merger probability on the two corresponding variables. In the top row, we show the dependence on different masses, the second row shows the dependence on proper motions, and the third row shows the dependence on distance moduli and line of sight velocities. The axes extend to ± 2 σ of the fiducial model. The merger probability is positively correlated with the mass of the MW, M31, and M33, and negatively correlated with the mass of the LMC. The merger probability strongly depends on µ δ ( M 31) and µ ∗ α ( M 31), but also varies significantly with µ δ ( M 33) and µ ∗ α ( M 33). The uncertainties in the distance moduli for M31 and M33 also contribute to the uncertainty of the outcome, while the effect of the line-of-sight velocities is small.  \n<!-- image -->", 'Methods': 'Our results are based on numerically integrated initial conditions, which are in turn based on Monte Carlo samples of the observational data. We describe here the generation of the Monte Carlo samples, the numerical integration, and our treatment of dynamical friction and galaxy mergers. We also demonstrate the robustness of our findings to the particular choices made and show that our fiducial model and the results are conservative in predicting an uncertain future and relatively low merger probability. We also present results when the SMC is included in the analysis, in addition to the four main galaxies. To facilitate the reproduction of our results and allow future work to easily incorporate new observational data, our analysis code is flexible and public.', 'Monte Carlo Samples': 'Our fiducial model consists of the Milky Way and three additional galaxies, M31, M33, and the LMC. We assume that the sky coordinates (RA & Dec) for the centres of M31, M33, and the LMC are known. We furthermore assume that the position of the Galactic centre with respect to the sun is fixed at ( Ra,Dec ) = (266.4051, 28.936175) [30], d = 8 . 122 kpc [31], and that the galactocentric velocity of the Sun is (12 . 9 , 245 . 6 , 7 . 78) km -1 [32, 31, 30].  \nWe create Monte Carlo samples for the remaining 20 variables, the four halo masses, M MW , M M 31 , M M 33 , M LMC , the four halo concentration parameters, the three distance moduli, µ , the three sets of proper motions, µ δ and µ ∗ α , and three lineof-sight velocities, v los . Sampling directly in the space of the observables rather than sampling in a space of derived variables such as Cartesian coordinates or velocities minimises the effect of possible correlations.  \nTo allow reproducibility and identification of individual orbits across our figures when changing parameters of the model such as the set of galaxies included and numerical parameters such as the gravitational softening length and merger threshold, we re-initialise the pseudo-random number generator with the same values for every new set of Monte Carlo samples. While we generally only show the first 100 orbits on each plot, all quoted probabilities are computed from at least 2500 samples each, so that statistical errors are less than 1%.  \nThe most likely values of each variable together with the ± 1 σ uncertainty are either taken directly from single sources in the literature or in the case of masses and concentrations (see below), estimated by us based on multiple sources. In creating our Monte Carlo samples, we assume that each variable follows a Gaussian probability distribution. However, in our fiducial model, we truncate all distributions at ± 2 σ , corresponding to the central ∼ 95% of values. While Gaussian distributions are a natural assumption for measurement errors, this is not always explicitly stated, and it may not reflect the true probability distribution, especially at some distance from the most likely values. Indeed, for some variables, untruncated distributions extend to unphysical values with finite probabilities. A truncation at 2 σ also ensures that all variables in the samples are physical.  \nGiven that our central claim is uncertainty in the future evolution of the LG, truncating the probability distributions of the observables is a conservative assumption that leads to lower uncertainty about the outcome. However, as shown in Extended Data Figure 2, our results are not sensitive to the truncation at ± 2 σ , and the merger probability is only slightly lower when the distributions are not truncated. On the other hand, even a truncation at only ± 1 σ leaves approximately 1 / 4 of orbits that do not merge within 10 Gyr.', 'Numerical Integration': 'The orbits are integrated using a symplectic leapfrog algorithm in the centre-of-mass frame, with a step size of 1 Myr, approximately the time it takes for a galaxy moving at 100 kms -1 to travel 0.1 kpc - 1/200 th of our merger threshold (see below). Our results are not affected by the finite time step.  \nTo account for the fact that the haloes are extended objects, the gravitational force between the haloes is softened with a constant softening length of 20 kpc, similar to the scale radius of an NFW halo in the mass range we consider here. We also consider different choices for the softening and show in Extended Data Figure 4 results with softening lengths of 10, 20, or 30 kpc, respectively. A softening length that is too small can lead to unphysical hard scattering events during close encounters while a softening length that is too large artificially reduces the gravitational force. In the context of the Local Group, both of these effects could reduce the merger probability. However, we find no strong dependence of the merger probability on the softening length, with our adopted fiducial value of 20 kpc resulting in the highest merger probability.', 'Dynamical Friction': "To estimate the effect of dynamical friction, we use a modified Chandrasekhar formula, similar to that used in [24]. The classic Chandrasekhar formula assumes a point mass orbiting in the potential of a much more massive host halo, which is in turn composed of much less massive particles. This approach has been expanded to account for extended satellites [46], and we use the following expression to calculate the acceleration of a satellite due to dynamical friction [16] once inside r 200 of the host halo:  \nd v d t = -4 πG 2 Mρ ln Λ v 2 [ erf( X ) -2 X √ π e -X 2 ] v v , (1)  \nwhere G is the gravitational constant, M is the mass of the satellite, v is the velocity of the satellite relative to the host, ρ is the density of the host at the position of the satellite, X = v/ (2 σ ) is the ratio between the orbital speed of the satellite and the 1D-velocity dispersion, σ , of the host at the location of the satellite, and Λ is the Coulomb factor expressed as r//epsilon1 , with /epsilon1 a scale length that depends on the density of the satellite. To determine /epsilon1 , we adopt an empirical expression derived from N-body simulations [47]:  \n/epsilon1 = { 2 . 2 r s -14 kpc if r s ≥ 8kpc 0 . 45 r s if r s < 8kpc (2)  \nwhere r s is the scale radius of the satellite's NFW halo. Finally, to approximate the velocity dispersion of the host at the location of the satellite, we use the expression derived in [23] for NFW haloes.  \nStandard dynamical friction schemes assume a clear hierarchy between the (much more massive) host and the (much less massive) satellite. According to the assumptions underlying Equation 1, the satellite and host enter the calculation in clearly defined and distinct roles, with the dynamical friction force applied only on the satellite, while the host remains unaffected.  \nHowever, in the Local Group context where galaxies and haloes of similar mass are interacting, this introduces an inconsistency. In particular, in the (relatively likely) scenario that M31 has a similar mass to the MW, the roles of the satellite and host are unclear, but their assignment changes the result of the calculation. For example, if the MW is considered the satellite, it would be accelerated by its interaction with M31 while M31 would remain unaffected, and only the motion of the MW with respect to the other galaxies would be affected while that of M31 would remain unchanged. If M31 is considered the satellite, the roles would be reversed. Accelerating only one galaxy also violates momentum conservation.  \nTo make the calculation more symmetrical, in our dynamical friction calculation we distribute the dynamical friction force proportionally between the satellite ( s ) and host ( h ), conserving the total momentum:  \nd v s d t = a DF M h M s + M h , (3)  \nd v h d t = -a DF M s M s + M h , (4)  \nwhere v s and M s are the velocity and mass of the satellite, v h and M h are the mass and velocity of the host, and a DF is the acceleration computed using Equation 1. In the limit that the satellite is much less massive than the host, the standard hierarchical scheme is recovered and only the satellite is accelerated, while in the limit that both galaxies have equal mass, both receive equal and opposite accelerations.  \nA small inconsistency remains in that even when the differences in mass are small, we still assign the more massive galaxy as the host and the less massive galaxy as the satellite when calculating the magnitude of the dynamical friction, where the velocity dispersion of the host, σ , but not that of the satellite, and through the Coulomb factor, the scale radius of the satellite, but not that of the host, are considered. In our spherical halo models, both the velocity dispersion and Coulomb factor depend only on the assumed masses and concentrations, and as we discuss below, the concentration of the satellite has a greater impact on the dynamical friction force. For two halos of significantly different masses, the assignment of host and satellite is clear. For an individual case of two haloes of nearly identical masses but different (randomly  \nassigned) concentrations, the choice of calculating the dynamical friction force by treating either halo as the satellite seems arbitrary. However, for a large number of samples of nearly equal-mass interactions with randomly assigned concentrations, the dynamical friction calculations are not biased. We also assume identical distributions of concentration parameters for all galaxies.  \nIn Extended Data Figure 6, we compare the orbits of the fiducial system using no dynamical friction, 'hierarchical' dynamical friction (only from the more massive host to the less massive satellite), and our default 'proportional' dynamical friction. It is clear that without dynamical friction a MW-M31 merger is highly unlikely. In fact, the finite merger rate without dynamical friction depends strongly on our default impact parameter threshold of 20 kpc, with a lower threshold, the merger rate can become arbitrarily small. On the other hand, when dynamical friction is included, the evolution of each orbit is quite similar in the 'hierarchical' and 'proportional' schemes, and the merger rate is not significantly affected by the exact choice of scheme.  \nIt is worth noting that our semi-analytical approach to dynamical friction is still quite simplistic, and while the average behaviour of N-body simulations has been used to calibrate parameters, numerical simulations also show that individual systems with non-zero internal angular momenta and substructures can have different merger times than predicted by these simple equations [48]. A precise prediction of the MW-M31 orbit will likely require full N-body simulations. On the other hand, we show that even a simple dynamical model that assumes no spin, no triaxiality, and no substructure results in considerable uncertainty in the future evolution of the MW-M31 system.", 'Mergers': 'Below a certain distance, the interactions of the gas and stellar components become significant, and our simple approach is no longer appropriate for predicting the remaining orbital evolution. Our aim here is not to predict the precise time of the merger (in fact, we argue that such a prediction is futile based on the current data), so we simply assume that any system that passes below a threshold distance will eventually merge, and identify this time as a lower limit for the time of the merger.  \nIn our fiducial model, we have adopted a value of 20 kpc as the merger threshold for all galaxy interactions. When such a merger occurs between two galaxies, we combine the masses and momenta of the two galaxies at their common centre of mass and continue the integration. The concentration parameter is set to that of the most massive galaxy.  \nIn Extended Data Figure 5, we compare the sensitivity of our results to adopting merger thresholds of 10, 20, and 30 kpc. With a threshold of 10 kpc, we find a slightly reduced MW-M31 merger probability. However, raising the threshold from 20 to 30 kpc has no significant impact on our results. This confirms our assertion in the main section: due to the effect of dynamical friction, orbits either inspiral and eventually merge, or do not come close enough for dynamical friction to become effective and hence do not merge.  \nThe small fraction of MW-M31 orbits that merge with a threshold of 20 kpc, but do not when the threshold is set to 10 kpc, approach with a small impact parameter and high velocity. This reduces the effect of dynamical friction and also allows them  \nto escape to a large apocentre. While this may be a possible scenario for the Local Group, our methods are not adequate for studying such close interactions between galaxies. To be conservative in our prediction of a low merger probability, by setting the merger threshold at 20 kpc we assume that these orbits also merge, and with an even larger threshold of 30 kpc, we would still predict a similar merger rate.  \nFor comparison, the best-studied merger of two galaxies with similar stellar masses to the MW and M31 are the Antenna galaxies. Their past evolution is reproduced with an orbit with pericentre ∼ 10 kpc [49], predicted to lead to coalescence ∼ 1 . 3 Gyr later [50].', 'Fate of the LMC and M33': 'While our main focus here is on the future evolution of the MW-M31 orbit, we naturally also make predictions for the evolution of the LMC and M33.  \nIn the fiducial model, we find that the LMC is certain to merge with the MW before any eventual MW-M31 merger. With a merger threshold of 20 kpc, we find a median time of the LMC-MW merger of 1.3 Gyr, while with a threshold of 10 kpc, we find a median time of 1.9 Gyr.  \nFor M33, with a merger threshold of 20 kpc, we find an ∼ 86% chance of a merger with M31 and a median time of 3.3 Gyr, while with a merger threshold of 10 kpc, we find an ∼ 83% probability of a merger with M31 and a median time of 3.9 Gyr. In both cases, we also find a small probability of ∼ 1 -2% for a merger of M33 with the MW-M31 remnant after a MW-M31 merger within the next 10 Gyr.  \nIt must be noted that our simulations are not designed to study these mergers in detail and ignore, for example, the impact of the disk of the Milky Way. Nevertheless, they broadly agree with the results of [24], who have previously studied the LMC-MW encounter in the presence of M31.', 'Additional galaxies': 'The next most massive LG member galaxy for which proper motion data [28] is available is the Small Magellanic Cloud (SMC), a satellite galaxy of the MW that has likely been accreted together with the LMC. We repeat our analysis including the SMC as a fifth system, and show results in Extended Data Figure 7. Adding the SMC, whose mass is approximately 10% of that of the LMC, has no significant effect on the merger rate. The SMC properties used are listed along with those of the other galaxies in Extended Data Table 1.', 'DR2 + HST proper motions': 'Due to the fact that the M31 proper motions have the largest impact on the probability of a MW-M31 orbit, and for easier comparison with earlier works, particularly [13], we also repeat our analysis adopting the HST+ Gaia DR2 M31 proper motions of [13]. In Extended Data Figure 1, we show the corresponding evolution of the MWM31 distance in the same two-body MW-M31, three-body MW-M31-M33 and MWM31-LMC, and four-body MW-M31-M33-LMC systems. The results can be directly compared to Figure 2 which shows the same quantities in our fiducial model that used  \nGaia DR3 proper motions. In each case, the merger probability is slightly lower using the HST+ Gaia DR2 proper motions: 40% instead of 44% for the MW-M31 system, 56% instead of 63% for the MW-M31-M33 system, 34% instead of 37% for the MWM31-LMC system, and 48% instead of 54% for the full MW-M31-M33-LMC system. However, both sources of proper motions predict similar distributions of outcomes and a similar uncertainty about the MW-M31 merger. This difference can be attributed to the slightly lower precision of the HST+ Gaia DR2 proper motions.', 'Sources for the masses': 'As discussed in Figure 6, the assumed masses of all four galaxies and their associated uncertainties have a strong impact on the likely evolution of the MW-M31 orbit in our fiducial model. Here, we review the mass measurements and show the implications for some alternative scenarios.  \nMW The total mass of the MW has been extensively studied with different methods and tracers, and the accurate astrometry of the Gaia space telescope has brought a flurry of recent measurements. Estimates for the total MW mass based on Gaia DR2 or DR3 satellite dynamics include 1 . 17 +0 . 21 -0 . 15 × 10 12 M /circledot [51], 1 . 51 0 . 45 -0 . 40 × 10 12 M /circledot [52], 1 . 23 +0 . 21 -0 . 18 × 10 12 M /circledot [53], and 1 . 1 +0 . 1 -0 . 1 × 10 12 M /circledot [54] (where the latter two works also use a simulation based prior).  \nEstimates using rotation curves include 1 . 08 +0 . 20 -0 . 14 × 10 12 M /circledot based on Gaia DR2 [55], 0 . 89 +0 . 1 -0 . 08 × 10 12 M /circledot using stars in the galkin catalogue [56], 0 . 822 ± 0 . 052 × 10 12 M /circledot using classical Cepheids [57], and 1 . 08 +0 . 12 -0 . 11 × 10 12 M /circledot from the H3 survey and Gaia DR3 [58]. Other recent measurements include 1 . 54 +0 . 75 -0 . 44 × 10 12 M /circledot using combined Gaia DR2 and HST kinematics of globular clusters [59], 1 . 16 ± 0 . 24 × 10 12 M /circledot (including the mass of the LMC) using the kinematics of halo stars, 1 . 26 +0 . 40 -0 . 22 × 10 12 M /circledot from high velocity RR-Lyrae stars [60], and 1 . 19 +0 . 49 -0 . 32 × 10 12 M /circledot from a Bayesian estimate using dwarf galaxy kinematics from multiple sources [61].  \n[36] contains an overview of recent measurements and [62] includes a comprehensive review of earlier results.  \nAnalysis from several different tracers, and particularly from the latest studies using the latest Gaia observations consistently point towards a Milky Way total mass that is close to 10 12 M /circledot . Both the simple mean and median of the above measurements are 1 . 16 × 10 12 M /circledot , but it is worth noting that most methods of measuring the total mass of the MW include that of the LMC (at a galactocentric distance of ∼ 50 kpc, well within r vir or r 200 ), which is not always made explicit. When the LMC is excluded, the mass of the MW is reduced by ∼ 0 . 15 × 10 12 M /circledot (see our discussion of the LMC mass below). We adopt M MW = 1 . 0 ± 0 . 2 × 10 12 M /circledot excluding the LMC in our fiducial model, reflecting the consensus of recent observations.  \nM31 While M31 does not enjoy the benefit of accurate Gaia proper motions, there nevertheless exist a number of studies estimating its mass. Most recently, [63] measured a total mass of 1 . 14 +0 . 51 -0 . 35 × 10 12 M /circledot using rotation curves based on LAMOST data release 9 and DESI. [64] measured 1 . 4 ± 0 . 4 × 10 12 M /circledot derived from satellite kinematics, [65] derived 1 . 2 +0 . 9 -0 . 7 × 10 12 M /circledot from kinematics of M31 dwarf spheroidals, and [66] measured a total mass of 1 . 05 +0 . 15 -0 . 15 × 10 12 M /circledot from SED fitting together with  \nthe rotation curve and the kinematics of outer globular clusters and satellite galaxies. [67] measured 1 . 0 × 10 12 M /circledot from the HI rotation curve, [68] measured 2 . 0 +0 . 4 -0 . 3 × 10 12 M /circledot from kinematics of the Giant southern stream. [69] measured 1 . 35 +0 . 15 -0 . 15 × 10 12 M /circledot and [70] found 1 . 4 +0 . 2 -0 . 2 × 10 12 M /circledot both using outer halo globular clusters, while [71] found 0 . 8 ± 0 . 1 × 10 12 M /circledot from high-velocity planetary nebulae, and [72] found 1 . 39 ± 0 . 26 × 10 12 for the total mass combining disk rotation velocities and radial velocities of satellite galaxies and globular clusters. Both the simple mean and median of the above values are 1 . 27 × 10 12 M /circledot . Our adopted mass of M M 31 = 1 . 3 ± 0 . 4 × 10 12 M /circledot in the fiducial model reflects the broad consensus of M31 mass estimates using different methods, but also the considerable remaining uncertainty.  \nM33 Mass estimates of M33 are much more sparse. More than twenty years ago, [73, 74] measured a dark matter mass of 5 × 10 10 M /circledot , extrapolated out to a virial mass of 5 × 10 11 M /circledot from the measured HI rotation curve, but noted that this results in a very low baryon fraction. More recently, [75] obtained a similar result using the Hα rotation curve, but noted that because the measurements only extend to a few percent of the virial radius, there are no strong constraints on the total dark matter halo. On the other hand, abundance matching based on the observed stellar mass results in a significantly lower total mass of 1 . 7 ± 0 . 55 × 10 11 M /circledot [76, 77]. Citing both the direct measurements and abundance matching, [77] adopt a mass range of 0 . 8 -3 . 2 × 10 11 M /circledot in their dynamical models of the M33-M31 interaction, while [78] and [13] assume a total mass of 2 . 5 × 10 11 M /circledot .  \nWe adopt M M 33 = 3 ± 1 × 10 11 M /circledot in the fiducial model, marginally compatible with both the extrapolated masses from rotation curve measurements and the results of abundance matching and in line with previous studies. However, the results from the two methods are certainly in tension.  \nLMC For the mass of the LMC, recent mass estimates based on its effect on Galactic stellar streams include 1 . 38 +0 . 37 -0 . 24 × 10 11 M /circledot [79], 1 . 30 ± 0 . 3 × 10 11 M /circledot [80], 1 . 88 +0 . 4 -0 . 35 × 10 11 M /circledot [81] and 1 . 29 +0 . 28 -0 . 23 × 10 11 M /circledot . [82] obtain good agreement between the perturbations of Milky Way halo stars with an LMC mass of 1 . 5 × 10 11 M /circledot . Using the abundance of likely satellites of the LMC, [83] obtain a lower limit of 1 . 24 × 10 11 M /circledot , and using kinematics of satellites associated with the LMC, [84] find 1 . 65 +0 . 47 -0 . 49 × 10 11 M /circledot . Most recently, [85] used 30 LMC globular clusters to infer a total mass of 1 . 80 +1 . 05 -0 . 54 × 10 11 M /circledot while [35] find an even higher value of 2 . 5 +0 . 9 -0 . 8 × 10 11 M /circledot using a timing argument and to be most consistent with the Hubble Flow around the Local Group. The simple mean of the above values is 1 . 6 × 10 11 M /circledot while the median is 1 . 5 × 10 11 M /circledot . A comprehensive recent review on the effect of the LMC on the Milky Way, including a discussion of mass measurements, is given by [86], who conclude that the LMC mass is likely in the range 1 -2 × 10 11 M /circledot , which matches the choice of M LMC = 1 . 5 ± 0 . 5 × 10 11 M /circledot in our fiducial model.', 'Halo Concentrations': "As explained above, we assume that for the purpose of integrating their orbits, all galaxies are represented by NFW haloes [22] with the total masses M 200 defined above. We assume concentration parameters of 10 ± 2 for all galaxies, consistent with results  \nof cosmological simulations [87, 88, 89, 90, 91] in this mass range. While cosmological simulations show the average concentration parameter to be mass-dependent, the halo-to-halo scatter is significantly larger than the change in the mean concentration over this narrow mass range [87].  \nThe concentration parameter does not affect the orbital calculation, except for the dynamical friction, where (for a given mass), it sets the scale radius, r s = r 200 /c , of the 'satellite' that enters in Equation 2.  \nIn Extended Data Figure 3, we show the dependence of the merger probability on the concentration. Except for a very low M31 mass, there is only a weak dependence of the merger probability on the concentration of M31, which is more likely to be the more massive 'host' galaxy in the MW-M31 encounter. For the concentration of the MW, which, in the MW-M31 interaction is more likely to be the satellite, we find that the merger probability is reduced if the concentration is below ∼ 8, i.e. below -1 σ of our fiducial value. We find no significant dependence on the merger probability on the concentration assumed for M33 or the LMC.  \nAcknowledgments. We thank Marius Cautun for his generous help, and Joonas Nattila for helpful suggestions. TS and JD are supported by the Research Council of Finland grant 354905, and TS and PHJ are also supported by the Research Council of Finland grant 339127. JD is supported by an Erasmus+ grant. TS and CSF are supported by the European Research Council (ERC) Advanced Investigator grant DMIDAS (GA 786910) and the STFC Consolidated Grant ST/T000244/1. PHJ, AK, AR and RW also acknowledge support from the European Research Council (ERC) Consolidator Grant KETJU (no. 818930). This work used facilities hosted by the CSC-IT Centre for Science, Finland, and the DiRAC@Durham facility, managed by the Institute for Computational Cosmology on behalf of the STFC DiRAC HPC Facility (www.dirac.ac.uk) and funded by BEIS capital funding via STFC capital grants ST/K00042X/1, ST/P002293/1, ST/R002371/1 and ST/S002502/1, Durham University and STFC operations grant ST/R000832/1. DiRAC is part of the UK National e-Infrastructure. We thank the authors of the open source software listed below.  \nAuthors' contributions. TS planned the project and performed the analysis together with JD. TS, AJD, CSF and PHJ planned the paper, and TS wrote the analysis code, created the figures and wrote the first draft. AK and AR contributed to the dynamical friction calculations. TS, JD, AJD, CSF, PHJ, AK, AR and RW jointly discussed and edited the manuscript.", 'Data availability': 'All data used in this work is publicly available and provided as part of the analysis code listed below.', 'Code availability': "The analysis in this paper was performed using Python 3.10, and makes extensive use of the following open-source libraries: Astropy 6.0.1 [42], Matplotlib 3.8.3, NumPy 1.26.4 [43], SciPy 1.13.0 [44] and Colossus 1.3.5 [45]. A documented Juypter notebook containing the code to produce all figures in this paper is available at: https://github. com/TillSawala/MW-M31.  \nConflict of interest. The authors declare that they have no conflict of interest.  \nExtended Data Table 1 Parameters of the Fiducial model  \n<!-- image -->  \n<!-- image -->  \n<!-- image -->  \nExtended Data Figure 1 Distance between the MW and M31, using HST+ Gaia DR2 proper motions for M31 [13], analogous to Figure 2. As in Figure 2, solid and dashed white lines denote the orbits using the most likely values using either Gaia DR3 proper motions [29] or HST+ Gaia DR2 proper motions. White markers denote MW-M31 mergers, percentages indicate the fraction of orbits that merge within 10 Gyr. The merger rate is slightly lower in all cases when compared to the fiducual model that uses the more precise Gaia DR3 proper motions, but the results are qualitatively similar.  \n<!-- image -->  \nExtended Data Figure 2 Effect of truncating the assumed probability distributions of observables. The distance between the MW and M31 (analogous to Figure 2) in the four-body MW-M31-M33LMC system for different truncations of the observables. From left to right, we show results where the probability distribution for each observable in the fiducial model is truncated to ± 1 σ , ± 2 σ (our default model, same as in Figure 2), or left untruncated. The fraction of systems that merge is only slightly increased when the distributions are clipped at ± 2 σ or above.  \n<!-- image -->  \n<!-- image -->  \n<!-- image -->  \nc  \n<!-- image -->  \nc  \nc  \nExtended Data Figure 3 Effect of the concentration parameters on the merger probability in the four-body MW-M31-M33-LMC system, similar to Figure 6. From left to right, in the top row, we show the dependencies on mass and concentration of the MW and M31, while in the bottom row, we show those of M33 and the LMC, respectively. The concentration parameter affects the merger rate only for the lower mass system in the MW-M31 encounter, which is the MW in most cases. A concentration parameter of the MW below ∼ 7 ( -1 . 5 σ ), particularly in combination with a low MW mass results in a significantly lower merger rate. Unlike their masses, the concentration parameters for M33 and the LMC have no discernible effects on the MW-M31 merger probability.  \n<!-- image -->  \n600  \n400  \n200  \n0  \n200  \n400  \n600  \n<!-- image -->  \nx  \nxMW  \n[kpc]  \n600  \n400  \n200  \n0  \n200  \n400  \n600  \n<!-- image -->  \n<!-- image -->  \n1000  \n800  \n600  \n400  \n200  \n0  \n<!-- image -->  \n<!-- image -->  \n<!-- image -->  \nx  \nxMW  \n[kpc]  \n600  \n400  \n200  \n0  \n200  \n400  \n600  \n<!-- image -->  \nx  \nxMW  \n[kpc]  \n600  \n400  \n200  \n0  \n200  \n400  \n600  \n<!-- image -->  \n1000  \n800  \n600  \n400  \n200  \n0  \nExtended Data Figure 4 Effect of the gravitational softening. MW-M31 orbits (analogous to Figure 1) and distance between the MW and M31 (analogous to Figure 2) in the four-body MW-M31M33-LMC system for different softening lengths. From left to right, we show results with a softening length of 10 kpc, 20 kpc (our default value), and 30 kpc. A softening length that is too small can lead to some unrealistically strong kicks in close encounters, while a softening length that is too large weakens the overall gravitational attraction. However, the merger fraction is not significantly affected by the choice of softening length within this range.  \n<!-- image -->  \nExtended Data Figure 5 Effect of the merger threshold on the MW-M31 distance evolution and merger rate. The distance between the MW and M31 (analogous to Figure 2) in the four-body MWM31-M33-LMC system for different merger thresholds. From left to right, we show results with a threshold of 10 kpc, 20 kpc or 30 kpc for all mergers. The MW-M31 merger probability is not very sensitive to the assumed merger threshold, and we obtain almost the same merger probability even with a (very generous) threshold of 30 kpc.  \n<!-- image -->  \n600  \n400  \n200  \n0  \n200  \n400  \n600  \n<!-- image -->  \n<!-- image -->  \n600  \n400  \n200  \n0  \n200  \n400  \n600  \n<!-- image -->  \n<!-- image -->  \n1000  \n800  \n600  \n400  \n200  \n0  \n<!-- image -->  \n<!-- image -->  \n<!-- image -->  \n<!-- image -->  \n600  \n400  \n200  \n0  \n200  \n400  \n600  \n600  \n400  \n200  \n0  \n200  \n400  \n600  \n1000  \n800  \n600  \n400  \n200  \n0  \nExtended Data Figure 6 Effect of different schemes of dynamical friction. MW-M31 orbits (topt two rows, analogous to Figure 1) and distance between the MW and M31 (bottom row, analogous to Figure 2) in the four-body MW-M31-M33-LMC system for different softening lengths. The left column assumes no dynamical friction, the middle column uses our default 'proportional' scheme where the dynamical friction force is divided such that equal and opposite dynamical forces are applied to both host and satellite, the right column uses the 'hierarchical' scheme where dynamical friction is only applied to the less massive object. Dynamical friction is essential for orbits to decay and for the MWM31 merger to occur, but the probability of an MW-M31 merger is not sensitive to the scheme used. However, due to its momentum-conserving property, mergers in the proportional scheme are more likely to occur close to the original orbital plane compared to the hierarchical scheme.  \n<!-- image -->  \n<!-- image -->  \n<!-- image -->  \n<!-- image -->  \nExtended Data Figure 7 The effect of including the SMC on the MW-M31 distance and merger rate. In the top row, we show the MW-M31 two-body system, and the MW-M31-M33-LMC fourbody system in the fiducial model, corresponding to the top row of Figure 2. In the bottom row, we add the SMC to both systems, i.e. we show the MW-M31-SMC three-body system and the MWM31-M33-LMC-SMC five-body system. The inclusion of the SMC has only a small effect on most MW-M31 orbits and does not significantly affect the total merger rate.  \n<!-- image -->", 'Methods References': "- [46] White, S. D. M. A note on the minimum impact parameter for dynamical friction involving spherical clusters. MNRAS 174 , 467-470 (1976).\n- [47] Jethwa, P., Erkal, D. & Belokurov, V. A Magellanic origin of the DES dwarfs. MNRAS 461 , 2212-2233 (2016). 1603.04420.\n- [48] Boylan-Kolchin, M., Ma, C.-P. & Quataert, E. Dynamical friction and galaxy merging time-scales. MNRAS 383 , 93-101 (2008). 0707.2960.\n- [49] Karl, S. J. et al. One Moment in Time-Modeling Star Formation in the Antennae. ApJ 715 , L88-L93 (2010). 1003.0685.\n- [50] Lah'en, N., Johansson, P. H., Rantala, A., Naab, T. & Frigo, M. The fate of the Antennae galaxies. MNRAS 475 , 3934-3958 (2018). 1709.00010.\n- [51] Callingham, T. M. et al. The mass of the Milky Way from satellite dynamics. MNRAS 484 , 5453-5467 (2019). 1808.10456.\n- [52] Fritz, T. K., Di Cintio, A., Battaglia, G., Brook, C. & Taibi, S. The mass of our Galaxy from satellite proper motions in the Gaia era. MNRAS 494 , 5178-5193 (2020). 2001.02651.\n- [53] Li, Z.-Z. et al. Constraining the Milky Way Mass Profile with Phase-space Distribution of Satellite Galaxies. ApJ 894 , 10 (2020). 1912.02086.\n- [54] Rodriguez Wimberly, M. K. et al. Sizing from the smallest scales: the mass of the Milky Way. MNRAS 513 , 4968-4982 (2022). 2109.00633.\n- [55] Cautun, M. et al. The milky way total mass profile as inferred from Gaia DR2. MNRAS 494 , 4291-4313 (2020). 1911.04557.\n- [56] Karukes, E. V., Benito, M., Iocco, F., Trotta, R. & Geringer-Sameth, A. A robust estimate of the Milky Way mass from rotation curve data. 'Journal of Cosmology and Astroparticle Physics' 2020 , 033 (2020). 1912. 04296.\n- [57] Ablimit, I., Zhao, G., Flynn, C. & Bird, S. A. The Rotation Curve, Mass Distribution, and Dark Matter Content of the Milky Way from Classical Cepheids. ApJ 895 , L12 (2020). 2004.13768.\n- [58] Shen, J. et al. The Mass of the Milky Way from the H3 Survey. ApJ 925 , 1 (2022). 2111.09327.\n- [59] Watkins, L. L., van der Marel, R. P., Sohn, S. T. & Evans, N. W. Evidence for an Intermediate-mass Milky Way from Gaia DR2 Halo Globular Cluster Motions. ApJ 873 , 118 (2019). 1804.11348.\n- [60] Prudil, Z. et al. Milky Way archaeology using RR Lyrae and type II Cepheids II. High velocity RR Lyrae stars, and mass of the Milky Way. arXiv e-prints arXiv:2206.00417 (2022). 2206.00417.\n- [61] Slizewski, A. et al. Galactic Mass Estimates Using Dwarf Galaxies as Kinematic Tracers. ApJ 924 , 131 (2022). 2108.12474.\n- [62] Wang, W., Han, J., Cautun, M., Li, Z. & Ishigaki, M. N. The mass of our Milky Way. Science China Physics, Mechanics, and Astronomy 63 , 109801 (2020). 1912.02599.\n- [63] Zhang, X., Chen, B., Chen, P., Sun, J. & Tian, Z. The rotation curve and mass distribution of M31. MNRAS 528 , 2653-2666 (2024). 2401.01517.\n- [64] Watkins, L. L., Evans, N. W. & An, J. H. The masses of the Milky Way and  \nAndromeda galaxies. MNRAS 406 , 264-278 (2010). 1002.4565.  \n- [65] Tollerud, E. J. et al. The SPLASH Survey: Spectroscopy of 15 M31 Dwarf Spheroidal Satellite Galaxies. ApJ 752 , 45 (2012). 1112.1067.\n- [66] Tamm, A., Tempel, E., Tenjes, P., Tihhonova, O. & Tuvikene, T. Stellar mass map and dark matter distribution in M 31. A&A 546 , A4 (2012). 1208.5712.\n- [67] Chemin, L., Carignan, C. & Foster, T. H I Kinematics and Dynamics of Messier 31. ApJ 705 , 1395-1415 (2009). 0909.3846.\n- [68] Fardal, M. A. et al. Inferring the Andromeda Galaxy's mass from its giant southern stream with Bayesian simulation sampling. MNRAS 434 , 2779-2802 (2013). 1307.3219.\n- [69] Veljanoski, J. et al. Kinematics of Outer Halo Globular Clusters in M31. ApJ 768 , L33 (2013). 1303.7368.\n- [70] Veljanoski, J. et al. The outer halo globular cluster system of M31 - II. Kinematics. MNRAS 442 , 2929-2950 (2014). 1406.0186.\n- [71] Kafle, P. R., Sharma, S., Lewis, G. F., Robotham, A. S. G. & Driver, S. P. The need for speed: escape velocity and dynamical mass measurements of the Andromeda galaxy. MNRAS 475 , 4043-4054 (2018). 1801.03949.\n- [72] Sofue, Y. Dark halos of M 31 and the Milky Way. PASJ 67 , 75 (2015). 1504.05368.\n- [73] Corbelli, E. & Salucci, P. The extended rotation curve and the dark matter halo of M33. MNRAS 311 , 441-447 (2000). astro-ph/9909252.\n- [74] Corbelli, E. Dark matter and visible baryons in M33. MNRAS 342 , 199-207 (2003). astro-ph/0302318.\n- [75] Kam, S. Z. et al. H I Kinematics and Mass Distribution of Messier 33. AJ 154 , 41 (2017).\n- [76] Guo, Q., White, S., Li, C. & Boylan-Kolchin, M. How do galaxies populate dark matter haloes? MNRAS 404 , 1111-1120 (2010). 0909.4305.\n- [77] Patel, E., Besla, G. & Sohn, S. T. Orbits of massive satellite galaxies - I. A close look at the Large Magellanic Cloud and a new orbital history for M33. MNRAS 464 , 3825-3849 (2017). 1609.04823.\n- [78] Patel, E., Carlin, J. L., Tollerud, E. J., Collins, M. L. M. & Dooley, G. A. ΛCDM predictions for the satellite population of M33. MNRAS 480 , 1883-1897 (2018). 1807.05318.\n- [79] Erkal, D. et al. The total mass of the Large Magellanic Cloud from its perturbation on the Orphan stream. MNRAS 487 , 2685-2700 (2019). 1812.08192.\n- [80] Vasiliev, E., Belokurov, V. & Erkal, D. Tango for three: Sagittarius, LMC, and the Milky Way. MNRAS 501 , 2279-2304 (2021). 2009.10726.\n- [81] Shipp, N. et al. Measuring the Mass of the Large Magellanic Cloud with Stellar Streams Observed by S 5 . ApJ 923 , 149 (2021). 2107.13004.\n- [82] Erkal, D. et al. Detection of the LMC-induced sloshing of the Galactic halo. MNRAS 506 , 2677-2684 (2021). 2010.13789.\n- [83] Erkal, D. & Belokurov, V. A. Limit on the LMC mass from a census of its satellites. MNRAS 495 , 2554-2563 (2020). 1907.09484.\n- [84] Correa Magnus, L. & Vasiliev, E. Measuring the Milky Way mass distribution in the presence of the LMC. MNRAS 511 , 2610-2630 (2022). 2110.00018.\n- [85] Watkins, L. L., van der Marel, R. P. & Bennet, P. The Mass of the Large Magellanic Cloud from the Three-dimensional Kinematics of Its Globular Clusters.  \nApJ 963 , 84 (2024). 2401.14458.  \n- [86] Vasiliev, E. The Effect of the LMC on the Milky Way System. Galaxies 11 , 59 (2023). 2304.09136.\n- [87] Neto, A. F. et al. The statistics of Λ CDM halo concentrations. MNRAS 381 , 1450-1462 (2007). 0706.2919.\n- [88] Ludlow, A. D. et al. The mass-concentration-redshift relation of cold dark matter haloes. MNRAS 441 , 378-388 (2014). 1312.0945.\n- [89] Wang, W. et al. Estimating the dark matter halo mass of our Milky Way using dynamical tracers. MNRAS 453 , 377-400 (2015). 1502.03477.\n- [90] Correa, C. A., Wyithe, J. S. B., Schaye, J. & Duffy, A. R. The accretion history of dark matter haloes - III. A physical model for the concentration-mass relation. MNRAS 452 , 1217-1232 (2015). 1502.00391.\n- [91] Diemer, B. & Kravtsov, A. V. A Universal Model for Halo Concentrations. ApJ 799 , 108 (2015). 1407.4730.\n- [92] Watkins, L. L., Evans, N. W. & van de Ven, G. A census of orbital properties of the M31 satellites. MNRAS 430 , 971-985 (2013). 1211.2638.\n- [93] Graczyk, D. et al. The Araucaria Project. The Distance to the Small Magellanic Cloud from Late-type Eclipsing Binaries. ApJ 780 , 59 (2014). 1311.2340.\n- [94] Harris, J. & Zaritsky, D. Spectroscopic Survey of Red Giants in the Small Magellanic Cloud. I. Kinematics. AJ 131 , 2514-2524 (2006). astro-ph/0601025."}