Implications of the Milky Way Travel Velocity for Dynamical Mass Estimates of the Local Group

Following past work that utilizes the "Timing Argument," we assume that the orbital trajectories of the MW and M31—the LG system—over cosmic history are well described by Keplerian orbits (e.g., Kahn & Woltjer 1959; Lynden-Bell 1981; Kroeker & Carlberg 1991; Li & White 2008; van der Marel et al. 2012; Peñarrubia et al. 2016). By assuming that M31 and the MW are gravitationally bound and were last at closest approach in the early universe (i.e., the two galaxies have not yet strongly interacted), we can then use the present-day kinematics of M31 relative to the MW to estimate the total mass of the LG (i.e., using the Timing Argument).

In this work, we largely follow the methodology and notation defined in Peñarrubia et al. (2016). Briefly recapping the classical Timing Argument method, we assume that the dynamics of the MW and M31 pair is dominated by the local gravitational potential of the LG, and therefore the Hubble flow can be neglected for computing the relative orbits of the galaxies (see, e.g., Peñarrubia et al. 2014). Since we observe the relative position and motion between M31 and the MW, we reduce the dynamics of the galaxies in the LG system to a single Keplerian orbit that specifies the relative orbit between the galaxies and is completely determined by four model parameters: the total mass of the LG, M_LG, the semimajor axis, a, the eccentricity, e, and the present value of the eccentric anomaly, η.

In terms of these four model parameters, the closed-form equations for relevant two-body quantities that are closer to observables, like the separation between the masses, r, the elapsed time since last pericenter, t, and the radial and tangential velocity components, v_rad and ${v}_{\tan }$, are given by

$$\begin{eqnarray}&&r=a\,(1-e\,\cos \eta ),\end{eqnarray} \tag{ 1 }$$

$$\begin{eqnarray}&&t={\left(\displaystyle \frac{{a}^{3}}{{GM}}\right)}^{1/2}(\eta -e\,\sin \eta ),\end{eqnarray} \tag{ 2 }$$

$$\begin{eqnarray}&&{v}_{\mathrm{rad}}={\left(\displaystyle \frac{{GM}}{a}\right)}^{1/2}\displaystyle \frac{e\,\sin \eta }{1-e\,\cos \eta },\end{eqnarray} \tag{ 3 }$$

$$\begin{eqnarray}&&{v}_{\tan }={\left(\displaystyle \frac{{GM}}{a}\right)}^{1/2}\displaystyle \frac{\sqrt{1-{e}^{2}}}{1-e\,\cos \eta }.\end{eqnarray} \tag{ 4 }$$

In the expressions above, r is the separation between the centers of the MW and M31 halos, the time since last pericenter, t, is the age of the universe, and the velocity components, $({v}_{\mathrm{rad}},{v}_{\tan })$, express the radial and tangential velocity components of M31 relative to the center of the MW halo.

In a simpler universe where the MW and M31 are point masses and there are no other massive bodies in the LG system, we could transform the observed heliocentric sky position, distance, and velocity of M31 to a MW Galactocentric reference frame and combine these with an estimate of the Hubble time to obtain the four "observables" $(r,t,{v}_{\mathrm{rad}},{v}_{\tan })$. These observables would be enough to infer the four model parameters (M_LG, a, e, η) using Equations (1)–(4).

To describe this "classical" Timing Argument approach in more detail and set the stage for extending it, we adopt the notation of Peñarrubia et al. (2016) in which v _A→B represents the velocity vector of A as measured in the reference frame of B and x _A→B represents the position vector of A as measured from B. With this notation, v _A→B = − v _B→A and v _A→C = v _A→B + v _B→C.

In the classical Timing Argument, the MW disk and M31 are assumed to occupy the center of the potential well of their dark matter halos and have zero velocity with respect to the halos. We refer to this reference frame in the MW dark matter halo as MW_halo. Thus, the position and velocity of M31 with respect to the MW can be represented by ${{\boldsymbol{x}}}_{{\rm{M}}31\to {{\rm{MW}}}_{\mathrm{halo}}}$ and ${{\boldsymbol{v}}}_{{\rm{M}}31\to {{\rm{MW}}}_{\mathrm{halo}}}$, which are assumed to be equivalent to the position and velocity of M31 with respect to the center of the MW disk. Then, the observed position and velocity of M31, measured in a heliocentric reference frame, are given by

$$\begin{eqnarray}&&{{\boldsymbol{x}}}_{{\rm{M}}31\to \odot }={{\boldsymbol{x}}}_{{\rm{M}}31\to {{\rm{MW}}}_{\mathrm{halo}}}+{{\boldsymbol{x}}}_{{{\rm{MW}}}_{\mathrm{halo}}\to \odot },\end{eqnarray} \tag{ 5 }$$

$$\begin{eqnarray}&&{{\boldsymbol{v}}}_{{\rm{M}}31\to \odot }={{\boldsymbol{v}}}_{{\rm{M}}31\to {{\rm{MW}}}_{\mathrm{halo}}}+{{\boldsymbol{v}}}_{{{\rm{MW}}}_{\mathrm{halo}}\to \odot }.\end{eqnarray} \tag{ 6 }$$

Here $\left|{{\boldsymbol{x}}}_{{\rm{M}}31\to {{\rm{MW}}}_{\mathrm{halo}}}\right|=r$ as determined from Equation (1), ${{\boldsymbol{v}}}_{{\rm{M}}31\to {{\rm{MW}}}_{\mathrm{halo}}}$ is determined completely by the Keplerian model parameters (through v_rad and ${v}_{\tan }$), and the position and velocity of the center of the MW halo as measured from the Sun are ${{\boldsymbol{x}}}_{{{\rm{MW}}}_{\mathrm{halo}}\to \odot }$ and ${{\boldsymbol{v}}}_{{{\rm{MW}}}_{\mathrm{halo}}\to \odot }$.

However, the true dynamics of the MW–M31 system are not so simple. Perturbations introduced by interactions and mergers between the MW and its satellite galaxies, such as the merger of the Sagittarius dwarf galaxy or the recent pericentric passage of the LMC, break the assumption that the MW disk is stationary in the center of its dark matter halo. In fact, these interactions will introduce an additional reflex motion component in observations from the MW disk compared to the MW_halo reference frame. The LMC's impact on the dynamics of the MW disk and inner halo have been studied in detail by, e.g., Gómez et al. (2015), Garavito-Camargo et al. (2019, 2021), and Petersen & Peñarrubia (2020). These works imply that we must include additional terms in Equations (5) and (6) to account for the travel velocity of the MW disk with respect to the center of the halo in its unperturbed state. Thus, the observed position and velocity vectors of M31 from the solar reference frame become

$$\begin{eqnarray}&&{{\boldsymbol{x}}}_{{\rm{M}}31\to \odot }={{\boldsymbol{x}}}_{{\rm{M}}31\to {{\rm{MW}}}_{\mathrm{halo}}}+{{\boldsymbol{x}}}_{{{\rm{MW}}}_{\mathrm{halo}}\to {{\rm{MW}}}_{\mathrm{disk}}}+{{\boldsymbol{x}}}_{{{\rm{MW}}}_{\mathrm{disk}}\to \odot },\end{eqnarray} \tag{ 7 }$$

$$\begin{eqnarray}&&{{\boldsymbol{v}}}_{{\rm{M}}31\to \odot }={{\boldsymbol{v}}}_{{\rm{M}}31\to {\mathrm{MW}}_{\mathrm{halo}}}+{{\boldsymbol{v}}}_{{\mathrm{MW}}_{\mathrm{halo}}\to {\mathrm{MW}}_{\mathrm{disk}}}+{{\boldsymbol{v}}}_{{\mathrm{MW}}_{\mathrm{disk}}\to \odot },\end{eqnarray} \tag{ 8 }$$

where MW_halo refers to a reference frame centered at and moving with the center of mass of the outer MW halo, MW_disk refers to a reference frame centered at and moving with the center of the MW disk, and ${{\boldsymbol{x}}}_{{{\rm{MW}}}_{\mathrm{disk}}\to \odot }$ and ${{\boldsymbol{v}}}_{{{\rm{MW}}}_{\mathrm{disk}}\to \odot }$, respectively, are the adopted solar position and velocity in the Galaxy. The values we adopt for ${{\boldsymbol{x}}}_{{{\rm{MW}}}_{\mathrm{disk}}\to \odot }$ and ${{\boldsymbol{v}}}_{{{\rm{MW}}}_{\mathrm{disk}}\to \odot }$ (shortened to x _⊙ and v _⊙) are given in Table 1 below. ⁵

Table 1. Observational Datasets Used for Comparison Throughout Analysis and Their References

	vdMG08 Dist. + HST PM	Cepheid Dist. + Gaia PM	Cepheid Dist. + HST PM
D (kpc)	770 ± 40a	761 ± 11 kpcf	761 ± 11f
vrad (km s−1)	−301 ± 1b	−301 ± 1b	−301 ± 1b
${\mu }_{{\alpha }^{\ast }}\,(\mu {\rm{as}}\,{{\rm{yr}}}^{-1})$	34.30 ± 8.25c	48.98 ± 10.47g	34.30 ± 8.25c
μδ (μas yr−1)	−20.22 ± 7.71c	−36.85 ± 8.03g	−20.22 ± 7.71c
x ⊙ (kpc)	${\left(-\mathrm{8.29,0,0}\right)}^{{\rm{d}}}$	(−8.122, 0, 20.8)h	(−8.122, 0, 20.8)h
v ⊙ (km s−1)	(11.1, 251.54, 7.25)d	(12.9, 245.6, 7.78)i	(12.9, 245.6, 7.78)i
tperi (Gyr)	13.75 ± 0.11e	13.801 ± 0.024j	13.801 ± 0.024j

Note. Each value is measured for M31 with respect to the Sun. D is the distance, v_rad is the radial velocity, and $({\mu }_{\alpha }^{* },{\mu }_{\delta })$ are proper motions in R.A. cosdec and Decl. x _⊙ = (x, y, z) and v _⊙ = (U_pec, V_pec + V₀, W_pec) are the the position of the Sun and the solar motion with respect to the Galactic center, with the x-axis pointing from the projection of the Sun on the disk toward the Galactic center, and the z-axis pointing in the direction of the North Galactic pole. t_peri is the time elapsed since the last pericenter of the M31 Keplerian orbit, which in this case is the age of the universe.

References. (a) van der Marel & Guhathakurta (2008); (b) Courteau & van den Bergh (1999); (c) van der Marel et al. (2012); (d) Schönrich et al. (2010), McMillan (2011); (e) Jarosik et al. (2011); (f) Li et al. (2021); (g) Salomon et al. (2021); (h) Gravity Collaboration et al. (2018), Bennett & Bovy (2019); (i) Drimmel & Poggio (2018); (j) Planck Collaboration et al. (2020).

Download table as:  ASCIITypeset image

Observationally, the reflex motion of the disk imprints itself on velocity measurements as an instantaneous velocity shift. Recently, Petersen & Peñarrubia (2021) used tracers in the outer MW stellar halo to measure this instantaneous travel velocity ${{\boldsymbol{v}}}_{\mathrm{travel}}={{\boldsymbol{v}}}_{{{\rm{MW}}}_{\mathrm{halo}}\to {{\rm{MW}}}_{\mathrm{disk}}}$. They find ∣ v _travel∣ = 32 ± 4 km s⁻¹ with a highest likelihood apex direction in Galactocentric coordinates of ${({\ell },b)}_{\mathrm{apex}}=({56}_{-9}^{+9},-{34}_{-9}^{+10})$ degrees. We assume that ${{\boldsymbol{x}}}_{{{\rm{MW}}}_{\mathrm{halo}}\to {{\rm{MW}}}_{\mathrm{disk}}}\approx 0$ motivated by the fact that this displacement is likely much smaller than the distance between the MW and M31 ${{\boldsymbol{x}}}_{{{\rm{MW}}}_{\mathrm{halo}}\to {{\rm{MW}}}_{\mathrm{disk}}}\ll r$ (as expected from simulations; e.g., Garavito-Camargo et al. 2021). However, it is important to note that this displacement is still significant on scales relevant for many other MW studies (see Section 4.3 for more details).

Figure 1 shows a schematic of these different vectors—all drawn in a frame that is comoving with the MW_halo frame—and a rough illustration of the geometry we assume. For clarity, we show ${{\boldsymbol{x}}}_{{{\rm{MW}}}_{\mathrm{halo}}\to {{\rm{MW}}}_{\mathrm{disk}}}={{\boldsymbol{x}}}_{\mathrm{travel}}$, ${{\boldsymbol{v}}}_{{{\rm{MW}}}_{\mathrm{halo}}\to {{\rm{MW}}}_{\mathrm{disk}}}={{\boldsymbol{v}}}_{\mathrm{travel}}$, and v _M31→⊙ = v _obs.

Zoom In Zoom Out Reset image size

The present distance and relative velocity of M31, as well as the age of the universe (used in Equation (2)), are key observables that are used to constrain our Timing Argument model. In this paper, we consider three different compilations of data to understand how different measurements might affect the Timing Argument model with the addition of the travel velocity. In particular, we consider two different M31 distance measures: an approximated distance measure from van der Marel & Guhathakurta (2008), and a more accurate Cepheid-based distance measure from Li et al. (2021). We also consider two different M31 proper-motion measurements: Hubble Space Telescope (HST)-based proper motions from van der Marel et al. (2012), and a more recent Gaia early Data Release 3 (eDR3)-based proper-motion measure from Salomon et al. (2021).

We have split these measurements into three compiled data sets:

• 1.  
  vdMG08 Dist. + HST PM: the M31 distance measure from van der Marel & Guhathakurta (2008) and HST proper motions from van der Marel et al. (2012), the same data set as used to constrain the LG mass via the Timing Argument in van der Marel et al. (2012).
• 2.  
  Cepheid Dist. + Gaia PM: a compilation of more recent M31 kinematic measurements, including a more precise Cepheid-based distance measure to M31 (Li et al. 2021) and updated Gaia eDR3 proper motions from Salomon et al. (2021).
• 3.  
  Cepheid Dist. + HST PM: a hybrid data set with the Cepheid-based distance measure to M31 and the HST-based proper-motion measurement (Li et al. 2021; van der Marel et al. 2012).

The HST-based proper motions were originally presented in Sohn et al. (2012), and were then corrected for the internal kinematics and space motion of M31 in van der Marel et al. (2012), from which we used the "weighted average" heliocentric velocities in Table 4. ⁶ The Gaia-based M31 proper-motion measurement is slightly larger than the HST proper motion of M31, leading to an increased implied transverse velocity that, a priori, should lead to a higher inferred LG mass compared to the more radial orbit implied by the HST proper motions. See Table 1 for numerical values used in each of these data sets.

We construct a likelihood function, ${ \mathcal L }$, to quantify the probability of measuring the observed quantities ${\boldsymbol{y}}\,=\{D,{v}_{\mathrm{rad}},{\mu }_{\alpha }^{* },{\mu }_{\delta },{t}_{\mathrm{peri}}\}$ given a timing argument model with parameters θ = (M_LG, a, e, η, α). Here, D is the distance to M31, v_rad is the radial velocity of M31, ${\mu }_{\alpha }^{* }$ and μ_δ are the proper-motion components, and t_peri is the time since last pericenter. The parameter vector θ contains M_LG, the total mass of the LG, a, the semimajor-axis, e, the eccentricity, η, the present value of the eccentric anomaly, and α, a nuisance parameter discussed in detail later in this section. We assume that the measurements are independent and have Gaussian uncertainties such that the likelihood function is a product:

$$\begin{eqnarray}&&{ \mathcal L }=p({\boldsymbol{y}}| {\boldsymbol{\theta }})=\prod ^{n}\displaystyle \frac{1}{\sqrt{2\pi }{\sigma }_{n}}\exp \left[-\displaystyle \frac{1}{2}\,\displaystyle \frac{{\left({y}_{n}-{\tilde{y}}_{n}({\boldsymbol{\theta }})\right)}^{2}}{2{\sigma }_{n}^{2}}\right],\end{eqnarray} \tag{ 9 }$$

where n indexes the elements of the data vector, σ_n is the corresponding uncertainty for the nth data element, and ${\tilde{y}}_{n}({\boldsymbol{\theta }})$ is the model-predicted value for a given data component.

We then adopt prior probability distribution functions (pdfs) for the parameters and use these pdfs to compute the posterior pdf over the parameters given the data in order to generate samples from the posterior pdf using a Markov Chain Monte Carlo (MCMC) method. In order to recover the estimated mass distribution, we marginalize over all other model parameters.

In detail, we first use the four Timing Argument parameters to compute the present-day separation between the MW and M31 halos and their relative radial and tangential velocities as defined in Equations (1)–(4). These velocity components represent the relative velocity M31 would have as observed from the center of an unperturbed MW halo. We then use the measured "travel velocity" of the MW disk, v _travel, to find the relative velocity of M31 with respect to the center of the moving MW disk (i.e., a moving MW Galactocentric frame). We finally transform from this Galactocentric frame to a heliocentric reference frame moving with the solar system barycenter (i.e., ICRS coordinates). At this final stage, we must introduce an additional nuisance parameter, α, that represents the orientation of the MW–M31 orbital plane as it intersects the tangent plane located at the sky position of M31 as viewed from the MW disk center. This parameter is needed to convert from the two-dimensional velocity components given by Equations (3)–(4) to the three-dimensional velocity components represented by the two proper-motion components and the radial velocity of M31. However, we stress that this position angle has no impact on the fundamental dynamical parameters and is only used for coordinate transformations.

We specify this model using the Python probabilistic programming package pymc3 (Salvatier et al. 2016) and use the No-U-Turn Sampler (Homan & Gelman 2014) implemented in pymc3 to generate samples from this posterior pdf, given data from each of the data sets defined in Table 1. We sample over the parameters LG mass M_LG, the present-day MW–M31 halo separation r, log eccentricity $\mathrm{ln}\left(1-e\right)$, eccentric anomaly η, and the orbital plane orientation nuisance parameter α. Our adopted prior pdfs are defined in Table 2. For each data set, we run the sampler with four chains for 4000 tuning steps and 40,000 draws.

Table 2. A Description of Our Adopted Prior Probability Distribution Functions Over the Timing Argument Model Parameters

Prior	Description
MLG: ${{ \mathcal N }}_{T}(4.5,3)\times {10}^{12}{M}_{\odot }$	Mass of the Local Group
r: ${{ \mathcal N }}_{T}(700,100)$ kpc	Distance from M31 to MWdisk
$\mathrm{ln}(1-e)$: ${ \mathcal U }(-10,0)$	Eccentricity (close to 1)
η: ${ \mathcal U }(-\pi ,\pi )$	Eccentric anomaly
α: ${ \mathcal U }(-\pi ,\pi )$	Position angle of M31 orbital
	plane from MW disk center

Note. Here, ${ \mathcal U }(a,b)$ represents a uniform distribution over the domain (a, b), and ${{ \mathcal N }}_{T}(\mu ,\sigma )$ represents a truncated normal distribution with mean μ and standard deviation σ. We truncate the mass prior pdf to the range (0.5, 20) × 10¹² M_⊙ and the distance prior pdf to the range (100, 10⁴) kpc.

Download table as:  ASCIITypeset image

There have been two primary pathways toward measuring the mass of the LG: measure the masses of MW and M31 individually and add them, or go after the total mass directly with LV dynamics, the Timing Argument, or cosmological simulations.

Historically, the total mass estimates have been much larger than the sum of the individual MW and M31 masses: typical total LG masses are upwards of 4 × 10¹² M_⊙, while the sum of independent MW+M31 mass measures result in a total mass closer to 2–2.5 × 10¹² M_⊙, as can be seen in the collection of previous mass estimates for the LG, M31, and the MW in Table 4. It is not surprising that there are many discrepancies between the total and summed values of the LG mass, since MW+M31 (i.e., individual summed) mass estimates often cannot measure the full extent of the distribution of dark matter within each galaxy, let alone in the LG, and thus must extrapolate in regions where there may be large uncertainties.

Table 4. A Collection of Discussed Previous Mass Measurements from LG Dynamics Focusing on Previous Timing Argument Results

Mass	Method	Result (1012 M⊙ )	Citation
MLG	TA	3.6	Lynden-Bell (1981)
	TA (radial + cosmo sim calibration)	5.27	Li & White (2008)
	TA only	4.27 ± 0.45	van der Marel et al. (2012)
	TA (3D + cosmic bias and scatter)	4.93 ± 1.63	van der Marel et al. 2012
	Local Group (LG) dynamics	2.5 ± 0.4	Diaz et al. (2014)
	Local Volume galaxies + TA + Λ	2.64±0.4	Peñarrubia et al. (2016)
	Machine learning	4.9 ± 0.8	McLeod et al. (2017)
	Machine learning (+large M31 transverse motion)	3.6 ± 0.3	McLeod et al. (2017)
	Cosmological sims	${4.4}_{-1.5}^{+2.4}$	Zhai et al. (2020)
	Cosmological sims (likelihood-free inference)	${4.6}_{-1.8}^{+2.3}$	Lemos et al. (2021)
	TA + Λ + cosmological sims	${4.75}_{-2.41}^{+2.22}$	Hartl & Strigari (2022)
	TA + Λ + LMC	${5.6}_{-1.2}^{+1.6}$	Benisty et al. (2022)
	TA + Λ + cosmic biad + LMC	${3.4}_{-1.1}^{+1.4}$	Benisty et al. (2022)
	TA + v travel (Cepheid + HST)	${4.0}_{-0.3}^{+0.5}$	Chamberlain et al. (2022; this work)
	TA + v travel (Cepheid + Gaia)	${4.5}_{-0.6}^{+0.8}$	Chamberlain et al. (2022; this work)
MM31	Kinematics of M31 sats	1.4 ± 0.4 (<300 kpc)	Watkins et al. (2010)
	Giant Stellar Stream	${2.00}_{-0.41}^{+0.52}$	Fardal et al. (2013)
	LG dynamics	1.7 ± 0.3	Diaz et al. (2014)
	Local Hubble Flow	1.33 ± 0.4	Peñarrubia et al. (2016)
	M31 orbital ang. mom.	${1.37}_{-0.75}^{+1.39}$	Patel et al. (2017a)
	Cosmological sims	1.0 − 2.0	Carlesi et al. (2017)
	Cosmological sims	${2.5}_{-1.1}^{+1.3}$	Zhai et al. (2020)
	Machine learning	2.3 − 2.5	Villanueva-Domingo et al. (2021)
	Machine learning (+ velocity information)	2.2–2.5	Villanueva-Domingo et al. (2021)
MMW	Kinematics of LG sats	1.4 ± 0.3 (<300 kpc)	Watkins et al. (2010)
	LG dynamics	0.8 ± 0.5	Diaz et al. (2014)
	Local Hubble Flow	1.04 ± 0.26	Peñarrubia et al. (2016)
	LMC orbital ang. mom.	${1.02}_{-0.55}^{+0.77}$	Patel et al. (2017a)
	Cosmological sims	0.6 − 0.8	Carlesi et al. (2017)
	MW sats	${0.96}_{-0.28}^{+0.29}$	Patel et al. (2018)
	Cosmological sims	${1.5}_{-0.7}^{+1.4}$	Zhai et al. (2020)
	Machine learning	1.0–1.3	Villanueva-Domingo et al. (2021)
	Machine learning (+ velocity information)	2.3–2.6	Villanueva-Domingo et al. (2021)

Download table as:  ASCIITypeset image

There have been a few exceptions in the trend of high total mass measurements, namely in Diaz et al. (2014) and Peñarrubia et al. (2016). Diaz et al. (2014) utilize the fact that the LG momentum should balance to zero in the frame of the LG barycenter to determine the total mass of the LG, as well as the mass ratio between M31 and the MW. Using the LG barycenter, indicated by a set of 17 LG satellites at >350 kpc from the MW and M31, and the velocities of M31 and the MW with respect to the barycenter, they found a LG mass of 2.5 ± 0.4 × 10¹² M_⊙, and a mass ratio M_M31/M_MW > 2.29. The impact of the LMC and M33 are absorbed by assuming they contribute to the masses of their host galaxies. However, the recent measurement of the MW travel velocity will change the measured barycenter and velocities of the MW, M31, and the LG satellites, though it is unclear how this would affect the total mass in their analysis.

Additionally, Peñarrubia et al. (2016) found a LG mass of ∼2.64 ± 0.4 × 10¹² M_⊙, which is significantly lower than our findings, although this constraint combines the Timing Argument dynamics of the M31–MW system in addition to the observed kinematics of 35 LV galaxies. They parameterize the offset of the LMC+MW barycenter from a MW-only barycenter as a function of the mass ratio between the LMC and MW, and find that the LMC likely has ∼25% of the mass of the MW halo, resulting in a large shift in the barycenter of the LMC+MW. The recovered mass using only the dynamics of the LV galaxies was quite low (∼2 × 10¹² M_⊙), though including the Timing Argument dynamics of the MW–M31 system increased their recovered total mass to ∼2.64 ±0.4 × 10¹² M_⊙.

Our TA+ v _travel LG mass estimates are consistent with a number of other recent studies that estimate the mass of the LG through dynamical methods. For example, we find agreement with the previous Timing Argument model of van der Marel et al. (2012), which found a total mass of 4.27 ± 0.45 × 10¹² M_⊙ (neglecting for cosmic bias and scatter) using the same values for the distance and velocity of M31 as in our vdMG08 Dist. + HST PM data set.

In two N-body cosmological simulations, Millenium-WMAP7 and MilleniumII, Zhai et al. (2020) identified pairs of stellar analogs to the MW and M31, then applied a series of kinematic cuts on the separation, isolation, and velocities of the pair to determine LG analogs. They find stellar and dynamical LG analogs on mostly radial orbits have total masses of ${4.4}_{-1.5}^{+2.4}\times {10}^{12}\,{M}_{\odot }$, which is consistent with our findings for each data set. They also find that low-ellipticity orbits (where ${v}_{\mathrm{rad}}\sim {v}_{\tan }$), result in a higher LG mass, M31 mass, and MW mass, reporting ${6.6}_{-1.5}^{+2.7}\times {10}^{12}\,{M}_{\odot }$, ${3.8}_{-1.8}^{+2.8}\times {10}^{12}\,{M}_{\odot }$, and ${2.5}_{-1.4}^{+2.2}\times {10}^{12}\,{M}_{\odot }$, respectively. We find that the ellipticity of the orbit decreases dramatically with larger travel velocities; however, in contradiction to Zhai et al. (2020), we find that this leads to a lower total LG mass.

Additionally, Benisty et al. (2022) recently modeled the contribution of the LMC-induced shift in the barycenter of the MW and used the Timing Argument to place constraints on the LG mass while removing the impact of the LMC. More specifically, they estimate the contribution of the reflex motion of the MW disk to the observed velocity vector of M31 by modeling the orbital history of the MW–M31 system and the MW+LMC–M31 system. The impact of the LMC on the mass measurements of the LG thus depend on the orbital and mass models of the LMC about the MW. This is in contrast to our work, which need make no assumption about the mass of the LMC, its orbital history, or the merger history of the MW. Upon implementing the Timing Argument, including a cosmological constant, and removing the LMC-induced reflex motion, they find a LG mass of ${5.6}_{-1.2}^{+1.6}\times {10}^{12}\,{M}_{\odot }$, which is roughly 25% larger than our findings, and that accounting for cosmic bias and scatter lowers the mass by an additional 40% to ${3.4}_{-1.1}^{+1.4}\times {10}^{12}\,{M}_{\odot }$, 25% lower than our findings. However, they find that, in general, inclusion of the motion of the MW disk due to the interaction with the LMC lowers the LG mass by ∼10%, which is consistent with our finding of a reduction in the total LG mass by ∼10%–12%.

A notable difference between this study and Benisty et al. (2022) is the method by which the reflex motion of the MW disk is accounted for. Rather than relying on accurately simulating the reflex motion of the disk, we let the travel velocity of the MW disk introduce a coordinate transformation (boost) of the measured velocity vectors of M31, and fit for the model parameters given the observable data. This method allows us to avoid model uncertainties in the mass profiles of each galaxy in the orbital models of the interaction between the MW, M31, and the LMC. Additionally, using the measured travel velocity also allows us to innately account for possible additional contributions to the present-day, instantaneous travel velocity induced by the extensive merger history of the MW in a robust way, without the need to simulate the entire interaction history.

Given the simplicity of the Timing Argument dynamical model—in particular, the assumption that the MW and M31 are point masses with constant masses—it is reasonable to wonder whether this methodology provides unbiased estimates of the true LG mass. An early study of a dark-matter-only cosmological simulation found that the Timing Argument applied to pairs of galaxies did provide unbiased estimates of the sum of masses of the pairs (Li & White 2008). However, more recently it was found that conditioning on LG analogs with similar radial and tangential velocities to the MW and M31 leads to slightly biased (overestimated) inferred total masses of those systems (González et al. 2014; Hartl & Strigari 2022). In this work, we do not attempt to "correct" our inferred LG masses for this cosmic bias effect, because it is unclear whether cosmological simulations accurately reproduce the detailed properties of LG systems. Accounting for this effect would likely lower our reported LG mass measurements. However, as we have shown, the existence of a reflex motion of the MW disk as a response to the MW's interaction with its satellite population decreases is an additional perturbation to the TA that must be considered in future studies alongside cosmic bias and a cosmological constant.

The existence of the travel velocity of the MW disk results in measurable differences in the estimated mass of the LG through the Timing Argument. These results are consistent with Erkal et al. (2020), who find that neglecting the LMC-induced reflex motion of the MW can result in masses that are overestimated by up to 50%. As shown in our results above, neglecting this motion at its currently measured value of v _travel = 32 ± 4 km s⁻¹ (Petersen & Peñarrubia 2021) leads to LG masses that are overestimated by ∼30%. However, both the magnitude and direction of the travel velocity are directly tied to the inferred mass of the LG. As it is currently measured, the (highest likelihood) direction of v _travel is ∼60° from the sky position of M31, meaning that the travel velocity impacts the conversion of both M31's observed proper motion and radial velocity from a heliocentric reference frame to the "outer halo" reference frame used above. At fixed magnitude, if the true apex of the travel velocity motion is closer (farther) to M31's sky position, it would primarily affect the radial velocity (proper motions).

There is reason to believe that the recently measured MW disk travel velocity could be a lower bound on the true value, which could be up to a factor of ∼2–3 higher than the currently measured value. Using an idealized simulation of an equilibrium dark matter halo that has a recent merger with an LMC-like halo, Garavito-Camargo et al. (2021) showed that stellar halo tracers at different distances from the MW disk center may result in different measured travel velocities. While this simulation does not span previous mergers in the MW's history, it gives a good first-order approximation of what we may expect to observe. At fixed apex direction, a larger travel velocity would correspond to a lower inferred LG mass. Figure 3 shows the effect of increasing the measured travel velocity magnitude up to these predicted values and the impact on the inferred mass of the LG for each of the data sets used in this work (see Section 2.2). The vertical lines in this figure show the LMC-induced travel velocities that are predicted for three tracer distances in simulations from Garavito-Camargo et al. (2021). Thus, future measurements of the travel velocity of the disk that use tracers at larger distance around the MW stellar halo will likely lead to a lower inferred LG mass. We note again that the value of the eccentricity is derived assuming a matter-only universe (i.e., we neglect other cosmological effects in the orbit computation, as discussed above).

Zoom In Zoom Out Reset image size

The biggest source of uncertainty in our empirically inferred LG mass, M_LG, currently comes from the proper-motion measurements, which have signal-to-noise ratios of just 3–4. Future data releases from the Gaia Mission (Gaia Collaboration et al. 2016) will lead to more precise mean proper motions of M31 and thus more precise Timing Argument constraints on the LG mass. For example, between Gaia eDR3 and the end of the extended (10 yr) mission, the expected individual-source proper-motion precision improvement for a G = 20 source (i.e., an upper-giant-branch star in M31) is a factor of ∼6. Naïvely scaling the proper-motion uncertainties of M31 as measured with Gaia (Salomon et al. 2021) by a factor of 6 leads to a ∼2×improvement in the M_LG precision. Of course, the true improvement of the mean M31 proper motion with improved individual source kinematics could be even better than linear because more sources will be detected and usable in the measurement.

While M31 has a massive satellite (M33) of comparable mass ratio to the MW–LMC system, we do not expect there to be a significant reflex velocity of M31's disk relative to its equivalent outer halo reference frame. Recent work predicts that M33 is likely on first infall into the M31 halo and has a much larger orbital pericenter than the MCs (e.g., Patel et al. 2017b). Additionally, M31 has likely experienced other significant mergers, as evidenced by the double nucleus and Giant Southern Stream (e.g., Ibata et al. 2001; Font et al. 2006), but these were likely lower mass-ratio mergers (e.g., Gilbert et al. 2019; Milošević et al. 2022) and thus will have less of an impact on the bulk motion of the M31 disk. Given current knowledge of the M31 system and uncertainties in the orbital histories of its most massive satellites, here we neglect the reflex motion of the M31 disk. However, a measurement or upper limit on the M31 disk travel velocity would enable further unbiased constraints on the LG mass.

Reconciling techniques that compute the LG mass from the summed MW+M31 mass and from the Timing Argument is not straightforward, but the two approaches are complimentary. Future measurements of the travel velocity at large Galactocentric distances will likely exceed current measurements (see Section 4.3), which directly implies a lower LG mass and may improve agreement between these two general methods for estimating the LG mass. Since constraints on the MW mass consistently find a mass of ∼10¹² M_⊙, if not slightly higher (see, for reference, Table 4), our mass limits from the Timing Argument may begin to place meaningful upper limits on the mass of the M31 system.

The Timing Argument is highly sensitive to the tangential motion of M31, and larger proper motions will generically lead to lower eccentricities and higher inferred LG masses. Recent Gaia proper motions of M31 suggest the orbit of MW–M31 is less radial than previously believed (van der Marel et al. 2019; Salomon et al. 2021). Neglecting the travel velocity ( v _travel = 0 in the above plots), we find that proper motions from HST (van der Marel et al. 2012) are consistent with a highly radial orbit e ∼ 0.96, while proper motions from Gaia (Salomon et al. 2021) result in a slightly lower eccentricity of e ∼ 0.88. As seen in Figure 3, we find that as the travel velocity increases, the inferred eccentricity of the MW–M31 orbit decreases. However, contrary to expectations, we find that the inferred LG mass also decreases. This implies that the velocity contribution to the relative and transverse velocity of M31 is dominant to the change in eccentricity.

Additionally, studies of the Bolshoi N-body cosmological simulation by Forero-Romero et al. (2013) find that typical LG analogs, when selected via mass and isolation criteria, do not have completely radial orbits. As with studies that measure a higher transverse velocity for M31, and thus a slightly less radial orbit (van der Marel et al. 2019; Salomon et al. 2021), we find that the decrease in the eccentricity due to the measured travel velocity makes the LG less eccentric and, thus, more cosmologically typical.