Main

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 galaxies1. In addition, it probably hosts other galaxies yet to be discovered2, and, according to the standard cosmological model, a vast number of completely dark substructures3. The negative radial velocity of M31 towards the MW has been known for over a century4, even before its distance was first accurately measured5. However, although indirect methods have since been used to constrain the transverse components of M31’s velocity vector6,7,8, direct measurements of the minute proper motions have been achieved only much more recently with the Hubble Space Telescope (HST)9. The first numerical studies10 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 merger11,12,13. This scenario has since become the prevalent narrative14,15 and textbook knowledge16,17.

Predicting the future of the LG

The MW and M31 both contain remnants of past mergers and interactions with other galaxies18,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, 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) profile22 with an isotropic velocity distribution, and we used standard assumptions about their concentrations, c, and velocity dispersion profiles23. 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 model24.

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 M200, the mass enclosed within a sphere whose mean density is 200 × 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 1012 M⊙. We adopted M200 = 1 ± 0.2 × 1012 M⊙ (excluding the mass of the LMC, which we treated separately). For all other galaxies, there is considerably more uncertainty. For M31, we adopted M200 = 1.3 ± 0.4 × 1012 M⊙. For M33, we assumed M200 = 3 ± 1 × 1011 M⊙, and for the LMC, we assumed M200 = 1.5 ± 0.5 × 1011 M⊙. See Methods for a review of mass estimates.

The line-of-sight velocities (vlos) are known well, and we adopted the values given by ref. 1. For the distance moduli μ, 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 results27.

We use the notation μδ for the proper motion in declination (dec.) and μα* = μα cos δ 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) astrometry29. However, we found very similar results using the combined Gaia DR2 and HST proper motions of ref. 13, as shown in Extended Data Fig. 1. All assumed values are listed in Table 1. We assumed galactocentric coordinates RAGC = 266.41°, dec.GC = −28.94° (ref. 30), dGC = 8.122 kpc (ref. 31), and a velocity of (12.9, 245.6, 7.78) km s−1 with respect to the Sun30,31,32.

Table 1 Parameters of the fiducial model
Full size table

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 ±2σ. 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 of orbits and, additionally, for the MW-M31-M33 and MW-M31-LMC three-body systems.

Fig. 1: Possible future MW-M31 orbits.
figure 1

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.

Full size image
Fig. 2: Distance between the MW and M31.
figure 2

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 (fM). 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 motions29 of the fiducial model (solid) or HST + Gaia DR2 proper motions13 (dashed). ML, fiducial model.

Full size image

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 galaxy10,33,34. 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.

Fig. 3: Distributions of the merger time tm, the minimum distance min(dM31–MW) and the survival rate of the MW in the MW-M31-M33-LMC system.
figure 3

Blue and red colours show results using a distance threshold (dm) 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, 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.

Full size image

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.

Fig. 4: Effects of M33 and the LMC on the trajectory of the MW.
figure 4

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.

Full size image
Fig. 5: Effects of M33 and the LMC on the trajectory of M31.
figure 5

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

Full size image

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.

Fig. 6: Dependence of the merger probability on observables of the MW-M31-M33-LMC system.
figure 6

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 μδ(M31) and μα*(M31) but also varies significantly with μδ(M33) and μα*(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.

Full size image

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}_\mathrm{t}=7{5}_{-40}^{+65}\) km s−1 (ref. 36), which is compatible with recent measurements13,29 but not with a perfectly radial orbit12. There is, thus, no reason to assume that the transverse velocity measured using Gaia DR3 (most probable value vt = 76 km s−1) is overestimated or to expect 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 (vt = 0, μα* = 38.03 mas yr−1 and μδ = −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 haloes37. The next most massive individual LG galaxies that could impact the MW-M31 orbit are M32 (an M31 satellite and possible merger remnant38) 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 simulations39,40,41 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 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−1 (refs. 30,31,32).

We created Monte Carlo samples for the remaining 20 variables: the four halo masses (MMW, MM31, MM33 and MLMC), the four halo concentration parameters c, the three distance moduli μ, the three sets of proper motions μδ and μα*, and three line-of-sight velocities vlos. 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−1 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 satellites42, and we use the following expression to calculate the acceleration of a satellite due to dynamical friction16 once inside r200 (the radius of a sphere that encloses a mass M200) of the host halo:

$$\frac{{\rm{d}}{\bf{v}}}{{\rm{d}}t}=-\frac{4\uppi {G}^{2}M\rho \ln \Lambda }{{v}^{2}}\left[{\operatorname{erf}}(X)-\frac{2X}{\sqrt{\uppi }}\operatorname{e}^{-{X}^{2}}\right]\frac{{\bf{v}}}{v},$$
(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 simulations43:

$$\epsilon =\begin{cases}2.2{r}_\mathrm{s}-14\,{\rm{kpc}},&\text{if}\;{r}_\mathrm{s}\ge 8\,{\rm{kpc}},\\ 0.45{r}_\mathrm{s},&\text{if}\;{r}_\mathrm{s} < 8\,{\rm{kpc}},\end{cases}$$
(2)

where rs 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{{\rm{d}}{{\bf{v}}}_\mathrm{s}}{{\rm{d}}t}={{\bf{a}}}_\mathrm{DF}\,\frac{{M}_\mathrm{h}}{{M}_\mathrm{s}+{M}_\mathrm{h}},$$
(3)
$$\frac{{\rm{d}}{{\bf{v}}}_\mathrm{h}}{{\rm{d}}t}=-{{\bf{a}}}_\mathrm{DF}\,\frac{{M}_\mathrm{s}}{{M}_\mathrm{s}+{M}_\mathrm{h}},$$
(4)

where vs and Ms are the velocity and mass of the satellite, vh and Mh are the velocity and mass of the host, and aDF 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 equations44. 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 apocentre. 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 later46.

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 data28 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 r200, with a sigmoid function. We used Monte Carlo sampling for the distribution function \({f}_{i}({\mathcal{E}})\) for the NFW density profile ρNFW with discrete particles using the Eddington formula:

$${f}_{i}({\mathcal{E}})=\frac{1}{2\sqrt{2}{\uppi}^{2}}\int_{{\varPhi}_{{\rm{T}}} = 0}^{{\varPhi}_{{\rm{T}}} = {\mathcal{E}}}\frac{{{\rm{d}}}^{2}{\rho}_{{\rm{NFW}},i}}{{\rm{d}}{\varPhi}_{{\rm{T}}}^{2}}\frac{{\rm{d}}{\varPhi}_{{\rm{T}}}}{\sqrt{{\mathcal{E}}-{\varPhi}_{{\rm{T}}}}},$$
(5)

where \({\mathcal{E}}\) is the relative energy and ΦT is the total gravitational potential.

The simulations were performed with the Gadget-4 code47 using a particle mass resolution of 106.5 M⊙ (resulting in ~106 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 works23,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.1{7}_{-0.15}^{+0.21}\times 1{0}^{12}\,{{{M}}}_{\odot }\) (ref. 48), \(1.5{1}_{-0.40}^{0.45}\times 1{0}^{12}\,{{{M}}}_{\odot }\) (ref. 49), \(1.2{3}_{-0.18}^{+0.21}\times 1{0}^{12}\,{{{M}}}_{\odot }\) (ref. 50) and \(1.{1}_{-0.1}^{+0.1}\times 1{0}^{12}\,{{{M}}}_{\odot }\) (ref. 51) (where the latter two works also used a simulation-based prior).

Estimates made using rotation curves include \(1.0{8}_{-0.14}^{+0.20}\times 1{0}^{12}\,{{{M}}}_{\odot }\) based on Gaia DR2 (ref. 52), \(0.8{9}_{-0.08}^{+0.1}\times 1{0}^{12}\,{{{M}}}_{\odot }\) using stars in the galkin catalogue53, 0.822 ± 0.052 × 1012 M⊙ using classical cepheids54 and \(1.0{8}_{-0.11}^{+0.12}\times 1{0}^{12}\,{{{M}}}_{\odot }\) from the H3 survey and Gaia DR3 (ref. 55). Other recent measurements include \(1.5{4}_{-0.44}^{+0.75}\times 1{0}^{12}\,{{{M}}}_{\odot }\) using combined Gaia DR2 and HST kinematics of globular clusters56, 1.16 ± 0.24 × 1012 M⊙ (including the mass of the LMC) using the kinematics of halo stars, \(1.2{6}_{-0.22}^{+0.40}\times 1{0}^{12}\,{{{M}}}_{\odot }\) from high-velocity RR Lyrae stars57 and \(1.1{9}_{-0.32}^{+0.49}\times 1{0}^{12}\,{{{M}}}_{\odot }\) from a Bayesian estimate using dwarf galaxy kinematics from several sources58. 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 1012 M⊙. Both the simple mean and median of the above measurements are 1.16 × 1012 M⊙, 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, rvir, or r200), which is not always made explicit. When the LMC is excluded, the mass of the MW is reduced by ~0.15 × 1012 M⊙ (see our discussion on the LMC mass below). We adopted MMW = 1.0 ± 0.2 × 1012 M⊙ 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.1{4}_{-0.35}^{+0.51}\times 1{0}^{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 ± 0.4 × 1012 M⊙ from satellite kinematics, ref. 62 derived \(1.{2}_{-0.7}^{+0.9}\times 1{0}^{12}\,{{{M}}}_{\odot }\) from the kinematics of M31 dwarf spheroids and ref. 63 measured a total mass of \(1.0{5}_{-0.15}^{+0.15}\times 1{0}^{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 × 1012 M⊙ from the H i rotation curve, and ref. 65 measured \(2.{0}_{-0.3}^{+0.4}\times 1{0}^{12}\,{{{M}}}_{\odot }\) from the kinematics of the Giant Southern Stream. Reference 66 measured \(1.3{5}_{-0.15}^{+0.15}\times 1{0}^{12}\,{{{M}}}_{\odot }\), and ref. 67 found \(1.{4}_{-0.2}^{+0.2}\times 1{0}^{12}\,{{{M}}}_{\odot }\), both using outer halo globular clusters, whereas ref. 68 found 0.8 ± 0.1 × 1012 M⊙ from high-velocity planetary nebulae, and ref. 69 found 1.39 ± 0.26 × 1012 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 × 1012 M⊙. Our adopted mass of MM31 = 1.3 ± 0.4 × 1012 M⊙ 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 × 1010 M⊙, extrapolated out to a virial mass of 5 × 1011 M⊙ 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 ± 0.55 × 1011 M⊙ (refs. 73,74). Citing both the direct measurements and abundance-matching, ref. 74 adopted a mass range of 0.8–3.2 × 1011 M⊙ in their dynamical models of the M33-M31 interaction, whereas refs. 13,75 assumed a total mass of 2.5 × 1011 M⊙. We adopted MM33 = 3 ± 1 × 1011 M⊙ 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.3{8}_{-0.24}^{+0.37}\times 1{0}^{11}\,{{{M}}}_{\odot }\) (ref. 76), 1.30 ± 0.3 × 1011 M⊙ (ref. 77), \(1.8{8}_{-0.35}^{+0.4}\times 1{0}^{11}\,{{{M}}}_{\odot }\) (ref. 78) and \(1.2{9}_{-0.23}^{+0.28}\times 1{0}^{11}\,{{{M}}}_{\odot }\). Reference 79 obtained good agreement between the perturbations of MW halo stars with an LMC mass of 1.5 × 1011 M⊙. Using the abundance of probable satellites of the LMC, ref. 80 obtained a lower limit of 1.24 × 1011 M⊙, and using the kinematics of satellites associated with the LMC, ref. 81 found \(1.6{5}_{-0.49}^{+0.47}\times 1{0}^{11}\,{{{M}}}_{\odot }\). Most recently, ref. 82 used 30 LMC globular clusters to infer a total mass of \(1.8{0}_{-0.54}^{+1.05}\times 1{0}^{11}\,{{{M}}}_{\odot }\), whereas ref. 35 found an even higher value of \(2.{5}_{-0.8}^{+0.9}\times 1{0}^{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 × 1011 M⊙, and the median is 1.5 × 1011 M⊙. 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 × 1011 M⊙, which matches the choice MLMC = 1.5 ± 0.5 × 1011 M⊙ in our fiducial model.

Other mass

Cosmological simulations show that the total mass of LG analogues is \({\sim} 2{5}_{-8}^{+9}\%\) larger than the sum of the M200 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 haloes22 with the total masses M200 defined above. We assumed concentration parameters of 10 ± 2 for all galaxies, consistent with the results of cosmological simulations84,85,86,87,88 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 range84. The concentration parameter does not affect the orbital calculation, except for the dynamical friction, where (for a given mass), it sets the scale radius, rs = r200/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.