The Milky Way and M31 orbital history: did the Local Group evolve in isolation?

ABSTRACT

We use new measurements of the M31 proper motion to examine the Milky Way (MW) – M31 orbit and angular momentum. For Local Group (LG) mass consistent with measured values, and assuming the system evolves in isolation, we show a wide range of orbits is possible. We compare to a sample of LG-like systems in the Illustris simulation, and find that ∼13 per cent of these pairs have undergone a pericentric passage. Using the simulated sample, we examine how accurately an isolated, two-body model describes the MW–M31 orbit, and show that ∼10 per cent of the analogues in the simulation are well-modelled by such an orbit. Systems that evolve in isolation by this definition are found to have a lower rate of major mergers, and in particular have no major mergers since |$z \approx 0.3$|⁠. For all systems, we find an increase in the orbital angular momentum which is fairly independent of the merger rate, and is possibly explained by the influence of tidal torques on the LG. Given the likely quiet recent major merger history of the MW, it is plausible that the isolated two-body model appropriately describes the orbit, though recent evidence for a major merger in M31 may complicate this interpretation.

1 INTRODUCTION

The Milky Way (MW) and M31 are the dominant constituents of the Local Group (LG). The sum of the masses of the MW and M31 dark matter haloes appears to be consistent with the mass of LG as a whole (Lemos et al. 2021; Benisty et al. 2022; Sawala, Teeriaho & Johansson 2023a), though systematics such as the non-equilibrium nature of the MW and M31 may affect the interpretation of the mass estimators (Patel & Mandel 2023); see Wang et al. (2020) and Bhattacharya (2023) for recent reviews. A possible difference between the sum of the MW and M31 masses and the estimated LG mass as a whole may hint at additional, dark mass contained within the LG.

An important, related property associated with the LG, MW, and M31 is the nature of its orbit. While it has long been assumed that the MW and M31 are on a radial-approaching orbit (Kahn & Woltjer 1959; Gunn 1974; Lynden-Bell 1981), recent measurements of the 3D kinematics of M31 have challenged this standard assumption. In particular, proper motion measurements from Gaia DR2, EDR3, and the Hubble Space Telescope (HST) suggest a non-zero value for the MW–M31 relative tangential velocity (van der Marel et al. 2012, 2019; Salomon et al. 2021; Rusterucci et al. 2024). Though it still likely sub-dominant relative to the radial component, any tangential component of motion must be taken into account for mass estimates of the LG, and to understand the formation of the LG.

In addition to these kinematic measurements, there has been a new understanding that galaxies within the LG, such as the Large Magellanic Cloud (LMC) or possibly even M33, have a significant effect on the kinematics of LG galaxies. In particular, the LMC may induce a relative motion between the MW disc and the extended MW stellar and dark matter halo, resulting in a displacement in the velocities between the disc and the halo components (Gómez et al. 2015; Garavito-Camargo et al. 2019; Conroy et al. 2021; Erkal et al. 2021). Petersen & Peñarrubia (2021) use halo stars to measure this effect, and the resulting velocity shift adjusts the relative velocity between the MW and M31. This is in addition to the kinematic effect that changes the location of the MW–M31 barycenter due to a massive LMC (Peñarrubia et al. 2016).

Though it is a fundamental means to characterize the LG, and dictates the future fate of the MW and M31 (Cox & Loeb 2008; Sawala et al. 2024), the above emphasizes how little is still known on the true nature of the MW–M31 orbit. The simplest model for the MW–M31 orbit appeals to two-body kinematics combined with cosmology. This model originally motivated Kahn & Woltjer (1959) to examine the Keplerian equation for a MW–M31 two-body orbit and use this to estimate the mass of the LG. The simple input assumptions were that the LG mass is dominated by the MW and M31, and that they can be modelled as spatially coincident near the time of the Big Bang. While the original analysis assumed a radial orbit, more generally the model may include tangential motion. Applying these assumptions along with the present day measured values for the relative separation and velocity, they reported a lower limit on the LG mass of |$\gtrsim 10^{12} \mathrm{{\rm M}_{\odot }}$|⁠. This so-called timing model was calibrated to cosmological simulations which identify LG-like systems (Kroeker & Carlberg 1991; Li & White 2008). These analyses found that generally the timing mass estimator is reasonable, though biases may be incurred when accounting for a large tangential velocity and the cosmological constant (González, Kravtsov & Gnedin 2014; Hartl & Strigari 2022), and the aforementioned affect of the LMC (Benisty et al. 2022). An appropriate definition of the halo mass may help to mitigate these biases (Kroeker & Carlberg 1991; Sawala et al. 2023b).

While the assumption of a Keplerian two-body orbit in which the LG evolves in isolation provides a reasonable starting point, it has yet to be established that it accurately describes the LG. For example major mergers, either within the MW or M31, may affect the LG kinematics. Though it appears the MW has had a relatively quiet recent major merger history, in which the last significant merger occurred |$\gtrsim 8$| Gyr ago (Helmi et al. 2018; Belokurov et al. 2020), there is evidence for a more recent merger in M31 (D’Souza & Bell 2018; Hammer et al. 2018). The merger histories of LG-like systems may also be extracted from cosmological simulations (Carlesi et al. 2020). Aside from mergers, tidal interactions between the LG and the local environment may influence the LG, for example via the influence of external torques that may alter the angular momentum. Though there has been progress in characterizing the matter distribution in the local volume of the universe (Hoffman et al. 2012), its gravitational influence on the LG is still largely uncertain (McCall 2014).

In this paper, we undertake a critical examination of the MW–M31 orbit, in particular the assumption of isolation, the impact of mergers, and the angular momentum evolution of the orbit. Using the measured M31 kinematics and the estimates of the LG mass, we determine the possible orbits for the MW and M31 using the simple assumption of an isolated two-body orbit. Then, we determine how well this model describes the orbital history between the MW and M31 by comparing it to analogues in cosmological simulations, and we define a criteria for systems that are well-described by the two-body model.

For the LG-analogue systems identified in the simulations, we determine the merger rates, and examine how this merger rate affects the orbital history of the LG. We additionally extract the angular momentum for our simulated sample, and compare to the observed orbital angular momentum of the LG. For the simulated sample, we examine the evolution of the angular momentum over cosmic time.

This paper is organized as follows. In Section 2, we examine the possible orbits given the LG mass and kinematics. In Section 3, we discuss the simulation sample used in our analysis, and the calculation of the merger rate in the simulations. In Section 4, we present the results of our analysis, and then in the final section we move on to the conclusions and discussion.

2 LOCAL GROUP KINEMATIC PROPERTIES AND ORBITAL CHARACTERISTICS

In this section, we review the kinematics properties of the LG, and use these to estimate the LG orbital angular momentum. We then move on to calculate possible LG orbits given the kinematic constraints.

2.1 LG mass and angular momentum

The distance from the Sun to M31 is |$770 \pm 40$| kpc, and the heliocentric velocity is |$-301$| km s⁻¹ (van der Marel & Guhathakurta 2008). For best-fitting values of the MW circular and Local Standard Rest velocities (McMillan 2011), this corresponds to a MW-centric radial velocity of |$v_{rad} = -110$| km s⁻¹. The results from Gaia EDR3 find that the tangential velocity is |$v_{tan} = 82.4 \pm 31.2 \mathrm{\, km \, s^{-1}}$| (Salomon et al. 2021), which significantly rules out a purely radial orbit and is larger than previous measurements from HST (Marel et al. 2012). Earlier results from Gaia DR2 are consistent with an even larger tangential velocity, |$v_{tan} = 170 \pm 51 \mathrm{\, km \, s^{-1}}$| (der Marel et al. 2019). Systematics on these results have been considered in Rusterucci et al. (2024). Large tangential velocities may also be motivated by indirect measurements of the M31–MW relative motion (Diaz et al. 2014; Salomon et al. 2016).

The deduced radial and tangential velocities may be corrected to account for the impact of the LMC, which may induce a significant relative motion between the MW disc and the frame defined by the extended halo (Gómez et al. 2015; Garavito-Camargo et al. 2019; Conroy et al. 2021; Erkal et al. 2021). Using simulations of LG-like systems, Benisty et al. (2022) found that correcting for the LMC may increase the tangential velocity by 25–30 km s⁻¹, increase the radial velocity by up to 50 |$\mathrm{\, km \, s^{-1}}$|⁠, and increase the relative separation up to 40 kpc. This is consistent with the results from Chamberlain et al. (2023), who use the measured ‘travel velocity’ of the MW disc with respect to the outer halo to find an upwards radial velocity correction of approximately 20 |$\mathrm{\, km \, s^{-1}}$| and a corresponding increase in the tangential velocity by a few |$\mathrm{\, km \, s^{-1}}$|⁠.

We consider a range of mass estimates for the LG, motivated by the results from several independent measurements. Fitting to the measurements of the Local Hubble Flow, Peñarrubia et al. (2014) obtain a LG mass of |$(2.3 \pm 0.7) \times 10^{12}$| M|$_\odot$| [which Peñarrubia & Fattahi (2017) find is larger than the sum of the MW and M31 virial masses]. Using a numerical least action method and fitting to orbits of LG galaxies, Phelps, Nusser & Desjacques (2013) obtain a LG mass of |$(2 \pm 1) \times 10^{12}$|  M|$_\odot$|⁠. Diaz et al. (2014) apply the virial theorem in combination with momentum conservation and obtain a LG mass of |$(2.5 \pm 0.4) \times 10^{12} \, {\rm M}_\odot$|⁠. These results may be compared to a pure application of the timing argument combined with the measured tangential velocity from Gaia EDR3, which gives an estimated LG mass of |$(5.8 \pm 1.0) \times 10^{12}$| M|$_\odot$|⁠. These results are broadly consistent with measured values for the sum of the masses of the MW and M31 (Sawala et al. 2023a), and may be also compared to results based on cosmological simulations (González et al. 2014; McLeod & Lahav 2020; Lemos et al. 2021; Carlesi, Hoffman & Libeskind 2022). Taking into account dynamical corrections due to the LMC, Benisty et al. (2022) find a total LG mass of |$(3.7 \pm 0.5) \times 10^{12}$| M|$_\odot$|⁠. Changing the embedding from |$\lambda$| to |$\ddot{a}/a$|⁠, Benisty (2024) found the lowest timing argument mass estimate of |$~ 2 \times 10^{12}$| M|$_{\odot }$|⁠.

The orbital angular momentum per reduced mass for the MW/M31 system is |$\vec{\boldsymbol{L}}_{\mathrm{ orb}} = \vec{\boldsymbol{r}} \times \vec{\boldsymbol{v}}$|⁠, where |$\vec{\boldsymbol{r}} = \vec{\boldsymbol{r}}_{\mathrm{ MW}}-\vec{\boldsymbol{r}}_{\mathrm{ M}31}$| represents the separation vector, and |$\vec{\boldsymbol{v}} = \vec{\boldsymbol{v}}_{\mathrm{ MW}}-\vec{\boldsymbol{v}}_{M31}$| is the relative velocity. The magnitude of the orbital angular momentum is simply given by |$|\vec{\boldsymbol{L}}_{\mathrm{ orb}}| = rv_{tan}$|⁠. Using the EDR3 measurement for the tangential velocity, the blue curve in Fig. 1 depicts the derived probability distribution for the magnitude of the orbital angular momentum. Here, the uncertainties on the proper motion are propagated to the orbital angular momentum, and the mean and the 90 per cent confidence interval of the distribution is given by

2.2 Analytic orbits

From the measured LG mass and angular momentum, we are in position to examine possible orbital trajectories. We consider a simplified analysis assuming that the MW and M31 evolve in isolation, and are unaffected by mergers or interactions with nearby structures. This serves as a simplified model which we can compare to the results from simulations in the following sections.

Defining the total mass as |$M_{\mathrm{ LG}}=M_{\mathrm{ MW}}+ M_{\mathrm{ M}31}$|⁠, along with a polar coordinate system in which the separation r defines the orbital plane, the equation of motion is

where |$H_0 = 70$| km s⁻¹ Mpc⁻¹ is the Hubble constant, and |$\Omega _{\Lambda } = 0.7$|⁠. Note that here we interpret |$M_{\mathrm{ LG}}$| as the sum of the virial masses of the MW and M31, and as discussed below there is some ambiguity as to its precise definition. For example Sawala et al. (2023b) advocate for a definition of |$M_{\mathrm{ LG}}$| as the mass contained within a sphere with radius defined by the present day separation between the two galaxies. We use the range of measured angular momentum and LG mass and solve equation (2) and plot resulting example two-body orbits in Fig. 2. In all cases, we wind the orbit back to |$z=20$|⁠, which matches the earliest snapshot of the simulation sample we compare to in Section 3. The first row uses the canonical case in which the relative radial velocity is |$v_{rad} = -110 \mathrm{\, km \, s^{-1}}$|⁠. In this case, for the lower mass assumptions, the physical separation increases back to high redshift, implying that the galaxies are currently on their first approach. In this case, only for masses on the larger end of the allowed range do the two galaxies have a well-defined apocenter within the redshift range considered. For high mass and low tangential velocity, we see that a pericentric passage is even plausible at early times.

The bottom row of Fig. 2 show cases for which the radial velocity is increased (so the modulus of the velocity is decreased), which is motivated by the possible shift due to the LMC. Increasing the radial velocity to as large as |$-60 \mathrm{\, km \, s^{-1}}$|⁠, plotted in the bottom row, for over half of the allowed mass range there is a pericentric passage for the range of the tangential velocities. Interestingly, even the median estimated mass shows a pericentric passage for the low tangential velocity system. This suggests that if the LMC significantly increases the relative radial velocity between the MW and M31, the pair is more likely to have undergone a pericentric passage.

The above analysis is interesting when considered in the context of the classical timing model for the formation of the LG. This model relies on two assumptions: (1) |$r \longrightarrow 0$| as |$t \longrightarrow 0$|⁠, and (2) the LG is on its first infall. For the given set of present day kinematic measurements, the LG mass is then adjusted to satisfy the above boundary condition in the equation of motion. Focusing on the (red) set of curves in the first row of Fig. 2 with the lowest tangential velocity of |$\sim 51$| km s⁻¹, which shows a pericentric passage in the rightmost panel, if we were to fit |$M_{\mathrm{ LG}}$| to satisfy the timing mass conditions, we would recover a LG mass that is less than the true mass, as shown in the panel in the second column. The timing mass estimate would then systematically underestimate the true mass of the LG. The opposite extreme is illustrated for the case of the same tangential velocity and a lower LG mass, in which the two galaxies are on first infall, but had a larger separation at early times. If we fit the mass by enforcing the timing conditions, we would have estimated a larger LG mass, for example as given in the panel in the second column. The timing mass estimate would then overestimate the true mass of the LG. Note that this was the case found in Hartl & Strigari (2022) for the simulated sample discussed in the section below. Thus, we can conclude that it is likely the case that the two galaxies are on first infall and evolved from physical separations comparable to or larger than their present day separation.

Regardless if the pair is on first infall or not, when we wind the orbits back we notice two populations of MW/M31 pairs: (1) those which evolved from separations smaller than their current separation, and (2) those which evolved from separations larger than their current separation. The first population reflects the canonical timing argument model, while the second population breaks with this traditional |$r \longrightarrow 0$| as |$t \longrightarrow 0$| assumption, in particular for systems with higher tangential velocity. Below we study this classification in more detail in the context of simulations.

Fig. 3 illustrates the general dynamics of two-body motion in terms of tangential and radial velocities. We consider systems with a current separation of 770 kpc and masses of |$2 \times 10^{12} {\rm M}_{\odot }$| (top row) and |$4 \times 10^{12} {\rm M}_{\odot }$| (bottom row), and compute orbits for tangential velocities ranging from 0 to 130 |$\mathrm{km \, s^{-1}}$| and radial velocities from −150 to 0 |$\mathrm{km \, s^{-1}}$|⁠. The left plots are colour-coded to display the maximum separation calculated for each scenario. The analysis indicates that the maximum separation is primarily influenced by the radial velocity. However, as the radial velocity decreases, the influence of tangential velocity on maximum separation becomes more pronounced. The right plots show a subset of the left plots, but focus on orbits that experience a pericentric passage, with colour-coding representing the distance of these passages. Notably, the distance of the pericentric passage is driven by the tangential velocity and appears independent of radial velocity, where a larger tangential velocity corresponds to a larger distance for these passages. A comparison between the top and bottom rows reveals that for a higher present day LG mass, a broader range of radial and tangential velocities corresponds to systems that exhibit pericentric passages.

To summarize the above analysis, due to the uncertainty in the LG kinematics, we deduce a wide range of possible orbits. As a next step, we will explore how well equation (2) predicts the orbits when known kinematic parameters are provided. We aim to answer these questions by analysing a large sample of MW/M31 analogue pairs in IllustrisTNG.

3 LG-ANALOGUES IN SIMULATIONS

In this section, we describe the simulation sample used for our analysis, followed by a brief description of how their merger histories are calculated.

3.1 Simulation

We identify a sample of LG-like systems in the Illustris TNG300-1 simulation (Marinacci et al. 2018; Naiman et al. 2018; Nelson et al. 2018, 2019; Pillepich et al. 2018; Springel et al. 2018). The TNG300-1 simulation has a comoving volume of |$302.6 \, {\rm Mpc}$| and |$N=2500^{3}$| particles. This simulation includes baryonic physics with a particle mass of |$m_{\mathrm{ baryon}}=7.6 \times 10^{6}$| M|$_\odot h^{-1}$| and a dark matter particle mass of |$m_{\mathrm{ dm}}=4.0 \times 10^{7}$| M|$_\odot h^{-1}$|⁠. The simulation assumes cosmological parameters measured by Planck (Planck Collaboration I 2016).

The Illustris data base stores galaxy properties at 100 snapshots. We utilize the position of each galaxy, defined as position of the particle with the lowest gravitational potential, and the velocity of each galaxy which is calculated by summing over all particles and weighting by the mass. The cuts imposed to obtain the sample of 597 LG-analogues, described in detail in Hartl & Strigari (2022), were chosen to match the observed kinematics and stellar mass of the MW and M31. Briefly, we consider MW/M31 type galaxies as those with a B-band magnitude within |$-22.3 \lt M_{B} \lt -19.3$|⁠, and identify approaching (⁠|$V_r \lt 0 \, \mathrm{km \, s^{-1}}$|⁠) MW–M31 analogue pairs with a separation 500 kpc |$\lt |\vec{r}| \lt 1000$| kpc, motivated by the measured distance between the MW and M31. For the full motivation for the cuts used to obtain the sample, we refer to Hartl & Strigari (2022).

For the halo masses extracted from the simulations, we use the subfind masses as provided in the Illustris data base. This tends to be on average larger than that |$M_{200}$| masses, defined as the mass within a sphere with a radius, |$r_{200}$| that encompasses 200 times the critical density (Li & White 2008). This is because the subfind algorithm includes particles outside of |$r_{200}$|⁠, though it excludes particles within subhaloes of the main halo.

3.2 Mergers

We calculate the merger history of each galaxy in our sample with the merger catalogues constructed by the sublink algorithm, see Rodriguez-Gomez et al. (2015) for a detailed description of this process. Note, these authors use the merger trees following the stellar particles, and we use the merger trees following the dark matter particles. Additionally, they provide the details for calculating the merger rate. In determining the merger history, we classify each merger as either a major or a minor merger. The type of merger is determined by the merger mass ratio, which is specifically the ratio of the mass between the main halo and the halo which merged into the main halo. Major mergers are classified by a merger mass ratio |$\mu$| of |$0.3 \lt \mu \le 1.0$|⁠, and minor mergers with |$0.05 \lt \mu \le 0.3$|⁠. Note that, as discussed in more detail below, the LMC impact into the MW is classified as a minor merger.

Briefly, the merger rate, |$\mathrm{ d}N_m/\mathrm{ d}z$|⁠, is calculated in a few straightforward steps, where |$N_{m}$| is the number of mergers and |$\mathrm{ d}z$| is the redshift bin width. The steps are as follows:

• Determine the range of merger mass ratios: |$0.05\lt \mu _{\mathrm{ DM}} \le 0.3$| (minor), |$0.3\lt \mu _{\mathrm{ DM}} \le 1.0$| (major).

• Define redshift bins.

• Count the number of mergers in each redshift bin.

• Divide by the redshift bin width.

• Divide by the number of galaxies in the sample at |$z=0 (1194)$|⁠.

Fakhouri, Ma & Boylan-Kolchin (2010) have explored the merger rate of all galaxies in the Millennium simulation for a halo of mass M, and reported a fitting function for the merger rate, given by

The best fit parameters were found to be |$(\alpha _{1}, \beta _{1}, \gamma ,\eta) = (0.133, -1.995, 0.263, 0.0993)$| and |$(A, \bar{\mu }) = 0.0104, 9.72 \times 10^{-3})$|⁠. Note that the halo mass, M, drives the behavior of this fitting function and depends on z. The redshift dependence of the halo mass is given by (Carlesi et al. 2020),

with the best fit parameters of |$\alpha$| and |$\beta$| obtained for MW/M31 like galaxies. For our results, we use the M31(CS) value of |$\alpha = 6.01$| and |$\beta = 2.66$| from Carlesi et al. (2020).

We compare the merger rate calculated from our sample of analogue pairs to the prediction from the above fitting function. In performing this comparison, we can determine if MW/M31 type galaxies in pairs have a similar merger history to what is predicted for a typical MW/M31 mass galaxy.

In Fig. 4, we plot the cumulative form of the differential merger rate |$\mathrm{ d}N_{m}/\mathrm{ d}z$|⁠, which represents the total number of mergers per halo from redshift |$z=0$| out to redshift z. Fig. 4 shows that the majority of MW and M31-like galaxies in LG pairs have each had more than 10 minor mergers and a few major mergers since |$z \sim 1$|⁠. For comparison, we plot the Fakhouri et al. (2010) fitting function using the median, 5 per cent, and 95 per cent points on the total mass distribution. Focussing on the minor merger history for MW/M31 in LGs, the blue curve in Fig. 4 suggests our sample is in agreement with the Fakhouri et al. (2010) prediction. Note that, while the normalization of the major mergers is below the prediction, this is most likely explained by the scatter in the major merger rate between different simulations (Hopkins et al. 2010). Importantly, the shape of the major merger is very similar to that of the Fakhouri et al. (2010) prediction.

4 RESULTS

We now move on to examine the orbital histories, angular momentum, and merger histories for our LG analogues. We determine how well the orbital histories match the simple, isolated two-body model discussed above, and examine how the properties of the analogues depend on their orbital histories.

4.1 MW-M31 orbital history in Illustris

We calculate the orbital history of the MW–M31 analogue pairs identified from the simulation by calculating the relative separation between the MW and M31 at each snapshot, and examining the separation as a function of time and redshift. In order to highlight the diversity of orbital histories, the top panel in Fig. 5 shows orbits for two example LG analogues. The orbital history plotted for the analogue pair on the left is an example of an orbit which is perfectly described in the timing model, in which |$r \rightarrow 0$| as |$t \rightarrow 0$|⁠, and the pairs are approaching the pericenter for the first time today. This system has a total mass of |$5.3 \times 10^{12} \, \mathrm{{\rm M}_{\odot }}$|⁠, a separation of 970 kpc and the radial and tangential velocities are |$-65 \, \mathrm{ km \, s^{-1}}$| and |$37 \, \mathrm{ km \, s^{-1}}$|⁠, respectively. For comparison, the right plot shows a pair that has undergone a pericentric passage and is now approaching its pericenter for the second time today. This pair that has undergone a percentric passage has a total mass of |$5.6 \times 10^{12} \, \mathrm{{\rm M}_{\odot }}$|⁠, a separation of 651 kpc and the respective radial and tangential velocities are calculated to be |$-40 \, \mathrm{ km \, s^{-1}}$| and |$144 \, \mathrm{ km \, s^{-1}}$|⁠.

Fig. 5 indicates that the individual LG-analogue pairs can have vastly differing orbital histories. To get a sense of the average orbital history for our entire sample, we plot the cumulative distribution of the relative separation for all our pairs in Fig. 6 (top). This shows that the mean orbit reached an apocenter several Gyr ago, and is approaching the apocenter for the first time. This result was also found in Sawala et al. (2023b), and shows that even though an individual orbit may be significantly different, the mean orbit is well-described by the classical timing orbit model.

To further highlight the range of orbital histories for our entire sample, we note that 77 have undergone a pericentric passage, and are approaching their pericenters for the second time. There are no systems that have had two pericenter passages. Systems that have undergone pericentric passages have a wide range of masses (⁠|$\mathrm{\sim 0.8 \times 10^{12} {\rm M}_{\odot }-1.0 \times 10^{13} {\rm M}_{\odot }}$|⁠), though they tend to have lower modulus for their radial velocities, |$\gtrapprox -60 \, \mathrm{km \, s^{-1}}$|⁠. This, perhaps surprisingly, suggests that it is still possible for a MW/M31 pair to have already had a previous passage. Further, through a visual inspection and anticipating the analysis below, we define ‘isolated’ pairs as those with a smooth trajectory. We find 15 of 77 pairs with pericentric passages satisfy this smooth criteria, with pericenters ranging from 60-700kpc.

4.2 Orbital angular momentum

We now move on to consider the orbital angular momentum. We determine the orbital angular momentum from the separation and tangential velocity of the LG-analogues. Fig. 1 shows the probability density of the resulting magnitude of the orbital angular momentum calculated for all 597 pairs (red curve) at |$z=0$|⁠. We compare this to the blue curve which shows the uncertainty in the angular momentum of the LG calculated from the EDR3 proper motions, and find these two calculations are in good agreement. So, in this sense the observed orbital angular momentum is consistent with our LG-analogue sample.

We then calculate the orbital angular momentum for all snapshots to study how this quantity evolves for LG-like systems. The bottom panel in Fig. 5 shows the example orbital angular momentum growth over cosmic history for the two analogues shown previously. We see both of these pairs, regardless of whether they are on first or second infall, show a steady increase in angular momentum up to the present day snapshot. To see if this angular momentum growth is a general trend across all LG pairs, we plot the cumulative orbital angular momentum at each snapshot in Fig. 6 (bottom). It is interesting to note that, when examining the 95  per cent containment for the orbital angular momentum, the width of this region is fairly constant across time, indicating that this trend is similar for most LG analogues.

4.2.1 Origin of |$\vec{L}_{\mathrm{ orb}}$|

Referring back to Fig. 2, haloes are characterized by a wide range of physical separations, up to as large as |$\sim 2400$| kpc. Pairs with a larger magnitude of approach velocity, i.e. more negative values of radial velocity, tend to have larger separations at identification. However, a more significant trend is noted when considering the tangential velocity, which varies with separation at identification in all cases. We examine this trend further with our simulation sample, where we define |$r_{i}$| as the separation at the earliest snapshot at which both of the main haloes are identified in the system. This can be roughly thought of as the distance the two galaxies were when the LG ‘formed’ as a system, and is related to the size of the overdensity that collapsed to form the LG (Boylan-Kolchin et al. 2016).

We split the sample of LG-analogues into two, those with physical separations |$r_i\lt 400$| kpc and those with |$r_i \gt 400$| kpc, resulting in 135 and 434 pairs, respectively. We found no preference for pairs to be in either of the initial separation bins based on whether they were on their first infall or second infall. However, Fig. 7 plots the orbital angular momentum of these two initial separation samples, where we see a very notable effect that the sample with larger such formation distances (red curve) have larger orbital angular momentum than those pairs that formed more in closer proximity (blue curve). Thus, the present day tangential velocity appears to be a reflection of the earliest time in the formation of the LG system.

To investigate this issue further, we calculate the angular momentum at the earliest snapshot the pair is identified and plot it against the separation at this snapshot. As seen in Fig. 8 (left), the initial orbital angular momentum increases with larger formation distances. The curves on the plot are constructed by binning in formation distance and computing the median, 68 per cent and 95 per cent containment regions shown as the solid black, dashed-black, and dashed-grey lines, respectively. We also examine the orbital angular momentum at |$z=0$|⁠, and show in Fig. 8 (right) that this relationship to the initial formation distance is preserved. Note that we investigated the present day orbital angular momentum against the present day separation, and we found no such relationship, which suggests that the tangential velocity drives the relationship shown in Fig. 8.

Finally, we determine the maximum separation of the pair at any point in the simulation, |$D_{\mathrm{ max}}$|⁠, and plot this value against the present day orbital angular momentum in Fig. 9. Fig. 9 shows there is also a correlation between the present day orbital angular momentum and the maximum separation. For pairs that follow the canonical timing argument criteria, their maximum separation is the distance of the pair at apocenter: |$D_{\mathrm{ max}}$| = |$r_{\mathrm{ apo}}$|⁠. Addditionaly, note that since the MW and M31 are approaching today, their maximum separation must be less than their current separation: |$D_{\mathrm{ max}} \lt r_{0}$|⁠. According to Fig. 9, apocenters of about 950 kpc correspond to the most likely value of the present day angular momentum |$({\rm in} \, {\rm km} \, {\rm s}^{-1} \, {\rm Mpc})$| of 64.4 with the respective upper and lower bounds of 113.8 and 26.3. We can divide by the present day separation of 770 kpc to get the median tangential velocity of 83.7 |$\mathrm{km \, s^{-1}}$|⁠, with 147.8 |$\mathrm{km \, s^{-1}}$| and 34.2 |$\mathrm{km \, s^{-1}}$| as the respective upper and lower bounds. It is interesting to note that the median estimated value of the tangential velocity using this method is |$\lt 2 \mathrm{\, km \, s^{-1}}$| from the most recent measured value.

4.3 Orbital matching

With an understanding of typical orbits and angular momenta for our LG analogues, we now move on to examining how well a simple, isolated two-body model describes the orbits. We start from equation (2), and compare the implied trajectory to the true orbit from the simulation, and thereby determine how well equation (2) predicts the orbital history of a typical MW/M31 pair.

Starting from the present day snapshot, we compare the predicted separation from the analytic orbit, |$|\vec{r}_a|$|⁠, to the true separation in the simulation |$|\vec{r}_s|$| by defining the quantity |$C = (|\vec{r}_a|-|\vec{r}_s|/)|\vec{r}_s|$|⁠. We track each of the orbits from the present day back to when the universe was 3.2 Gyr old (⁠|$z=2$|⁠). Rather than tracing back to the earliest snapshot, |$z\approx 20$|⁠, we choose 3.2 Gyr in order to remove noise that is introduced due to different halo identification times. To identify matched orbits we consider a strict criteria for |$|C| \lt 0.1$| yielding 51 pairs. If we instead loosen this criteria to |$|C| \lt 0.2$| the sample increases to 169 pairs. However, for this paper we will focus our analysis on the stricter sample of 51 pairs.

We briefly investigate the significance of incorporating the |$\Lambda$| term in our orbital predictions. If we instead set |$\Lambda = 0$|⁠, our analysis reveals only one instance where |$|C| \lt 0.1$|⁠. This finding stresses the critical necessity of including the cosmological constant in our models to accurately reconstruct and predict the orbital history of the MW and M31.

Fig. 10 compares the angular momentum and the merger rate for the sample of 51 pairs, which closely matched the predicted analytical orbit, to the remainder of the sample. In the left panel of Fig. 10, we see a slight downward shift in the angular momentum with a median shift from |$\mathrm{\log _{10}|} {\vec{\boldsymbol{L}}}_{orb} (\, {\rm km} \, {\rm s}^{-1} \, {\rm Mpc})| = 1.70$| to 1.60 for the 51 systems that match the analytic orbit. As compared to the measured orbital angular momentum derived from the Gaia EDR3 proper motions, the calculated orbital angular momentum of the real LG is consistent with both samples, though the median is in slightly better agreement with that of the non-isolated LG pairs.

A more interesting notable result is shown in the right plot in Fig. 10, which plots the merger history for each of these samples. Here, we show that pairs with |$|C| \lt 0.1$| have had much more quiet recent major merger histories as compared to the sample with |$|C| \gt 0.1$|⁠. There appears to be little difference in the minor merger histories for the two samples. This implies that the analytic two-body orbit discussed in Section 2 is a good proxy for those systems that evolved in isolation, at least in the sense that systems that evolved in isolation are those that have had a quiet major merger history.

In our analysis, we confirmed that there were no close encounters among the 51 systems that align with the analytic orbit criteria. Out of the 77 pairs identified, only 4 were found to have experienced ambiguous pericentric passages that fell within the subset of 51 systems. Furthermore, all 4 systems were characterized by significant pericentric distances (⁠|$\gt 700 \, \mathrm{kpc}$|⁠). This leads us to conclude that accurately predicting the orbits of systems with very close passages remains a challenge.

A question may arise as to whether the MW and M31 are isolated due to the proximity of the massive LMC satellite. Current estimates of the LMC’s mass faction to the MW range from 10 per cent to 20 per cent (Peñarrubia et al. 2016; Erkal et al. 2019; Shipp et al. 2021), which is considered a minor merger by our definition. The right plot of Fig. 10 shows that minor mergers do not affect our ability to predict the orbit; only the recent (⁠|$z\lt 0.25$|⁠) major (⁠|$\mu _{\mathrm{ DM}} \ge 0.3$|⁠) mergers do. By this definition, the LMC is included in this isolation criteria and thus may not significantly affect our ability to predict the orbital history of the MW/M31 using this two-body approximation.

This conclusion remains valid even if the simulation pairs do not represent systems with an LMC-like satellite at present day. This is because our criteria for isolated pairs is to predict the orbit within 10 per cent accuracy from present day to a universe age of |$t=3.2$| Gyr. This traces the orbit far enough back in time to where there must have been an LMC- like object for most of these systems at some point. Hartl & Strigari (2022) studied these pairs and could not find LMC-like candidates with present day distances comparable to the observed distance between the LMC and the MW. However despite the simulation sample lacking present day LMC satellites, due to the time frame covered by our analysis to predict the orbits we can still conclude that the system may still be in isolation with the existence of the LMC, as long as there are no major mergers after |$z=0.25$|⁠.

5 DISCUSSION AND CONCLUSION

In recent years, there has been significant progress in the measurements of the relative motion of the MW and M31. Most notably, it is no longer valid to assume that the MW and M31 are on a purely radial orbit, as new proper motions measurements suggest a non-negligible tangential velocity. Furthermore, it has been realized that the LMC may introduce an additional correction to the relative velocity, possibly the nature of the orbit. Accounting for these new measurements, in this paper we study the orbit, angular momentum, and merger rate for MW–M31 like systems, using a combination of an analytic orbit model and simulations of LG-like systems.

We find that the uncertainty on the tangential velocity, combined with the uncertainty in the total mass of the MW and M31, results in a wide possible range of orbital histories. Additionally, correcting for the LMC influence may alter the approach velocity significantly, which we showed when taken into account, results in a higher probability that the MW and M31 have undergone a pericentric passage and are approaching pericenter for the second time today (see Benisty (2021) for previous discussion of possible past encounters for the MW and M31).

We use a large sample of |$\sim 600$| MW/M31 analogue pairs in the IllustrisTNG 300–1 simulation to study their orbital history and orbital angular momentum. We look at all previous simulation snapshots and calculate each pair’s orbital history. While we show the average orbit of all pairs implies that the MW and M31 are on first infall, we find 77 (13 per cent) of these pairs have already undergone a pericentric passage. Out of these 77 pairs, we identified 15 that were isolated by satisfying a smooth criterion, with orbits showing clear pericentric passages with pericenters ranging from 60 to 700 kpc. Although, we did not observe a difference in the present day orbital angular momentum between this sample and the sample of pairs at first infall, we did find a slight correlation with the merger rate. The pairs that had undergone a pericentric passage generally had a more active merger history.

We use these analogue pairs to determine how well the two-body approximation predicts the orbital history between the pair when known kinematic parameters are provided. We find we are only able to predict the orbit within 10 per cent accuracy for |$\sim 10~{{\ \rm per\ cent}}$| of the sample, and that we can predict the orbit best for those pairs in which the MW and M31 have not undergone a major merger since |$z \approx 0.25$|⁠. So, if we can conclude that the MW and M31 have not undergone a merger since |$z=0.25$|⁠, then we can be more confident in the orbital history we recover using this two-body model. While this is likely true for the merger history of the MW, M31 may have had a more active recent major merger history.

As a result of our analysis we derive the LG orbital angular momentum to be |$\mathrm{\log _{10}}|\vec{\boldsymbol{L}}_{\mathrm{ orb}} \mathrm{(\, {\rm km} \, {\rm s}^{-1} \, {\rm Mpc})| }= {1.80_{-0.21}^{+0.14}}$|⁠, and check its consistency with our sample from the Lambda cold dark matter (⁠|$\Lambda \mathrm{ CDM}$|⁠) simulation. We find that the magnitude of the angular momentum as derived from the methods is in good agreement and that the range of the magnitude of the angular momentum is consistent with what is seen in various simulations which identify LG candidates (Forero-Romero et al. 2013; Carlesi et al. 2016; Hartl & Strigari 2022). We then used these pairs to study the evolution of the orbital angular momentum. We calculate the orbital angular momentum for each pair at each snapshot, and we find that each the orbital angular momentum of the pairs steadily increase up to the present day (Fig. 6 bottom). Though we do not specifically address the origin of this angular momentum, one possibility is tidal torques acting to the LG during the course of its evolution (White 1984; Porciani, Dekel & Hoffman 2002).

We found a wide range of initial separations between the pair at the snapshot of their first identification in the simulation. We split our sample based on those with initial distance at the time of identification as |$\gt 400$| kpc and those with |$\lt 400$| kpc. While we did not find a correlation between the present day orbital angular momentum and the present day separation, we did interestingly find a correlation such that pairs with a larger initial separation have a larger orbital angular momentum today. The fact this correlation is preserved through cosmic time strongly suggests that the magnitude of the orbital angular momentum is an imprint of the conditions on the system at the time of formation.

The finding that pairs that evolved from larger separations have a larger orbital angular momentum may also reflect the influence of tidal torques on these systems. Most simply, the magnitude of the tidal torque is proportional to the quadrupole moment of the proto-LG mass distribution, in addition to the quadrupole moment of the external gravitational potential. The fact that all systems indicate a steady increase in orbital angular momentum suggests that a build of tidal torques occurs for all systems. But the fact the initial and final angular momentum are both correlated to the formation distance, as seen in Fig. 8, indicates the magnitude may be set by the initial formation conditions of the LG.

Our results have highlighted that the uncertainty in tangential velocity measurement and the total mass of the MW and M31 preclude, at this stage, definitive answers to the questions of the LG formation history. The kinematics and mass measurements may improve with future measurements from Gaia, which are expected to improve in precision due to single-epoch astrometry. Further, improved measurements of the LG kinematics would provide information on the direction of the LG orbital angular momentum vector, which may be compared to the matter distribution in the local volume. These topics will be addressed in a forthcoming study.

ACKNOWLEDGEMENTS

We are grateful to Jorge Penarrubia and Carton Zeng for discussions on this paper. We acknowledge support from DOE Grant de-sc0010813, and by the Texas A&M University System National Laboratories Office and Los Alamos National Laboratory.

DATA AVAILABILITY

The data underlying this article will be shared on reasonable request to the corresponding author.

REFERENCES