Local Group timing argument and virial theorem mass estimators from cosmological simulations

IllustrisTNG (Nelson et al., 2019, 2018; Pillepich et al., 2018; Springel et al., 2018; Naiman et al., 2018; Marinacci et al., 2018) is a suite of cosmological gravo-magnetohydrodynamic simulations which follows the evolution of dark matter, cosmic gas, luminous stars, and supermassive black holes within the centers of galaxies. Within this suite, the TNG300-1 is the simulation with the largest number of particles, $N=2500^{3}$, with a dark matter particle mass of $m_{dm}=4.0\times 10^{7}$ M${}_{\odot}h^{-1}$ and a baryonic particle mass of $m_{baryon}=7.6\times 10^{6}$ M${}_{\odot}h^{-1}$. The TNG300-1 simulation is evolved in the large comoving volume of $205\,h^{-1}{\rm Mpc}$. The initial conditions, in agreement with Planck 2015 results (Aghanim et al., 2016), are determined by the cosmological parameters: $\Omega_{(\Lambda,0)}=0.6911$, $\Omega_{(m,0)}=0.3089$, $\Omega_{(b,0)}=0.0486$, $\sigma_{8}=0.8159$, $n_{s}=0.9667$, and $h=0.6774.$ This yields an age of the Universe of 13.8 Gyr.

The Illustris database provides the properties of the galaxies and dark matter halos identified in 100 snapshots each associated with a unique redshift. Here we use the SubLink algorithm (Rodriguez-Gomez et al., 2015) to access the data associated with different merger trees. To identify LG-like systems, we are interested in the following kinematic galaxy and halo properties: 1) The position of the center of mass of each halo, which Illustris reports as the spatial position within the halo volume for the particle with the lowest gravitational potential; 2) the velocity, which is calculated as the sum of all particles velocity weighted by mass; 3) the maximum circular velocity, $V_{max}$, which is calculated from the spherically-averaged rotation curve; and 4) the mass $M_{halo}$ for each halo, which is the sum of all particle masses bound to the subhalo. Note, $M_{halo}$ does not include particle masses bound to subhalos of this halo, and differs from mass definitions used in previous related studies (Peñarrubia & Fattahi, 2017).

In addition to these kinematic properties, we extract the absolute magnitudes for each galaxy. IllustrisTNG calculates the absolute magnitude as the sum of the luminosity for each stellar particle of the halo. IllustrisTNG reports these magnitudes in eight bands; here we utilize the B-band magnitude and denote this as $M_{B}$.

In the $z=0$ snapshot, the subfind algorithm identifies $14.5\times 10^{6}$ dark matter halos in TNG300-1, with the vast majority of these halos hosting galaxies. From this sample of halos, we choose systems that best resemble the LG, based on the following criteria. Our cuts and the corresponding sample sizes subsequent to each cut are summarized in Table 1.

We begin by identifying galaxies with B-band magnitude similar to that of the MW and M31. For the MW the estimated B-band magnitude is $M_{B}=-20.9$, while for M31 it is $M_{B}=-20.7$ (van den Bergh, 1999). To include a sufficiently large number of pairs in our sample, we expand this range and select galaxies with B-band magnitude within $-22.3<M_{B}<-19.3$. This cut results in a sample of 97,574 galaxies.

From this sample of galaxies chosen based on absolute magnitude, we add a cut to restrict to pairs of galaxies that have a separation similar to that of MW and M31. The distance between the MW and M31 is $D=783\pm 25$ kpc  (McConnachie, 2012). Cutting on galaxy pairs within this observed range would result in an unreasonably small sample of LG systems. In order to obtain a sufficiently large sample size of LG pairs, we take as our distance cut pairs of galaxies with a separation $500$ kpc $<|\vec{r}|<1000$ kpc. The lower limit of 500 kpc is motivated to ensure that the dark matter halos of the two pairs galaxies do not overlap, while the upper limit of $1000$ kpc is motivated to ensure that the galaxies are bound and are detached from the Hubble Flow.

Even after performing the above distance cut, we must be careful to ensure that there is not a third galaxy, with magnitude comparable (i.e. within $-22.3<M_{B}<-19.3$) to the MW or M31, in close proximity to the identified pair. To account for this possibility, we keep only those pairs that have two unique halo ids, so that each member of the pair is not also assigned to a third pair. With this additional cut, our combined sample after absolute magnitude and distance cuts contains 4133 pairs.

To compare our magnitude cut to previous cuts that consider stellar masses of galaxies,  Zhai et al. (2020) found 12 pairs in Illustris TNG100-1 using a stellar mass cut of $4.0\times 10^{10}\,\rm{M}_{\odot}<\rm{M}_{MW}<8.0\times 10^{10}\,\rm{M}_{\odot}$ for the Milky Way mass, and $8.0\times 10^{10}\,\rm{M}_{\odot}<\rm{M}_{M31}<13.0\times 10^{10}\,\rm{M}_{\odot}$ for M31 mass. In this work, they also use a slightly stricter distance cut of $600$ kpc $<|\vec{r}|<1000$ kpc. When we perform these stellar mass cuts and distance cuts to Illustris TNG300-1, we find 427 pairs. Accounting for the difference between the simulation volumes of TNG100-1 and TNG300-1, the sample of LGs that we obtain is consistent with that of  Zhai et al. (2020).

For our next cut, we ensure that each pair is sufficiently isolated from another massive system that would gravitationally perturb the LG pair. For this cut, we remove pairs if a galaxy with $M_{B}<-19.3$ is found within a distance of $3$ Mpc from the barycenter of the pair, or if a dark matter halo with $V_{max}>150$ km s^-1 is found within $1$ Mpc of the barycenter of the pair. This isolation criteria is similar to that of Fattahi et al. (2016), who use an isolation cut of 2.5 Mpc. These cuts are implemented to identify pairs with a similar local environment to the observed LG. Our magnitude cut is motivated by demanding that nearby galaxies are fainter than the most luminous nearby galaxy outside of the LG, Centaurus A, which is $\sim$ 3 Mpc from the LG. The $V_{max}$ cut ensures the pair dominates the gravitational potential within its local volume, and thus dominates the dynamics of the LG-like system. After this isolation cut we are left with 798 isolated pairs.

The MW and M31 are observed to approach one another with a velocity (in the Galactocentric frame) of $V_{r}=-109.3\pm 4.4\,{\rm km}\,{\rm s}^{-1}$. The bound on the tangential velocity from HST observations is $V_{t}<34$ km/s (van der Marel et al., 2012), while recent measurements from Gaia DR2 (van der Marel et al., 2019) and EDR3 (Fardal et al., 2021) indicate larger mean values of $V_{t}\simeq 79$ km/s. All these measurements are still consistent with the observed MW-M31 system being on a nearly radial orbit. To match this kinematic criteria with our sample of pairs, we calculate the total relative radial velocity, $V_{r}$, for the two main galaxies of each pair. Since the simulation provides peculiar velocities, we add in the expansion velocity to obtain the total relative radial velocity for the pair. If we cut our above sample to simply ensure that the galaxies are approaching with $V_{r}<0\,{\rm km}\,{\rm s}^{-1}$, on top of the above magnitude, distance, and isolation cuts we are left with a sample of 658 pairs. In the analysis below we examine the impact of restricting the range of allowed radial velocities even further.

Note that in choosing our sample based on absolute magnitude, we have neglected the impact of the dark matter halo mass, or similarly the maximum circular velocity, in our cuts. To better understand the range of maximum circular velocities that our absolute magnitude cuts correspond to, in Figure 1 we show the relation between absolute magnitude and maximum circular velocity. As indicated, there is a significant scatter in $V_{max}$ over the absolute magnitude range that we consider. To remove galaxies with $V_{max}$ much larger than that of the MW or M31 galaxies, we add the constraint $V_{max}<300$ km s^-1 on top of our above cuts, which results in a sample 613 galaxy pairs.

Our cut in $V_{max}$ is similar to that implemented in Li & White (2008), who used the Millenium simulation to identify LG-like systems based on the range $150$ km s${}^{-1}\leq V_{max}<300$ km s^-1 as their primary sample. Examining Figure 1 suggests that our initial magnitude cut above corresponds to a range in $V_{max}$ that is larger than that considered in Li & White (2008), and does indeed introduce unreasonably massive halos in our sample. An additional benefit of choosing our sample based on magnitude is that it allows us to compare to the results of Li & White (2008), who use halo mass as their primary criteria.

In addition to the two massive galaxies, we will be interested in the fainter galaxies that are gravitationallly-bound to the LG systems. In particular, we are interested in the subsample of galaxies that are bound to the LG system, but are not bound to either of the two massive galaxies. We refer to these galaxies as outer LG members (OLGM).

From our LG sample identified above, we consider two cuts to extract OLGMs: first a cut on the absolute magnitude, and second a cut on the dark matter mass. For an absolute magnitude limit, we take the limit on the B-band magnitude to be $M_{B}<-10$, which approximately corresponds to the faintest satellite in the observed LG that is not bound to either the MW or to M31 (McConnachie, 2012; Simon, 2019). For our dark matter halo cut, we take the lower bounds on the dark matter halo mass of $M_{dm}>10^{10}\,{\rm M}_{\odot}$ and $M_{dm}>10^{9}\,{\rm M}_{\odot}$. This corresponds to a cut on systems with dark matter mass approximately larger than the most massive luminous satellite galaxies of the MW and M31.

On top of these mass and luminosity cuts, we must be careful to identify OLGMs that are bound to either of the two massive galaxies or are too far out into the Hubble flow. For the former, we tag the specific OLGMs that are within $<350$ kpc from either of the two most massive galaxies. This cut tags galaxies that are within the approximate virial radius of the MW and M31. For the latter, we exclude OLGMs that are $>1.5$ Mpc from either of the two massive galaxies, thereby excluding systems that are strongly affected by the Hubble flow. Though we only utilize the $>1.5$ Mpc cut for our fiducial analysis below, we do compare to the results when the $<350$ kpc cut is also implemented and find a negligible difference. Note that the $>350$ kpc and $<1.5$ Mpc distance cut is similar to that used in the observational sample in Diaz et al. (2014).

For OLGMs identified with the distance cut of $D<1.5$ Mpc, Figure 2 shows the distribution of $N$. The left plot shows the distribution of $N$ for the OLGM luminosity cut, $M_{B}<-10$, and the strict OLGM dark matter mass cut of $M_{dm}>10^{10}$ M_⊙. The right plot shows the distribution of $N$ for the relaxed dark matter mass cut of $M_{dm}>10^{9}$ M_⊙. The two distributions on the left plot are similar, though the sample cut on dark matter halo mass $>10^{10}$ M_⊙ has a larger mean of $\langle N\rangle=12.4$, compared to the mean of $\langle N\rangle=11.4$ for the luminosity sample cut. Reducing the lower limit on dark matter halo mass to $>10^{9}$ M_⊙ yields a mean of $\langle N\rangle=67.5$ OLGMs per LG. When initially requiring each OLGM to be energetically bound, the number of OLGMs associated with each pair decreases. This reduces the LG sample size when we require $N>0$. This sample yields a similar behavior with a mean of $\langle N\rangle=8.3$ for the luminosity sample, $\langle N\rangle=7.9$ for the $>10^{10}$ M_⊙ dark matter halo mass cut, and $\langle N\rangle=56.4$ for the relaxed $>10^{9}$ M_⊙ dark matter halo mass cut.

We impose two final cuts to obtain our sample of mock-LG systems which we will use as our sample to study the VT mass estimator below. First, we require at least 10 OLGMs per LG, i.e., $N\geq 10$. Requiring a minimum number of OLGMs is motivated by the number of OLGMs observed in the LG, $N=17$. Second, we place a restriction on the mass ratio $m=M_{M31}/M_{MW}$. We take the constraint on this quantity from  Diaz et al. (2014), who used outer LG kinematics to estimate $m=2.30$. Since our choice of which galaxy in the main pair is the MW or M31 is arbitrary, we restrict this value to a range of $0.3\leq m\leq 3.3$.

In addition to their number, we are also interested in characterizing the kinematics of the OLGM sample. To examine the kinematics, we first identify the position of the barycenter, where the barycenter is calculated as the center of mass of the two massive galaxies in each system. We then determine the radial velocity, $V_{r,\imath}$, of each OLGM relative to the LG barycenter, excluding the two massive galaxies. The radial velocity dispersion of the OLGMs is then calculated as $\sigma_{r}^{2}=\frac{1}{N-1}\sum_{\imath=1}^{N}(V_{r,\imath}-\bar{V}_{r})^{2}$, where $\bar{V}_{r}$ is the mean radial velocity of the sample. This assumes that the radial velocity dispersion is constant throughout the LG, which is an appropriate assumption given the small sample sizes that we are considering.

Figure 3 (left) plots the distribution of the radial velocity dispersion for the sample with the distance cut of $D<1.5$ Mpc and the additional cuts $N\geq 10$ and $0.3\leq m\leq 3.3$ described above. Implementing these cuts yields 169 pairs, 234 pairs, 425 pairs for the absolute magnitude, $M_{dm}>10^{10}\,\rm{M}_{\odot}$, and $M_{dm}>10^{9}\,\rm{M}_{\odot}$ cuts respectively. The histogram distribution of the radial velocity dispersion for the three samples are similar. We see the sample of mock-LG’s with OLGMs identified with the luminosity cut has the largest mean radial velocity dispersion of $97.3\,{\rm km}\,{\rm s}^{-1}$. The mean radial velocity dispersions are $79.2\,{\rm km}\,{\rm s}^{-1}$ and $82.0\,{\rm km}\,{\rm s}^{-1}$ for $M_{dm}>10^{10}\,\rm{M}_{\odot}$, and $M_{dm}>10^{9}\,\rm{M}_{\odot}$ respectively.

The velocity dispersions we obtain for the three samples are similar, and further are consistent with the radial velocity dispersion measured from the observed LG (Diaz et al., 2014). Note that the observed velocities are heliocentric radial velocities, which must then be converted to a barycentric frame by measuring the location of the MW-M31 barycenter. While this conversion adds uncertainty to the velocity dispersion as deduced from observations, it does not add an uncertainty to our analysis given that we identify the barycenter directly from the simulation.

Figure 3 (right) shows the corresponding ratio of the radial to the tangential velocity dispersion, $\sigma_{t}^{2}$, for each LG in the sample. The tangential velocity is defined so that for an isotropic distribution, $\sigma_{r}^{2}/\sigma_{t}^{2}=0.5$, where this value is denoted by the vertical line. The peak at $\sigma_{r}^{2}/\sigma_{t}^{2}\simeq 2$ implies that OLGMs are on dominantly radial orbits, likely due to infall into the LG. Below we discuss the implications of $\sigma_{r}^{2}/\sigma_{t}^{2}$ on the equilibrium nature of the LG.

The TA has been long studied as a LG mass estimator (Kahn & Woltjer, 1959). The TA assumes that the MW and M31 can be described by a Keplerian orbit, with $z=0$ boundary conditions given by present-day separation, relative velocity, and the age of the Universe. From these assumptions and observations, the total LG mass $M_{LG}=M_{MW}+M_{M31}$ may be determined.

Implementing the TA involves solving for the evolution of the MW-M31 separation. In $\Lambda\mathrm{CDM}$, this relative separation, $\vec{r}$, is given by

Equation 1 is solved for $r(t)$ in 3D, given the above boundary conditions and the condition that $r\rightarrow 0$ as $t\rightarrow 0$.

In some prior implementations of the TA, two simplifications have been invoked. The first involves neglecting the vacuum energy term in Equation 1. A second simplification involves neglecting the tangential component of the relative motion between the MW and M31. This is a reasonable approximation, given that the magnitude of the tangential velocity is small relative to the radial velocity. In order to obtain the most accurate solution, we numerically solve Equation 1 and account for both the vacuum energy term and the non-zero tangential velocity.

For each mock LG in our sample, we determine the separation $|\vec{r}|$, as well as the relative velocity components $V_{r}$ and $V_{t}$. Then for a chosen $M_{LG}$, we take these as initial conditions to numerically solve Equation 1 by integrating backwards in time. The best solution for $M_{LG}$ is the one for which the separation between the two massive galaxies approaches zero as $t\rightarrow 0$ Gyr. Practically, the solution for $M_{LG}$ is obtained by extending the integration to times $t<0$, and generating a library of $M_{LG}$ and times at which the separation goes to zero. The best value for $M_{LG}$ is obtained by interpolating between the smallest positive $t$ and the smallest negative $t$ obtained. Then given this solution for the “timing argument mass," $M_{TA}$, we can compare to the true LG mass from the simulation to estimate the bias in the TA mass estimator.

Next we consider a LG mass estimator based on the VT. This estimator was first employed in Courteau & van den Bergh (1999) given the sample of OLGMs known at the time, and has been developed more recently with updated data on OLGMs in Diaz et al. (2014).

For a system in equilibrium, the VT relates the kinetic energy $(K)$ and the potential energy $(U)$ via $2K+U=0$. Assuming that the system is isotropic, so that the random motions are the same in each of the coordinate directions, the kinetic energy is $K=3/2NM_{0}\sigma_{r}^{2}$, where $N$ is the number of galaxies in the sample, $M_{0}$ is the mass of each galaxy, and $\sigma_{r}$ is the one-dimensional radial velocity dispersion of the system. The factor of 3 arises from the assumption of isotropy.

We take two contributions to the potential energy $U$: the gravitational contribution from the two massive galaxies, and the contribution from the cosmological constant. On the $\imath^{th}$ galaxy, the potential is then

where $\rho_{\Lambda}$ is the present value of the vacuum energy density, and

Here $a_{MW}$ is the scalelength of the MW’s dark matter halo, and $a_{M31}$ is the corresponding scalelength of the M31’s dark matter halo. To model the DM density profiles of the MW and M31, we have followed Diaz et al. (2014) and adopted a Hernquist profile, and for the scale radii we take $a_{MW}=a_{M31}=40$ kpc. We find that reasonable changes in the scale radius for each halo do not affect our results.

Taking the approximation that all the satellite galaxies have equal mass, the LG mass can be estimated as

where $\rho_{\Lambda}=0.69\times 1.5\times 10^{7}$ M_⊙ ${\rm Mpc}^{-3}$, and as defined above $m=M_{M31}/M_{MW}$ is the mass ratio between the MW and M31. Note that this formula is similar to that derived in Diaz et al. (2014), the only difference being the vacuum energy term which is added to the kinetic energy term.

With our LG sample, we implement Equation 5 by determining $\sigma_{r}$ and the mass ratio $m$ for each system, and from these determine $M_{LG}$. We refer to this solution as the virial theorem mass, $M_{VT}$. Note that since $m$ is not a priori known for the MW-M31 system, we must take realistic bounds on this quantity when comparing to observations. We discuss in detail below how the assumption for $m$ affects our results.

We begin with an analysis of the TA mass estimator. Our analysis follows a similar approach to that of Li & White (2008), who utilize the dark matter-only Millenium simulation. These authors determine the bias in the TA estimate, $M_{TA}$, as compared to the true LG mass as obtained from the simulation, where the true mass of the LG from the simulation, $M_{tr,halo}$, is defined as the sum of the MW/M31 halo analogs, $M_{tr,halo}=M_{halo,1}+M_{halo,2}$. These authors define $A\equiv M_{tr,halo}/M_{TA}$, and study the distribution of this quantity for their mock LG sample from Millenium. We use our sample of LGs to obtain a distribution of $A$, which provides an estimate of the bias in the TA mass estimator for the Illustris simulation.

To get a sense of the relative importance of the different velocity cuts, we consider three subsamples of our 613 LG pairs identified above. In the first subsample, we restrict the tangential velocity to be less than the absolute value of the radial velocity, i.e., $V_{t}<|V_{r}|$. This is our loosest cut, motivated to simply ensure that the dominant component of the velocity is in the radial direction. In the second subsample, we restrict the radial velocity range to $-195$ km s${}^{-1}\leq V_{r}\leq-65$ km s^-1, reducing our sample to 251 pairs. This cut is motivated to span the uncertainty in the observed radial velocity between the MW and M31, and was implemented with similar motivation by Li & White (2008). For our third subsample, we consider the union of the first and second cuts. This is our smallest subsample, resulting in a total of 177 pairs.

To compare our magnitude cut to previous cuts that consider $V_{max}$,  Li & White (2008) found 11838 pairs in Millennium simulation using $150$ km s${}^{-1}<V_{max}<300$ km s^-1 as their initial cut. For this sample of 11838 pairs, they use the same distance cut of $500$ kpc $<|\vec{r}|<1000$ kpc, and radial velocity cut of $-195$ km s ${}^{-1}\leq V_{r}\leq-65$ km s^-1. When we perform the $V_{max}$, distance and radial velocity cuts on Illustris TNG300-1, we find 1051 pairs. Accounting for the difference in the simulation volumes, $500\,h^{-1}{\rm Mpc}$ and $205\,h^{-1}{\rm Mpc}$ for the Millennium and IllustrisTNG simulations respectively, our final number of mock-LG pairs is consistent with the number of pairs found by  Li & White (2008).

Figure 4 plots the cumulative (left) and histogram (right) distributions of $A$ for all three velocity samples. In Figure 4 (left) also plots the best fit gaussian curve (dashed-black) for the preferred velocity sample with $-195$ km s${}^{-1}\leq V_{r}\leq-65$ km s${}^{-1},V_{t}<|V_{r}|$. From this we see that the distribution is very close to gaussian. In particular, comparing the mean of 0.75 to 0.77 and standard deviation of 0.22 to 0.27 for the curve fit and Gaussian fit respectively, we see our sample has a relatively normal distribution.

In Table 2 we summarize our statistical results, reporting the interquartile and the 90% range of the distributions for all three velocity cuts. The median value of $A$ and the interquartile ranges are less than unity, with $A\simeq 0.74-0.83$, depending on the specific velocity cut implemented. This implies that there is a slight bias in the TA mass estimator to predict a mass higher than the true mass obtained from the simulation. However, when considering the 90% range about the mean, this bias is removed, with the upper bounds on the distribution in the range $A\simeq 1.16-1.22$.

To understand this bias, we also calculate the timing argument mass assuming no cosmological constant, i.e. $\Lambda=0$. This calculation yields a median value of $A$ in the range $A\simeq 0.82-0.93$. Thus, including the cosmological term in calculating the TA mass estimator results in a shift to a median value of $A\simeq 0.74-0.83$, indicating a larger timing argument mass estimate. This particular shift in the TA mass estimator can be understood by considering the current distance and approach velocity between the pair. When we include the cosmological term in the calculation, we account for a repulsive term in the potential. Thus, a more massive system is estimated in order to account for the current distance and approach velocity. Recall these distributions record $M_{tr,halo}/M_{TA}$, thus a higher estimate in the timing argument mass lowers the distribution.

It is interesting to compare our median values of the distribution of $A$ to that of Li & White (2008). In addition to selecting LG systems based on maximum circular velocity, these authors do not include the effect of the cosmological constant term, nor do they include the tangential velocity. Assuming $V_{t}=0$ yields a lower bound to the TA mass estimate. In addition, excluding the cosmological term has the affect described above. Thus the combination of $V_{t}=0$ and $\Lambda=0$ results in their slight calculated bias in the median value of $A\simeq 0.99$. This is consistent with González et al. (2014), who find that including the tangential velocity overestimates the true mass by a factor $\sim 1.3-1.6$.

The different velocity cuts that we consider impact the tail of the $A$ distributions. In particular, constraining the range of the radial velocity to $-195$ km s${}^{-1}\leq V_{r}\leq-65$ km s^-1 reduces the width of the distribution. For different velocity cuts, the distribution widens because the sample includes systems with either very large or very small radial velocities. In addition requiring the tangential velocity to be less than the radial velocity, $V_{t}<|V_{r}|$, the width of the distribution is unchanged, though the median values shift up slightly. This upward shift in the median values reflects the smaller TA mass estimate for systems with small tangential velocity.

We are now in position to calculate the TA mass of the LG, and to apply a correction to this estimate based on our analysis above. We assume $\Omega_{(\Lambda,0)}=0.69$ and we utilize the observed Local Group kinematics of $V_{t}=79$ km s^-1, $V_{r}=-109.3$ km s^-1 and $r=785$ kpc. For these assumptions, we calculate a TA mass of

The observational uncertainty on the TA mass is primarily due to the measured uncertainties of the radial and tangential velocities. van der Marel et al. (2012) quote a statistical uncertainty of $0.45\times 10^{12}$ M_⊙. Noting 5% and 95% containment points of the distribution for our preferred sample are separated by a factor of $2.97$, we see that the systematic uncertainty of the TA mass is much larger than the measured uncertainty.

From Table 2, the median value for the distribution of $A$ for the preferred sample of mock-LGs is 0.79. Thus, the corrected value to the timing argument mass is

The 5% and 95% c.l. on the $A$ distribution are $0.39$ and $1.16$ respectively. This yields a range of $2.34\times 10^{12}\,\text{M}_{\odot}\leq M^{true}_{TA,LG}\leq 6.97\times 10^{12}\text{ M}_{\odot}$ with $90\%$ confidence.

For comparison, assuming that $\Lambda=0$, we obtain

For $\Lambda=0$, the median of the $A$ distribution is $0.87$, and the 5% and 95% range yield values of 0.45 and 1.29, respectively. Thus, using the median value to correct the TA mass yields

in a range of $2.44\times 10^{12}\text{ M}_{\odot}\leq M^{true}_{TA,LG}\leq 7.01\times 10^{12}\text{ M}_{\odot}$ with $90\%$ confidence. Note that the 90% confidence ranges for the $M_{true,LG}$ with and without the cosmological constant term are very similar.

Our finding of the mild impact of the cosmological constant in the TA mass is consistent with previous work (Partridge et al., 2013). This is because, when including the cosmological constant, the magnitude of the effective gravitational potential felt by the system is reduced relative to the case of a pure matter-only model. So for fixed observed parameters ($V_{r}$, $V_{tan}$, and the MW-M31 separation), the mass would have to increase to compensate for the reduced potential, and so the ratio $A$ is shifted to slightly smaller values.

We now move on to an analysis of the VT mass estimator, and compare this to the LG mass derived from the simulations. Here we define the VT mass estimate as $M_{VT}$, and the ratio of the VT mass estimate to the true mass as $B\equiv M_{tr,halo}/M_{VT}$. We consider the same 613 LG-analogs, and again define $M_{tr,halo}$ as the sum of the masses of the two most massive galaxies. For each sample, we use Eq. 5 to calculate the mass of each mock-LG, $M_{VT}$. To begin, we calculate the mass ratio $m$ for each mock-LG by setting the MW and M31 masses to their true values as extracted from the simulation.

Figure 5 shows the distributions of $\log_{10}B$. We plot the sample with the OLGM distance criteria of $D<1.5$ Mpc and the additional cuts, $N\geq 10$ and $0.3\leq m\leq 3.3$. The left plot of Figure 5 shows this distribution for all three samples in cumulative form. The dashed curves indicate the best fit gaussian for the three samples. The right plot of Figure 5 shows histogram distribution of $B$ for all three samples.

In Table 3 we summarize our statistical results, reporting the interquartile and the 90% points of the distributions for all three OLGM identification cuts. The median value of $B$ and the interquartile ranges are less than unity, with $B\simeq 0.80-1.07$, depending on the initial OLGM cut implemented. This implies that in general there is a slight bias in the VT mass estimator depending on the specific choice of OLGM sample. In particular in the case of our main sample, $M_{dm}>10^{9}$ M_⊙, we find a slight bias, $B=0.8$, with a 95% containment interval of $[0.20,1.63]$.

As noted above, current estimates of $m=M_{M31}/M_{MW}$ indicates that the mass of M31 is about a factor of two larger than that of the MW. This motivated an additional cut to our Local Group sample, keeping only those with $0.3\leq m\leq 3.3$. If we instead relax this restriction on $m$, but still use the true values of their masses in the VT, there is minimal change in the $B$ distribution, yielding median values of $B\simeq 0.79-1.01$

When we instead consider the sample with the initial OLGM distance criteria of $350$ kpc $<D<1.5$ Mpc and the additional cuts, $N\geq 10$ and $0.3\leq m\leq 3.3$, we find a larger range for the median value between cuts of $B\simeq 0.82-1.18$. This larger range may be in part because the stricter initial OLGM distance cut yields a smaller number of pairs in our sample.

For all of the above analysis, we chose $m$ to be the true value by extracting the true M31 and MW masses from the simulation. To check the results are not sensitive to the choice of $m$, we simply set $m=1$ in the VT, implying that the masses of the MW and M31 are equal. Note setting $m=1$ implies the barycenter is located at the midpoint of the pair, and this assumption propagates in calculating the velocity dispersion. Referring to our main sample in which we use the distance cut of $D<1.5$ Mpc and restrict the true value of $m$ to be within $0.3\leq m\leq 3.3$, we find the median value to be $B\simeq 0.80-1.09$. Comparing this to $B\simeq 0.80-1.07$ when we use the true value of $m$, we find the VT mass estimate is not sensitive to the value of $m$ assumed in the VT.

Diaz et al. (2014) have estimated the LG mass from OLGMs using the VT. These authors derive this result using radial velocity data on OLGMs within a distance of $350$ kpc $<D<1.5$ Mpc. The result of their analysis is an LG mass estimate of

Referring to our preferred sample with a median of $B=0.8$ and a 90% containment interval of 0.20-1.63, the corrected value for the LG mass estimate from the VT is then

This correction was performed with our preferred sample of all OLGMs $<1.5$ Mpc, though as noted above the $B$ distribution changes only slightly if OLGMs $<350$ kpc are excluded. The uncertainties in Equation 11 are just larger than the measured uncertainty of $0.4\times 10^{12}\text{ M}_{\odot}$ quoted in Diaz et al. (2014). The upper and lower limits of this estimate imply a range of is $0.5\times 10^{12}\text{ M}_{\odot}\leq M^{true}_{VT,LG}\leq 4.1\times 10^{12}\text{ M}_{\odot}$ with 90% confidence. This indicates, when estimating the mass using the VT from observations, the derived value is an overestimate of the true mass, though the correction is small.

The fact that the VT estimator is largely unbiased, i.e. $\langle B\rangle\simeq 1$, may come as a surprise, given the simplicity of this estimator and the assumptions made in deriving it. The two key assumptions in deriving the VT are: 1) that the LG is in dynamical equilibrium, and 2) that the orbits of the OLGMs are statistically isotropic. The former assumption is most likely violated simply because the MW and M31 are falling towards one another for the first time. The application of the VT is further complicated given the fact that a large fraction of LG satellites are bound to either the MW or M31. Regarding the second assumption, this is also manifestly violated in our LG sample since the OLGMs are on largely radial orbits (Fig. 3 right). If the system is assumed to be isotropic, but the underlying orbits are anisotropic, the result would be that the VT overestimates the true mass. This is indeed the trend that we observe in our distribution of $B$, albeit with a substantial scatter.

To further examine this point, in Figure 6 we plot the log of the true mock-LG mass vs the log of the mass derived from the VT. We plot the OLGM dark matter mass sample of $>10^{9}\,\rm{M}_{\odot}$, where the left plot considers OLGM’s with a distance of $D<1.5$ Mpc, and the right plot considers the distance of $350$ kpc $<D<1.5$ Mpc. In each plot the dashed-black line indicates the one-to-one line, where the true mass equals the VT mass. By taking bins of roughly $0.1$ in log of the true mass from 12.2-12.8, we calculate the median, 16% and 84% points of the VT mass in each of these bins. We show the median with the solid black line, and the 16% and 84% points with the gray lines. The color bar indicates the value of $(\sigma_{r}/\sigma_{t})^{2}$ where we have set an upper limit of 4 to cover the range on the majority of the sample. There are three points highlighted in these plots. First, in both plots, we generally see the trend of lighter points on the left and darker points on the right of the dashed-black line. This agrees with our conclusion that the anisotropic mock-LGs generally estimate a larger mass using the virial theorem. Second, the lower limit of $(\sigma_{r}/\sigma_{t})^{2}$ is 0.51 for the left plot and 0.26 for the right plot. This suggests that removing OLGM’s within $350$ kpc from one of the main pairs yields systems with larger tangential velocity and thus more circular orbits around the center of mass of the main pair. When we keep OLGM’s with a distance $<350$ kpc of either the mock-MW or mock-M31, we include satellites that follow the bulk motion of the mock-MW or mock-M31 toward the center of mass resulting in a large radial component. Third, in both plots we see the median (solid-black line) of the VT derived mass tracks the one-to-one (dashed-black) line. This suggests the VT gives a reasonable estimate of the true mass for systems with various true masses.

To construct our LG sample including LMC-like objects, we have considered two different cuts from our primary sample of LG-like systems. Our first cut selects LMCs as those galaxies found within 350 kpc from either of the main galaxies in the LG-like system, and requires a dark matter mass cut on LMC-like systems such that ratio of the LMC-like object relative to the main halo mass is $f_{c}>0.05$. Given this criteria, we consider this as our most “relaxed" cut. For our preferred TA cut, corresponding to $V_{t}<|V_{r}|,\,-195\,{\rm km\,s}^{-1}<V_{r}<-65\,{\rm km\,s}^{-1}$, we find 34 LG analogs with an LMC-like object. We refer to this as sample (a) in our discussion in the remainder of this section. Our second cut imposes a stricter requirement to identify LMCs. First, we choose the MW as the fainter of the two galaxies in the main pair and select LMCs as galaxies found within 350 kpc of the mock MW. In addition for this sample, we require $f_{c}>0.1$. Given these criteria, we consider this as our “strict" cut. Given our preferred TA sample as above, we find 12 LG analogs with an LMC-like object. We refer to this as sample (b) below. Note that for clarity sample (b) is a subset of sample (a).

For both of the above samples, we then solve the timing argument twice: once with $V_{r},V_{t}\,\mathrm{and}\,|\vec{r}|$ calculated when ignoring the LMC presence, and a second time applying the velocity and distance transformations described in Peñarrubia et al. (2016) to account for the LMC presence. As a result, we find that the mean and the 90% containment for the cumulative distributions are similar for both cases. For sample (a) including the effect of the LMS gives a mean of $0.70$ and the 90% containment regions are $[0.41,1.14]$, while for this same sample not accounting from the LMC gives a mean of $0.73$ and the 90% containment regions are $[0.38,1.12]$. For sample (b) when accounting for the mean is $0.70$ and the 90% containment regions are $[0.57,1.04]$, while for this same sample not accounting from the LMC gives a mean of $0.68$ and the 90% containment regions are $[0.54,1.03]$. These ranges are consistent with those reported in Table 2.

We now move on to discuss the impact of the LMC on our VT analysis. For this, we consider the same LMC selection cuts (a) and (b) described above. We apply these cuts to our preferred VT sample of OLGMs with $M_{dm}>10^{9}$ M_⊙ given in bottom panel of Table 3. For our cuts (a) defined above, we find 102 analog LGs, while for our cuts (b) defined above we find 26 analog LGs. Again we calculate $M_{VT}$ ignoring the LMC presence and taking into account the presence of the LMC. Here we take the LMC into account by calculation the barycenter of the MW-LMC-M31 and calculate the radial velocity of each OLGM with respect to this new barycenter. We note that we take the MW center of mass coordinate to be the barycenter of the LMC-MW system. As a result, we find that the mean and the 90% containment for the cumulative distributions are as follows. For sample (a) accounting for the presence of the LMC gives a mean is $0.71$ and the 90% containment regions are $[0.24,1.43]$, while not accounting for the LMC gives a mean is $0.72$ and the 90% containment regions are $[0.20,1.02]$. For sample (b) accounting for the effect of the LMC gives a mean of $0.71$ and the 90% containment regions are $[0.11,1.53]$, while not accounting for the LMC gives a mean of $0.65$ and the 90% containment regions are $[0.14,1.50]$.

We note that, given the relatively small numbers of LGs we find in Illustris, it is unlikely to find a “near exact" LMC-MW-M31 like LG system, indicating how rare the LMC-MW-M31 system is. Nonetheless, it is worthwhile to identify a “best" M31-MW-LMC candidate system from our cuts above, and examine the TA and the VT for this system. Out of our entire LMC samples (a) and (b) constructed above, we find the “best" candidate LG system has an LMC-like object 139 kpc from the fainter galaxy in the main pair, and has $f_{c}=0.26$. For this system we do not find a difference in the TA estimate when calculated with or without the presence of the LMC, with $A=0.68$ in both cases. When we calculate the mass using the VT we find a slight change when including the LMC, from $B=0.85$ to $B_{LMC}=0.73$. In both cases, these values are within 10% of the median value in Table 3, which is calculated for the entire sample of LG-like systems.