A dynamical model of the local cosmic expansion

Abstract

We combine the equations of motion that govern the dynamics of galaxies in the local volume with Bayesian techniques in order to fit orbits to published distances and velocities of galaxies within 3 Mpc. We find a Local Group (LG) mass 2.3 ± 0.7 × 10¹² M_⊙ that is consistent with the combined dynamical masses of M31 and the Milky Way, and a mass ratio |$0.54^{+0.23}_{-0.17}$| that rules out models where our Galaxy is more massive than M31 with ∼95 per cent confidence. The Milky Way's circular velocity at the solar radius is relatively high, 245 ± 23 km s⁻¹, which helps to reconcile the mass derived from the local Hubble flow with the larger value suggested by the ‘timing argument’. Adopting Planck's bounds on Ω_Λ yields a (local) Hubble constant H₀ = 67 ± 5 km s⁻¹ Mpc⁻¹ which is consistent with the value found on cosmological scales. Restricted N-body experiments show that substructures tend to fall on to the LG along the Milky Way–M31 axis, where the quadrupole attraction is maximum. Tests against mock data indicate that neglecting this effect slightly overestimates the LG mass without biasing the rest of model parameters. We also show that both the time dependence of the LG potential and the cosmological constant have little impact on the observed local Hubble flow.

Galaxy: fundamental parameters, Galaxy: kinematics and dynamics, Local Group, cosmological parameters, dark energy, dark matter

1 INTRODUCTION

The dynamics of galaxies in and around the Local Group provide strong support in favour of expansionary cosmological theories. Nowhere else in the Universe can we map in such detail the competition between the primordial expansionary velocities and the gravitational pull of galaxies. To appreciate the beauty of this contest it is worth adopting a Friedmann–Lemaître–Robertson–Walker (FLRW) model wherein the distribution of matter in the Universe is isotropic and homogeneous and gravity behaves as predicted by the theory of general relativity. Under these assumptions the relative motion between two massless particles (say A and G) in a flat Universe can be described by the Friedmann equations (e.g. Peacock 1999):

\begin{equation} \frac{\ddot{r}}{r}=-\frac{4 \pi G}{3}\rho + \frac{\Lambda c^2}{3}, \end{equation}

(1)

where |$r=|{\boldsymbol r}_{\rm A}-{\boldsymbol r}_{\rm G}|$|⁠, ρ is the density of the Universe, c is the speed of light and Λ is the cosmological constant. Defining the Hubble constant H₀ = (8πGρ/3)^1/2 and the fractional vacuum energy density at z = 0 as |$\Omega _{\Lambda }=\Lambda c^2/(3H_0^2)$| it is straightforward to show that the current constraints on these parameters, H₀ ≈ 70 km s⁻¹ Mpc⁻¹ and Ω_Λ ≈ 0.7, lead to a cosmological model in which the particles A and G must move away from each other regardless of the time at which they are observed. However, this is in stark contradiction with the large number of nearby galaxies that show blueshifted spectra. The ingredient missing in equation (1) is, of course, the gravity of the Local Group members.

Kahn & Woltjer (1959) realized that this equation can be modified in order to measure the Local Group mass (M), resulting in what is typically known as the ‘timing argument’. Under the assumption that the mass in the Local Group is approximately made up by our Galaxy and Andromeda (M31) so that M ≈ M_G + M_A, equation (1) can be re-written as

\begin{equation} {\ddot{r}}=-\frac{G M}{r^2} + H_0^2\Omega _{\Lambda }r. \end{equation}

(2)

Dropping the term with the cosmological constant in equation (2) reduces the relative motion between the particles A and G to a Keplerian (radial) orbit. If the age of the Universe is known, the mass M can be directly measured from the current separation and relative velocity between these galaxies. The original value obtained by Kahn & Woltjer (M ≥ 1.8 × 10¹² M_⊙) was considerably larger than the combined masses estimated from rotation curves of both galaxies (∼2 × 10¹¹ M_⊙), suggesting the presence of large amounts of ‘unseen intergalactic matter’, a result which holds true after five decades of active research.¹

Whether or not the cosmological constant plays an active role in the formation of the Local Group remains a matter of debate. For example, the velocity dispersion of nearby (<3 Mpc) galaxies about the Hubble flow is much lower than predicted by cold dark matter (CDM) simulations (Governato et al. 1997; Macciò, Governato & Horellou 2005; Karachentsev et al. 2008). Given the strong suppression of random motions in a fast expanding Universe (Baryshev, Chernin & Teerikorpi 2001; Chernin et al. 2001, 2004), this observation has been interpreted as a manifestation of vacuum energy on local scales (Teerikorpi, Chernin & Baryshev 2005). Such interpretation has been challenged by Hoffman et al. (2008), Peirani (2010) and Martinez-Vaquero et al. (2009), who find a marginal effect of the cosmological constant on local dynamics in N-body simulations with and without vacuum energy. Partridge, Lahav & Hoffman (2013) reach a similar conclusion through an analysis of the timing argument.

The local Hubble flow may be fairly sensitive to the density environment in the vicinity of the Local Group. For example, Aragon-Calvo, Silk & Szalay (2011) find that the velocity fluctuations of galaxies embedded in large-scale structures tend to be smaller than around field galaxies. In this scenario the location of the Milky Way within a vast ‘wall’ of structures connected to the Virgo Cluster (Tully & Fisher 1987) may be responsible for the ‘coldness’ of the local volume. Interestingly, these models also predict an enhanced value for the Hubble constant (≃77–113 km s⁻¹ Mpc⁻¹) on ∼3 Mpc scales. On the other hand, Wojtak et al. (2013) notice that the Hubble constant measured by observers embedded in low-density regions of the Universe tends to be systematically higher than the cosmological value.

Despite the neat theoretical background on which equation (2) rests, there remains a number of open questions regarding the dynamics of galaxies in the Local Group. For example, the mass returned from the timing argument when applied to the relative motion between the Milky Way and M31 (M ∼ 5 × 10¹² M_⊙; Li & White 2008; van der Marel et al. 2012a,b; Partridge et al. 2013, and references therein) is considerably larger than the combined dynamical masses of all galaxies in the Local Group (∼2 × 10¹² M_⊙). This result suggests the striking possibility that approximately half of the Local Group mass is missing. Van der Marel et al. (2012b) have inspected this long-standing issue in detail. First, using Hubble Space Telescope (HST) data these authors find that the transverse velocity of M31 is statistically negligible, supporting the assumption of radial motion in equation (2). In addition, their analysis also reveals a strong dependence between the Local Group mass derived from the timing argument and the azimuthal velocity of the Sun within the Milky Way plane. Unfortunately, both the solar distance to the Galaxy centre and the local circular velocity remain unsettled (e.g. Schönrich 2012, and references therein).

There is also mounting evidence that the orbits of satellite galaxies around the Milky Way and Andromeda deviate from the distribution of substructures found in numerical simulations of structure formation. In particular, both major galaxies in the Local Group are surrounded by vast planes of satellites exhibiting a coherent orbital motion. The disc of satellites around M31 appears to cover ∼400 kpc in diameter, but has a thickness of ∼20 kpc and lies perpendicularly aligned with the axis joining the Milky Way and M31 (Ibata et al. 2013). In the Milky Way most satellites appear distributed over a thicker plane which is inclined by 35° to this axis (Metz, Kroupa & Jerjen 2007; Pawlowski, Kroupa & Jerjen 2013), although both the internal location of the Sun and the lack of deep photometric surveys of the Southern hemisphere introduce severe uncertainties in this measurement. Recently, Tully (2013) has pointed out that the presence of discs of satellites may also extend to galaxies beyond the Local Group (M81 and Centaurus A). Accretion of structures along dark matter filaments (Libeskind et al. 2010; Lovell et al. 2011; Shaya & Tully 2013) and in groups (Lynden-Bell 1982; D'Onghia & Lake 2008; Deason et al. 2011) may be plausible mechanisms to explain the planar distributions of satellites. However, filaments tend to be as extended as the virial radius of the host dark matter halo (Vera-Ciro et al. 2011), whereas disrupted associations quickly thicken unless the progenitor's orbit is perfectly aligned with the long or short axis of a triaxial halo (Bowden, Evans & Belokurov 2013). Indeed, Bahl & Baumgardt (2014) and Ibata et al. (2014) find that satellite planes with the thinness, radial extent and coherent kinematics as the one observed in the Andromeda galaxy are extremely rare, occurring in ≲1 per cent of the halo analogues found in the Millennium II simulations.

The apparent mismatch between ΛCDM predictions and observations has motivated the exploration of alternative gravity theories (e.g. Kroupa et al. 2010). Recently, Zhao et al. (2013) have re-derived the timing argument under the assumption that dynamics on small scales are governed by the empirical modified Newtonian dynamics (MOND) law of Milgrom (1983). These authors show that as a result of the extended MONDian attraction the observed baryonic masses of the Milky Way and M31 imply a past close encounter between the two galaxies between 6 and 12 Gyr ago. Such interaction may explain the presence of extended satellite structures surrounding both galaxies (Pawlowski et al. 2013).

This paper approaches these issues from a ΛCDM framework. First, we inspect the Newtonian equations of motion that describe the perturbations in the Hubble flow by the Milky Way and M31 beyond the point mass approximation (Section 2). Section 3 presents constrained N-body experiments designed to check the assumptions on which the classical timing argument rests. Section 4 outlines the Bayesian methodology used to analyse the observations of the local Hubble flow presented in Section 5. Mock data are generated using the above N-body models in order to quantify the systematic behaviour of the errors identified in our method. Section 6 describes the results of our analysis. A discussion of these results is given in Section 7. Finally, Section 8 summarizes the main findings of this contribution.

2 THE LOCAL HUBBLE FLOW

2.1 Equations of motion

As a result of the combined gravitational pull of the Milky Way and M31, galaxies in the vicinity of the Local Group move on orbits that deviate from the FLRW model. At sufficiently large distances these modifications can be described as gravitational perturbations of the solutions to the Friedmann equation. To analyse the impact of the Local Group on the local Hubble flow it is useful to adopt a coordinate frame whose origin is located at the Local Group barycentre:

\begin{equation} M_{\rm G}{\boldsymbol r}_{\rm G}+ M_{\rm A}{\boldsymbol r}_{\rm A} = {\boldsymbol 0}, \end{equation}

(3)

so that |${\boldsymbol r}_{\rm A}=-M_{\rm G}/M_{\rm A}{\boldsymbol r}_{\rm G}= -f_{\rm m} {\boldsymbol r}_{\rm G}$|⁠, where f_m ≡ M_G/M_A is the Milky Way to M31 mass ratio. The current separation between both galaxies is |$d=|{\boldsymbol r}_{\rm A}-{\boldsymbol r}_{\rm G}|\approx 783\pm 25\,{\rm kpc}$| (McConnachie 2012).

The equations of motion of a massless particle orbiting in the Local Group are

\begin{equation} \ddot{\boldsymbol r} =-\frac{G M_{\rm G}}{({\boldsymbol r}-{\boldsymbol r}_{\rm G})^3} ({\boldsymbol r}-{\boldsymbol r}_{\rm G}) -\frac{G M_{\rm A}}{({\boldsymbol r}-{\boldsymbol r}_{\rm A})^3} ({\boldsymbol r}-{\boldsymbol r}_{\rm A}) + H_0^2\Omega _{\Lambda }{\boldsymbol r}. \end{equation}

(4)

At large distances (r ≫ d), the relative separation between A and G and the tracer galaxy can be approximated as |$|{\boldsymbol r}-{\boldsymbol r}_{\rm G}|\approx |{\boldsymbol r}-{\boldsymbol r}_{\rm A}|\approx r$|⁠. Thus equation (2), which faithfully describes the relative motion between M31 and the Milky Way, can be also used to compute the orbits of distant galaxies with respect to the Local Group barycentre.² Indeed, it was originally pointed out by Lynden-Bell (1981) and Sandage (1986) that the equations on which the ‘timing argument’ rests also govern the orbits of individual galaxies in the outskirts of the Local Group (see also Chernin et al. 2009).

When modelling the orbits of galaxies about the Local Group barycentre we integrate the equations of motion from a time close to the big bang (t_init ≈ 0) to the current epoch t = t₀. The initial conditions r(t_init) = r_ϵ ≪ d, and |$v_{\rm init}\equiv \dot{r}(t_{\rm init})$|⁠, are chosen to match the observed distance and radial velocity of nearby galaxies at t = t₀ (see Section 4 for further details).

Assuming null curvature the age of the Universe (t₀) follows from our choice of the cosmological parameters H₀ and |$\Omega _\Lambda$|⁠:

\begin{eqnarray} t_0=\int _0^{t_0}{\rm d}t=\int _0^1 \frac{{\rm d}\xi }{H \xi }= H_0^{-1}\int _0^1 \frac{{\rm d}\xi }{[(1-\Omega _\Lambda )\xi ^{-3}+\Omega _{\Lambda }]^{1/2} \xi },\nonumber \\ \end{eqnarray}

(5)

where |$H(t)=\dot{\xi }/\xi (t)$|⁠, ξ = (1 + z)⁻¹ and z is the redshift. Here we have used Friedmann's equation |$H(\xi )=H_0[(1-\Omega _\Lambda ) \xi ^{-3}+\Omega _{\Lambda }]^{1/2}$| with no radiation.

For illustrative purposes we shall adopt the following fiducial parameters: h ≡ H₀/(100 km s⁻¹ Mpc⁻¹) = 0.7 and |$\Omega _\Lambda =0.7$|⁠. Plugging these parameters in equation (5) yields t₀ ≃ 13.46 Gyr.

2.2 Can we measure the cosmological constant locally?

The motion of particles implied by equation (4) can be solved analytically if one sets |$\Omega _\Lambda =0$| and models the Local Group as a point mass. Under this approximation the orbits of galaxies reduce to Keplerian orbits, which can be expressed as

\begin{equation} r=a(1-\cos 2 \eta ), \end{equation}

(6)

where a = GM/(−2E) is the semimajor axis of the orbit, E is the orbital energy and η is an angle typically referred to as the eccentric anomaly, which can be calculated numerically from the following equation:

\begin{equation} 2\eta -\sin 2\eta = (GM/a^3)^{1/2}t. \end{equation}

(7)

Using equations (6) and (7), and defining the frequencies Ω = (GM/r³)^1/2 and ω = v/r, Lynden-Bell (1981) showed that in a matter-dominated Universe the current locations and velocities of galaxies around the Local Group follow a close-to-linear relation:

\begin{equation} \Omega t_0 + m \omega t_0 =n; \end{equation}

(8)

with m ≃ 0.9 and n = 2^{− 3/2}π ≃ 1.1. Equation (8) is exact for ω = v = 0, that is η = π/2, which corresponds to the radius r₀ where the expansion of the Universe is momentarily stopped due to the attraction of the Local Group (i.e. the so-called ‘turn-around radius’).

Fig. 1 shows that equation (8) also reproduces the relation between the position and velocity of the Local Group members in a Universe with vacuum energy. The bottom panel indicates that the slope m does not depend strongly on the pressure term of equation (2), while the abscissa n slightly increases from 1.1 to 1.38 depending whether the expansion of the Universe is matter or dark energy dominated.

Figure 1.

Upper panel: relation between the orbital frequency, ω = v/r, and the frequency associated with the dynamical time, Ω = (GM/r³)^1/2, for different values of |$\Omega _\Lambda =1-\Omega _{\rm m}$|⁠. These curves are derived integrating equation (2) numerically for orbits within −3 ≤ ωt₀ ≤ 1. Bottom left-hand panel: linear fits to the curves shown above for a wide range of vacuum energy densities, see equation (8). Note that the slope m ≃ 0.91 is barely sensitive to |$\Omega _\Lambda$|⁠, whereas the abscissa varies as |$n\simeq 1.1 + 0.28 \Omega _\Lambda$| (dotted line). Bottom right-hand panel: as in the left-hand panel with a Local Group mass that varies with time as M(t) = M₀[1 + (t − t₀)ϵ/t₀]. For ease of comparison we set |$\Omega _\Lambda =0$| in equation (2).

Combination of equation (8) with the linear relation |$n=n(\Omega _\Lambda )$| fitted in the lower panel of Fig. 1 provides an analytical description of the perturbed Hubble flow in a ΛCDM Universe:

\begin{equation} v\approx (n-\Omega t_0)\frac{r}{m t_0}\simeq 1.2 \frac{r}{t_0} - 1.1\bigg (\frac{GM}{r}\bigg )^{1/2} + 0.31 \Omega _\Lambda \frac{r}{t_0}. \end{equation}

(9)

It is worth highlighting the small impact of the cosmological constant on the dynamics of Local Group galaxies. From equation (9) the pressure term in the equation of motion leads to a velocity increase |$\delta _\Lambda \equiv v-v(\Omega _\Lambda =0)\simeq 0.31 \Omega _\Lambda r/t_0$| independently of the Local Group mass. For galaxies within r < 3 Mpc and t₀ = 13.46 Gyr this implies |$\delta _\Lambda \lesssim 46 \,{\rm km \, s^{-1}}$|⁠, which is of the same magnitude as the local velocity dispersion about the Hubble flow (see Section 6). The location of the zero-velocity radius is also scarcely sensitive to the pressure term, |$r_0=[GM t_0^2/n^2]^{1/3}\simeq r_0(\Omega _\Lambda =0)(1-0.17\Omega _\Lambda )$|⁠. Given that the observed value is r₀ ≈ 1 Mpc (McConnachie 2012), the cosmological constant shifts the location of r₀ by an amount that is comparable to the error in the distance of several galaxies in our sample (see Section 5). The small (albeit non-negligible) contribution of the cosmological constant to equation (9) suggests that the shape of the local Hubble flow is mainly constrained by the combined masses of the Milky Way and M31 and the age of the Universe, as discussed below.

2.3 Effects of a time-dependent Local Group potential

The analytical derivation of the timing argument given by equation (9) assumes that the mass of the Local Group remains constant throughout the expansion of the local Universe. Given that this is clearly at odds with hierarchical galaxy formation theories we explore below to what extent the local Hubble flow is sensitive to time variations of the Local Group potential.

To this end we use the method of Peñarrubia (2013) for constructing dynamical invariants in time-dependent gravitational potentials. The technique relies on a canonical transformation |${\boldsymbol r}\mapsto {\boldsymbol r}^{\prime }R$|⁠, and a time-coordinate transformation, dt↦dτR², that removes the explicit time dependence from the equations of motion. Here |${\boldsymbol r}^{\prime }(\tau )$| corresponds to orbits calculated in a static potential with M(t) = M₀. For power-law forces F(r, t) = −GM(t)rⁿ, the scale factor is approximately R(t) ≈ [M₀/M(t)]^1/(3+n) (see Peñarrubia 2013, for details).

In the adiabatic limit, i.e. |$(\dot{M}/M_0)^{-1}\gg T$|⁠, where M₀ = M(t₀) the current Local Group mass and T is the radial period of a galaxy, the evolution of the orbital energy is

\begin{eqnarray} E(t)= \frac{E_0}{R^2(t)} + {\cal O}\bigg (\frac{\dot{R}}{R}\bigg ), \end{eqnarray}

(10)

where E₀ = E(t₀) is the energy measured from the current position and velocity vectors of a particle. As expected equation (10) implies no energy variation, regardless of the intermediate evolution of the system, if the final potential is the same as the initial.

The frequencies Ω and ω appearing in equation (8) can be re-written in the generic case where the Local Group evolves adiabatically:

\begin{eqnarray} \Omega (t)\mapsto \frac{1}{R^2}\bigg [\frac{G M_0}{r^{\prime 3}}\bigg ]^{1/2}=\frac{\Omega ^{\prime }}{R^2(t)}, \nonumber \\ \omega (t)\mapsto \frac{1}{r^{\prime }}\frac{{\rm d}r^{\prime }}{{\rm d}\tau }= \frac{1}{R(t)r^{\prime }}\bigg [2\bigg (\frac{E_0}{R^2(t)} + \frac{GM_0}{r^{\prime }R^2(t)}\bigg )\bigg ]^{1/2}=\frac{\omega ^{\prime }}{R^2(t)},\nonumber \\ \end{eqnarray}

(11)

where Ω^′ and ω^′ denote frequencies measured in a fixed potential with M(t) = M₀. Given that at t = t₀ the scale factor is R(t₀) = 1, we find that in the adiabatic limit the measured Hubble flow is independent of the past evolution of the Local Group mass.

To obtain a first-order correction beyond the adiabatic approximation we notice that in the transformed coordinates the age of the Universe is

\begin{equation} \int _0^{\tau _0} {\rm d}\tau =\int _0^{t_0}\frac{{\rm d}t }{R^2(t)}. \end{equation}

(12)

By integrating both sides of equation (12) one can define the time shift Δt as τ₀ = t₀ + Δt.

In the limit Δt/t₀ ≪ 1 the velocity measured at a fixed radius r(t₀) = r′(τ₀) can be written as

\begin{equation} v=v^{\prime }(r^{\prime }[\tau _0])\simeq v^{\prime } (1 -\Omega \Delta t), \end{equation}

(13)

where v^′(r^′) again corresponds to the Hubble relation if the case M(t) = M₀. Equation (13) indicates that the time dependence of the gravitational force changes the radial dependence of ω, while leaving Ω unperturbed. It is straightforward to show that both the slope and abscissa defined in equation (8) vary by the same amount, namely m = m^′(1 − m^′w^′Δt) and n = n^′(1 − m^′w^′Δt).

For illustrative purposes let us consider the case where the Local Group mass varies linearly with time, M(t) = M₀[1 + ϵ(t − t₀)/t₀], with a mass growth rate 0 ≤ ϵ ≤ 1. Recall that for a Keplerian potential n = −2 and the scale factor is R(t) ≈ M₀/M(t), so that the time shift given by equation (12) is simply Δt = −ϵt₀/2. Given that the Local Group mass is expected to grow with time we set ϵ > 0, so that both the slope and abscissa increase by a factor m′(w′t₀)ϵ/2, as shown in the lower right-hand panel of Fig. 1.

As we shall see in Section 5, most galaxies in our sample are moving away from the Local Group barycentre (i.e. ωt₀ ≳ 0). According to Fig. 1 this implies Ωt₀ ≲ 1. Therefore, for a linear mass growth Δv/v = Ωt₀ϵ/2 ≪ 1, suggesting that local Hubble flow holds little information on the time dependence of the Local Group potential. For simplicity the dynamical models discussed below are built in static potentials.

2.4 The Local Group quadrupole

The equations of motion outlined in Section 2.2 assume that the Local Group can be modelled as a central point mass. However, a far more accurate representation of the Local Group corresponds to a gravitational field dominated by the Milky Way and M31 masses. To study the gravitational potential of such a system it is useful to choose a Cartesian system where the galaxies A and G with masses M_A = M/(1 + f_m) and M_G = Mf_m/(1 + f_m) are located at |${\boldsymbol r}_{\rm A}=(-d/[1+f_{\rm m}],0,0)$| and |${\boldsymbol r}_{\rm G}=(+d f_{\rm m}/[1+f_{\rm m}],0,0)$|⁠. In these coordinates the gravitational potential associated with equation (4) can be written as

\begin{eqnarray} \Phi (x,y,z)=-\frac{G M }{(1+f_{\rm m})\lbrace [x+d f_{\rm m}/(1+f_{\rm m})]^2+y^2+z^2\rbrace ^{1/2}} \nonumber\\ -\frac{G Mf_{\rm m} }{(1+f_{\rm m})\lbrace [x-d /(1+f_{\rm m})]^2+y^2+z^2 \rbrace ^{1/2}}+ \frac{1}{2}H_0^2\Omega _{\Lambda }r^2.\nonumber \\ \end{eqnarray}

(14)

Let us now consider galaxies at large distances from the Local Group barycentre, i.e. d/r ≲ 1. At large radii the potential can be approximated as a Taylor expansion of (14):

\begin{eqnarray} \Phi (r,\theta )\approx -\frac{G M }{r} + \frac{1}{2}H_0^2\Omega _{\Lambda }r^2 +\frac{GM f_{\rm m}}{2(1+f_{\rm m})^2}\frac{(1-3 \cos ^2\theta )d^2}{r^3},\nonumber \\ \end{eqnarray}

(15)

where cos θ ≡ x/r. Expressing the total potential as the sum of a spherical plus axisymmetric terms, Φ(r, θ) = Φ⁽⁰⁾(r) + Φ⁽²⁾(r, θ), we find that the right-hand term of equation (15) corresponds to a gravitational quadrupole whose strength decays as Φ⁽²⁾ ∼ 1/r³.

At small distances from the Local Group barycentre, |$r\lesssim [4GM/(H_0^2\Omega _\Lambda )]^{1/3}$|⁠, the contribution of the quadrupole to the potential, Φ⁽²⁾/Φ, does not depend on the total mass of the pair, M. It is, however, fairly sensitive to the mass ratio between the galaxies A and G. For example, the quadrupole term reaches its maximum if the galaxies have equal masses (i.e. f_m = 1), and approaches zero asymptotically if one of the galaxy pair dominates the total mass (i.e. either f_m → 0, or f_m → ∞).

Dashed lines in Fig. 2 follow the iso-potential contours of the quadrupole:

\begin{eqnarray} \Phi ^{(2)}=\Phi ^{(2)}_0 \frac{(1-3 \cos ^2\theta )d^3}{r^3}, \end{eqnarray}

(16)

where |$\Phi ^{(2)}_0 =GM f_{\rm m}/[2d(1+f_{\rm m})^2]$|⁠. Contours in this plot correspond to ρ = 1, 2, 4 and 8 in the equation r(ρ, θ) = d|1 − 3cos ²θ|^1/3ρ^1/3. Notice that the sign of Φ⁽²⁾ flips at |$\theta =\cos ^{-1}(1/\sqrt{3})\approx 54 {^{\circ}_{.}} 73$|⁠. For ease of reference contours are colour coded according to the quadrupole sign (blue/magenta for positive/negative values).

Figure 2.

Potential quadrupole of two equal-mass point masses (red dots) located at |${\boldsymbol r}_{\rm A}/d=(-0.5,0,0)$| and |${\boldsymbol r}_{\rm G}/d=(+0.5,0,0)$|⁠. Iso-potential contours are colour coded according to their sign (magenta/blue for an attractive/repulsive quadrupole). Force lines are plotted with black dotted lines, with arrows marking the direction of the force (see text).

Fig. 2 also shows Faraday's lines of force (black dotted lines). These lines provide a useful representation of the quadrupole, as the number of lines at a given point is related to the strength of the field, whereas the tangent of any curve at a particular point is oriented along the direction of the force (marked with arrows for reference). Each line of force corresponds to a solution to the differential equation rdθ/dr = F_θ/F_r, where F_r = ∂Φ⁽²⁾/∂r and F_θ = r^{− 1}∂Φ⁽²⁾/∂θ. From equation (16) we find that the line of forces follow the equation r/d = ρ′cos ^1/4θsin ^1/2θ.

Note that the gravitational quadrupole induces an attractive force along the axis joining galaxies A and G which becomes repulsive in any transverse direction. Given the azimuthal symmetry of the field it is convenient to write the force in polar coordinates aligned with this axis:

\begin{eqnarray} F_\parallel =F_x=-\frac{\mathrm{\partial} \Phi ^{(2)}}{\mathrm{\partial} x}= \frac{3\Phi ^{(2)}_0d^3}{r^4}(2\cos ^2\theta -3\sin ^2\theta )\cos \theta ,\nonumber\\ F_\perp =\sqrt{F_y^2+F_z^2}=\frac{3\Phi ^{(2)}_0d^3}{r^4}(4\cos ^2\theta -\sin ^2\theta )\sin \theta . \end{eqnarray}

(17)

The force along the symmetry axis is twice as large as in any perpendicular direction and has an opposite sign, i.e. |F_∥(θ = 0)| = −2|F_⊥(θ = π/2)|.

These results suggest that the Local Group quadrupole may have a sustained impact on the orbits of galaxies that define the local Hubble flow, a possibility that we inspect below with the aid of N-body experiments.

3 RESTRICTED N-BODY MODELS

In this section we carry a suite of N-body experiments that follow the expansion of the local (<4 Mpc) Universe from a time close to the big bang to the present. Although these experiments do not capture the complexity of the non-linear growth of structures in a ΛCDM cosmology, they do share essential features with the hierarchical formation of galaxies in an expanding Universe and provide useful insight into the perturbations in the Hubble flow by the Local Group.

The initial conditions are set up so that the final configuration of particles can be approximately described by a FLRW model, where all bodies move away from each other on radial orbits isotropically distributed in space. To this end we place all particles initially at a radius r = r_ϵ. Subsequently, orbital energies are randomly generated within the interval E_min < E_i < E_max, with the range chosen so that at t = t₀ test particles are homogeneously distributed within 0 ≲ r/Mpc ≲ 4. For each particle the velocity associated with E_i is v_{init, i} = [2(E_i + GM/r_ϵ)]^1/2. The directions of the velocity vectors are randomly distributed on the surface of a sphere.

Each combination of (r_ϵ, v_i) defines an orbit which is integrated from t_init,i = r_ϵ/v_init,i to the present, t = t₀. This is done through a leapfrog integration of equation (2), with a time-step chosen so that energy is conserved at a 10⁻³ accuracy level. Because of the central divergence of the Keplerian potential the value of r_ϵ cannot be arbitrarily close to zero. Yet, the choice of r_ϵ should not influence the properties of the Hubble flow at t = t₀. We find that r_ϵ = 0.04 Mpc is sufficiently small so that this condition is met with ease.

Upper panels in Fig. 3 show three snapshots of the expansion of an idealized (local) Universe. Small black dots correspond to massless (‘dust’) particles that move on a Keplerian ‘central’ potential with M = 5 × 10¹² M_⊙ (thick red dot). At early stages all particles move with very high velocities and occupy a densely packed volume. As the Universe expands the mean density of particles decreases in a monotonic fashion. The gravitational pull of the Local Group slows down the motion the particles nearest to the central galaxy, so that eventually a fraction of them reach a turn-around radius and start to fall back on to the central regions of the potential. As expected, our initial conditions lead to a local universe at t = t₀ that resembles a quasi-homogeneous sphere.

Figure 3.

Snapshots of the distribution of particles along an axis perpendicular to the relative motion between the galaxies A and G (red dots). The upper panels adopt a Local Group model where all the mass is located at the barycentre. In the lower panels the Local Group mass is made up by the galaxies A and G, which move on trajectories defined by equation (2). In both cases we adopt M = 2M_A = 2M_G = 5 × 10¹² M_⊙, h = 0.7 and |$\Omega _\Lambda =0.7$|⁠. The separation between A and G at t = t₀ is d ≈ 0.78 Mpc. In both experiments the initial velocity distribution of the tracer particles (black dots) is isotropic. Yet, in models where the potential is dominated by a galaxy pair the infall of particles occurs preferentially along the axis joining both galaxies. Comparison with Fig. 2 shows that the distribution of particles traces the lines of force associated with the potential quadrupole.

The lower panels of Fig. 3 show that a bipolar mass distribution in the Local Group breaks the underlying symmetry built in the initial conditions. To construct this experiment we replace the central point mass by a ‘pair’ of point masses.³ The relative distance between the galaxies A and G evolves according to equation (2). The initial separation corresponds to the initial distance of the tracer particles, i.e. d(t_init) = r_ϵ, and the relative velocity v_init has been chosen so that the final separation is d(t₀) = 0.78 Mpc.

Comparison with Fig. 2 shows that ‘dust’ particles behave in a manner akin to the alignment of iron fillings with a magnetic field, i.e. they distribute along the lines of force defined by the potential quadrupole. Hence, we find tracer particles preferentially along the axis that joins A and G, which is the direction where the gravitational quadrupole of the Local Group is strongest. In contrast, the density of particles drops in transverse directions to the axis defined by the main galaxies, where the quadrupole force has a positive (repulsive) sign. Notice that some of the bodies that fall back towards the Local Group barycentre become bound to either A or G, inducing a strong anisotropy in the spatial distribution of particles around the two main galaxies.

Fig. 4 shows the Hubble flows associated with the snapshots plotted in Fig. 3. By construction the phase-space distribution of particles in a ‘central’ model (upper panels) is an exact solution to equation (2) (cyan dashed curves). Equation (9) can be used to determine the radius at which the Universe expansion halts owing to the gravitational pull of the Local Group, i.e. |$r_0\approx (0.7 GM t_0^2)^{1/3}\simeq 1.42\,{\rm Mpc}$|⁠. In comparison the lower panels exhibit a remarkable contrast. In these models the Local Group mass is made up by the combined masses of the galaxy pair, whose relative distance and velocity is marked with red dots for ease of reference. We find that the bipolar mass distribution induces strong perturbations in the local flow. In particular, the scattered velocity distribution found at small radii, r ≲ d(t₀), results from particles that become bound to either of these galaxies as they fall back towards the Local Group barycentre. In spite of the visible impact of the potential quadrupole on the distribution of tracer particles, comparison with the cyan dashed line shows that equation (2) still provides a reasonable match to the kinematics of galaxies at r ≳ d(t₀). It is worth bearing in mind, however, that this equation systematically underestimates the radial velocity of galaxies at intermediate distances, d(t₀) ≲ r ≲ r₀. This results from the strong spatial anisotropy of tests particles, which tend to be found along the bipolar direction where the gravitational attraction is enhanced.

Figure 4.

Hubble flows associated with the models shown in Fig. 3. Red dots mark the relative distance and velocity between the galaxies A and G at each snapshot. Their separation at z = 0 is d(t₀) = 0.78 Mpc. By construction the Hubble flow of the ‘central’ Local Group model (upper panel) is an exact solution to equation (2) (cyan dashed curves). In contrast the lower panels show that the potential quadrupole causes strong perturbations on the kinematics of nearby galaxies. Some of the particles that fall back towards the Local Group barycentre become bound to the galaxies A and G, which leads to a large velocity scatter at r ≲ d(t₀). Note also that at intermediate distances (1 ≲ r/d(t₀) ≲ 3) the stronger gravitational pull along the axis joining A and G tends to increase the infall velocity of particles with respect to solutions to equation (2).

4 BAYESIAN ANALYSIS

In this section we describe the fundamentals of our Bayesian analysis of the local Hubble flow and perform a number of tests using the N-body experiments outlined in Section 3. Our models contain six free parameters that are fitted simultaneously to the data. These are the Local Group mass, M = M_G + M_A, the mass ratio between the Milky Way and Andromeda, f_m = M_G/M_A, the circular velocity of the Milky Way at the solar radius, V₀ = V_c(R_⊙), the reduced Hubble constant, h, and the fractional vacuum energy density, |$\Omega _\Lambda$|⁠.

Our choice of h and |$\Omega _\Lambda$| as independent quantities as motivated by recent papers which show that the value of the Hubble constant is sensitive to environment (Aragon-Calvo et al. 2011; Wojtak et al. 2013), whereas the age of the Universe and the cosmological constant are not. This approach deviates from the standard ‘timing argument’ described in Section 2, which typically assumes that the age of the Universe is a known quantity. In theory, whether we consider |$\Omega _\Lambda$| and h, or |$\Omega _\Lambda h^2$| and t₀, as free parameters is a subjective decision which should not have a measurable impact on our fits. In practice, the choice of priors can in some cases modify the posterior distributions. Following the suggestion of the anonymous referee we have explicitly checked that the bounds derived in Section 6 and summarized in Table 3 are independent of the combination of free cosmological parameters.

Table 1.

Boundaries of uniform priors for free parameters in our likelihood function.

Parameter	Min.	Max.	Description
M/(1012 M⊙)	0.1	20.0	M = MG + MA
log10fm	−1.0	1.0	fm = MG/MA
V0/ km s−1	100	400	MW's circular velocity at R⊙
|$\Omega _\Lambda$|	0.0	1.0	Fractional vacuum density
h	0.2	1.0	Reduced Hubble constant
σm/ km s− 1	1.0	200	Hyperparameter (no phys. meaning)

Table 2.

Heliocentric coordinates and radial velocities of dwarf galaxies within a radial range 0.8 ≤ r/Mpc ≤ 3 from the Local Group barycentre and with separations to IC 342, M81 and Centaurus A larger than 1 Mpc (see text). Data taken from (1) McConnachie (2012), (2) Karachentsev et al. (2002b), (3) Kirby, Cohen & Bellazzini (2012), (4) Giovanelli et al. (2013) and (5) McQuinn et al. (2013). Relative distances to the major associations in the vicinity of the Local Group are calculated using a fiducial f_m = 1 and the following heliocentric positions: (D, l, b)_IC 348 = (3.3 Mpc, 138| $_{.}^{\circ}$|17, +10| $_{.}^{\circ}$|58); (D, l, b)_M81 = (3.6 Mpc, 142| $_{.}^{\circ}$|09, +40| $_{.}^{\circ}$|90) and (D, l, b)_Cen A = (3.8 Mpc, 309| $_{.}^{\circ}$|52, +19| $_{.}^{\circ}$|42).

Name	l	b	D	ϵD	Vh	ϵV	dIC 342	dM81	dCen A	Ref.
	(°)	(°)	(Mpc)	(Mpc)	( km s−1)	( km s−1)	(Mpc)	(Mpc)	(Mpc)
Andromeda (M31)	121.2	−21.6	0.783	0.025	−300.0	4.0	2.71	3.35	4.57	(1)
Leo A	196.9	52.4	0.798	0.044	22.3	2.9	3.02	3.00	3.85	(1)
Tucana	322.9	−47.4	0.887	0.049	194.0	4.3	4.04	4.47	3.56	(1)
WLM	75.9	−73.6	0.933	0.034	−130.0	1.0	3.47	4.18	4.32	(1)
Sagittarius dIrr	21.1	−16.3	1.067	0.088	−78.5	1.0	3.92	4.29	3.74	(1)
Aquarius (DDO 210)	34.0	−31.3	1.072	0.039	−140.7	2.5	3.76	4.47	3.56	(1)
NGC 3109	262.1	23.1	1.300	0.048	403.0	2.0	4.03	3.94	3.00	(1)
Antlia	263.1	22.3	1.349	0.062	362.0	2.0	4.09	3.98	2.97	(1)
Andromeda XVIII	113.9	−16.9	1.355	0.081	−326.2	2.7	2.35	3.23	5.12	(1)
UGC 4879	164.7	42.9	1.361	0.025	−29.1	1.3	2.41	2.33	4.44	(1), (3)
Sextans B	233.2	43.8	1.426	0.020	304.0	1.0	3.51	3.23	3.50	(1)
Sextans A	246.1	39.9	1.432	0.053	324.0	2.0	3.75	3.48	3.26	(1)
HIZSS 3[A]	217.7	0.1	1.675	0.108	288.0	2.5	3.42	3.67	4.19	(1)
HIZSS 3[B]	217.7	0.1	1.675	0.108	322.6	1.4	3.42	3.67	4.19	(1)
Leo P	219.6	54.4	1.720	0.400	264.0	2.0	3.34	2.85	3.72	(4), (5)
KKR 25	83.9	44.4	1.905	0.061	−139.5	1.0	2.78	2.53	4.64	(1)
NGC 55	332.9	−75.7	1.932	0.107	129.0	2.0	4.46	5.30	4.45	(1)
IC 5152	343.9	−50.2	1.950	0.045	122.0	2.0	4.88	5.49	3.83	(1)
ESO 294−G010	320.4	−74.4	2.032	0.037	117.0	5.0	4.57	5.41	4.44	(1)
NGC 300	299.2	−79.4	2.080	0.057	146.0	2.0	4.48	5.37	4.61	(1)
GR 8	310.7	77.0	2.178	0.120	213.9	2.5	4.03	3.21	3.21	(1)
KKR 3 (KK 230)	63.7	72.0	2.188	0.121	63.3	1.8	3.44	2.69	3.98	(1)
UKS 2323−326 (UGCA 438)	11.9	−70.9	2.208	0.092	62.0	5.0	4.57	5.47	4.69	(1)
IC 3104	301.4	−17.0	2.270	0.188	429.0	4.0	5.51	5.68	2.42	(1)
UGC 9128 (DDO 187)	25.6	70.5	2.291	0.042	152.0	1.0	3.93	3.15	3.59	(1)
IC 4662	328.5	−17.8	2.443	0.191	302.0	3.0	5.71	5.92	2.56	(1)
KKH 98	109.1	−22.4	2.523	0.105	−136.9	1.0	2.28	3.63	6.24	(1)
DDO 125	137.8	72.9	2.582	0.059	194.9	0.2	3.11	1.97	4.50	(1)
UGC 8508	111.1	61.3	2.582	0.036	56.0	5.0	2.77	1.78	4.88	(1)
KKH 86	339.0	62.6	2.582	0.190	287.2	0.7	4.69	3.87	2.81	(1)
DDO 99	166.2	72.7	2.594	0.167	251.0	4.0	3.20	2.05	4.40	(1)
DDO 190	82.0	64.5	2.793	0.039	150.0	4.0	3.36	2.37	4.66	(1)
NGC 4163	163.2	77.7	2.858	0.039	165.0	5.0	3.48	2.21	4.38	(1)
NGC 404	127.1	−27.0	3.060	0.370	−48.0	9.0	2.14	3.81	6.83	(2)

Table 3.

Constraints on model parameters using flat priors (Table 1), and a Gaussian prior on the fractional vacuum energy density based on Planck data (⁠|$\Omega _\Lambda =0.686\pm 0.020$|⁠). Error bars enclose the central 68 per cent (95 per cent) of area under the marginalized 1D posterior probability distribution functions shown in Fig. 11.

Model parameters	Flat priors	Planck prior on |$\Omega _\Lambda$|
M/(1012 M⊙)	|$2.3_{-0.7(-1.2)}^{+0.7(+1.7)}$|	|$2.3_{-0.7(-1.2)}^{+0.7(+1.7)}$|
fm	|$0.54_{-0.17(-0.30)}^{+0.23(+0.60)}$|	|$0.54_{-0.16(-0.30)}^{+0.24(+0.60)}$|
V0/km s−1	|$245_{-23(-45)}^{+23(+47)}$|	|$245_{-23(-45)}^{+23(+51)}$|
σm/km s− 1	|$35_{-4(-8)}^{+5(+11)}$|	|$35_{-4(-8)}^{+6(+12)}$|
|$\Omega _\Lambda$|	|$0.54_{-0.35(-0.51)}^{+0.32(+0.44)}$|	|$0.69_{-0.02(-0.04)}^{+0.02(+0.04)}$|
h	|$0.64_{-0.07(-0.12)}^{+0.10(+0.20)}$|	|$0.67_{-0.04(-0.09)}^{+0.04(+0.09)}$|
Derived quantities
|v⊙|/km s− 1	|$312_{-11(-22)}^{+11(+23)}$|	|$313_{-11(-22)}^{+11(+23)}$|
l⊙ (°)	|$93.9_{-1.6(-3.0)}^{+1.8(+3.6)}$|	|$93.9_{-1.5(-2.9)}^{+1.7(+3.5)}$|
b⊙ (°)	|$-3.2_{-1.3(-2.7)}^{+1.2(+2.3)}$|	|$-3.2_{-1.3(-2.7)}^{+1.2(+2.2)}$|
MG/(1012 M⊙)	|$0.8_{-0.3(-0.5)}^{+0.4(+0.9)}$|	|$0.8_{-0.3(-0.5)}^{+0.4(+0.9)}$|
MA/(1012 M⊙)	|$1.5_{-0.4(-0.8)}^{+0.5(+1.2)}$|	|$1.5_{-0.4(-0.8)}^{+0.5(+1.1)}$|
t0/Gyr	|$13.2_{-1.4(-2.4)}^{+2.9(+8.5)}$|	|$13.8_{-0.8(-1.7)}^{+1.0(+2.1)}$|
σH/km s−1	|$50_{-4(-8)}^{+4(+8)}$|	|$50_{-4(-8)}^{+4(+8)}$|

4.1 Likelihood function

Consider the Gaussian likelihood function:

\begin{eqnarray} &&\mathcal {L}(\lbrace D_i,l_i,b_i,V_{{\rm h},i}\rbrace ^{N_{\rm sample}}_{i=1}|\boldsymbol {S})\nonumber\\ &&\quad=\prod _{i=1}^{N_{\rm sample}}\frac{1}{\sqrt{2\pi \sigma _i^2}}\exp \bigg [-\frac{(V_i-V_{{\rm h},i})^2}{2\sigma _i^2}\bigg], \end{eqnarray}

(18)

where |$\boldsymbol {S}=(M,f_{\rm m},V_0,h,\Omega _\Lambda ,\sigma _{\rm m})$| is a vector that comprises the model parameters; D and V_h are heliocentric distances and velocities, respectively, and (l, b) the Galactocentric coordinates.

For a given set of parameters, |$\boldsymbol {S}$|⁠, our model returns a heliocentric velocity V at the location (D, l, b). The velocity of the ith galaxy, V_i, is calculated according the following procedure: (i) first, the distance of a galaxy to the Local Group barycentre, r_i(t₀), is derived using the coordinate transformation of Appendix A. Note that this conversion is model dependent, as it requires setting the values of f_m and V₀. (ii) Subsequently, we choose an initial (small) radius r_ϵ and solve for the initial velocity v_init(t = t_init), with t_init = r_ϵ/v_init ≪ t₀, so that r(t₀) matches the value obtained in step (i). This is done through a leapfrog integration of equation (2) from t = t_init until t = t₀, with a time-step chosen so that energy is conserved at a 10⁻³ accuracy level (see Section 3). The model parameters that enter in this step are M, h and |$\Omega _\Lambda$|⁠. (iii) To derive the value of V_i that goes in the likelihood function we convert the Local Group-centric coordinates of the galaxy into heliocentric ones using the transformation outlined in Appendix A.

4.2 Peculiar motions

Equation (2) is built upon the assumption that all galaxies in the Local Group move on radial orbits. Deviations from this assumption contribute to the presence of a peculiar motion with respect to the Hubble flow.

It is useful to outline the main processes that may contribute to peculiar motions. First, one should consider statistical fluctuations in the data, such as the presence of observational errors and systematic biases, which typically arise from the complexities involved in the measurement of distances and velocities for such distant objects, as recently shown by Fraternali et al. (2009). Second, the presence of galaxy associations in our tracer population, which tend to be relatively common in the local volume (e.g. Tully et al. 2006; Bellazzini et al. 2013; Chapman et al. 2013; Fattahi et al. 2013), is a non-negligible source of peculiar motions. Galaxies belonging to an association move around a common barycentre, which tends to inflate the distance–velocity relation derived purely from the cosmological expansion. Also, nearby galaxy groups may perturb the motion of kinematic tracers in the periphery of the Local Group.

Unfortunately, it is difficult to determine on observational grounds which objects have been acted on by external tidal fields or by encounters with neighbour systems. Given this strong limitation, here we follow a statistical approach in order to incorporate peculiar motions in our fits. Our method relies on the assumption that the above perturbations are small and do not lead to systematic biases in the distance–velocity relation.

To account for the effect of peculiar motions we incorporate the hyperparameter σ_m in our analysis. Its squared value is added linearly to the observational variance in order to minimize covariance with the rest of model parameters. Hyperparameters provide a useful tool to assign weights to data sets beyond those derived from statistical errors (e.g. Hobson, Bridle & Lahav 2002). Indeed, by introducing σ_m in our likelihood function we allow that the observed scatter in the distance–velocity relation may not be fully accounted by random errors in the data set. The role of σ_m is thus analogous to a ‘nuisance’ parameter. As such, we marginalize over σ_m in order to obtain a joint estimation of the parameters of interest.

The two-dimensional variance of the ith measurement is calculated as follows:

\begin{equation} \sigma _i^2=\epsilon _{V,i}^2 + \epsilon _{D,i}^2 \bigg (\frac{{\rm d}V}{{\rm d}D}\bigg )^2_{(D_i,l_i,b_i)} + \sigma _m^2, \end{equation}

(19)

which includes velocity (ϵ_V) as well as distance (ϵ_D) errors (see Ma, Hinshaw & Scott 2013), plus the additional ‘freedom’ provided by the hyperparameter σ_m. Here dV/dD denotes the gradient of the Hubble flow at a given location (D, l, b) returned by our model.

We apply a nested sampling technique (Skilling 2004) in order to calculate posterior distributions for our parameters and the evidence of the model. In particular we use the code MultiNest, a Bayesian inference tool which also produces posterior samplings and returns error estimates of the evidence (see Feroz & Hobson 2008; Feroz, Hobson & Bridges 2009, for details). Unless otherwise indicated, all our measurements adopt flat priors over ranges that include reasonable parameter values (see Table 1).

4.3 Mock data

Prior to applying the technique outlined in the previous section to actual data, we must examine the reliability of our method when it operates on realistic data sets that violate our model assumptions. With this aim in mind we use the models outlined in Section 3 to construct synthetic data sets which shall help us to detect possible biases in the statistical inference of the model parameters. Below we describe how mock data sets are generated and the systematic behaviour of the errors identified in our method.

4.3.1 The effect of peculiar velocities

We consider two sources of peculiar velocities in our tests: those arising from observational errors, which are assumed to be Gaussian, plus a transverse component, |$v_{\rm t}=\sqrt{v_\theta ^2 + v_\phi ^2}$|⁠, where ϕ and θ are, respectively, the azimuthal and polar angles in spherical coordinates, which we add to the orbital velocity vector at a random direction. We explore models with transverse velocity components comparable to the radial one, v_t = 50 and 100 km s⁻¹. Subsequently, from the position and velocity vectors in a Local Group-centric frame we derive heliocentric distances and velocities (see Appendix A), which are then convolved with Gaussian errors. In particular, we set ϵ_D = 50 kpc and ϵ_V = 5 km s⁻¹, which are broadly consistent with the observational errors in the sample of galaxies gathered from the literature (see Section 5).

Fig. 5 shows that the main effect of peculiar velocities is to thicken the distance–velocity relation derived from the coordinate transformation of Appendix A, which rests upon the assumption that all galaxies in and around the Local Group move on radial orbits. This is because the contribution of peculiar motions to the heliocentric velocity corresponds to the projection of the transverse component on to the line-of-sight vector. The fact that our galaxies are distributed over a very large area of the sky leads to a large range of projection angles, and thus to a ‘scattered’ contribution of peculiar motions to the bulk flow. Owing to a well-known geometrical effect the apparent scatter wanes progressively as the distance to the barycentre increases.

Figure 5.

Mock Hubble diagrams generated from the models shown in Figs 3 and 4. Distances are given in units of the zero-velocity radius r₀ ≈ 1.42 Mpc. Black dotted–dashed and cyan dashed lines show isochrones derived from equation (2) for M = 0 (unperturbed Hubble flow) and M = 5 × 10¹² M_⊙, respectively, in a Universe with |$\Omega _\Lambda =0.7$| and h = 0.7. Dots correspond to estimates of the radial velocity in a Local Group-centric frame (see Appendix A) after introducing Gaussian errors in heliocentric distance (ϵ_D = 50 kpc), velocity (ϵ_V = 5 km s⁻¹), as well as a randomly oriented tangential velocity component with moduli v_t = 50 (blue open dots) and 100 km s⁻¹ (red filled dots). Note that the presence of peculiar velocities plus observational errors tends to increase the scatter in the observed Hubble flow. The derivation of radial velocities becomes progressively less accurate as the distance to the barycentre decreases. The large velocity scatter in the lower panel arises from particles bound to the main galaxies. Mock data are constructed with particles within the distance range marked with green arrows.

Fig. 6 illustrates how the hyperparameter σ_m helps to lessen the impact of peculiar motions on our analysis. Here we plot posterior samplings for 10 independent mock data sets, each comprising 35 galaxies within a radial range 1.2 ≤ r/Mpc ≤ 3.2 and a transverse velocity v_t = 100 km s⁻¹. In units of the zero-velocity radius the radial range corresponds to 0.85 ≤ r/r₀ ≤ 2.25. The left-hand panel shows that setting the parameter σ_m = 0 leads to uncertainties in M and f_m that are unrealistically small, as many of the probability clouds appear isolated and do not overlap with the true values. In contrast, fitting σ_m and marginalizing over it yields uncertainties that do comprise the true values and can be therefore deemed more realistic.

Figure 6.

Sampling of the posterior distributions of M and f_m calculated for 10 mock data sets built from the ’central’ model, each containing 35 galaxies located within a radial range 0.85 ≤ r/r₀ ≤ 2.25 (see Fig. 5). We consider errors in the heliocentric distances (ϵ_D = 50 kpc) and velocities (ϵ_V = 5 km s⁻¹), plus a randomly oriented tangential velocity component with a magnitude v_t = 100 km s⁻¹. Red dots mark true parameter values. Comparison between the left- and right-hand panels illustrates the effects of introducing the hyperparameter σ_m in the variance. Allowing the presence of a dispersion in our models yields more realistic uncertainties in the fitted parameters.

4.3.2 Milky Way–M31 quadrupole

Thus far we have assumed that the equations that define the motion of nearby galaxies around the Local Group are those on which the timing argument rests, i.e. equation (2). Recall, however, that this approximation is only accurate for distant galaxies, i.e. in the limit |$r\gg d\equiv |{\boldsymbol r}_{\rm A}-{\boldsymbol r}_{\rm G}|$|⁠. Unfortunately, we shall see in Section 5 that the majority of galaxies in the galaxy sample do not obey this condition, as their distances from the Local Group barycentre tend to be comparable to the current separation between the Milky Way and Andromeda. Following the results of the restricted N-body experiments outlined in Section 3, which suggest that a bipolar mass distribution in the Local Group may introduce visible effects on the kinematics of nearby galaxies, we explicitly explore here whether the point mass approximation may be a source of systematic biases in our analysis.

The lower panel of Fig. 5 shows mock Hubble diagrams associated with the ‘pair’ models for different values of v_t. The first point to notice is the large velocity scatter at r ≲ d(t₀), which results from the particle population bound to either of the two main galaxies. Second, we also find that the accelerated infall velocity of accreting particles along the axis joining the main galaxies leads to systematic deviations from the ‘timing–argument’ relation (cyan dashed line). In particular, equation (2) tends to underestimate the magnitude of infall velocities at intermediate distances d(t₀) ≲ r ≲ r₀, where r₀ ≈ 1.42 Mpc is the zero-velocity radius.

Fig. 7 shows the error distributions, E(x) = (x − x_true)/x_true, obtained from fitting 50 mock data sets generated from the models shown in Fig. 5. Each data set contains 35 galaxies within a distance range 0.85 ≤ r/r₀ ≤ 2.25. Focusing first on the ‘central’ model, we find that increasing the magnitude of the tangential velocity component tends to augment the uncertainty in M, f_m and V₀, without leading to any visible bias. In contrast, our analysis returns slightly overestimated masses when applied to the ‘pair’ models. This bias increases with the magnitude of the peculiar velocity component. For the range of transverse motions explored here, we find that M can be overestimated by 20–30 per cent. Remarkably, the parameters f_m and V₀ are measured with high precision in both sets of models. Such tight constraints result from the strong sensitivity of the solar apex to the parameter f_m and V₀, as discussed in Appendix B. Fig. 8 shows posterior samplings for the cosmological parameters |$\Omega _\Lambda$| and h calculated from the 50 mock data sets shown in Fig. 7. This plot includes a few results of note. First, the parameters |$\Omega _\Lambda$| and h exhibit a strong degeneracy which roughly runs parallel to the isochrone |$t_0=t_0(\Omega _\Lambda , h)=13.46$| Gyr (dashed line) derived from equation (5) in a flat Universe. Yet, the degenerate solutions are still consistent with the true values (red dots) for both ‘central’ and ‘pair’ N-body models. Note also that the potential quadrupole does not introduce a bias on the cosmological parameters in spite of the relatively large magnitude of the transverse velocity component added to the velocity vector of all test particles (v_t = 100 km s^{− 1}). Indeed, the most distant objects in the galaxy sample considered here (r ∼ 3 Mpc) put the strongest constraints on |$\Omega _\Lambda$| and h. As shown in Fig. 5, the effects of (unknown) peculiar velocities tend to be negligible on these scales.

Figure 7.

Errors in the inference of the model parameters M, f_m and V₀, where E(x) ≡ (x − x_true)/x_true and M_true = 5 × 10¹² M_⊙, f_{m, true} = 1 and V_{0, true} = 220 km s⁻¹ derived from fits of 50 independent sets of 35 mock galaxies within the distance interval 0.85 ≤ r/r₀ ≤ 2.25 (see Fig. 5). Note that although peculiar velocities tend to increase the uncertainty of the fits, they do not lead to significant biases in the inferred parameters. In contrast, neglecting the contribution of the potential quadrupole in equation (2) tends to overestimate the Local Group mass by 20–30 per cent.

Figure 8.

Sampling of the posterior distributions of |$\Omega _\Lambda$| and h derived from 50 mock data sets with v_t = 100 km s⁻¹ (see Fig. 7). Red dots mark the true parameter values |$h=\Omega _\Lambda =0.7$|⁠. Note that the strong degeneracy between |$\Omega _\Lambda$| and h roughly traces the isochrone curve |$t_0=t_0(h,\Omega _\Lambda )=13.46$| Gyr derived in a Universe with null curvature (cyan dashed line).

As expected from the analytical estimates obtained in Section 2.2, we find that the constraints on the cosmological constant are considerably weaker than those on h. These tests suggest that the kinematics of nearby (r ≲ 3 Mpc) galaxies can be used to put meaningful bounds on the Hubble constant only if prior bounds on |$\Omega _\Lambda$| are incorporated in the analysis. We shall return to this point in Section 6.

5 OBSERVATIONAL DATA

Comparison between model and observations is done through the likelihood function built in Section 4, which incorporates measurements of distances and radial velocities of individual galaxies. Recently, two major catalogues of nearby galaxies have been made publicly available: McConnachie (2012), which provides a detailed description of the properties of dwarfs at heliocentric distances D ≤ 3 Mpc, and Karachentsev, Makarov & Kaisina (2013) who have compiled data for 869 galaxies within 11 Mpc from the Sun. Both data sets provide equivalent information within the range of overlap.

Fig. 9 shows the spatial location of galaxies within a 5 Mpc volume. For simplicity we have aligned the coordinate system so that the axis joining the Milky Way and M31 is the X-axis. Note that the local Universe is considerably more rich in substructures than the restricted N-body experiments built in Section 3. Clearly, the idealized initial conditions of our N-body models, which are based on a homogeneous and isotropic Universe, cannot reproduce the complex spectrum of overdensities that drive the formation of structures in the local volume. In particular, three prominent associations (or clusters) of galaxies stand out at r ≳ 3 Mpc (outer green dotted line): IC 342/Maffei-I (Karachentsev et al. 2003), M81 (Karachentsev et al. 2002a; Chiboucas et al. 2013) and Centaurus A/M83 (Karachentsev et al. 2007), which we mark with blue asterisks, squares and triangles, respectively. The middle and right-hand panels show glimpses of an even higher level of the hierarchical galaxy formation ladder, as many of the visible structures lie on a vast plane that connects to the Virgo cluster and its filamentary network (Tully & Fisher 1987).

Figure 9.

Spatial distribution of galaxies within 5 Mpc. The Milky Way and M31 are marked with red dots. We use Cartesian coordinates where the axis joining the Milky Way and M31 is aligned with the X-axis. Dashed lines in the left-hand panel mark 0.75 and 3.0 Mpc radii from the Local Group barycentre, which is derived under the assumption that both galaxies have equal masses. For ease of reference, we highlight the main associations/clusters around the Local Group, namely IC 342 (blue asterisk), M81 (blue square) and Centaurus A (blue triangle). Fill dots denote galaxies within the dashed lines whose relative distance to any of the three major associations is larger than d_min = 1 Mpc (see text).

We identify the Milky Way and M31 (red dots) as well as the overdensities that surround them. These correspond to gravitationally bound ‘satellite’ galaxies, which tend to be located at barycentric distances r ∼ d (inner green dotted line). Galaxies that are gravitationally bound to our Galaxy, Andromeda or any of the external associations/clusters move on orbits that strongly deviate from the model assumptions on which the timing argument rests. To illustrate this point we plot in Fig. 10 the distance and radial velocity with respect to the Local Group barycentre of the galaxies shown in Fig. 9. Notice the large velocity scatter shown by satellites in the neighbourhood of the main galaxies (see Section 4). In particular the velocity dispersion about the bulk flow increases noticeably at r ≲ d ≈ 0.78 Mpc and r ≳ 3 Mpc.

Figure 10.

Phase-space location of the galaxies plotted in Fig. 9. Black solid dots denote galaxies incorporated in our Bayesian analysis. A red dot marks the current separation and relative velocity between the Milky Way and M31 using V₀ = 220 km s⁻¹. For reference we also plot two isochrone lines for M = 5 × 10¹¹ and 5 × 10¹² M_⊙. Note that the zero-velocity radius of the Local Group is located at r₀ ∼ 1 Mpc. Regions of large velocity scatter point towards the presence of a large number of satellite galaxies.

The impact of external perturbers such as IC 342, M81 and Centaurus A on our fits can be strongly suppressed by excluding galaxies in the vicinity of those systems.⁴ Accordingly, our data set only include galaxies whose distance to any of the three major associations is larger than a given d_min. In addition, we also impose a distance cut, r_max, in our selection criteria, which is motivated by the decreasing accuracy of the timing argument at distances where the contribution of the Local Group to the local gravitational field becomes negligible.

Black solid dots in Figs 9 and 10 show galaxies that obey the above criteria for d_min = 1 Mpc and r_max = 3 Mpc (see also Table 2). Although these hard cuts are put ad hoc, we have checked that any combination within the range 1.0 ≤ d_min/Mpc ≤ 2.0 and 2.7 ≤ r_max/Mpc ≤ 3.0 yields a similar fit to the model parameters, indicating that the results discussed in Section 6 are not overly sensitive to the choice of d_min and r_max.

In spite of the small number of dwarfs that match the above conditions (N_sample ∼ 30), the analysis of mock data in Section 4 suggests that the set size may be sufficiently large to provide meaningful constraints on the individual masses of the Milky Way and M31, the circular velocity of the Milky Way at the solar radius, as well as on the cosmological parameters |$\Omega _\Lambda$| and h. We explore this possibility below.

6 RESULTS

We now apply our method to galaxies located within the radial range 0.8 Mpc ≤ r ≤ r_max and a separation to IC 342, M81 and Centaurus A larger than d_min = 1 Mpc. The largest number of galaxies in the sample is N_sample = 30 for r_max = 3 Mpc (see Table 2). Fig. 11 displays posterior distributions for each model parameter as returned by MultiNest. Table 3 lists for each parameter the median value from the posterior probability distribution function (PDF), with error bars indicating the interval that encloses the central 68 and 95 per cent of values.

Figure 11.

Posterior distributions for our model parameters. Fits include galaxies within the radial range 0.8 ≤ r/Mpc ≤ 3.0 and a separation to IC 342, M81 and Centaurus A larger than d_min = 1 Mpc (see Table 2).

6.1 Solar apex

A correct estimation of the motion of the Sun relative to the other members of the Local Group is a key aspect of our study. The conversion between helio- and Local Group-centric coordinates (see Appendix A) requires an understanding of the relative motion between the Sun and the Milky Way, as well as between the Milky Way and the Local Group barycentre. Both steps introduce a non-negligible degree of uncertainty in our models.

Although the motion of the Sun has been studied for many decades, this is the first attempt to measure the location of its apex by modelling the kinematics of individual Local Group members while simultaneously dealing with uncertainties in the circular velocity of the Sun as well as on the relative mass between the Milky Way and M31.

It is thus reassuring that the coordinates listed in Table 3 agree very well with previous measurements. For example, Karachentsev & Makarov (1996) minimize scatter in the distribution of radial velocities with respect to the bulk flow, v = Hr, where H is a free parameter related to the (local) expansion of the Universe (see also Appendix B). These authors find (v, l, b)_⊙ = (316 ± 5 km s^{− 1}, 93° ± 2°, −4° ± 2°), which is consistent with our measurement within 1σ uncertainties. Similarly, Courteau & van den Bergh (1999) assume that the velocity distribution of nearby galaxies with respect to the Hubble flow is Maxwellian. These authors argue that if the velocity distribution is invariant under spatial translations, the true solar motion is the one that minimizes the velocity dispersion of their models.⁵ Their result (v, l, b)_⊙ = (306 ± 18 km s^{− 1}, 99° ± 5°, −3° ± 4°) also agrees with our measurement at 1σ confidence level.

6.2 The Sun's transverse motion

The circular velocity of the Milky Way at the solar radius, V₀ ≡ V_c(R_⊙), is a crucial parameter for the derivation of the solar apex. Unfortunately, the wide range of values reported in the literature (e.g. Bhattacharjee, Chaudhury & Kundu 2014) cannot be explained by the quoted errors of individual measurements, which may be indicative of systematic biases in some of the published methods. Given that a wrong choice of V₀ propagates through our whole analysis, here we opt for not imposing an external constraint on its value. In a Bayesian framework this is equivalent to adopting a ‘diffuse’ or ‘uninformative’ prior (see Table 1). Note that by fitting V₀ simultaneously with the rest of parameters we are effectively incorporating the uncertainty in the value of V₀ into the joint posterior distributions of all measured quantities.

Fig. 11 shows that the kinematics of nearby galaxies can be used to put meaningful constraints on the value of V₀ if we adopt prior information on the motion of the local standard of rest (LSR; see Appendix A). Indeed, the small covariance of V₀ with the rest of parameters is at the core of the relative narrowness of the confidence intervals given in Table 3.

Adopting the LSR velocity vector measured by Schönrich, Binney & Dehnen (2010), the median value returned from our Bayesian fits is V₀ = 245 ± 23 km s⁻¹ at a 68 per cent confidence level. Although these bounds lie above the value of 220 km s⁻¹ adopted by the IAU (Kerr & Lynden-Bell 1986), our measurement appears in excellent agreement with the recent estimates of McMillan (2011), who finds V₀ = 239 ± 5 km s⁻¹ via modelling the kinematics of stars in the Milky Way disc, and with Bovy, Hogg & Rix (2009) who obtain V₀ = 246 ± 30 km s⁻¹ from trigonometric parallaxes. It is interesting to note that Arp (1986) found V₀ = 239 ± 17 km s⁻¹ by minimizing the velocity dispersion of nearby galaxies about the bulk flow, an argument very similar to the one exposed in Appendix B. Recently, Schönrich (2012) has devised a model-independent method for measuring V₀ using the position and velocities of kinematically hot stars in the solar neighbourhood, which yields V₀ = 238 ± 9 km s⁻¹, which also falls within the confidence interval of our fits. In contrast, Bovy et al. (2011) find V₀ = 218 ± 6 km s⁻¹ using data from the APO Galactic Evolution Experiment (APOGEE) spectroscopic survey. However, these authors also detect an offset between the Sun's rotational velocity with respect to the LSR of ≈22 km s⁻¹, which is a factor of ∼2 larger than the one measured by Schönrich et al. (2010). Accounting for this offset seems to reconcile their value with the one listed in Table 3. Finally, if we add the LSR azimuthal component to the value of V₀ we find that the transverse velocity component of the Sun with respect to the Milky Way is V_ϕ ≈ 257 km s⁻¹, which is in good agreement with the value derived by Reid et al. (2009), V_ϕ = 254 ± 16 km s⁻¹, from proper motions of masers in the Milky Way.

6.3 The (very) local value of |$\boldsymbol {H_0}$|

Fig. 11 shows a strong covariance between the cosmological parameters of our analysis, |$\Omega _\Lambda$| and h. The tests run in Section 3 suggest that the correlation follows approximately the isochrone |$t_0=t_0(\Omega _\Lambda ,h)$|⁠, where t₀ is the age of the Universe. The strong degeneracy arises from the scant sensitivity of the kinematics of nearby galaxies to the vacuum energy term in the equations of motion (see Section 2.2).

However, the shape of the covariance is such that even a modest prior on |$\Omega _\Lambda$| may be sufficient to put a tight bound on the value of the Hubble constant. Indeed, Fig. 12 shows that the posterior distribution function of H₀ returned by our analysis becomes relatively narrow once we incorporate the bounds on |$\Omega _\Lambda$| obtained from the spectrum of fluctuations in the cosmic microwave background (CMB) as priors in our analysis. For simplicity, we adopt a Gaussian prior on the value of the fractional vacuum density that follows the posterior distribution function derived from the analysis of the Planck data, i.e. |$\Omega _\Lambda =0.686\pm 0.020$| (Planck collaboration, paper XVI). Imposing distance cuts to the galaxy sample in the range 2.7 ≤ r_max/Mpc ≤ 3.0 does not significantly change our constraints.

Figure 12.

Posterior distribution functions for the Hubble constant H₀ derived from four different samples of nearby galaxies with distances to the Local Group barycentre r < r_max. In order to break the degeneracy between H₀ and |$\Omega _\Lambda$| shown in Fig. 11 we adopt Planck's measurement of the fractional vacuum density, |$\Omega _\Lambda =0.686\pm 0.020$|⁠, as a Gaussian prior in our fits. Note that our constraints on H₀ are insensitive to the value of r_max ≲3 Mpc. Beyond ∼3 Mpc the motion of galaxies appears to be strongly perturbed by galaxy clusters in the Local Group vicinity (see Figs 9 and 10).

Combination of CMB data and the dynamics of Local Group galaxies yields a local Hubble constant H₀ = 67 ± 5 km s⁻¹ Mpc⁻¹ at a 68 per cent confidence level. This result agrees very well with Planck's constraint, H₀ = 67.4 ± 1.4 km s⁻¹ Mpc⁻¹, and is also consistent at a 1σ level with the value derived from Cepheid data, H₀ = 72.5 ± 2.5 km s⁻¹ Mpc⁻¹ (Efstathiou 2014 and references therein), as well as with fits to the Hubble diagram over cosmological scales, e.g. H₀ = 74.4 ± 3.0 km s⁻¹ Mpc⁻¹ (Tully et al. 2013) and H₀ = 74 ± 4 km s⁻¹ Mpc⁻¹ (Sorce et al. 2013).

The velocity dispersion of the flow in the radial range sampled by our data (σ_H = 50 ± 4 km s⁻¹) is consistent with that expected galaxies embedded in large-scale walls (Aragon-Calvo et al. 2011).

Combination of Planck's constraint on the vacuum energy density with our measurement of the Hubble constant gives an estimate of the Universe's age, |$t_0=13.8_{-0.9}^{+1.0}$| Gyr (see Table 3).

Alternatively, one can use the age of the Universe obtained by Planck, t₀ = 13.813 ± 0.058 Gyr, as a prior in our models. In this case the parameter |$\Omega _\Lambda h^2$| is entirely constrained from the effects of the pressure term on the dynamics of nearby galaxies. Although note shown here, we obtain bounds on |$\Omega _\Lambda h^2$| that are in agreement with those given in Table 3, highlighting the consistency between the Hubble constant derived from local dynamics and its cosmological value.

6.4 Milky Way and Andromeda masses

The masses of the Milky Way and Andromeda follow directly from the bounds on M and f_m. The measurement rests on the assumption that both systems account for the entire mass of the Local Group, i.e. M = M_G + M_A, which appears reasonable given that the third brightest galaxy in our vicinity, M33, rotates with a relatively low speed,⁶V_rot,M33 ≈ 120 km s⁻¹ (Corbelli & Salucci 2000), whereas both the Milky Way and M31 have rotation velocities peaking at ∼250 km s⁻¹. Indeed, combination of the Tully–Fisher relation (Tully & Fisher 1977) with the rotational velocity for our own Galaxy indicates that M33 accounts for a tiny fraction of the Local Group mass, M_M33/M ∼ (V_rot,M33/V₀)⁴f_m/(1 + f_m) ≃ 0.028.

The median values of the Milky Way and Andromeda masses are |$M_{\rm G}=0.8_{-0.3}^{+0.4}\times 10^{12}\,{\rm M}_{\odot }$| and |$M_{\rm A}=1.5_{-0.4}^{+0.5}\times 10^{12}\,{\rm M}_{\odot }$| at a 68 per cent level (see Table 3). Models where the Milky Way is more massive than Andromeda (i.e. f_m ≥ 1) are ruled out with high significance (∼95 per cent).

It is worth stressing that, by virtue of the large distance range covered by the sample of Local Group galaxies, the quoted values correspond to the total mass of both galaxies. In contrast, most dynamical models in the literature use kinematic tracers (e.g. gas, stars, planetary nebulae and/or globular clusters) that only populate the inner regions of galactic haloes. From Newton's theorem any mass distribution outside the limiting radius of the data has no observational effect in a spherical or elliptical system. Hence, these models must assume a density profile in order to extrapolate the inner mass bounds to radii devoid of visible tracers (e.g. the virial radius, r_vir ∼ 260 kpc). Unfortunately, the outer mass profile of the Milky Way remains unknown, increasing the uncertainty of the extrapolation. For example, Smith et al. (2007) use high-velocity stars from the Radial Velocity Experiment (RAVE) survey (Steinmetz et al. 2006; Zwitter et al. 2008) to measure the local escape speed of our Galaxy. The range of values found by these authors (498 < v_esc/km s⁻¹ < 608) leads to a virial mass |$M_{\rm G,vir}=0.85_{-0.29}^{+0.55}\times 10^{12}\,{\rm M}_{\odot }$| under the assumption that the Milky Way halo follows a Navarro–Frenk–White (NFW) profile (Navarro, Frenk & White 1997). However, if the dark matter halo contracts as a result of dissipative processes (Mo, Mao & White 1998) the mass estimate increases to |$1.42_{-0.54}^{+1.14} \times 10^{12}\,{\rm M}_{\odot }$|⁠. Xue et al. (2008) extend the mass constraints to D ≲ 60 kpc by modelling the kinematics of blue horizontal branch stars (BHBs). Depending on whether or not the halo models are contracted the mass estimate varies between |$M_{\rm G, vir}=1.2_{-0.3}^{+0.4}\times 10^{12}$| and |$0.8_{-0.2}^{+0.2}\times 10^{12}\,{\rm M}_{\odot }$|⁠, respectively. This range is broadly consistent with Sakamoto, Chiba & Beers (2003) and Battaglia et al. (2005) who find |$0.8_{-0.2}^{+1.2}\times 10^{12}$| and |$2.5_{-1.0}^{+0.5}\times 10^{12}\,{\rm M}_{\odot }$|⁠, respectively, from a mix sample of globular cluster, giant stars and satellite galaxies; as well as with the recent estimates of Deason et al. (2012) based on Jeans modelling a sample of BHB stars out to D ∼ 90 kpc.

In principle, satellite galaxies sample a much larger volume of the galactic halo and may provide a tighter constraint on the virial mass of the host galaxy. In practice, their poorly known orbital distribution introduces a considerable uncertainty in this type of analysis. Using the locations and kinematics of 26 Milky Way satellites, Watkins, Evans & An (2010) constrain the Milky Way mass to lie within 0.7–3.4 × 10¹² M_⊙ depending on the assumed velocity anisotropy. Incorporating the proper motions of six satellites into the analysis narrows the mass range to 1.4 ± 0.3 × 10¹² M_⊙. Recently, Barber et al. (2014) have carried a comparison between the orbits of Milky Way satellites with known proper motions and the eccentricity distribution observed in the Aquarius N-body simulations.⁷ From this exercise they find that the Milky Way mass lies 0.6–3.1 × 10¹² M_⊙ with a best-fitting value of ∼1.1 × 10¹² M_⊙.

Our measurement of M31 mass, M_A = (1.5 ± 0.3) × 10¹² M_⊙, is in agreement with existing estimates. Seigar, Barth & Bullock (2008), who improve on Klypin, Zhao & Somerville (2002) analysis of the H i rotation curve using Spitzer 3.6-μm data, adopt an adiabatically contracted NFW halo profile, finding M_A,vir = (0.82 ± 0.02) × 10¹² M_⊙, which is broadly consistent with the value of ∼0.77 × 10¹² M_⊙ inferred by Geehan et al. (2006) and with the estimate of (0.37–2.1) × 10¹² M_⊙ derived from satellite kinematics (Côté et al. 2000). It is also very close to the lower limit of 0.9 × 10¹² M_⊙ obtained from kinematics of halo stars (Chapman et al. 2006). From the kinematics of 23 satellites Watkins et al. (2010) find that a mass in the range 0.85–1.6 × 10¹² M_⊙, concluding that the large uncertainty arising from the unknown orbital anisotropy prevents them to determine which of the two galaxies is actually the more massive. Using line-of-sight velocities for a sample of globular clusters in the stellar halo of M31 Veljanoski et al. (2013) estimate a mass 1.2–1.5 × 10¹² M_⊙, although the apparent association of many of these clusters to stellar streams (Mackey et al. 2010) may have some impact on the result. Stellar streams provide an independent constraint on the potential of M31. Fitting the kinematics the giant stream Ibata et al. (2004) and Fardal et al. (2013) derive a virial mass of (1.0 ± 0.5) × 10¹² and (2.0 ± 0.5) × 10¹² M_⊙, respectively.

Overall we find that the combined dynamical masses of the Milky Way and Andromeda published in the literature roughly match the Local Group mass derived from our analysis of the local Hubble flow. We discuss the implications of this result below.

7 DISCUSSION: MISSING MASS IN THE LOCAL GROUP?

7.1 Timing argument versus Hubble flow

The Local Group mass derived from the timing argument (∼5 × 10¹² M_⊙, see Fig. 10; also Li & White 2008; van der Marel et al. 2012a,b; Partridge et al. 2013; Yepes, Gottloeber & Hoffman 2014) is considerably larger than the combined masses of the Milky Way and M31 (∼2 × 10¹² M_⊙). The existence of large amounts of ‘missing’ mass in the Local Group with no visible counterpart poses a difficult problem to current galaxy formation models.

Interestingly, the mass derived from the local flow, M = 2.3 ± 0.7 × 10¹² M_⊙, is also a factor of ∼2–3 lower than the value suggested by the timing argument (see Section 6.4). This mismatch becomes the more intriguing if we take into account that both estimates rest on equation (2). It is worth discussing a number of mechanisms that could potentially bias the above measurements.

Let us start with the mass estimate derived from the local Hubble flow. In Section 4.3.2 we show that neglecting the quadrupole term in the equations of motion tends to overestimate the Local Group mass, a bias that grows under the presence of large peculiar motions in the synthetic data sets. Hence, correcting for this effect would widen the discrepancy with the timing argument even further. Also, the low velocity dispersion of the local Hubble flow, σ_H = 50 ± 4 km s⁻¹ (see Table 3), suggests that the overall magnitude of the peculiar motions is small. In this case our tests indicate that neglecting the quadrupole term in equation (2) leads to a minor bias in the mass estimate (see Section 4.3.2; ‘pair’ models).

The hierarchical mass growth of the Local Group also has little impact on the observed Hubble flow. In Section 2.3 we show that ignoring the time dependence of the Local Group potential in equation (2) tends to underestimate the Local Group mass, but the bias is so small that only a very rapid growth of the Milky Way and Andromeda masses would have a measurable effect on the observed kinematics of nearby galaxies.

It thus appears more simple to envision mechanisms that could potentially affect the mass obtained from the timing argument. In particular, it is worth following up the results of Van der Marel et al. (2012b), who observe a strong dependence between the mass derived from the relative motion between the Milky Way and M31 and the circular velocity of the Milky Way at the solar radius (V₀). Combination of equations (9) and (A2) yields

\begin{eqnarray} M\approx \frac{0.83d}{G}\bigg [(1.2+0.31\Omega _\Lambda )\frac{d}{t_0}-V_{{\rm h},{\rm A}} - {\boldsymbol v}_{\rm LSR}\cdot \hat{\boldsymbol r}_{\rm A} \nonumber\\ -V_0\sin (l_{\rm A})\cos (b_{\rm A})\bigg ]^2, \end{eqnarray}

(20)

where |${\boldsymbol v}_{\rm LSR}$| is the velocity vector of the local standard of rest (see Appendix A).

Equation (20) shows two points of interest. First, ignoring the cosmological constant term in the equations of motion lowers the Local Group mass by a small factor |${\sim } 0.55\Omega _\Lambda d^{3/2}/(\sqrt{GM}t_0)\sim 0.1$|⁠, which is in good agreement with the findings of Partridge et al. (2013). Second, the mass suggested by the timing argument is strongly sensitive to the value of V₀. In particular, the minus sign in front of this parameter and the fact that sin (l_A)cos (b_A) > 0 imply that the estimated Local Group mass drops if the circular velocity at the solar radius lies above the standard IAU value.

Fig. 13 illustrates this point in more detail. The blue dashed line shows the relation implied by equation (20). Colour-coded dots sample the posterior distributions on M and V₀, which reflect the uncertainty in the value of d and V_{h, A}, as well as in the cosmological parameters. For a better comparison with the results of previous sections we adopt Gaussian priors on these parameters, with H₀ = 67 ± 5 km s⁻¹ Mpc⁻¹ and |$\Omega _\Lambda =0.686\pm 0.020$|⁠. Comparison with the constraints derived in previous sections indicates that the local Hubble flow provides a much tighter bound on the mass of the Local Group (M = 2.3 ± 0.7 × 10¹² M_⊙; marked with dotted lines) than the timing argument, and that the amount of ‘missing’ mass is not statistically significant once we take into account the uncertain value of V₀ (see Section 6.2).

Figure 13.

Constraints on the Local Group mass as a function of the circular velocity of the Milky Way at the solar radius (V₀) using the relative motion between the Milky Way and M31 (the so-called ‘timing argument’) and the kinematics of nearby galaxies (‘Hubble flow’). Colour-coded dots sample the posterior distribution derived from the relative motion between the Milky Way and M31 using h = 0.67 ± 0.05 and Planck's prior on the fractional vacuum density, |$\Omega _\Lambda =0.686\pm 0.020$|⁠. For reference, a blue dashed line marks the solution to equation (20) with t₀ = 13.8 Gyr, r = d(t₀) = 0.783 Mpc and V_{h, A} = −300 km s^{− 1}. Dotted lines mark 68 and 95 per cent confidence intervals derived from our fits to the local Hubble flow (see Fig. 11) using uninformative priors on the cosmological parameters. The discrepancy between both methods eases by setting V₀ above the IAU concordance value.

7.2 The Local Group mass in a cosmological context

The dynamical models outlined in Section 2 assume that the Milky Way and M31 can be treated as point masses, thus neglecting the internal distribution of matter in these galaxies. In a ΛCDM Universe galactic haloes are expected to follow a close-to-universal density profile that falls with radius as ρ ∼ r⁻³ (e.g. Navarro et al. 1997). It is trivial to show that the mass profile |$M(r)=4\pi \int _0^r \rho (r^{\prime }) r^{\prime 2}{\rm d}r^{\prime }$| diverges at large radii, which complicates the interpretation of the mass estimates derived from the timing argument.

Li & White (2008) examined this problem and found that if the ‘true’ Local Group mass is taken to be the sum of the ‘virial’ masses (M₂₀₀) of the two dominant galaxies, where M₂₀₀ is defined as the mass within a sphere of mean density 200 times the critical value, the ratio of the true virial mass to that estimated from the ‘classical’ (⁠|$\Omega _\Lambda =0$| in equation 2) timing argument is ∼1.5 (see also Yepes et al. 2014). Adding a cosmological constant term in the equations of motion reduces the mismatch to a factor of ∼1.15 (Partridge et al 2013), with a considerable scatter.

Recent HST measurements of the transverse velocity vector between the Milky Way and M31 (van der Marel et al. 2012a) suggest that the galaxies are approaching each other on a head-on trajectory. Gonzalez, Kravtsov & Gnedin (2013) have surveyed the Bolshoi simulations (Klypin, Trujillo-Gomez & Primack 2011) in an attempt to find Local Group analogues on similar orbital configurations. When comparing the timing argument mass against the true virial masses these authors find that for systems moving on nearly radial orbits the timing argument overestimates the combined mass by a factor of ∼1.3–1.6. Fig. 13 shows than correcting for this bias and adopting a circular velocity V₀ ∼ 245 km s⁻¹ would bring the masses derived from the timing argument and the local Hubble flow into good agreement. The result M ∼ 2 × 10¹² M_⊙ would be compatible with the combined dynamical masses of the Milky Way and Andromeda, thus removing the need for ‘missing’ mass in the Local Group.

8 SUMMARY

We have inspected the equations of motion that govern the dynamics of nearby galaxies within the ΛCDM paradigm. Using analytical arguments we show that the time dependence of the Local Group potential has a small impact on the observed local Hubble flow. In contrast, our analysis indicates that the orbits of galaxies in the local volume can be strongly perturbed by the gravitational quadrupole that arises from by the dominant masses of the Milky Way and Andromeda. In particular, the quadrupole induces an attractive force along the axis joining the main galaxies which becomes repulsive in any transverse direction. Restricted N-body experiments of the local cosmic expansion show that the infall of galaxies on to the Local Group occurs preferentially along the axis joining the Milky Way and Andromeda, leading a highly anisotropic distribution of galaxies within several Mpc of the Local Group barycentre.

We devise a Bayesian method for analysing observations of the local Hubble flow using the timing-argument equations with a vacuum energy term. Our model fits simultaneously the Local Group mass, M = M_G + M_A, the mass ratio between our Galaxy and Andromeda, f_m = M_G/M_A, the circular velocity of the Milky Way at the solar radius, V₀ = V_c(R_⊙), the reduced Hubble constant, h, and the fractional vacuum density, |$\Omega _\Lambda$|⁠. Tests with synthetic data drawn from the restricted N-body models indicate that neglecting the potential quadrupole leads to Local Group mass estimates that can be overestimated up to ∼30 per cent, but does not affect the constraints on the rest of model parameters.

Applying our method to published locations and radial velocities of nearby galaxies returns a mass M = 2.3 ± 0.7 × 10¹² M_⊙, which is consistent with the combined dynamical masses of the Milky Way and M31 and does not require the presence of extra ‘missing’ mass between both galaxies, and a mass ratio |$f_{\rm m}=0.54^{+0.23}_{-0.17}$| which rules out dynamical models where the Milky Way is more massive than Andromeda with ∼95 per cent confidence. Both estimates are in very good agreement with those recently found by Diaz et al. (2014) who implement their own Bayesian modelling of the kinematics of Local Group galaxies (r ≲ 1.5 Mpc). The Milky Way's circular velocity at the solar radius, V₀ = 245 ± 23 km s⁻¹, is slightly lower than the peak velocity of M31 (Klypin et al. 2002), lending additional support to dynamical models where f_m < 1.

The cosmological parameters |$\Omega _\Lambda$| and h are strongly covariant. The shape of the covariance is such that introducing the CMB prior on |$\Omega _\Lambda$| is sufficient to put a tight bound on the value of the local Hubble constant, H₀ = 67 ± 5 km s⁻¹ Mpc⁻¹, which is consistent with that derived on cosmological scales and does not show evidence for a local ‘super-Hubble’ flow.

Overall we find that studying the dynamics of nearby galaxies in a broader cosmological context gives us clues on the hierarchical formation of the Local Group constituents as well as meaningful constraints on key cosmological parameters.

We are indebted to Jonathan Diaz who helped us to spot an error in the apex calculation that led to biased constraints on the mass ratio between Andromeda and the Milky Way. We thank Chervin Laporte and Alex Mead for contrasting our numerical set-up against the results of self-consistent cosmological N-body simulations. We also thank the anonymous referee for very useful comments. MGW is supported by NSF grant AST-1313045.

1

At the time the distance between our Galaxy and M31 was thought to be d ∼ 600 kpc, approximately 3/4 of the current value (McConnachie 2012).

2

The sample we have gathered from the literature in Section 5 includes a large number of galaxies at distances comparable to the separation between the Milky Way and Andromeda. Section 3 explicitly tests the validity of the point mass approximation with aid of N-body experiments.

3

Our force calculation includes a softening length ϵ = 5 kpc to avoid a divergence of the force during close encounters. We have checked that the choice of ϵ does not lead to qualitative changes in the results.

4

An alternative route can be pursued by adding additional terms in the equations of motion to account for the gravitational perturbations induced by the structures surrounding the Local Group (e.g. Mohayaee & Tully 2005; Courtois et al. 2012). Such undertaking, however, introduces a number of complexities that go beyond the scope of this paper.

5

Interested readers can find a formal proof of this argument in Peñarrubia, Koposov & Walker (2012).

6

Deep photometric surveys of the M31–M33 system show that M33 acted on by tides (McConnachie et al. 2009), which may alter the outer rotation curve of this galaxy.

7

The Aquarius project consists of six dark matter-only realizations of a ∼10¹² M_⊙ halo which do not account for the enhanced disruption rate of satellites moving on highly eccentric orbits under the presence of a disc component (D'Onghia et al. 2010; Peñarrubia et al. 2010).

8

Note that Karachentsev et al. (2009) find f_m ∼ 1 using very similar arguments. The slightly discrepant result may be due to revised measurements of distances and velocities (e.g. Tucana) and/or the addition of newly discovered systems (e.g. Leo P).

APPENDIX A: COORDINATE FRAME CONVERSION

In order to calculate the velocity of galaxy tracers with respect to the Local Group we follow Karachentsev & Makarov (1996) method. First, using equation (3) we calculate the solar velocity vector in a Local Group-centric frame (the so-called apex) as

\begin{equation} {\boldsymbol v}_{\odot }= {\boldsymbol V}_0+ {\boldsymbol v}_{\rm LSR} - \frac{v_{\rm A}}{1+f_{\rm m}} \hat{\boldsymbol r}_{\rm A}, \end{equation}

(A1)

where v_A is the velocity of M31 in the Galactic standard of rest (GSR) and f_m = M_G/M_A is the Milky Way-to-M31 mass ratio. The quantity |${\boldsymbol V}_0$| denotes the rotational velocity vector at the solar radius R₀, |${\boldsymbol v}_{\rm LSR}$| is the solar motion with respect to the GSR and |$\hat{\boldsymbol r}_{\rm A}$| is the unit vector pointing towards the centre of M31. We choose a right-handed Galactocentric coordinate system wherein the sun is located at (−R₀, 0, 0) and moves with a velocity |${\boldsymbol V}_0=(0,V_0,0)$| with respect to the MW centre. Following Schönrich et al. (2010) the LSR vector is fixed to |${\boldsymbol v}_{\rm LSR}=(11.1, 12.2,7.2)\,{\rm km \, s^{-1}}$|⁠. While we find that the current uncertainty in the value of R₀ (of the order of a kpc) has no visible impact on our results owing to the large heliocentric distances of the sample galaxies gathered from the literature, the value of V₀ does introduce a significant element of uncertainty in our measurements.

The parameter v_A in equation (A1) corresponds to the radial velocity of M31 in the GSR, that is

\begin{equation} v_{\rm A}=V_{{\rm h},{\rm A}} + ({\boldsymbol V}_0+ {\boldsymbol v}_{\rm LSR})\cdot \hat{\boldsymbol r}_{\rm A}, \end{equation}

(A2)

with |$\hat{\boldsymbol r}_{\rm A}=(\cos [l_{\rm A}]\cos [b_{\rm A}], \sin [l_{\rm A}]\cos [b_{\rm A}],\sin [b_{\rm A}])$|⁠. Following McConnachie (2012) we adopt the following Galactocentric coordinates for M31: (l, b)_A = ( − 121| $_{.}^{\circ}$|2, −21| $_{.}^{\circ}$|6), d = 0.783 Mpc and V_{h, A} = −300 km s^{− 1}.

For a given solar apex the radial velocity of a tracer galaxy with respect to the Local Group centre can be calculated as

\begin{equation} v=V_{{\rm h}} + \Delta v, \end{equation}

(A3)

where V_h is the heliocentric radial velocity, and Δv is the projection of the galaxy position vector on to the solar apex, i.e. |$\Delta v={\boldsymbol v}_{\odot }\cdot \hat{\boldsymbol r}_{\rm g}$| and |$\hat{\boldsymbol r}_{\rm g}={\boldsymbol r}_{\rm g}/|{\boldsymbol r}_{\rm g}|$|⁠.

APPENDIX B: THE MASS RATIO BETWEEN OUR GALAXY AND ANDROMEDA

For a given set of solar parameters the main uncertainty in the determination of the apex vector reduces to the mass ratio between our Galaxy and M31 (f_m). In this work we have explored two methods for constraining the value of the mass ratio between our Galaxy and M31 (f_m). In Section 4 this parameter is implemented in the likelihood function that fits orbits to the location and velocities of Local Group neighbour galaxies.

Here we also explore a geometrical method to constrain both parameters which does not require solutions to the equations of motion (2), but whose results turn out to be in excellent agreement with the convolved Bayesian fits explored above. A similar approach was followed by Arp (1986) to measure the circular velocity of the Milky Way at the solar radius (V₀), by Karachentsev et al. (2009) to pin down the location of the Local Group barycentre and by Karachentsev & Makarov (1996) and Courteau & van den Bergh (1999) to constrain the solar apex.

The method rests upon the assumption that that the distribution of peculiar velocity about the Hubble flow is cold and invariant under spatial translations. Hence, the fact that neighbour galaxies are distributed over a large area of the sky implies that a biased choice of f_m and/or V₀ must necessarily introduce a scatter in the distribution of radial velocities derived from equations (A1) and (A3). It follows that the proper choice of these parameters must be that that minimizes the scatter of the distance–velocity relation of Local Group galaxies when expressed in a Local Group-centric coordinate frame.

Fig. B1 illustrates the dependence of the radial phase-space location of the galaxy sample (see Section 5 for details). From left to right the first three panels adopt f_m = 0.5, 1.0 and 2.0. By eye it is clear that the values f_m ≲ 1 yield narrower distributions than f_m > 1. To measure the scatter in these distributions we fit straight lines y = a + br (dotted lines). The rightmost panel shows the reduced χ² values as a function of f_m. Here χ² is defined as

\begin{equation} \chi ^2=\sum _{i=1}^{N_{\rm g}}\frac{(v_{{\rm rad},i}-y_i)^2}{\sigma _i^2}, \end{equation}

(B1)

where |$\sigma _i^2$| is the variance associated with the measurements of distance (ϵ_D) and heliocentric velocity (ϵ_V) for the ith galaxy in the sample,

\begin{equation} \sigma _i^2= \epsilon _{V,i}^2 + (b \epsilon _{D,i})^2 + \sigma _{\rm m}^2. \end{equation}

(B2)

In Section 6 we find that setting the parameter σ_m ≈ 38 km s^{− 1} accounts for the presence of an intrinsic dispersion in the distance–velocity relation that goes beyond that introduced by observational errors.

Figure B1.

Distribution of radii and velocities of the galaxies in Table 1 for different values of Milky Way to M31 mass ratio (f_m). From left to right blue dotted lines show best-fitting straight lines with a = −78, −101, −101 km s⁻¹ and b = 87.4, 108.2, 111.7 km s⁻¹ Mpc⁻¹, respectively. The fourth panel on the right shows the reduced χ² value as a function of f_m. Notice that models in which the M31 is more massive than the Milky Way significantly increase the scatter of the distribution.

The right-hand panel of Fig. B1 shows that choosing Local Group models where M31 is more massive than our Galaxy leads to strongly scatted distributions of radii and velocities. The scatter minimizes at f_m ≈ 0.5 if we adopt a circular velocity V₀ = 220 km s⁻¹, and at f_m ≈ 0.55 for V₀ = 240 km s⁻¹. This measurement is fully consistent with the constraints on f_m derived in Section 6, and suggests that our Galaxy may be a factor of ∼2 less massive than M31.⁸

Author notes

†

Ramón y Cajal fellow.