# No certainty of a Milky Way-Andromeda collision

Till Sawala, Jehanne Delhomelle, Alis J Deason, Carlos S Frenk, Jenni Häkkinen, Peter H Johansson, Atte Keitaanranta, Alexander Rawlings, Ruby Wright

► To cite this version:

Till Sawala, Jehanne Delhomelle, Alis J Deason, Carlos S Frenk, Jenni Häkkinen, et al. No certainty of a Milky Way-Andromeda collision. Nature Astronomy, 2025, 22 p. 10.1038/s41550-025-02563-1. hal-05154780

HAL Id: hal-05154780

https://ut3-toulouseinp.hal.science/hal-05154780v1

Submitted on 9 Jul 2025

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# No certainty of a Milky Way–Andromeda collision

Received: 21 June 2024

Accepted: 22 April 2025

Published online: 02 June 2025

<figure><img src="image_2.png" /></figure>

Check for updates

Till Sawala Ⓔ<fnref n="1" />,²✉, Jehanne Delhomelle<fnref n="1" />,³, Alis J. Deason Ⓔ², Carlos S. Frenk Ⓔ², Jenni Häkkinen Ⓔ<fnref n="1" />, Peter H. Johansson Ⓔ<fnref n="1" />, Atte Keitaanranta Ⓔ<fnref n="1" />, Alexander Rawlings<fnref n="1" /> & Ruby Wright Ⓔ<fnref n="1" />,⁴

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 5 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 main sources of uncertainty in the evolution of the Local Group over the next 10 billion years. We found 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. Although 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 probable. In the full system, we found 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 will be no Milky Way–Andromeda merger during the next 10 billion years. Based on the best available data, the fate of our Galaxy is still completely open.

The Local Group (LG) contains two large spiral galaxies, our own Milky Way (MW) and the Andromeda galaxy (M31), along with approximately 100 known smaller galaxies<fnref n="1" />. In addition, it probably hosts other galaxies yet to be discovered², and, according to the standard cosmological model, a vast number of completely dark substructures³. The negative radial velocity of M31 towards the MW has been known for over a century⁴, even before its distance was first accurately measured⁵. However, although indirect methods have since been used to constrain the transverse components of M31's velocity vector⁶⁻⁸, direct measurements of the minute proper motions have been achieved only much more recently with the Hubble Space Telescope (HST)⁹. The first numerical studies¹⁰ 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 probable

future collision and merger¹¹⁻¹³. This scenario has since become the prevalent narrative¹⁴,¹⁵ and textbook knowledge¹⁶,¹⁷.

## Predicting the future of the LG

The MW and M31 both contain remnants of past mergers and interactions with other galaxies¹⁸⁻²¹. 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, as it transfers orbital kinetic energy to internal energy of the objects involved and, consequently, leads to the decay of galactic orbits.

In this study, we parameterized the density profile of each halo as a spherical Navarro–Frenk–White (NFW) profile²² with an isotropic

<figure><img src="image_3.png" /></figure>

<footnote n="1">¹Department of Physics, University of Helsinki, Helsinki, Finland. ²Institute for Computational Cosmology, Durham University, Durham, UK. ³Université de Toulouse, Toulouse, France. ⁴International Centre for Radio Astronomy Research, University of Western Australia, Crawley, Western Australia, Australia.</footnote>

<footnote>e-mail: till.sawala@helsinki.fi</footnote>
<table><caption>Table 1 | Parameters of the fiducial model</caption><thead><tr><th>Galaxy</th><th><i>M</i><sub>200</sub> (10<sup>10</sup> <i>M</i><sub>☉</sub>)</th><th>c</th><th>μ</th><th>μ<sub>a</sub><sup>*</sup> (μas yr<sup>-1</sup>)</th><th>μ<sub>b</sub> (μas yr<sup>-1</sup>)</th><th><i>v</i><sub>los</sub> (kms<sup>-1</sup>)</th></tr></thead><tbody><tr><td>MW</td><td>100±20</td><td>10±2</td><td>-</td><td>-</td><td>-</td><td>-</td></tr><tr><td>M31</td><td>130±40</td><td>10±2</td><td>24.407±0.032 (ref. 27)</td><td>48.9±10.5 (ref. 29)</td><td>-36.9±8.1 (ref. 29)</td><td>-301±1 (ref. 94)</td></tr><tr><td>M33</td><td>30±10</td><td>10±2</td><td>24.67±0.07 (ref. 26)</td><td>31±19 (ref. 13)</td><td>-29±16 (ref. 13)</td><td>-179.2±1.8 (ref. 1)</td></tr><tr><td>LMC</td><td>15±5</td><td>10±2</td><td>18.477±0.026 (ref. 25)</td><td>1,910±20 (ref. 28)</td><td>229±47 (ref. 28)</td><td>262.2±3.4 (ref. 1)</td></tr><tr><td>SMC</td><td>1.5±0.5</td><td>10±2</td><td>18.99±0.03 (ref. 95)</td><td>722±63 (ref. 28)</td><td>-1,117±61 (ref. 28)</td><td>145.6±0.6 (ref. 96)</td></tr></tbody></table>

velocity distribution, and we used standard assumptions about their concentrations, c, and velocity dispersion profiles$^{23}$. We calculated the dynamical friction using an analytic formalism (see Methods for a detailed description). Although non-gravitational processes (such as gas drag and star formation leading to increased central densities and so on) 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}$.

The simplest model for the MW-M31 orbit, as considered by ref. 11, contains only the two main galaxies. Because of the planar symmetry, only five parameters are required: the two masses, their initial separation and the two velocity components. In general, an orbit with $N > 2$ galaxies requires $3(N-1)$ coordinates and $3(N-1)$ velocity components, along with the $N$ masses. More recently, refs. 12,13 considered three-body orbits by including M33, the third most massive LG galaxy. Reference 13 also considered the Large Magellanic Cloud (LMC), the fourth most massive LG galaxy. They still concluded that a merger is certain. Here, we 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 LG, 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 used Monte Carlo 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.

Our 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. We express all masses in terms of $M_{200}$, the mass enclosed within a sphere whose mean density is $200 \times$ the critical density. 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_{☉}$. We adopted $M_{200} = 1 \pm 0.2 \times 10^{12} M_{☉}$ (excluding the mass of the LMC, which we treated separately). For all other galaxies, there is considerably more uncertainty. For M31, we adopted $M_{200} = 1.3 \pm 0.4 \times 10^{12} M_{☉}$. For M33, we assumed $M_{200} = 3 \pm 1 \times 10^{11} M_{☉}$, and for the LMC, we assumed $M_{200} = 1.5 \pm 0.5 \times 10^{11} M_{☉}$. See Methods for a review of mass estimates.

The line-of-sight velocities ($v_{los}$) are known well, and we adopted the values given by ref. 1. For the distance moduli $\mu$, we chose the most accurate and precise recent values in the literature: from ref. 25 for the LMC and from ref. 26 for M33. We used for M31 recent HST cepheid results$^{27}$.

We use the notation $\mu_b$ for the proper motion in declination (dec.) and $\mu_a^*$ = $\mu_a \cos \delta$ for the proper motion in right ascension (RA). We used the HST proper motions of ref. 28 for the LMC and the combined HST and Gaia Data Release 2 (DR2) proper motions of ref. 13 for M33. Our fiducial M31 proper motions are the most precise values in the literature based on Gaia Data Release 3 (DR3) astrometry$^{29}$. However, we found very similar results using the combined Gaia DR2 and HST

<figure><img src="image_4.png" /><figcaption>**Fig. 1 | Possible future MW-M31 orbits.** Coloured lines show probability densities for the positions of the MW (bright yellow to blue) and M31 (bright orange to purple) in 100 Monte Carlo samples of the fiducial model, integrated over 10 Gyr or until the merger. On the top and bottom rows, respectively, trajectories are projected in the orbital plane or perpendicular to the orbital plane defined by the initial positions and velocities of the MW and M31. White markers denote MW-M31 mergers. Colours were assigned according to the probability of each galaxy being at the projected location during the next 10 Gyr for the 100 orbits shown. More orbits crossing a given location and with lower projected speed resulted in brighter colours. In the left-hand column, we show the MW-M31 two-body system, whereas in the right, we show the four-body system with the MW, M31, M33 and the LMC.</figcaption></figure>

proper motions of ref. 13, as shown in Extended Data Fig. 1. All assumed values are listed in Table 1. We assumed galactocentric coordinates $RA_{GC} = 266.41^{\circ}$, $dec_{GC} = -28.94^{\circ}$ (ref. 30), $d_{GC} = 8.122$ kpc (ref. 31), and a velocity of (12.9, 245.6, 7.78) km s$^{-1}$ with respect to the Sun$^{30-32}$.

To 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 truncated all distributions at $\pm 2\sigma$. This truncation increased the probability of a MW-M31 merger by only -10% (Extended Data Fig. 2). For the fiducial model, we use 50,000 Monte Carlo samples, whereas for all other variants, we used 2,500 samples to ensure that all statistical errors on the merger rate were below 1%.

## Evolution of the MW-M31-M33-LMC system

The Monte Carlo initial conditions were integrated numerically using a symplectic direct scheme (see Methods for details). Figure 1 shows 100 realizations 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
<figure><img src="image_5.png" /></figure>

**Fig. 2 | Distance between the MW and M31.** On each panel, we show 100 realizations of our fiducial Monte Carlo 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 Fig. 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 ($f_M$). 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 probable values of every variable, either assuming the Gaia DR3 proper motions<sup>29</sup> of the fiducial model (solid) or HST + Gaia DR2 proper motions<sup>13</sup> (dashed). ML, fiducial model.

of orbits and, additionally, for the MW-M31-M33 and MW-M31-LMC three-body systems.

The MW-M31 two-body orbit evolves in a plane and leads to a merger in slightly less than half of cases, most 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 MW-M31-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%.

For each system, we also show the single orbit obtained by adopting the most probable value of each observable, either in the fiducial model or assuming HST + Gaia DR2 proper motions for M31 (ref. 13). For the MW-M31-M33 system, which is the one considered by refs. 12,13, and using the same parameter values as considered in ref. 13, corresponding to the solid white line in the bottom left panel of Fig. 2, we reproduced a very similar time for the first pericentre, ~4.5 Gyr, and a similar first pericentre distance, ~100 kpc, as reported by ref. 13.

However, despite using newer and even more precise measurements, when we performed a Monte Carlo analysis, we found 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, which considered only the most probable values without accounting for the numerous and significant uncertainties.

Figure 3 shows 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, which is the probability that no merger with M31 occurs. In the fiducial model, we consider that a merger has occurred when the distance between any two galaxies is below 20 kpc, but our results are not sensitive to this particular choice (see Methods for details). Because of the effect of dynamical friction, we found that there are two distinct possibilities for the eventual fate of the MW and M31: (1) Orbits that come within less than ~200 kpc eventually merge, which would probably lead to the formation of an intermediate-mass elliptical galaxy<sup>10,33,34</sup>. For systems that merge, we found a median time of 7.6 Gyr in the fiducial model, or 8.0 Gyr adopting a 10-kpc threshold. (2) 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 probable.

## The roles of M33 and the LMC

The distinct effects of each satellite on the MW-M31 orbit are illustrated in Figs. 4 and 5, which 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 extra 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 ref. 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, whereas the inclusion of the LMC also provides important momentum perpendicular to the MW-M31 plane. In our analysis, the LMC is certain to merge with the MW, and M33 is highly probable (~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, whereas the net effect of adding the LMC is to decrease it.

## Sources of uncertainty

Figure 6 shows 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 Fig. 3. We found 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 MW and 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 probable proper motions (assuming no errors) lead to a merger in only ~2/3 of cases.

Future, more precise, proper motion measurements may significantly change the expected outcome: they could make the merger either more or less probable. However, if they fall within ±1σ of the currently most probable values, even precise M31 proper motion measurements alone will not suffice to determine the outcome. Even the comparatively high precision of the line-of-sight velocity and distance moduli for M33 and M31 contribute substantial uncertainty, with the probability of a merger varying between 40% and over 60% for different combinations in the ±2σ ranges around the most probable values.

Given the considerable measurement errors, note that cosmological simulations result in a present-day MW-M31 transverse velocity prior of $v_t = 75_{-40}^{+65}$ km s⁻¹ (ref. 36), which is compatible with recent measurements<sup>13,29</sup> but not with a perfectly radial orbit<sup>12</sup>. There is, thus, no reason to assume that the transverse velocity measured using Gaia DR3 (most probable value $v_t = 76$ km s⁻¹) is overestimated or to expect
<figure><img src="image_6.png" /><figcaption>**Fig. 3 | Distributions of the merger time $t_m$, the minimum distance $min(d_{M31-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 ($d_m$) 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</figcaption></figure>

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, whereas the vast majority of the remaining systems do not approach closer than 200 kpc. The survival rate of the MW drops sharply between -6 and 9 Gyr and levels off afterwards.

that more precise measurements will result in a lower value. Also note that, as shown in Fig. 6, a perfectly radial present-day M31-MW motion ($v_t = 0, \mu_\alpha^* = 38.03$ mas yr$^{-1}$ and $\mu_\beta = -21.37$ mas yr$^{-1}$) does not 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 LG is uncertain. Intriguingly, we found an almost equal probability for the widely publicized merger scenario (albeit with a later median time to merger) and one where the MW and M31 survive unscathed. We reached this conclusion by including the LMC and, importantly, considering the relevant uncertainties in the observables.

Our results are not sensitive to the necessary choices of gravitational softening (Extended Data Fig. 4), merger threshold (Extended Data Fig. 5) or dynamical friction scheme (Extended Data Fig. 6). We used N-body simulations to validate the procedure for several examples (Extended Data Figs. 7–9). Importantly, our results were obtained using methods very much like those previously used to predict the MW-M31 merger, and when we applied our methods to the same restricted initial conditions, we obtained very similar results.

Although we have shown that considerable uncertainty results from the proper motion measurements of M31, we also found 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 indicate that ~25% of the bound mass of the LG is outside

the two main haloes³⁷. The next most massive individual LG galaxies that could impact the MW-M31 orbit are M32 (an M31 satellite and possible merger remnant³⁸) 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 found that including the SMC, for which proper motion measurements are available, does not significantly change the MW-M31 merger probability (Extended Data Fig. 10).

In our calculations, we made several simplifying assumptions. We considered isotropic haloes with no substructures, constant masses and constant concentrations, and we employed idealized treatments of dynamical friction and mergers. The unaccounted effects of substructures as well as those 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 cosmologically constrained simulations³⁹⁻⁴¹ that could account for these effects. Moreover, the assumptions and simplifications we have made here are probably conservative regarding our central claim, that there is considerable uncertainty about the MW-M31 merger.

Upcoming Gaia data releases will improve the proper motion constraints. Moreover, the mass models are continually being refined. However, it is clear that Galactic eschatology is still in its infancy, and substantial work is required before the eventual fate of the LG can be predicted with any certainty. As it stands, proclamations of the impending demise of our Galaxy seem greatly exaggerated.

## Methods

Our results are based on numerically integrated initial conditions, which are, in turn, based on Monte Carlo samples of the observational
<figure><img src="image_7.png" /><figcaption>**Fig. 4 | Effects of M33 and the LMC on the trajectory of the MW.**</figcaption></figure>

As in Figs. 1 and 5, panels in the top row are projected in the orbital plane defined by the initial MW and M31 positions and velocities, whereas 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. In the centre and right, we show the orbit of the LMC and that of the MW before

(centre) and after (right) 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.

<figure><img src="image_8.png" /><figcaption>**Fig. 5 | Effects of M33 and the LMC on the trajectory of M31.** As in Figs. 1 and 4, panels in the top row are projected in the orbital plane defined by the initial MW and M31 positions and velocities, whereas 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. In the centre and on the left, we show the orbit of M33 and that of M31 before and after a</figcaption></figure>

possible merger with M33 in the MW-M31-M33 system. Compared to the two-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 results in very little motion perpendicular to the plane of the two-body orbit and reduces the transverse velocity of M33 with respect to the MW.
<figure><img src="image_9.png" /></figure>

**Fig. 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 in each panel, which indicates the sensitivity of the merger probability on the two corresponding variables. Top row, dependence on different masses. Second row, dependence on proper motions. Third row, 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 most strongly depends on $\mu_\delta(M31)$ and $\mu_\alpha^*(M31)$ but also varies significantly with $\mu_\delta(M33)$ and $\mu_\alpha^*(M33)$. The uncertainties in the distance moduli for M31 and M33 contribute to the uncertainty of the outcome, whereas the effect of the line-of-sight velocities is small.

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 MW and three other galaxies: M31, M33 and the LMC. We assumed that the sky coordinates (RA and dec.) for the centres of M31, M33 and the LMC are known. Furthermore, we assumed that the position of the Galactic centre with respect to the Sun is fixed at (RA, dec.) = (266.4051, -28.936175) (ref. 30) and d = 8.122 kpc (ref. 31) and that the galactocentric velocity of the Sun is (12.9, 245.6, 7.78) km s⁻¹ (refs. 30–32).

We created Monte Carlo samples for the remaining 20 variables: the four halo masses ($M_{MW}$, $M_{M31}$, $M_{M33}$ and $M_{LMC}$), the four halo concentration parameters c, the three distance moduli $\mu$, the three sets of proper motions $\mu_\delta$ and $\mu_\alpha^*$, and three line-of-sight velocities $v_{105}$. Sampling directly in the space of the observables rather than sampling in a space of derived variables, such as Cartesian coordinates or velocities, minimizes the effect of possible correlations.

To allow reproducibility and identification of individual orbits across our figures when changing the parameters of the model, such as the set of galaxies included, or the numerical parameters, such as the gravitational softening length and merger threshold, we re-initialized the pseudorandom number generator with the same values for each new set of Monte Carlo samples. Although generally we show only the first 100 orbits on each plot, all quoted probabilities are computed from at least 2,500 samples each, so that the statistical errors are less than 1%.

The most probable values of each variable together with the ±1σ uncertainty were either taken directly from single sources in the literature or, for masses and concentrations (see 'Halo concentrations'), estimated by us based on several sources. In creating our Monte Carlo samples, we assumed that each variable follows a Gaussian probability distribution. However, in our fiducial model, we truncated all distributions at ±2σ, corresponding to the central -95% of values. Although 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 probable values. Indeed, for some variables, untruncated distributions extend to unphysical values with finite probabilities. A truncation at 2σ ensured that all variables in the samples were physical.

Given 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 Fig. 2, our results are not sensitive to the truncation at ±2σ, and the merger probability was
only slightly lower when the distributions were not truncated. On the other hand, even a truncation at only ±1σ means that approximately 1/4 of orbits do not merge within 10 Gyr.

## Numerical integration

The orbits were 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 km s⁻¹ to travel 0.1 kpc, which is 1/200th of our merger threshold (see ‘Mergers’). Our results were not affected by the finite time step.

To account for the fact that the haloes are extended objects, the gravitational force between the haloes was softened with a constant softening length of 20 kpc, which is like the scale radius of an NFW halo in the mass range we considered. We also considered different choices for the softening, and we show in Extended Data Fig. 4 results with softening lengths of 10, 20 or 30 kpc. Using a softening length that was too small led to unphysical hard scattering events during close encounters, whereas a softening length that was too large artificially reduced the gravitational force. In the context of the LG, both of these effects could reduce the merger probability. However, we found 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 like that used in ref. 24. The classic Chandrasekhar formula assumes a point mass orbiting in the potential of a much more massive host halo, which is itself composed of much less massive particles. This approach has been expanded to account for extended satellites⁴², and we use the following expression to calculate the acceleration of a satellite due to dynamical friction¹⁶ once inside $r_{200}$ (the radius of a sphere that encloses a mass $M_{200}$) of the host halo:

$$ \frac{dv}{dt} = -\frac{4\pi G^2 M \rho \ln \Lambda}{v^2} \left[ \text{erf}(X) - \frac{2X}{\sqrt{\pi}} e^{-x^2} \right] \frac{v}{v}, \quad (1) $$

where 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 v and the one-dimensional-velocity dispersion σ of the host at the location of the satellite, and Λ is the Coulomb factor expressed as r/ε, with ε a scale length that depends on the density of the satellite. To determine ε, we adopted an empirical expression derived from N-body simulations⁴³:

$$ \epsilon = \begin{cases} 2.2r_s - 14 \text{ kpc}, & \text{if } r_s \ge 8 \text{ kpc}, \\ 0.45r_s, & \text{if } r_s < 8 \text{ kpc}, \end{cases} \quad (2) $$

where r_s is the scale radius of the NFW halo of the satellite. Finally, to approximate the velocity dispersion of the host at the location of the satellite, we used the expression derived in ref. 23 for NFW haloes.

Standard 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.

However, in the LG context where galaxies and haloes of similar mass are interacting, this introduces an inconsistency. In particular, in the (relatively probable) scenario that M31 has a mass close to that of 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.

To make the calculation more symmetrical, in our dynamical friction calculation we distributed the dynamical friction force proportionally between the satellite (s) and host (h), thus conserving the total momentum:

$$ \frac{dv_s}{dt} = a_{DF} \frac{M_h}{M_s + M_h}, \quad (3) $$

$$ \frac{dv_h}{dt} = -a_{DF} \frac{M_s}{M_s + M_h}, \quad (4) $$

where v_s and M_s are the velocity and mass of the satellite, v_h and M_h are the velocity and mass 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, whereas in the limit that the galaxies have equal masses, both receive equal and opposite accelerations. In effect, although the standard hierarchical scheme assumes an infinite mass ratio so that the centre of mass of the system is identical to that of the host, the ‘proportional’ scheme applies dynamical friction in the centre-of-mass frame of the system. The magnitude of the dynamical friction force remains the same.

A 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 in ‘Halo concentrations’, the concentration of the satellite has a greater impact on the dynamical friction force. For two haloes with 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 assigned) concentrations, the choice of whether to calculate 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 assumed identical distributions of concentration parameters for all galaxies.

Extended Data Fig. 6 compares 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 a MW-M31 merger is highly unlikely without dynamical friction. 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 was included, the evolution of each orbit was very similar in the hierarchical and proportional schemes, and the merger rate was not significantly affected by the exact choice of scheme.

Note that our semi-analytical approach to dynamical friction is still very simplistic, and although 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¹⁴. A precise prediction of the MW-M31 orbit will probably 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 important, and our simple approach is no longer appropriate for predicting the remaining orbital evolution. Our aim was 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 assumed that any system that passes below a threshold distance will eventually merge, and we identified this time as a lower limit for the time of the merger.

In our fiducial model, we adopted a value of 20 kpc as the merger threshold for all galaxy interactions. When such a merger occurred between two galaxies, we combined the masses and momenta of the two galaxies at their common centre of mass and continued the integration. The concentration parameter was then set to that of the more massive galaxy.

Extended Data Fig. 5 compares the sensitivity of our results to adopting merger thresholds of 10, 20 and 30 kpc. With a threshold of 10 kpc, we found a slightly reduced MW-M31 merger probability. However, raising the threshold from 20 to 30 kpc had no significant impact on our results. This confirms our assertion in the main text: 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.

The 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 to escape to a large apogentre. Although this may be a possible scenario for the LG, 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 assumed that these orbits also merge, and with an even larger threshold of 30 kpc, we would still predict a similar merger rate.

For comparison, the best-studied merger of two galaxies with stellar masses like those of the MW and M31 are the Antennae galaxies. Their past evolution is reproduced with an orbit with pericentre -10 kpc (ref. 45), which is predicted to lead to coalescence -1.3 Gyr later$^{46}$.

## Fate of the LMC and M33

Although our main focus was on the future evolution of the MW-M31 orbit, we naturally also make predictions for the evolution of the LMC and M33. We found with the fiducial model that the LMC is certain to merge with the MW before any eventual MW-M31 merger. With a merger threshold of 20 kpc, we found a median time for the LMC-MW merger of 1.3 Gyr, whereas with a threshold of 10 kpc, we found a median time of 1.9 Gyr. For M33, with a merger threshold of 20 kpc, we found an -86% chance of a merger with M31 and a median time of 3.3 Gyr, whereas with a merger threshold of 10 kpc, we found an -83% probability of a merger with M31 and a median time of 3.9 Gyr. In both cases, we also found 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. Note that our simulations were not designed to study these mergers in detail and ignore, for example, the impact of the disk of the MW. Nevertheless, they broadly agree with the results of ref. 24, who previously studied the LMC-MW encounter in the presence of M31.

## Other galaxies

The next most massive LG member galaxy for which proper motion data$^{28}$ are available is the SMC, a satellite galaxy of the MW that was probably accreted together with the LMC. We repeated our analysis including the SMC as a fifth system and show the results in Extended Data Fig. 10. Adding the SMC, whose mass is approximately 10% of that of the LMC, had no significant effect on the merger rate. The SMC properties used are listed along with those of the other galaxies in Table 1.

## HST + DR2 proper motions

Because the M31 proper motions have the largest impact on the probability of a MW-M31 orbit and for easier comparison with earlier works, particularly ref. 13, we repeated our analysis adopting the HST + Gaia DR2 M31 proper motions of ref. 13. Extended Data Fig. 1 shows the corresponding evolution of the MW-M31 distance in the same two-body MW-M31, three-body MW-M31-M33 and MW-M31-LMC, and four-body MW-M31-M33-LMC systems. The results can be directly compared to Fig. 2, which shows the same quantities in our fiducial model that used Gaia 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 MW-M31-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.

## Comparison to N-body simulations

To validate the semi-analytical orbital integration method and our assumptions about dynamical friction, gravitational softening and mergers used throughout this paper, we also performed N-body simulations for five of the four-body MW-M31-M33-LMC orbits. The initial masses, concentrations, positions and velocities were chosen from among the same sample of fiducial systems shown in Figs. 1 and 2. The five examples were selected based on the outcome predicted by the semi-analytical model: one system with a predicted early merger at first pericentre ('early merger'), two with typical merger times ('typical merger 1' and 'typical merger 2') and two that do not merge within 10 Gyr ('no merger 1' and 'no merger 2').

In each case, we created a set of NFW haloes, truncated around $r_{200}$, with a sigmoid function. We used Monte Carlo sampling for the distribution function $f_i(\epsilon)$ for the NFW density profile $\rho_{\text{NFW}}$ with discrete particles using the Eddington formula:

$$f_i(\epsilon) = \frac{1}{2\sqrt{2\pi^2}} \int_{\Phi_T=0}^{\Phi_T=\epsilon} \frac{d^2 \rho_{\text{NFW},i}}{d\Phi_T^2} \frac{d\Phi_T}{\sqrt{\epsilon - \Phi_T}}, \quad (5)$$

where $\epsilon$ is the relative energy and $\Phi_T$ is the total gravitational potential.

The simulations were performed with the Gadget-4 code$^{47}$ using a particle mass resolution of $10^{6.5} M_{\odot}$ (resulting in $-10^6$ particles per simulation) and a gravitational softening of 200 pc. The simulations were performed for the same 10-Gyr period considered throughout this work, with 100 outputs at intervals of 100 Myr. We determined the centres of the systems in each output using the shrinking spheres method implemented in Pynbody. For comparison to the orbits predicted by the semi-analytical integration, these positions were linearly interpolated between snapshots. As in the fiducial model of the semi-analytical integration, we identified the time of a possible MW-M31 merger as the time when the two halo centres were separated by less than 20 kpc.

Extended Data Fig. 7 shows the projected matter distribution in the initial conditions in each of the five systems in the centre-of-mass frame. Reflecting the high accuracy of the positions, the five initial conditions are visually very similar. However, as shown in Extended Data Fig. 8, after 10 Gyr, the density distributions have diverged. In particular, in agreement with the prediction of our semi-analytical integration, the MW and M31 have merged in the first three cases, whereas in the last two cases, the two large haloes are still clearly separated by several hundreds of megaparsecs. For the early merger, the newly formed halo looks spherical, whereas in the two mergers predicted to happen later, the newly formed halo still contains substructures.

In a more detailed comparison, Extended Data Fig. 9 compares the evolution of the MW-M31 distances in the N-body simulation and in
the semi-analytical integration, as in Fig. 2. The left-hand panel shows the three orbits that lead to a MW-M31 merger, and the right-hand panel shows the two orbits with no merger within 10 Gyr, according to the semi-analytical prediction. Each pair of solid and dashed lines shows the evolution using the *N*-body and semi-analytical methods, respectively. The early merger case proceeds almost identically using both methods. The typical merger 1 proceeds nearly identically until the first pericentre, but the apocentre in the semi-analytical case is slightly larger, leading to a somewhat later second pericentre, when the merger occurs in both cases. Conversely, in the typical merger 2, the merger is delayed in the *N*-body case, following a slightly larger first pericentre and a larger apocentre than predicted by the semi-analytical method. However, in all three cases, the galaxies merge within 10 Gyr using both methods.

The two no-merger cases, shown in the right-hand panel, also evolve in a similar way in both methods. Both reach somewhat smaller pericentres in the *N*-body simulation compared to the semi-analytical integration, but in both cases, dynamical friction is not strong enough to allow the orbit to decay within 10 Gyr using either method, and the outcome is the same.

As expected, and consistent with previous works$^{23,44}$, the evolution in the *N*-body simulations and the semi-analytical methods are not identical. However, in all five cases, they are qualitatively similar. In particular, the agreement is sufficient to show that the uncertainty about the outcome is driven by the uncertainty of the initial conditions, rather than caused by the necessary assumptions and simplifications inherent to the semi-analytical integration.

## Sources for the masses

As discussed in Fig. 6, the assumed masses of all four galaxies and their associated uncertainties have a strong impact on the probable evolution of the MW-M31 orbit in our fiducial model. Here, we review the mass measurements and show the implications for some alternative scenarios.

**Milky Way.** 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} \times 10^{12} M_{\odot}$ (ref. 48), $1.51^{+0.45}_{-0.40} \times 10^{12} M_{\odot}$ (ref. 49), $1.23^{+0.21}_{-0.18} \times 10^{12} M_{\odot}$ (ref. 50) and $1.11^{+0.1}_{-0.1} \times 10^{12} M_{\odot}$ (ref. 51) (where the latter two works also used a simulation-based prior).

Estimates made using rotation curves include $1.08^{+0.20}_{-0.14} \times 10^{12} M_{\odot}$ based on Gaia DR2 (ref. 52), $0.89^{+0.1}_{-0.08} \times 10^{12} M_{\odot}$ using stars in the galkin catalogue$^{23}$, $0.822 \pm 0.052 \times 10^{12} M_{\odot}$ using classical cepheids$^{54}$ and $1.08^{+0.12}_{-0.11} \times 10^{12} M_{\odot}$ from the H3 survey and Gaia DR3 (ref. 55). Other recent measurements include $1.54^{+0.75}_{-0.44} \times 10^{12} M_{\odot}$ using combined Gaia DR2 and HST kinematics of globular clusters$^{36}$, $1.16 \pm 0.24 \times 10^{12} M_{\odot}$ (including the mass of the LMC) using the kinematics of halo stars, $1.26^{+0.40}_{-0.22} \times 10^{12} M_{\odot}$ from high-velocity RR Lyrae stars$^{57}$ and $1.19^{+0.49}_{-0.32} \times 10^{12} M_{\odot}$ from a Bayesian estimate using dwarf galaxy kinematics from several sources$^{58}$. Reference 36 contains an overview of recent measurements, and ref. 59 includes a comprehensive review of earlier results.

An analysis from several different tracers, particularly from the latest studies using the latest Gaia observations, consistently point towards a MW total mass that is close to $10^{12} M_{\odot}$. Both the simple mean and median of the above measurements are $1.16 \times 10^{12} M_{\odot}$, but note that most methods of measuring the total mass of the MW include the mass of the LMC (at a galactocentric distance of ~50 kpc, well within the virial radius, $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_{\odot}$ (see our discussion on the LMC mass below). We adopted $M_{MW} = 1.0 \pm 0.2 \times 10^{12} M_{\odot}$ excluding the LMC in our fiducial model, reflecting the consensus of recent observations.

## M31

Although, there are no accurate Gaia proper motions for M31, several studies have estimated its mass. Most recently, ref. 60 measured a total mass of $1.14^{+0.51}_{-0.35} \times 10^{12} M_{\odot}$ using rotation curves based on LAMOST data release 9 and data provided by the Dark Energy Spectroscopic Instrument. Reference 61 measured $1.4 \pm 0.4 \times 10^{12} M_{\odot}$ from satellite kinematics, ref. 62 derived $1.2^{+0.9}_{-0.7} \times 10^{12} M_{\odot}$ from the kinematics of M31 dwarf spheroids and ref. 63 measured a total mass of $1.05^{+0.15}_{-0.15} \times 10^{12} M_{\odot}$ from fitting the spectral energy distribution and from the rotation curve and the kinematics of outer globular clusters and satellite galaxies. Reference 64 measured $1.0 \times 10^{12} M_{\odot}$ from the H I rotation curve, and ref. 65 measured $2.0^{+0.4}_{-0.3} \times 10^{12} M_{\odot}$ from the kinematics of the Giant Southern Stream. Reference 66 measured $1.35^{+0.15}_{-0.15} \times 10^{12} M_{\odot}$, and ref. 67 found $1.4^{+0.2}_{-0.2} \times 10^{12} M_{\odot}$, both using outer halo globular clusters, whereas ref. 68 found $0.8 \pm 0.1 \times 10^{12} M_{\odot}$ from high-velocity planetary nebulae, and ref. 69 found $1.39 \pm 0.26 \times 10^{12}$ for the total mass by 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 \times 10^{12} M_{\odot}$. Our adopted mass of $M_{M31} = 1.3 \pm 0.4 \times 10^{12} M_{\odot}$ in the fiducial model reflects the broad consensus of M31 mass estimates using different methods but also the considerable remaining uncertainty.

**M33.** Mass estimates of M33 are much more sparse. More than 20 years ago, refs. 70,71 measured a dark matter mass of $5 \times 10^{10} M_{\odot}$, extrapolated out to a virial mass of $5 \times 10^{11} M_{\odot}$ from the measured H I rotation curve, but they noted that this results in a very low baryon fraction. More recently, ref. 72 obtained a similar result using the Hα rotation curve, but noted that because the measurements extend only 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 resulted in a substantially lower total mass of $1.7 \pm 0.55 \times 10^{11} M_{\odot}$ (refs. 73,74). Citing both the direct measurements and abundance-matching, ref. 74 adopted a mass range of $0.8–3.2 \times 10^{11} M_{\odot}$ in their dynamical models of the M33-M31 interaction, whereas refs. 13,75 assumed a total mass of $2.5 \times 10^{11} M_{\odot}$. We adopted $M_{M33} = 3 \pm 1 \times 10^{11} M_{\odot}$ in the fiducial model, which is marginally compatible with both the extrapolated masses from rotation curve measurements and the results of abundance-matching and is in line with previous studies. However, the results from the two methods are certainly in tension.

**LMC.** Recent mass estimates of the LMC based on its effect on Galactic stellar streams include $1.38^{+0.37}_{-0.24} \times 10^{11} M_{\odot}$ (ref. 76), $1.30 \pm 0.3 \times 10^{11} M_{\odot}$ (ref. 77), $1.88^{+0.4}_{-0.35} \times 10^{11} M_{\odot}$ (ref. 78) and $1.29^{+0.28}_{-0.23} \times 10^{11} M_{\odot}$. Reference 79 obtained good agreement between the perturbations of MW halo stars with an LMC mass of $1.5 \times 10^{11} M_{\odot}$. Using the abundance of probable satellites of the LMC, ref. 80 obtained a lower limit of $1.24 \times 10^{11} M_{\odot}$, and using the kinematics of satellites associated with the LMC, ref. 81 found $1.65^{+0.47}_{-0.49} \times 10^{11} M_{\odot}$. Most recently, ref. 82 used 30 LMC globular clusters to infer a total mass of $1.80^{+1.05}_{-0.54} \times 10^{11} M_{\odot}$, whereas ref. 35 found an even higher value of $2.5^{+0.9}_{-0.8} \times 10^{11} M_{\odot}$ using a timing argument, which is most consistent with the Hubble flow around the LG. The simple mean of the above values is $1.6 \times 10^{11} M_{\odot}$, and the median is $1.5 \times 10^{11} M_{\odot}$. A comprehensive recent review on the effect of the LMC on the MW, including a discussion of mass measurements, was given by ref. 83, who concluded that the LMC mass is probably in the range 1–2 × $10^{11} M_{\odot}$, which matches the choice $M_{LMC} = 1.5 \pm 0.5 \times 10^{11} M_{\odot}$ in our fiducial model.

**Other mass.** Cosmological simulations show that the total mass of LG analogues is ~25+9-6% larger than the sum of the $M_{200}$ masses of the MW and M31 (ref. 37). The amount of extra matter (in addition to that of M33, which adds around 10–15%) and its distribution within the LG is not known. This forms part of the extra uncertainty that we cannot address here.
## Halo concentrations

As explained above, we assume that for the purpose of integrating their orbits, all galaxies are represented by NFW haloes²² with the total masses M₂₀₀ defined above. We assumed concentration parameters of 10 ± 2 for all galaxies, consistent with the results of cosmological simulations⁸⁴⁻⁸⁸ in this mass range. Although cosmological simulations show that the average concentration parameter depends on mass, the halo-to-halo scatter is significantly larger than the change in the mean concentration over this narrow mass range⁸⁴. The 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₂₀₀/c, of the 'satellite' that in equation (2). Extended Data Fig. 3 shows 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 found that the merger probability is reduced if the concentration is below -8, that is below -1σ of our fiducial value. We found no significant dependence on the merger probability on the concentration assumed for M33 or the LMC.

## Data availability

All data used in this work are publicly available and provided as part of the analysis code listed under 'Code Availability'.

## Code availability

The analysis in this paper was performed using Python 3.10, and made extensive use of the following open-source libraries: Astropy v.6.0.1 (ref. 89), Matplotlib v.3.8.3, NumPy v.1.26.4 (ref. 90), SciPy v.1.13.0 (ref. 91), Colossus v.1.3.5 (ref. 92) and PyNbody v.2.0.0 (ref. 93). A documented Jupyter notebook containing the code used to produce all the figures in this paper is available via https://github.com/TillSawala/MW-M31.

## References

1. McConnachie, A. W. The observed properties of dwarf galaxies in and around the Local Group. *Astron. J.* **144**, 4 (2012).

2. Newton, O. et al. The undiscovered ultradiffuse galaxies of the Local Group. *Astrophys. J. Lett.* **946**, L37 (2023).

3. Sawala, T. et al. The APOSTLE simulations: solutions to the Local Group's cosmic puzzles. *Mon. Not. R. Astron. Soc.* **457**, 1931–1943 (2016).

4. Slipher, V. M. The radial velocity of the Andromeda Nebula. *Lowell Obs. Bull.* **2**, 56–57 (1913).

5. Hubble, E. P. A spiral nebula as a stellar system, Messier 31. *Astrophys. J.* **69**, 103–158 (1929).

6. Peebles, P. J. E., Phelps, S. D., Shaya, E. J. & Tully, R. B. Radial and transverse velocities of nearby galaxies. *Astrophys. J.* **554**, 104–113 (2001).

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. *Astrophys. J.* **633**, 894–898 (2005).

8. van der Marel, R. P. & Guhathakurta, P. M31 transverse velocity and Local Group mass from satellite kinematics. *Astrophys. J.* **678**, 187–199 (2008).

9. Sohn, S. T., Anderson, J. & van der Marel, R. P. The M31 velocity vector. I. Hubble Space Telescope proper-motion measurements. *Astrophys. J.* **753**, 7 (2012).

10. Dubinski, J., Mihos, J. C. & Hernquist, L. Using tidal tails to probe dark matter halos. *Astrophys. J.* **462**, 576 (1996).

11. Cox, T. J. & Loeb, A. The collision between the Milky Way and Andromeda. *Mon. Not. R. Astron. Soc.* **386**, 461–474 (2008).

12. van der Marel, R. P. et al. The M31 velocity vector. II. Radial orbit toward the Milky Way and implied Local Group mass. *Astrophys. J.* **753**, 8 (2012).

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. *Astrophys. J.* **872**, 24 (2019).

14. Cowen, R. Andromeda on collision course with the Milky Way. *Nature* https://doi.org/10.1038/nature.2012.10765 (2012).

15. Harvey-Smith, L. *When Galaxies Collide* (Melbourne Univ. Press, 2020).

16. Binney, J. & Tremaine, S. Galactic Dynamics 2nd edn (Princeton Univ. Press, 2008).

17. Eicher, D. J. *The New Cosmos: Answering Astronomy's Big Questions* (Cambridge Univ., 2015).

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

19. Belokurov, V., Erkal, D., Evans, N. W., Koposov, S. E. & Deason, A. J. Co-formation of the disc and the stellar halo. *Mon. Not. R. Astron. Soc.* **478**, 611–619 (2018).

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. *Nat. Astron.* **4**, 965–973 (2020).

21. Naidu, R. P. et al. Reconstructing the last major merger of the Milky Way with the H3 survey. *Astrophys. J.* **923**, 92 (2021).

22. Navarro, J. F., Frenk, C. S. & White, S. D. M. A universal density profile from hierarchical clustering. *Astrophys. J.* **490**, 493–508 (1997).

23. Zentner, A. R. & Bullock, J. S. Halo substructure and the power spectrum. *Astrophys. J.* **598**, 49–72 (2003).

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. *Mon. Not. R. Astron. Soc.* **483**, 2185–2196 (2019).

25. Pietrzyński, G. et al. A distance to the Large Magellanic Cloud that is precise to one per cent. *Nature* **567**, 200–203 (2019).

26. Ou, J.-Y. et al. A distance measurement to M33 using optical photometry of Mira variables. *Astron. J.* **165**, 137 (2023).

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. *Astrophys. J.* **920**, 84 (2021).

28. Kallivayalil, N., van der Marel, R. P., Besla, G., Anderson, J. & Alcock, C. Third-epoch Magellanic Cloud proper motions. I. Hubble Space Telescope/WFC3 data and orbit implications. *Astrophys. J.* **764**, 161 (2013).

29. Salomon, J. B. et al. The proper motion of Andromeda from Gaia EDR3: confirming a nearly radial orbit. *Mon. Not. R. Astron. Soc.* **507**, 2592–2601 (2021).

30. Reid, M. J. & Brunthaler, A. The proper motion of Sagittarius A*. II. The mass of Sagittarius A*. *Astrophys. J.* **616**, 872–884 (2004).

31. GRAVITY Collaboration et al. Detection of the gravitational redshift in the orbit of the star S2 near the Galactic centre massive black hole. *Astron. Astrophys.* **615**, L15 (2018).

32. Drimmel, R. & Poggio, E. On the solar velocity. *Res. Notes AAS* **2**, 210 (2018).

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. *Astron. Astrophys.* **358**, 471–480 (2000).

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. *Proc. Int. Astron. Union* **14**, 161–164 (2020).

35. Peñarrubia, J., Gómez, F. A., Besla, G., Erkal, D. & Ma, Y.-Z. A timing constraint on the (total) mass of the Large Magellanic Cloud. *Mon. Not. R. Astron. Soc.* **456**, L54–L58 (2016).

36. Sawala, T., Teeriaho, M. & Johansson, P. H. The Local Group's mass: probably no more than the sum of its parts. *Mon. Not. R. Astron. Soc.* **521**, 4863–4877 (2023).
37. Sawala, T., Peñarrubia, J., Liao, S. & Johansson, P. H. The timeless timing argument and the total mass of the Local Group. *Mon. Not. R. Astron. Soc.* **526**, L77–L82 (2023).

38. D'Souza, R. & Bell, E. F. The Andromeda galaxy's most important merger about 2 billion years ago as M32's likely progenitor. *Nat. Astron.* **2**, 737–743 (2018).

39. Libeskind, N. I. et al. The HESTIA project: simulations of the Local Group. *Mon. Not. R. Astron. Soc.* **498**, 2968–2983 (2020).

40. Sawala, T. et al. The SIBELIUS project: a pluribus unum. *Mon. Not. R. Astron. Soc.* **509**, 1432–1446 (2022).

41. Wempe, E. et al. Constrained cosmological simulations of the Local Group using Bayesian hierarchical field-level inference. *Astron. Astrophys.* **691**, A348 (2024).

42. White, S. D. M. A note on the minimum impact parameter for dynamical friction involving spherical clusters. *Mon. Not. R. Astron. Soc.* **174**, 467–470 (1976).

43. Jethwa, P., Erkal, D. & Belokurov, V. A Magellanic origin of the DES dwarfs. *Mon. Not. R. Astron. Soc.* **461**, 2212–2233 (2016).

44. Boylan-Kolchin, M., Ma, C.-P. & Quataert, E. Dynamical friction and galaxy merging time-scales. *Mon. Not. R. Astron. Soc.* **383**, 93–101 (2008).

45. Karl, S. J. et al. One moment in time—modeling star formation in the Antennae. *Astrophys. J. Lett.* **715**, L88–L93 (2010).

46. Lahén, N., Johansson, P. H., Rantala, A., Naab, T. & Frigo, M. The fate of the Antennae galaxies. *Mon. Not. R. Astron. Soc.* **475**, 3934–3958 (2018).

47. Springel, V., Pakmor, R., Zier, O. & Reinecke, M. Simulating cosmic structure formation with the GADGET-4 code. *Mon. Not. R. Astron. Soc.* **506**, 2871–2949 (2021).

48. Callingham, T. M. et al. The mass of the Milky Way from satellite dynamics. *Mon. Not. R. Astron. Soc.* **484**, 5453–5467 (2019).

49. 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. *Mon. Not. R. Astron. Soc.* **494**, 5178–5193 (2020).

50. Li, Z.-Z. et al. Constraining the Milky Way mass profile with phase-space distribution of satellite galaxies. *Astrophys. J.* **894**, 10 (2020).

51. Rodriguez Wimberly, M. K. et al. Sizing from the smallest scales: the mass of the Milky Way. *Mon. Not. R. Astron. Soc.* **513**, 4968–4982 (2022).

52. Cautun, M. et al. The Milky Way total mass profile as inferred from Gaia DR2. *Mon. Not. R. Astron. Soc.* **494**, 4291–4313 (2020).

53. Karukes, E. V., Benito, M., Iocco, F., Trotta, R. & Geringer-Sameth, A. A robust estimate of the Milky Way mass from rotation curve data. *J. Cosmol. Astropart. Phys.* **2020**, 033 (2020).

54. 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. *Astrophys. J. Lett.* **895**, L12 (2020).

55. Shen, J. et al. The mass of the Milky Way from the H3 survey. *Astrophys. J.* **925**, 1 (2022).

56. 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. *Astrophys. J.* **873**, 118 (2019).

57. Prudil, Z. et al. Milky Way archaeology using RR Lyrae and type II cepheids. II. High-velocity RR Lyrae stars and Milky Way mass. *Astron. Astrophys.* **664**, A148 (2022).

58. Slizewski, A. et al. Galactic mass estimates using dwarf galaxies as kinematic tracers. *Astrophys. J.* **924**, 131 (2022).

59. Wang, W., Han, J., Cautun, M., Li, Z. & Ishigaki, M. N. The mass of our Milky Way. *Sci. China Phys. Mech. Astron.* **63**, 109801 (2020).

60. Zhang, X., Chen, B., Chen, P., Sun, J. & Tian, Z. The rotation curve and mass distribution of M31. *Mon. Not. R. Astron. Soc.* **528**, 2653–2666 (2024).

61. Watkins, L. L., Evans, N. W. & An, J. H. The masses of the Milky Way and Andromeda galaxies. *Mon. Not. R. Astron. Soc.* **406**, 264–278 (2010).

62. Tollerud, E. J. et al. The SPLASH survey: spectroscopy of 15 M31 dwarf spheroidal satellite galaxies. *Astrophys. J.* **752**, 45 (2012).

63. Tamm, A., Tempel, E., Tenjes, P., Tihhonova, O. & Tuvikene, T. Stellar mass map and dark matter distribution in M 31. *Astron. Astrophys.* **546**, A4 (2012).

64. Chemin, L., Carignan, C. & Foster, T. H. Kinematics and dynamics of Messier 31. *Astrophys. J.* **705**, 1395–1415 (2009).

65. Fardal, M. A. et al. Inferring the Andromeda Galaxy's mass from its giant southern stream with Bayesian simulation sampling. *Mon. Not. R. Astron. Soc.* **434**, 2779–2802 (2013).

66. Veljanoski, J. et al. Kinematics of outer halo globular clusters in M31. *Astrophys. J. Lett.* **768**, L33 (2013).

67. Veljanoski, J. et al. The outer halo globular cluster system of M31. II. Kinematics. *Mon. Not. R. Astron. Soc.* **442**, 2929–2950 (2014).

68. 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. *Mon. Not. R. Astron. Soc.* **475**, 4043–4054 (2018).

69. Sofue, Y. Dark halos of M 31 and the Milky Way. *Publ. Astron. Soc. Jpn* **67**, 75 (2015).

70. Corbelli, E. & Salucci, P. The extended rotation curve and the dark matter halo of M33. *Mon. Not. R. Astron. Soc.* **311**, 441–447 (2000).

71. Corbelli, E. Dark matter and visible baryons in M33. *Mon. Not. R. Astron. Soc.* **342**, 199–207 (2003).

72. Kam, S. Z. et al. H I kinematics and mass distribution of Messier 33. *Astron. J.* **154**, 41 (2017).

73. Guo, Q., White, S., Li, C. & Boylan-Kolchin, M. How do galaxies populate dark matter haloes? *Mon. Not. R. Astron. Soc.* **404**, 1111–1120 (2010).

74. 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. *Mon. Not. R. Astron. Soc.* **464**, 3825–3849 (2017).

75. Patel, E., Carlin, J. L., Tollerud, E. J., Collins, M. L. M. & Dooley, G. A. ΛCDM predictions for the satellite population of M33. *Mon. Not. R. Astron. Soc.* **480**, 1883–1897 (2018).

76. Erkal, D. et al. The total mass of the Large Magellanic Cloud from its perturbation on the Orphan stream. *Mon. Not. R. Astron. Soc.* **487**, 2685–2700 (2019).

77. Vasiliev, E., Belokurov, V. & Erkal, D. Tango for three: Sagittarius, LMC, and the Milky Way. *Mon. Not. R. Astron. Soc.* **501**, 2279–2304 (2021).

78. Shipp, N. et al. Measuring the mass of the Large Magellanic Cloud with stellar streams observed by S⁵. *Astrophys. J.* **923**, 149 (2021).

79. Erkal, D. et al. Detection of the LMC-induced sloshing of the Galactic halo. *Mon. Not. R. Astron. Soc.* **506**, 2677–2684 (2021).

80. Erkal, D. & Belokurov, V. A. Limit on the LMC mass from a census of its satellites. *Mon. Not. R. Astron. Soc.* **495**, 2554–2563 (2020).

81. Correa Magnus, L. & Vasiliev, E. Measuring the Milky Way mass distribution in the presence of the LMC. *Mon. Not. R. Astron. Soc.* **511**, 2610–2630 (2022).

82. 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. *Astrophys. J.* **963**, 84 (2024).

83. Vasiliev, E. The effect of the LMC on the Milky Way system. *Galaxies* **11**, 59 (2023).

84. Neto, A. F. et al. The statistics of Λ CDM halo concentrations. *Mon. Not. R. Astron. Soc.* **381**, 1450–1462 (2007).

85. Ludlow, A. D. et al. The mass-concentration-redshift relation of cold dark matter haloes. *Mon. Not. R. Astron. Soc.* **441**, 378–388 (2014).

86. Wang, W. et al. Estimating the dark matter halo mass of our Milky Way using dynamical tracers. *Mon. Not. R. Astron. Soc.* **453**, 377–400 (2015).
87. 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. *Mon. Not. R. Astron. Soc.* **452**, 1217–1232 (2015).

88. Diemer, B. & Kravtsov, A. V. A universal model for halo concentrations. *Astrophys. J.* **799**, 108 (2015).

89. 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. *Astrophys. J.* **935**, 167 (2022).

90. Harris, C. R. et al. Array programming with NumPy. *Nature* **585**, 357–362 (2020).

91. Virtanen, P. et al. SciPy 1.0: fundamental algorithms for scientific computing in Python. *Nat. Methods* **17**, 261–272 (2020).

92. Diemer, B. COLOSSUS: a Python toolkit for cosmology, large-scale structure, and dark matter halos. *J. Suppl. Ser.* **239**, 35 (2018).

93. Pontzen, A., Roškar, R., Stinson, G. & Woods, R. pynbody: N-body/SPH analysis for Python. *Astrophysics Source Code Library* ascl:1305.002 (2013).

94. Watkins, L. L., Evans, N. W. & van de Ven, G. A census of orbital properties of the M31 satellites. *Mon. Not. R. Astron. Soc.* **430**, 971–985 (2013).

95. Graczyk, D. et al. The Araucaria Project. The distance to the Small Magellanic Cloud from late-type eclipsing binaries. *Astrophys. J.* **780**, 59 (2014).

96. Harris, J. & Zaritsky, D. Spectroscopic survey of red giants in the Small Magellanic Cloud. I. Kinematics. *Astron. J.* **131**, 2514–2524 (2006).

## Acknowledgements

We thank M. Cautun for his generous help and J. Nättilä for helpful suggestions. T.S., J.D. and J.H. are supported by the Research Council of Finland (Grant No. 354905). T.S., J.H. and P.H.J. are also supported by the Research Council of Finland (Grant No. 339127). J.D. is supported by an Erasmus+ grant. T.S. and C.S.F. are supported by the European Research Council through the Advanced Investigator Grant DMIDAS (Grant No. GA 786910) and an STFC Consolidated Grant (Grant No. ST/T000244/1). P.H.J., A.K., A.R. and R.W. also acknowledge support from the European Research Council (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 (STFC Capital Grant Nos. ST/K00042X/1, ST/PO02293/1, ST/RO02371/1 and ST/S002502/1), Durham University and the STFC (Operations Grant No. ST/RO00832/1). DiRAC is part of the UK national e-Infrastructure. We thank the authors of the open-source software listed in 'Code Availability'.

## Author contributions

T.S. planned the project and performed the analysis together with J.D. T.S., A.J.D., C.S.F. and P.H.J. planned the paper, and T.S. wrote the analysis code, created the figures and wrote the first draft. J.H. performed and analysed the N-body simulations based on initial conditions prepared by A.R. A.K. and A.R. contributed to the dynamical friction calculations. T.S., J.D., A.J.D., C.S.F., J.H., P.H.J., A.K., A.R. and R.W. jointly discussed and edited the paper.

## Funding

Open Access funding provided by University of Helsinki (including Helsinki University Central Hospital).

## Competing interests

The authors declare no competing interests.

## Additional information

Extended data is available for this paper at
https://doi.org/10.1038/s41550-025-02563-1.

Correspondence and requests for materials
should be addressed to Till Sawala.

Peer review information *Nature Astronomy* thanks Andreea Font, Matthias Steinmetz and the other, anonymous, reviewer(s) for their contribution to the peer review of this work.

Reprints and permissions information is available at
www.nature.com/reprints.

**Publisher's note** Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.

**Open Access** This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons.org/licenses/by/4.0/.

© The Author(s) 2025
<figure><img src="image_10.png" /></figure>

<figure><img src="image_11.png" /></figure>

<figure><img src="image_12.png" /></figure>

<figure><img src="image_13.png" /></figure>

**Extended Data Fig. 1 | Effect of using Gaia DR2+HST proper motions for M31.** Each panel shows the distance between the MW and M31, using HST+Gaia DR2 proper motions for M31 (ref. 13), analogous to Fig. 2 which uses the DR3 proper motions of the fiducial model. As in Fig. 2, solid and dashed white lines denote the orbits using the most likely values using either Gaia DR3 proper

motions²⁹ 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 fiducial model that uses the more precise Gaia DR3 proper motions, but the results are qualitatively similar.
<figure><img src="image_14.png" /></figure>

**Extended Data Fig. 2** | **Effect of truncating the assumed probability distributions of observables.** The distance between the MW and M31 (analogous to Fig. 2) in the four-body MW-M31-M33-LMC 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 Fig. 2), or left untruncated. The fraction of
systems that merge is only slightly increased when the distributions are clipped
at ±2σ or above.
<figure><img src="image_15.png" /></figure>

**Extended Data Fig. 3 | Effect of the concentration parameters on the merger probability.** In each case, we consider the four-body MW-M31-M33-LMC system, and plot the merger probability similar to Fig. 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

<figure><img src="image_16.png" /></figure>

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.
<figure><img src="image_17.png" /></figure>

**Extended Data Fig. 4 | Effect of the gravitational softening.** MW-M31 orbits (analogous to Fig. 1) and distance between the MW and M31 (analogous to Fig. 2) in the four-body MW-M31-M33-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.
<figure><img src="image_18.png" /></figure>

**Extended Data Fig. 5** | Effect of the merger threshold on the MW-M31 distance evolution and merger rate. The distance between the MW and M31 (analogous to Fig. 2) in the four-body MW-M31-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.
<figure><img src="image_19.png" /><figcaption>**Extended Data Fig. 6 | Effect of different schemes of dynamical friction.**</figcaption></figure>

MW-M31 orbits (top two rows, analogous to Fig. 1) and distance between the MW and M31 (bottom row, analogous to Fig. 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 MW-M31 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.
<figure><img src="image_20.png" /></figure>

**Extended Data Fig. 7** | Projected matter density in the present-day initial conditions for the five N-body simulations in the centre-of-mass frame. In each case, the MW is located at the bottom right, with the LMC perceptible as a small

substructure close to the centre of the MW. M31 and M33 are located in the top left. Differences between the five systems are subtle, reflecting the (small) differences in the present positions, and the differences in masses.
<figure><img src="image_21.png" /></figure>

**Extended Data Fig. 8** | **Projected matter density 10 Gyr in the future in the centre-of-mass frame.** The five systems are the same as those shown in Fig. 7. In the three cases where the semi-analytical integration predicted a merger, the MW-M31 haloes have merged into a single halo, with residual

substructures visible in the “Typical merger 2” scenario. In the two cases where the semi-analytical integration predicted no merger, the MW-M31 haloes remain separated by several hundred Mpc.
<figure><img src="image_22.png" /></figure>

**Extended Data Fig. 9** | Comparison between the N-body simulations and the semi-analytical integration. Solid and dashed lines show the distance between the MW and M31 in the N-body simulations and in the semi-analytical numerical integration, respectively, for five example systems. The left panel shows the distance for three cases where the MW and M31 were predicted to merge by the

<figure><img src="image_23.png" /></figure>

semi-analytical model, while the right panel shows two cases where no merger was predicted to occur within 10 Gyr. The dotted lines indicate the 20 kpc merger threshold of our fiducial model. In all five cases, the N-body simulations and the semi-analytical model predict the same outcome, although the merger time of the intermediate cases differ.
<figure><img src="image_24.png" /></figure>

<figure><img src="image_25.png" /></figure>

<figure><img src="image_26.png" /></figure>

<figure><img src="image_27.png" /></figure>

**Extended Data Fig. 10 | 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 four-body system in the fiducial model, corresponding to the top row of Fig. 2. In the bottom row, we add the SMC to

both systems, that is we show the MW-M31-SMC three-body system and the MW-M31-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.