The Local Group Mass in the light of Gaia

High accuracy proper motions (PMs) of M31 and other Local Group satellites have now been provided by the {\it Gaia} satellite. We revisit the Timing Argument to compute the total mass $M$ of the Local Group from the orbit of the Milky Way and M31, allowing for the Cosmological Constant. We rectify for a systematic effect caused by the presence of the Large Magellanic Cloud (LMC). The interaction of the LMC with the Milky Way induces a motion towards the LMC. This contribution to the measured velocity of approach of the Milky Way and M31 must be removed. We allow for cosmic bias and scatter by extracting correction factors tailored to the accretion history of the Local Group. The distribution of correction factors is centered around $0.63$ with a scatter $\pm 0.2$, indicating that the Timing Argument significantly overestimates the true mass. Adjusting for all these effects, the estimated mass of the Local Group is $ M = 3.4^{+1.4}_{-1.1} \times 10^{12} M_{\odot}$ (68 % CL) when using the M31 tangential velocity $ 82^{+38}_{-35}$ km/s. Lower tangential velocity models with $59^{+42}_{-38}$ km/s (derived from the same PM data with a flat prior on the tangential velocity) lead to an estimated mass of $ M = 3.1^{+1.3}_{-1.0} \times 10^{12} M_{\odot}$ (68 % CL). By making an inventory of the total mass associated with the 4 most substantial LG members (the Milky Way, M31, M33 and the LMC), we estimate the known mass is in the range $3.7^{+0.5}_{-0.5} \times 10^{12} \, M_{\odot}$.


INTRODUCTION
The Timing Argument (TA) is a simple way of working out the mass M of the Local Group (LG). In its earliest manifestation (Kahn & Woltjer 1959), the Milky Way and M31 proto-galaxies are assumed to have a small separation at the time of the Big Bang. They travel away from each other in the Hubble flow (see Banik & Zhao 2016, for a general relativistic derivation). If there is enough mass present, their expansion is reversed. Given an estimate of the age of the universe, together with their present separation and velocity of approach, then the equations of motion can be solved to give the mass of the Local Group. Kahn & Woltjer (1959) did this assuming the Milky Way and M31 are on an exactly radial orbit. Einasto & Lynden-Bell (1982) showed that the problem retained an analytic solution, even if M31 has some tangential motion.
Since then, many elaborations of the TA have been proposed, including: (i) the effects of the tidal influ-ence of galaxies outside the Local Group, which is minor (Raychaudhury & Lynden-Bell 1989); (ii) the introduction of a Cosmological Constant, which manifests itself as an additional expansion term and increases the LG mass by ≈ 10 % (Partridge et al. 2013); (iii) corrections for the effects of hierarchical growth of the two galaxies by comparisons with cosmological simulations (Kroeker & Carlberg 1991;Li & White 2008), which shows the TA is (mostly) unbiased though suffers from cosmic scatter; (iv) successive refinements of the M31 orbit in view of the improving observational accuracy of M31's tangential motion (van der Marel et al. 2012;van der Marel et al. 2019).
In this Letter, we revisit the TA. Our motivations are twofold. First, we wish to exploit the new measurement of the proper motion of M31 provided by the Gaia satellite Salomon et al. 2021). Second, we identify sources of systematic error in applications of the TA which needs correction. Recent work, again driven by data from the Gaia satellite, has shown that the Large Magellanic Cloud is much more massive than originally envisaged. It is pulling the central parts of the Milky Way towards it Petersen & Peñarrubia 2021;Garavito-Camargo et al. 2021), and so the measured line of sight velocity of M31 needs correction before the TA can be applied. This effect on the TA was noted before by Penarrubia et al. (2016), though their analysis method differs from that presented here. Equally important, we include the corrections for the TA-derived mass calibrated on pairs of galaxies extracted from cosmological simulations (Hartl & Strigari 2022), which have orbits similar to those of M31 and the Milky Way. In this procedure, it is of crucial importance to ensure that the mock pairs match the LG as closely as possible.
The material is arranged as follows. In § 2, we estimate the known mass in the LG from stellar/satellite kinematics. § 3 implements our corrections to the TA, while § 4 provides a discussion of the results.

THE LG MASS BUDGET FROM KINEMATICS
The LG mass can be estimated by modelling the kinematics of tracers (halo stars, satellite galaxies and globular clusters or HI gas) around its prominent members. This assumes that the dark matter is clustered around the major galaxies, and not distributed throughout the LG (as originally envisaged in Kahn & Woltjer 1959).
The advent of Gaia data has substantially reduced the uncertainty on the virial mass of the Milky Way M MW , with most recent measurements satisfying 1.17 +0.21 −0.15 × 10 12 M (Callingham et al. 2019, see also Watkins et al. (2019) and Fritz et al. (2020) for similar results). M31 is more massive than the Milky Way with dynamical arguments suggesting M M31 = 1.8 ± 0.5 (Shull 2014, see also Diaz et al. (2014) and Karachentsev & Kudrya (2014) for similar values). This suggests that the mass associated with the two largest galaxies in the LG is 3.0 +0.5 −0.5 × 10 12 M . Note that the virial mass of M31 is more uncertain than that of the Milky Way, and it remains (just) possible that the mass ratio is close to unity (e.g., Evans & Wilkinson 2000;Fardal et al. 2013;Kafle et al. 2018), which gives us a lower bound.
The next most massive members of the LG in decreasing order are: M33 with M M33 = 5.0 ± 1.0 × 10 11 M (Corbelli et al. 2014;Kam et al. 2017) and the Large Magellanic Cloud with M LMC = 1.8 ± 0.4 × 10 11 M (Erkal et al. 2019;Shipp et al. 2021). M32 is a compact dwarf elliptical with current mass at least an order of magnitude less than M33. Its progenitor may once have been more massive than M33, though much of its tidally stripped material is now in the halo of M31 (D'Souza & Bell 2018), and so already accounted for in our inventory. We surmise that the minor members of the LG contribute about 0.7 × 10 12 M to the mass budget. We conclude that the total mass in the LG -as judged from kinematics of tracers -is at least 3.0×10 12 M and most likely in the range 3.7 ± 0.5 × 10 12 M .

The Data
We take the current values of the separation between the Milky Way and M31 as r = 770 ± 40 kpc and the heliocentric line-of-sight velocity as v los = −301 ± 1 km s −1 (van der Marel et al. 2012). The measurement of the proper motion (PM) of M31 has been refined over the last decade. van der Marel et al. (2012) used Hubble Space Telescope (HST) observations in three small off-centered fields in conjunction with a model for its internal kinematics, deriving the mean PM of M31 to be µ α = 45 ± 13 µas yr −1 , µ δ = −32 ± 12 µas yr −1 . This value is dominated by the Solar velocity with respect to the Milky Way centre, which corresponds to µ reflex {α,δ} = {38, −22} µas yr −1 at the distance of M31 (770 ± 40 kpc). Thus the reflex-corrected PM is consistent with zero, and they gave a 1σ upper limit on the tangential velocity at 34 km s −1 . More recently, van der Marel et al. (2019) computed an independent estimate of M31's absolute PM from Gaia Data Release 2, which had twice larger uncertainties than the HSTbased value, and is also larger in absolute sense: µ α,δ = {65 ± 18, −57 ± 15} µas yr −1 , with an additional systematic uncertainty of 16 µas yr −1 in each component. The error-weighted average of the two independent measurements is µ α,δ = {49 ± 11, −38 ± 11} µas yr −1 and corresponds to a reflex-corrected tangential velocity of 57 +35 −31 km s −1 . Salomon et al. (2021) used the updated Gaia Early Data Release 3 astrometry to measure µ α,δ = {49 ± 11, −37 ± 8} µas yr −1 , which is very close to the weighted average derived in van der Marel et al. (2019); however, they reported the reflex-corrected tangential velocity to be 82 ± 31 km s −1 . The discrepancy between the two studies stems from the way the distribution of tangential velocities is derived from the distribution of PM. van der Marel et al. (2019), following their earlier work (section 3.1 in van der Marel & Guhathakurta 2008), derive the posterior distribution of the magnitude of the two-dimensional tangential velocity vector v tan by convolving the observed Gaussian PM distribution with a prior on the tangential velocity P |v tan | , which they take to be flat in |v tan |, favouring smaller values. By contrast, Salomon et al. (2021) simply convert both reflex-corrected PM components into velocity and sum them in quadrature, which effectively means using a flat prior on each component of v tan , i.e., P |v tan | ∝ |v tan |.
As there is no convincing reason in favour or against the use of van der Marel & Guhathakurta (2008)'s prior, we consider both alternatives, drawing Monte Carlo samples from the distribution of observed PM values (taken from Salomon et al. 2021 in both cases) with or without the additional reweighting by the prior. It is intuitively clear that increasing the M31 tangential velocity makes the M31 orbit less eccentric and therefore the LG mass must increase (for fixed age of the Universe).

The TA Algorithm with Cosmological Constant
The center of mass coordinate system is defined by the relative distance r = | r M31 − r MW | and the relative velocity v = d r/dt. The masses are replaced by the total mass M := m MW + m M31 . In polar coordinates (r, ϕ), the relative distance variation now reads (Emelyanov et al. 2015;Carrera & Giulini 2006;Emelyanov & Kovalyov 2013 where l is the conserved angular momentum per mass (l = r 2φ = r v tan   (2008) results in lower values (short dashes), while using the raw PM measurements without any reweighting produces higher values (long dashes). In both cases, the radial velocity is −114 ± 1 km s −1 . After compensating for the LMC perturbation as described in § 3.3, we find that vtan is increased by 25-30 km s −1 , and v rad is shifted to −75 ± 15 km s −1 . This distribution is shown by short and long dot-dashed contours and shaded in green. These distributions are compared to the radial and tangential velocities of galaxy pairs selected from the IllustrisTNG cosmological simulation, in which the mass correction factor is computed as the ratio of the actual combined mass of both galaxies to the mass obtained from TA, as described in Hartl & Strigari (2022)  and PM of M31 therefore include a contribution due to this downward motion, which we would like to remove. Penarrubia et al. (2016) accounted for the presence of the LMC by assuming that it forms a two-point-mass system with the Milky Way, and that M31 moves around the barycentre of this combined system. The displacement and velocity shift of the Milky Way relative to the barycentre are obtained by multiplying the relative position and velocity of the LMC in the Milky Way-centered frame by the mass ratio M LMC /(M MW + M LMC ). These offsets need to be subtracted from the current position and velocity of M31 in the Milky Way-centered frame prior to computing its trajectory in the barycentric system. This argument is qualitatively correct, but ignores the fact that the LMC is currently only ∼ 50 kpc from the Milky Way centre, and that the enclosed mass of the Milky Way within this radius is substantially smaller than its total mass -in other words, it underestimates the actual displacement of the central region of the Milky Way caused by the LMC.
A more sophisticated technique to compensate for the LMC perturbation was recently introduced by Correa Magnus & Vasiliev (2022). It starts by computing the past trajectory of the Milky Way and the LMC under their mutual gravitational attraction, using the actual (distance-dependent) force from each galaxy rather than the point-mass approximation implied by Penarrubia et al. (2016)'s method. Once this has been done for a given choice of Milky Way and LMC potentials, we inte-grate the orbit of M31 in this time-dependent potential of both galaxies backward in time until the LMC perturbation is negligible, and the integrate it forward without the LMC to the present epoch. For simplicity, the Milky Way potential is fixed to an NFW halo with virial mass M vir = 1.1 × 10 12 M , virial radius r vir = 270 kpc and a concentration c = 13.5, but we take into account the uncertainty on the LMC mass by sampling it from a lognormal distribution centered on log 10 (M LMC ) = 11.15 with width 0.15 dex and repeating the orbit rewinding step for each choice of LMC mass.
The left-hand panel of Figure 2 shows the posterior distribution of M31's Galactocentric radial and tangential velocity components v rad , v tan , with or without the LMC correction. The use of a prior from van der Marel & Guhathakurta (2008) results in a lower tangential velocity, v tan = 59 ± 34 km s −1 , while the use of raw PM measurements produces a higher v tan = 78 ± 32 km s −1 . In both cases, the compensation of the LMC perturbation increases v tan by 25-30 km s −1 , changes v rad from −114 ± 1 km s −1 to −75 ± 15 km s −1 , and increases the distance by ∼ 40 kpc. The marginalized posterior distributions of v tan for all four cases are also shown in the left-hand panel of Figure 3. Compared to the simpler model for the barycentric motion used in Penarrubia et al. (2016), our velocity correction is roughly twice higher for the given LMC mass, but since that paper used a significantly larger range of LMC masses with a median at 2.5 × 10 11 M , the velocity correction is quite similar in absolute terms.

Cosmic Bias and Scatter
Owing to simplifications in the TA, the mass estimate may suffer from systematic bias and scatter. Li & White (2008, see their Figure 1) found that the TA mass is unbiased, though with some scatter, using analogues of the Milky Way-M31 pair extracted from the dissipationless Millenium simulation. However, González et al. (2014) noted the TA mass is only unbiased on average, and can be an overestimate if the pairs are restricted to to have similar radial and tangential velocities as the true Milky Way and M31. The matter has been reinvestigated recently by Hartl & Strigari (2022), who used the IllustrisTNG N-body and hydrodynamical simulations. They also found a tendency of the TA mass to be overestimated. Specifically, Hartl & Strigari (2022) identify 580 bound analogues of the LG by a series of cuts on B band magnitude, separation, velocity of approach and total velocity, computing distributions of P (A), where A is the ratio of true mass to mass predicted by the TA. We tailor the Hartl & Strigari (2022) sample by imposing three new cuts: (i) a separation between 650 and 950 kpc; (ii) a mass ratio within a factor of 4; (iii) −150 < v rad < v tan − 100 km s −1 so it resembles the actual distribution of LMC-corrected velocities as shown in Figure 2. This retains 160 galaxy pairs, with the distribution P (A) shown the right panel of that figure centered around A = 0.63 with a scatter ±0.2.
To incorporate uncertainties, we use the Markov chain Monte Carlo method. We sample 10 4 values for the initial conditions vector {r,ṽ rad ,ṽ tan , t 0 , Λ} and calculate the corresponding predicted mass, and then convolve this TA-predicted mass distribution with the distribution of correction factors P (A). We assume the initial conditions for t 0 and Λ have a normal distribution with the mean and dispersion given by the estimate from Planck and its reported uncertainty.

RESULTS AND DISCUSSION
The right panel of Fig. 3 shows the posterior distributions of the the total LG mass obtained for the high M31 tangential velocity case, including or not the correction for the LMC perturbation and cosmic scatter. The central panel shows the two-dimensional posterior distributions. The results are summarized in Table 1, together with the low tangential velocity case. We also record the posterior probability P 1 that M lies in the range 3.7 +0.5 −0.5 × 10 12 M . This is informed by our inventory of the LG mass in § 2.
The effect of correcting for the LMC's perturbation is to decrease the inferred LG mass by ∼ 10%. Although  the correction increases the tangential velocity of M31, which normally would lead to a higher LG mass from TA, it also affects the radial velocity and the distance, so the net result is the opposite (a downward shift). To lowest order, this counteracts the effect of the Cosmological Constant, which acts in the opposite sense by a similar amount (e.g., Partridge et al. 2013;Benisty 2021;Benisty & Guendelman 2020). Note that further corrections to the infall velocity are probably also needed because of the effects of M33 and M32. M33 and M31 came within ∼ 50 kpc of each other in the past, ∼ 6.5 Gyr ago (Tepper-García et al. 2020), whilst M32 may even have been more massive than M33 before its catastrophic encounter with M31. Although the effects of these interactions require detailed modelling, they are probably 10% (as less important than the LMC).
The effect of correction for cosmic bias and scatter is substantial, reducing the median LG mass by a factor ∼ 1.5 and increasing its relative uncertainty by a similar amount. In constructing our distribution of correction factors A, we ensured as much as possible that our mock LGs match the distribution of LMC-corrected velocities of infall. There is a wide range of accretion histories in any mock LGs extracted from simulations. It is important to condition distributions on the true environment of the LG as much as possible, as first clearly realised by González et al. (2014). The pure TA mass -in our case, 6.0 +1.3 −0.9 M -can then be a serious overestimate of the true mass. From Table 1, we see that the mass of the LG is 3.4 +1.4 −1.1 × 10 12 M with the raw data (the high v tan case). If the van der Marel & Guhathakurta (2008) prior is used, then the mass is 3.1 +1.3 −1.0 × 10 12 M (the low v tan case). Both are in reasonable accord with the dynamically estimated LG mass as quoted in § 2.
Since the LG is assumed to be a closed system, the total energy is negative E < 0, otherwise M31 could approach infinity. From this limit, we get an expression for the minimal mass (cf Chernin et al. 2009): GM min = rv 2 /2 − Λc 2 r 3 /6. For the low tangential velocity, the minimal mass is 1.27±0.11×10 12 M , and for the high tangential velocity, the minimal mass is a bit larger: 1.61 ± 0.24 × 10 12 M . These numbers may be compared with the observationally derived lower limit to the LG mass of ≈ 3 × 10 12 M in § 2. Fig. 4 compares the value obtained in this paper with other recent measurements. Notice that our masses are somewhat lower than (though still consistent with) a number of other recent estimates, such as numerical im-plementation of least action (Phelps et al. 2013;Banik & Zhao 2017), simulations (González et al. 2014;Carlesi et al. 2017), neural networks (McLeod et al. 2017) and likelihood-free density estimation (Lemos et al. 2021). However, estimates in Fig. 4 based on the virial theorem (e.g., Diaz et al. 2014;Hartl & Strigari 2022) or on the assumption of pure radial orbits (Penarrubia et al. 2016) are systematically lower than our values. Whilst the outer parts of the Local Group are not virialized, Hartl & Strigari (2022) showed that that virial mass estimator is unbiased, but the scatter around the true value is much larger for virial mass estimators than for the TA, rendering it a much less satisfactory method. The high radial velocities of some outlying satellites like NGC 3109 have been suggested as evidence for a past encounter between the Milky Way and M31 (Banik & Zhao 2017). This is possible in some modified gravity theories, though not in the ΛCDM paradigm used in this paper.
Overall, we conclude that the LG mass derived via the TA -if calibrated against realistic analogues in simulations -is in reasonable agreement with the mass known to be associated with the Milky Way, the LMC, M31 and M33. Gaia has improved the accuracy with which the first two are known. The main observational uncertainty that remains is the virial mass of M31, on which future work could usefully be concentrated. But, the Timing Argument works better than we have a right to expect such a simple argument to do -once corrected from the effects of the LMC and Cosmic Bias! EV is supported by the Consolidated Grant, whilst DB thanks the Blavatnik and the Rothschild Foundations for support and partial support from European COST actions CA15117 and CA18108.