1. Introduction
The dynamics of galaxies constitute one of the most profound challenges faced by astrophysics. This is a long-standing problem: suggestions that the gravitational dynamics of the Milky Way galaxy could be dominated by nonluminous matter predate confirmation that the Milky Way is itself one of many galaxies in the universe (Appendix D, Kelvin Reference Kelvin1904). Conversely, appeals to ‘modified gravity’ at galactic length scales undermine the bedrock assumption that the laws of physics revealed by terrestrial experimentation apply equally on astrophysical scales—a paradigm with its origins in Newton’s explanation of the Keplerian dynamics of the solar system. It is thus hard to overstate the potential impact of this topic.
On scales of light years and above, however, we can never directly recover the local gravitational force; observables are positions, velocities (from redshifts or stellar proper motions) and integrated quantities such as photon trajectories, the outcome of structure formation, or microwave background temperature anisotropies. Forces are proportional to accelerations and measuring these directly requires the ability to discern changes in the velocities of objects over time. Given that stellar velocities are hundreds of kilometres a second in the galactic rest frame at radii of 
$\sim 10\ \text{kpc}$, we immediately estimate that typical stellar accelerations are 
$\sim10^{-10}~{\rm m\,s}^{-2}$, corresponding to a change in velocity of a few centimetres per second per decade.
Spectrographs now under construction aim to measure stellar radial velocities at precisions of 
$\mathcal{O}(10~{\rm cm\,s}^{-1})$ (Fischer et al. Reference Fischer, Jurgenson, McCracken, Sawyer, Blackman and Szymkowiak2017). Exoplanet searches provide the primary motivation for these ultra-high resolution instruments; an Earth-like planet orbiting a Sun-like star induces a velocity of roughly 
$10\ {\rm cm\,s}^{-1}$ in the star. Given that Doppler shifts are proportional to v/c, this is a parts-in-a-billion measurement of the stellar spectrum—a stunning technical achievement. However, it is also the latest step in a decades-long sequence of instruments. The challenges of continuing this rate of progress should not be underestimated but it is plausible that ‘next-to-next’ generation instruments will be capable of measuring stellar accelerations due to the galactic gravitational field over a baseline of several years.
The present paper explores this possibility and examines possible issues with its implementation. Such measurements would constitute a galactic analogue to the Sandage–Loeb effect (Sandage Reference Sandage1962; Loeb Reference Loeb1998). Just as resolving the changing expansion rate of the universe would provide novel constraints on the dynamics and composition of the universe as a whole, direct measurements of stellar accelerations open new lines of enquiry into the structure of the galaxy and the fundamental properties of gravity. In particular, these would include direct tests of the validity of Newtonian dynamics on kiloparsec length scales.
Clearly, any feasible timeline for the full realisation of this strategy is measured in decades rather than years. However, direct measurements of the acceleration field can map the galactic mass distribution and probe the substructure of the Milky Way in novel ways. Moreover, this strategy provides an new approach to disentangling modified gravity from dark matter, a problem over a century old. Given the vast range of possible dark matter scenarios—each with its own parameter space—it is entirely conceivable (and perhaps even likely from some standpoints) that dark matter is effectively undetectable by any plausible instrumentation or experiment, even if the concordance cosmology is fundamentally correct.Footnote a
There is no fundamental connection between the radial acceleration induced by terrestrial planets orbiting sunlike stars and galactic dynamics, so it is perhaps unsurprising that little attention has been paid to the fact that exoplanet searches based on Doppler spectroscopy are approaching the regime where these measurements become possible. Moreover, as modified gravity theories typically contain a fundamental acceleration parameter 
$a_0$, which is necessarily of the same order as the accelerations experienced by stars moving in a galactic potential (Milgrom Reference Milgrom1983), this target is one that has been clearly visible in the literature. That said, Quercellini, Amendola & Balbi (Reference Quercellini, Amendola and Balbi2008) proposes mapping the Milky Way potential via the peculiar accelerations of globular clusters, but its authors are less than optimistic about the immediate applicability of the technique. We also draw attention to Ravi et al. (Reference Ravi, Langellier, Phillips, Buschmann, Safdi and Walsworth2018), which provides a complementary analysis of a similar proposal to the one presented here, and which was completed contemporaneously with this work.
Computing both line-of-sight (LoS) and transverse accelerations shows that direct measures of transverse accelerations would require a substantially greater advances in astrometric precision than the corresponding improvement needed in spectroscopic measurements of radial velocities. However, these two effects are coupled via perspective acceleration, as the changing viewing angle associated with stellar proper motion changes the component of the velocity vector parallel to the LoS. The magnitude of this effect decreases with distance but can easily swamp the gravitational acceleration. Fortunately, this effect can be estimated and accounted for with suitable precision by next generation astrometric measurements of proper motion and parallax.
Having argued that it will apparently be feasible to directly measure stellar accelerations, we consider two specific physical scenarios we could investigate with this data. Significant, localised dark matter over-densities or other substructure can be revealed by direct measurements of local accelerations. Furthermore, assuming that gravity (outside of lensing effects) is well-described by Newtonian dynamics, stellar acceleration information could provide a direct, tomographic measurement of the mass distribution in the Milky Way galaxy. Secondly, this work facilitates detailed comparisons of modified Newtonian dynamics (MOND) and Newtonian dynamics. Their predictions of standard and modified gravity overlap in the plane of the galaxy (by design, since modified gravity was built to account for these observations) but predictions for motion in the vertical direction generically differ.
In what follows we explore strategies to obtain and make use of acceleration information. In Section 2, we map out the acceleration field of the Milky Way with a conventional disc+halo model and compute expected radial and transverse accelerations. In Section 3, we consider the practicalities of implementing this proposal in terms of stellar magnitudes, spectral types, and the likely road map for future instrumentation. In Section 4, we show that perspective accelerations from changing lines of sight resulting from proper motion are similar to or larger than that derived from the Milky Way potential, but that high-precision astrometric measurements of parallax and proper motion can correct for this effect. In Section 5, we explore strategies for extracting the galactic signal from accelerations induced by stellar planetary systems. Finally, in Section 6 we discuss the ability of these measurements to detect large-scale inhomogeneities in the dark matter distribution and to probe large-scale gravitational physics. We discuss our results in Section 7.
2. Milky Way acceleration field
The Poisson equation links the density 
$\rho$, the Newtonian potential 
$\Phi$, and local gravitational acceleration
g
:
\begin{equation}4 \pi G \rho = \nabla^2 \Phi = - \nabla \cdot {\boldsymbol g},\end{equation}
where G is the gravitational constant. We can investigate the mass distribution and potential of the Milky Way if we can measure the accelerations of a set of tracer stars at positions 
$\boldsymbol{x}_i$ given by 
$\boldsymbol{{g}}(\boldsymbol{{x}}_i)$. Currently, only positions and velocities 
$\boldsymbol{{v}}_i$ are measured. The direct use of accelerations—or even a single component of the full 3D acceleration—would remove the need for many of the assumptions and approximations used when recovering the potential with methods such as the collisionless Boltzmann equation in combination with distribution function modelling (e.g. Binney & Piffl Reference Binney and Piffl2015) or Jeans equation modelling (Jeans Reference Jeans1922; Silverwood et al. Reference Silverwood, Sivertsson, Steger, Read and Bertone2016, 2018; see, e.g., Read Reference Read2014 for an overview), and reveal a more direct and unbiased view of the Milky Way’s potential and mass distribution.
To begin our discussion, we map the acceleration field of the Milky Way, taking into account the acceleration of the Sun. We work with the MilkyWayPotential2014 model in the galpy package (Bovy Reference Bovy2015). Radial and vertical accelerations at each point in the Rz plane are converted into heliocentric values by subtracting the acceleration at the solar position [
$R=8.2\ \text{kpc}$, 
$z=25\ \text{pc}$; Bland-Hawthorn & Gerhard (Reference Bland-Hawthorn and Gerhard2016)]. The galactocentric radial and vertical accelerations are projected along the LoS between the Sun and a given position to provide the LoS acceleration 
$a_{\rm LoS}$. The vector rejection (the acceleration orthogonal to the LoS) then gives the perpendicular acceleration 
$a_{\rm perp}$, which is converted to the proper motion acceleration via 
$\dot{\mu} = \tan^{-1} \left(a_{\rm perp} /d \right)$, where d is the distance to the Sun. The LoS and proper motion accelerations are shown in Figure 1. For convenience, we also express the accelerations as a velocity change over a fiducial interval of decade.
Figure 1. These plots show the LoS (top) and proper motion (bottom) accelerations relative to the Sun for stars across the 
$(R, \phi = \phi_\odot, z)$ plane, assuming a potential given by MilkyWayPotential2014 of Bovy (Reference Bovy2015). Accelerations are indicated by the left-hand scales on the colour bars; total change over a decade is shown on the right-hand scale. In the top plot, dotted and solid lines indicate the 10-yr 
$\Delta V_{\rm LoS} = 1{\rm cm\,s}^{-1}$ and 
$\Delta V_{\rm LoS} = 2{\rm cm\,s}^{-1}$ contours, respectively. The latter 
$\Delta V_{\rm LoS}$ is the sensitivity goal of the planned CODEX spectrograph (Pasquini et al. Reference Pasquini, Cristiani, Garcia-Lopez, Haehnelt and Mayor2010a,b, see Section 3).
For typical stellar separations, the impact of neighbouring stars is negligible compared to that of the overall galactic gravitational field. The Newtonian gravitational acceleration a on a star from another star of mass m at distance r is 
$a = Gm/r^2$. For a 
$1~{\rm M}_\odot$ star at 4 light years (or 1.3 pc, roughly the distance to Alpha Centauri), the acceleration is 
$\sim 10^{-13}~{\rm m\,s}^{-2}$, far smaller than the signal we would be seeking to detect.
Within a kiloparsec of the Sun, we find regions where the 10-yr change in velocity, 
$\Delta V_{\rm LoS}$, is 
$\sim\!1~{\rm cm\,s}^{-1}$. Within this range, stars in the vertical direction have higher accelerations, reflecting the steeper potential gradient moving away from the baryonic disc. At greater distances towards the Galactic Centre, there are regions where 
$\Delta V_{\rm LoS} \sim 10~{\rm cm\,s}^{-1}$ at galactocentric radii of 
$\lesssim 3\ \text{kpc}$.
The changes in proper motions over a 10-yr period are 
$\mathcal{O}({\rm nas\,yr}^{-1})$Footnote b. To put this into context, the Gaia mission has at best a proper motion sensitivity on the order of 
$10~\mu {\rm as\,yr}^{-1}$ (Gaia Collaboration et al. 2018), along with parallax precision of 
$\mathcal{O}(10~\mu {\rm as})$. Looking ahead, Gaia-NIR, a proposed successor mission with a similar scanning strategy to Gaia but observing in the near-infrared, would have comparable sensitivity, and when combined with data from Gaia could yield a factor 14 improvement in proper motion precision due to the 
$\sim\!20$-yr time baseline between the missions (Hobbs et al. Reference Hobbs2016; Hobbs & Høg Reference Hobbs, Høg, Recio-Blanco, de Laverny, Brown and Prusti2018), and double the parallax precision (Høg Reference Høg2014). The proposed Theia mission would have at best 
$\mathcal{O}$(100 nas) and 
$\mathcal{O}(100~{\rm nas\,yr}^{-1})$ parallax and proper motion precision, respectively (Boehm et al. Reference Boehm2017).
Sub-
$\mu {\rm as}$ and nas precision astrometry was discussed as a goal for a large L-class mission (Brown et al. Reference Brown2013; Brown Reference Brown2014). Optical and infrared astrometric precision scales as 
$\sigma \propto \lambda / (B \sqrt{N})$, where 
$\lambda$ is the wavelength of light observed, B is the mirror aperture size of a single telescope system or the baseline distance of an interferometer system, and N is the number of photons collected (Lindegren Reference Lindegren, Turon, O’Flaherty and Perryman2005). With current photon collection efficiency already near 100% and most stars emitting in the optical or infrared bands, the only avenue for increasing precision is the mirror aperture or interferometer baseline B. Brown et al. (Reference Brown2013) and Brown (Reference Brown2014) explored an interferometer mission with a collecting area of a few square metres and a baseline of 100–1 000 m, requiring ‘formation flying’ and considerable advances in global astrometry. Thus 
${\rm nas\,yr}^{-1}$ proper motion precisions are imaginable, but it will not be achieved in the near or even medium term. By contrast, there is a concrete path towards 
$\mathcal{O}({\rm cm\,s}^{-1})$ precision radial velocities, so we focus on this approach in what follows.
Separately from precision stellar astrometry and spectroscopy, pulsar timing is an increasingly important tool for high-precision astrophysical measurements that extend beyond the properties of the pulsars themselves. Pulsars come in a wide variety of types, with millisecond pulsars providing the greatest stability, rivalling that of simple atomic clocks (Kiziltan & Thorsett Reference Kiziltan and Thorsett2009). Individual binary systems allow direct tests of relativistic effects, including orbital decay via the emission of gravitational radiation (Hulse & Taylor Reference Hulse and Taylor1975; Weisberg & Huang Reference Weisberg and Huang2016) and pulsar timing arrays are expected to be sensitive to both point source and stochastic nanohertz gravitational wave backgrounds in the coming decade (Burke-Spolaor et al. Reference Burke-Spolaor2018).
To our knowledge, the use of a pulsar timing array as a probe of the overall galactic gravitational field has not been investigated in detail. However, timing data is certainly sensitive to local accelerations and has been used for exoplanet detection (e.g. Wolszczan & Frail Reference Wolszczan and Frail1992). Likewise, pulsar timing measurements have already been used to measure the high accelerations found in extreme locations such as globular clusters or the vicinity of large black holes (e.g. Peuten et al. Reference Peuten, Brockamp, Küpper and Kroupa2014; Perera et al. Reference Perera2017a,b; Gieles et al. Reference Gieles, Balbinot, Yaaqib, Hénault-Brunet, Zocchi, Peuten and Jonker2018). Pulsars are subject to spin-down which will produce frequency changes several orders of magnitude greater than those produced by Milky Way acceleration field, so reliably removing this will be critical to the detection of any linear change in velocity (Phinney Reference Phinney1992; Johnston & Karastergiou Reference Johnston and Karastergiou2017).
As we discuss in Section 4, the changing LoS induced by relative motion between the pulsar and the Sun will typically need to be accounted for in order to reveal the galactic acceleration; pulsar timing data has been used to obtain astrometric data from the parallax acceleration, but requires a model of the gravitational field of the galaxy as input (Matthews et al. Reference Matthews2016). This degeneracy can be broken by direct Very Long Baseline Interferometry (VLBI) measurements of pulsar parallax and proper motion and analogous measurements to those discussed here are potential targets for future pulsar timing studies.
3. Spectroscopic performance
Hunting exoplanets is a key motivation for improving spectrographic precision. Modern precision spectrographs are vibrationally isolated, temperature stabilised, and in vacuum chambers. Stable reference sources are needed to provide a time-independent calibration to track changes over long time intervalsFootnote c, for which the current state of the art is laser frequency combs linked to an atomic clock; see, e.g., Wright (Reference Wright2017) and Fischer et al. (Reference Fischer2016) for reviews of spectrographic measurements.
Figure 2 shows the roughly log-linear progress in the precision of planet hunting spectrographs since 1980. Such Moore’s Law-like progressions are not automatic, but reflect a deliberate, challenging, decades-long campaign to realise a series of significant technical improvements in spectrographic technology. However, it is not unreasonable to project that the 
$\mathcal{O}({\rm cm\,s}^{-1})$ velocity changes resulting from the galactic gravitational field will soon be within reach.
Figure 2. Spectrographic precision since 1981, and an extrapolation to 2040. Drawn from data presented in Wright & Gaudi (Reference Wright and Gaudi2013) and Fischer et al. (Reference Fischer2016). Full references for each spectrograph are listed in Appendix A.
Not all stars are suitable for high-precision radial velocity measurements. If a star is too hot (
$T_{\rm eff} > 10\,000\ \text{K}$), the absorption lines are smeared into a continuum, and if a star is too cold (
$T_{\rm eff} <\ 3\,500 \text{K}$), the lines can become too densely packed and overlapping to be distinguishable. Colder stars are also typically fainter, making it harder to achieve a specified signal-to-noise ratio (Lovis & Fischer Reference Lovis and Fischer2010).
A stellar spectrograph measures the integrated output from the stellar surface, and stars with rapidly moving features in their photospheres exhibit jitter, creating variation in their apparent Doppler velocities (e.g. Queloz et al. Reference Queloz2001). The impact of this jitter can be reduced by selecting stars for which these effects are typically small; ideal candidates are cool, old G- and K-dwarf stars (Hekker et al. Reference Hekker, Reffert, Quirrenbach, Mitchell, Fischer, Marcy and Butler2006; Wright & Gaudi Reference Wright and Gaudi2013). G- and K- dwarves show variations of 
$\mathcal{O}(100~{\rm cm\,s}^{-1})$ (Wright Reference Wright2005), while K giants have variations of 
$\mathcal{O}(1\,000~{\rm cm\,s}^{-1})$ (Hekker et al. Reference Hekker, Reffert, Quirrenbach, Mitchell, Fischer, Marcy and Butler2006). Careful modelling can separate jitter from true stellar motion as the physical mechanisms that produce jitter are correlated with variations in the detailed line shapes and profiles (see, e.g., Wright Reference Wright2017; Davis et al. Reference Davis, Cisewski, Dumusque, Fischer and Ford2017, for a full overview).
The same photospheric features that cause jitter in radial velocities can also generate astrometric jitter, impacting precision parallax and proper motion measurements (see, e.g., Eriksson & Lindegren Reference Eriksson and Lindegren2007, and references therein). For instance, a dark spot on the surface of a rotating star will shift the observed radial velocity up or down, depending on whether it is moving towards or away from us. Astrometric methods rely on measuring the location of the photocentre of a star, and the presence of a dark spot will shift this photocentre towards the opposite side of the star’s visible disc. As dark spots form, rotate, and/or dissolve, the photocentre will shift, perhaps mimicking accelerations acting on the star. Due to their common origin, radial velocity and astrometric jitter are correlated (Eriksson & Lindegren Reference Eriksson and Lindegren2007), opening a further area of complementarity between high-precision radial velocity and astrometric measurements that might be exploited in future work.
The cost of using low-jitter dwarf stars is their intrinsic faintness, which limits the distances at which high quality spectra can be obtained. Figure 3 shows the apparent V-band magnitude 
$m_V$ of dwarves, giants, and supergiants of the G5 spectral class over a range of distances r, neglecting extinction and reddening, and assuming absolute magnitudes of 
$M_V = 5.1,\ M_V = 0.9$, and 
$M_V = -4.6$, respectively (Binney & Merrifield Reference Binney and Merrifield1998). In Figure 1, we saw that to access regions of 
$\Delta V_{\rm LoS} \sim 1~{\rm cm\,s}^{-1}$ we need to measure stars at distances of around 1 kpc. To derive useful data from a G5 dwarf at this distance, the spectrograph and telescope would need to have a limiting magnitude of 
$\gtrsim\!15$. For reference, the EXPRES spectrograph that will be installed on the 4.3-m Discovery Channel Telescope (Jurgenson et al. Reference Jurgenson2016; Fischer et al. Reference Fischer, Jurgenson, McCracken, Sawyer, Blackman and Szymkowiak2017) has a limiting magnitude of 
$m_V \sim 7$, while the ESPRESSO spectrograph, which combines light from the four 8.2-m telescopes of the Very Large Telescope (VLT), has a limiting magnitude of 
$m_V \sim 17$ (Pepe et al. Reference Pepe2014; ESPRESSO Collaboration 2018). Consequently, implementing this technique will require the largest present-day telescopes or next generation instruments.
Figure 3. Apparent magnitude of G5 dwarfs, giants, and supergiants at various distances from the observer, assuming absolute magnitudes of 
$M_V = 5.1$, 
$M_V = 0.9$, and 
$M_V = -4.6$, respectively (Binney & Merrifield Reference Binney and Merrifield1998), and no reddening or extinction. For reference, we also show the approximate limiting magnitude of the EXPRES (Jurgenson et al. Reference Jurgenson2016; Fischer et al. Reference Fischer, Jurgenson, McCracken, Sawyer, Blackman and Szymkowiak2017) and EXPRESSO (Pepe et al. Reference Pepe2014; ESPRESSO Collaboration 2018) spectrographs.
On the other hand, giant and supergiant stars are intrinsically brighter and so can be seen at greater distances; ignoring extinction and reddening, a telescope with limiting magnitude 
$m_V = 15$ would be able to see a giant star out to 
$\sim\! 7$ kpc, reaching into high-acceleration regions near the galactic centre. In this region, 
$\Delta V_{\rm LoS}$ values can approach 
$15~{\rm cm\,s}^{-1}$ over a decade, thus requiring less stringent jitter modelling compared to that needed to extract the 
$\mathcal{O}({\rm cm\,s}^{-1})$ velocity changes expected for stars closer to Sun. Note too that high-cadence measurements can reduce the impact of jitter in a long time series, but this will require a trade-off with telescope resources.
With Figure 3, we can also assess the feasibility of using globular clusters to probe the Milky Way’s potential as suggested by Quercellini, Amendola & Balbi (Reference Quercellini, Amendola and Balbi2008). Without extinction the maximum distance that EXPRES and EXPRESSO can observe a G5 dwarf star are 
$\lesssim\,25$ pc and 
$\lesssim\,2.5$ kpc, respectively. Globular clusters are distant objects, and only two are closer than 2.5 kpc: NGC 6121 and NGC 6397 at 2.2 and 2.3 kpc, respectively (Harris Reference Harris1996). This rules out the use of EXPRES level spectrographs, and makes the use of ESPRESSO marginal. Using globular clusters as tracers would thus require the use of giants or super giants and suitable jitter modelling for such stars.
4. Perspective acceleration
The apparent radial velocities and proper motions of stars are projections of their velocity vector along the LoS and the plane of the sky. These projected quantities will change due to relative motion, independently of the acceleration induced by the galactic gravitational field. For a star travelling in a straight line at constant velocity V, the instantaneous change in the radial velocity 
$V_R$ from this effect is (van de Kamp Reference van de Kamp1977)
\begin{align}a_{\rm persp} = \frac{{\rm d} V_R}{{\rm d} t} = \frac{V_T^2}{r} = \frac{V^2}{r_P} \cos^3\theta ,\end{align}
where 
$V_T$ is the transverse velocity relative to the Sun, r is the distance to the star, 
$r_P$ is the perihelion distance,Footnote d and 
$\theta$ is the angle to the star measured from that point. Clearly the maximum acceleration occurs at perihelion, i.e. when 
$\theta = 0$, 
$V_T = V$, and 
$V_R = 0$. The change in radial velocity from this effect is always positive (when radial velocities away from us are set as positive): stars approaching perihelion have negative 
$V_R$ that is increasing towards zero, while those moving away from perihelion have positive and increasing 
$V_R$.
In the top panel of Figure 4, we plot the perspective radial acceleration 
$a_{\rm pers}$ as function of tangential velocity 
$V_T$ and distance r, showing the contour 
$a_{\rm pers} = 10^{-10}~{\rm m\,s}^{-2}$, which roughly matches the expected gravitational accelerations shown in the top panel of Figure 1. The effect diminishes with distance, but any star within a few kiloparsecs of the Sun with a relative transverse velocity above 
$\sim\! 100~{\rm km\,s}^{-1}$ has a perspective acceleration comparable to or greater than that expected from the Milky Way potential.
Figure 4. Top: Apparent radial acceleration 
$a_{\rm pers}$ caused by perspective acceleration for given transverse velocities and distances. Bottom: Maximum 
$V_T$ compatible with a perspective acceleration absolute error of 
$\sigma(a_{\rm pers}) = 10^{-12}~{\rm m\,s}^{-2}$, for varying parallax and proper motion precisions.
On the face of it, this limits us to stars with relatively low transverse velocities, but sufficiently accurate astrometric measurements will let us correct for this effect. From a rudimentary error analysis, neglecting correlations, the uncertainty on 
$a_{\rm pers}$ is
\begin{align}\sigma(a_{\rm persp}) = V_T^2 \sqrt{\frac{ \sigma_p^2}{1 {\rm AU}^2} + \frac{4 \sigma_\mu^2}{V_T^2}},\end{align}
where 
$\sigma_p$ is the uncertainty in parallax (in radians) and 
$\sigma_\mu$ is the uncertainty in proper motion (in radians per secondFootnote e). At this level the uncertainty on 
$a_{\rm pers}$ is independent of distance r. Rearranging Equation (3) gives an expression for 
$V_T$ in terms of 
$\sigma_p$, 
$\sigma_\mu$, and the desired level of precision for 
$\sigma(a_{\rm pers})$
\begin{align}V_T^2 = \frac{1 {\rm AU}^2}{\sigma_p^2} \left[-2 \sigma_\mu^2 \pm \frac{1}{2}\sqrt{16\sigma_\mu^4 + \frac{4\sigma_p^2}{1 {\rm AU}^2} \sigma(a_{\rm persp})^2}\right].\end{align}
As 
$\sigma(a_{\rm pers}) \propto V_T^2$, we can find an upper limit on the transverse velocities for which we can deduce the perspective acceleration to a desired precision, given 
$\sigma_p$ and 
$\sigma_\mu$.
In the bottom panel of Figure 4, we plot the maximum transverse velocity compatible with 
$\sigma(a_{\rm pers}) = 10^{-12}~{\rm m\,s}^{-2}$ for a range of parallax and proper motion uncertainties. This is a stringent requirement, given that it is 1% of the typical Milky Way acceleration, but it is not an unattainable one. The best Gaia precisions are 
$\sigma_p \sim 10~\mu {\rm as}$ and 
$\sigma_\mu \sim 10~\mu {\rm as\,yr}^{-1}$, yielding a maximum 
$V_T = 55~{\rm km\,s}^{-1}$. Searching the Gaia Data Release 2 catalogue (Gaia Collaboration et al. 2016, Reference Collaboration2018) shows that there are 
$\sim \!36$ million stars within 1 kpc of the Sun and with parallax errors less than 10%, and 
$\sim\! 55\%$ of these have transverse velocities below this threshold. An order of magnitude improvement to 
$1~\mu {\rm as}$ parallax and 
$1~\mu {\rm as\,yr}^{-1}$ proper motion precision pushes the 
$V_T$ limit to 
$175~{\rm km\,s}^{-1}$ and the percentage of our Gaia DR2 subsample under the 
$V_T$ threshold to 99%. Achieving Theia-level precisions of 
$\mathcal{O}$(100 nas) and 
$\mathcal{O}(100~{\rm nas\,yr}^{-1}$) would increase the 
$V_T$ limit beyond 
$500~{\rm km\,s}^{-1}$.
5. The exoplanetary background
5.1. Signals from exoplanetary systems
High-precision radial velocity measurements are driven by searches for exoplanetary systems. Ironically, these same exoplanets will be a key source of ‘noise’ that must be accounted for and subtracted in order to reveal possible gravitational accelerations. Radial velocity exoplanet searches look for accelerations induced in stars by their orbiting planets, as the star and planet(s) move about their common centre of gravity. For a single planet with period P, the induced velocity oscillation has the semiamplitude
\begin{align}K = \left( \frac{P}{2\pi G} \right) ^{-\nicefrac{1}{3}} \frac{M_p \sin i}{M_*^{\nicefrac{2}{3}}} \left(1-e^2\right)^{-\nicefrac{1}{2}},\end{align}
where 
$M_*$ and 
$M_p$ are the masses of the star and planet, respectively, e is the eccentricity, and i is the inclination of the orbital plane normal vector with respect to the LoS, with 
$i=0$ indicating the system is face on (Wright & Gaudi Reference Wright and Gaudi2013). Radial velocity measurements determine the combined term 
$M_p \sin i$, and thus there is a degeneracy between the mass of the planet and the orbital inclination.
In Figure 5, we show the values of K for a range of planetary orbital radii and masses, assuming a simplified circular (
$e=0$) and edge-on (
$i=90$) orbit around a host star of mass 
$M_* = 1~{\rm M}_\odot$, along with the orbital radii and masses of the eight planets and one dwarf planet of our own solar system for comparison.
Figure 5. Semiamplitudes K (Equation (5)) for a range of exoplanet orbital radii and masses, assuming the host star mass is 
$M_* = 1~{\rm M}_\odot$, the inclination of the orbital plane is 
$i = 90^\circ$, and the eccentricity 
$e = 0$. Also shown are the orbital radii and masses for the eight planets of our solar system and the dwarf planet Pluto, along with a region of radii and masses for the so-called ‘Hot Jupiters’—large gas giants orbiting close to their host stars.
Assuming this same simplified circular and edge-on orbit, the reflex velocity of a star due to an orbiting planet is then 
$V_{\rm LoS} = K \cos(\nu(t))$, where 
$\nu \in [0,2\pi]$ is the angle on the orbital plane, with 
$\nu = 0$ corresponding to the 3 o’clock position, and increasing counterclockwise. Over a time interval 
$t_1-t_0$, the change in velocity of the star is
\begin{align}\Delta V_{\rm LoS} = K \left[\cos\left(\nu(t_1)\right) - \cos\left(\nu(t_0)\right)\right].\end{align}
To give a qualitative picture of the impact of various planetary configurations, Figure 6 shows the cumulative change in velocity over 10 yr of a 
$1~{\rm M}_\odot$ star induced by planets of given mass and orbital radius (linked to the period via 
$P=2\pi \sqrt{R^3/G M_*}$), assuming the planet is transiting from 
$\nu(t_0) = \pi/2 - \omega \times 10{\rm yr}/2$ to 
$\nu(t_1) = \pi/2 + \omega \times 10{\rm yr}/2 $, e.g. directly behind the star which maximises the LoS velocity changeFootnote f. For reference, we also show the four gas giant planets and the dwarf planet Pluto from our solar system, along with lines indicating the mass of Earth and the mass range for brown dwarfs. We show contours in 
$\Delta V_{\rm LoS}$ at 
$0.1~{\rm cm\,s}^{-1}$, 
$1.0~{\rm cm\,s}^{-1}$, and 
$10.0~{\rm cm\,s}^{-1}$, roughly bracketing the typical 
$\Delta V_{\rm LoS}$ values generated by the Milky Way potential (see Figure 1). We exclude planets with periods of less than 10 yr, assuming that planets completing a full orbit in the observation period can be detected, modelled, and marginalised out.
Figure 6. Cumulative change in velocity of a 
$1~{\rm M}_\odot$ star by a single planet of a given mass and orbital radius over a 10-yr span. This assumes a circular planetary orbit (
$e=0$) that is edge-on to our LoS (
$i=0$). Also shown are the orbital radii and masses for Jupiter, Saturn, Uranus, Neptune, and the dwarf planet Pluto.
Gas giants with radii of 
$\mathcal{O}(1-10~{\rm AU})$, such as Jupiter and Saturn, induce large, periodic (or at least time-varying) velocity changes that would be easily detectable. However, large planets with orbital radii in 
$\mathcal{O}(10-100 ~{\rm AU})$ range (‘Neptunes’) will be more challenging, as they can yield a 
$\Delta V_{\rm LoS}$ of a similar magnitude to the acceleration to the Milky Way, which will be roughly linear with time on intervals much less than their orbital periods. Likewise Earth mass planets with orbital radii of 
$\mathcal{O}(10 {\rm AU})$ will contribute similar signals, albeit with a distinct time derivative.
One might hope to reduce the exoplanetary background by selecting systems with small inclination angle i. However the orbital plane would need to be determined with exquisite accuracy: the 10-yr 
$\Delta V_{\rm LoS}$ for a Jupiter mass and period planet at 
$i=90^{\circ}$ is 
$\sim\! 640~ {\rm cm\,s}^{-1}$, and only for 
$i = 0.1^\circ$ does this drop to 
$\sim\! 1~ {\rm cm\,s}^{-1}$. Furthermore, studies of our own solar system and exoplanetary systems found with Kepler suggest mutual inclinations between the planets of 
$\sim\! 1^{\circ}-3^{\circ}$ (Fabrycky et al. Reference Fabrycky2014). Again using the solar system as an example, even if we had the Jupiter-type planet’s orbital plane at or below the required 
$i = 0.1^\circ$ precision, a Saturn mass and period planet with a mutual inclination of 
$2^\circ$ would produce a 
$\Delta V_{\rm LoS}$ of 
$\sim 10~ {\rm cm\,s}^{-1}$, overwhelming the Milky Way acceleration.
5.2. Reducing the exoplanetary background
We focus on broad strategies that have the potential to reduce or exclude the exoplanetary signal, leaving the detailed data analysis challenge associated with extracting the galactic acceleration signal from radial velocity data to future workFootnote g. We can see three specific approaches: identifying systems for which this background is likely to be intrinsically small; combining radial velocity information with other data on planetary systems to improve the fitting process; and combining spatially adjacent stars to separate the correlated acceleration arising from the galactic potential from the largely random residual accelerations caused by exoplanets.
5.2.1. Single-system measurements
If we can identify classes of stars which are substantially less likely to have exoplanets (or at least those in the most problematic categories identified above), the 
$\Delta V_{\rm LoS}$ due to the Milky Way potential can be more easily and reliably extracted. Given the timeframe in which this experiment is likely to be conducted, it is possible that signatures (based on detailed elemental abundances, spectral type, and dynamical properties) of stars with complex planetary systems will be more clearly quantified than is the case at present.
Stars moving with very high speeds relative to the Milky Way rest frame may be a fruitful category to consider, as they may be born without planets, or be stripped of their planets by the mechanism responsible for accelerating them. These stars come in several subcategories, as delineated by Boubert et al. (Reference Boubert, Erkal, Evans and Izzard2017). Hypervelocity stars (HVSs) are stars with sufficient velocity to be unbound from the Milky Way, regardless of production mechanism; this is roughly 
$\sim\! 500~{\rm km\,s}^{-1}$ in the solar neighbourhood, rising to 
$\sim\! 600 ~{\rm km\,s}^{-1}$ in the galactic centre (Williams et al. Reference Williams, Belokurov, Casey and Evans2017; Monari et al. Reference Monari2018). Hills stars are accelerated via a three-body interaction of a tight binary with the central black hole of either the Milky Way or the Large Magellanic Cloud (LMC), a process called the Hills mechanism (Hills Reference Hills1988), while runaway stars are accelerated by the supernova explosion of their binary companion (Blaauw Reference Blaauw1961).
Runaway stars (Blaauw Reference Blaauw1961) are posited to form from binary system where the larger star sheds much of its mass in a supernova explosion, releasing the less massive star at its (potentially high) orbital velocity. Velocities can be up to 
$\sim\! 200\,{\rm km\,s}^{-1}$ in the galactic rest frame and point back to where the supernova explosion took place. This ejection mechanism may not strip the star of its exoplanets, but stars that form in relatively tight binaries are possibly less likely to possess long-period planets to begin with.
Stars can also be accelerated to high velocities through dynamical encounters in dense clusters—these fall outside the categorisation of Boubert et al. (Reference Boubert, Erkal, Evans and Izzard2017) but can be referred to as dynamical ejection scenario (DES) runaways (Hoogerwerf et al. Reference Hoogerwerf, de Bruijne and de Zeeuw2001). The bulk of DES runaways are actually walkaway stars, with velocities 
$< 30~{\rm km\,s}^{-1}$ (Eldridge, Langer & Tout Reference Eldridge, Langer and Tout2011; Renzo et al. Reference Renzo2018). Whether the DES process can remove any planets requires further investigation, but would likely depend on the proximity of the encounters.
The Hills mechanism accelerates stars via the interaction of a tight binary with a massive black hole, with one star being captured by the black hole and the other ejected at velocities of 
$\mathcal{O}(1\,000~{\rm km\,s}^{-1})$. For Hills stars in the Milky Way, the black hole responsible would be Sagittarius 
$\text{A}^*$ (Hills Reference Hills1988), or possibly a massive black hole in the LMC (Boubert & Evans Reference Boubert and Evans2016).
The fate of planetary systems around Hills stars was analysed by Ginsburg, Loeb & Wegner (Reference Ginsburg, Loeb and Wegner2012), who performed simulations where both members of the binary have exoplanets, finding that hypervelocity Hills stars can retain planets. However, these were typically very tightly bound to the parent star with orbital radii in the range of 0.02–0.05 AU. At smaller radii planets will interact with the star itself; at larger radii, dynamical instabilities eject planets or cause them collide with one of the stars before the final interaction with the black hole. Thus, the surviving planets have very short periods making their radial velocity signal easier to extract.
Initially theorised in 1988, the first hypervelocity Hills star was discovered in 2005 (Brown et al. Reference Brown, Geller, Kenyon and Kurtz2005). This star, HVS1, is a B-type star with an effective temperature of 
$T_{\rm eff} \sim 10\,500$ K, making it far too hot for precise radial velocity measurements (see Section 3). Other candidates are likewise hot, massive, short-lived stellar types, as they are easily spotted in a background of old halo stars (Brown et al. Reference Brown, Geller, Kenyon and Kurtz2006, Reference Brown, Geller, Kenyon, Kurtz and Bromley2007). However, nothing appears to preclude older hypervelocity Hills stars—they could in fact be more numerous but simply harder to find in the halo (Kollmeier & Gould Reference Kollmeier and Gould2007; Kollmeier et al. Reference Kollmeier2010).
Several features allow Hills stars to be selected from stellar catalogues. Their velocity vectors should (roughly) point away from their origin, the Galactic Centre. Secondly, since Hills stars come from the Galactic Centre, they should have the high metallicity characteristic of their origin. For instance, Zhang, Smith & Carlin (Reference Zhang, Smith and Carlin2016) found 29 F- and G-type low-mass stars with orbits and high metallicities consistent with Galactic Centre origins, though none had velocities above the escape speed and so did not qualify as HVSs. However, such stars may still be suitable for our purposes if they are above the threshold required to remove problematic long-period exoplanets.
As discussed in Section 4, HVSs necessarily have large perspective accelerations; using these stars as probes of the galactic magnetic field will need either increased astrometric precision to calculate their proper motion and parallax and thus their radial perspective acceleration, or a further cut to remove stars that were not moving in directions closely aligned to the LoS, reducing their transverse velocity.
Finally, a separate approach would be to focus on relatively tight binaries—in such systems planetary systems are more likely to be disrupted, reducing the likelihood of finding the sort of ‘Neptunes’ that would potentially contribute the most challenging radial velocity changes. However, these systems would likely be spectroscopic binaries, increasing the difficulty of extracting the underlying radial velocity.
5.2.2. Combining exoplanet detection techniques
A different strategy would be to target stars for which complementary exoplanet detection methods had (at least partially) mapped their planetary systems, providing additional constraints on the fit to the radial velocity data and so facilitating the extraction of the Milky Way acceleration signal. Critically, this includes null results that would exclude possible planetary candidates. In particular, it would be valuable to rule out the presence of long-period/large-radius Neptunes (
$\gtrsim 10\ \text{yr}$, 
$\gtrsim 5\ \text{AU}$), and to break the 
$M_p \sin i$ degeneracy between exoplanetary masses and radial accelerations.
Direct imaging searches are ideally suited for detecting large planets at large radii, and thus long periods, for the obvious reason that this maximises the contrast and angular separation on the sky between star and planet (Wright & Gaudi Reference Wright and Gaudi2013). Direct imaging at multiple epochs can also break the 
$M_p \sin i$ degeneracy (Mede & Brandt Reference Mede and Brandt2017; Hełminiak et al. Reference Hełminiak2016). Several direct imaging instruments are currently operational, such as the Gemini Planet Imager (Macintosh et al. Reference Macintosh2006), VLT-SPHERE (Beuzit et al. Reference Beuzit2008), and Subaru Coronagraphic Extreme Adaptive Optics (SCExAO) (Jovanovic et al. Reference Jovanovic2015) which is also serving as a technology test bed for exoplanet imaging with the Extremely Large Telescope (ELT) (Jovanovic et al. Reference Jovanovic, Esposito and Fini2013; Chauvin Reference Chauvin2018). There are also several planned and proposed space missions to conduct direct imaging, such as WFIRST (Collaboration Reference Collaboration2018; Noecker et al. Reference Noecker, Zhao, Demers, Trauger, Guyon and Jeremy Kasdin2016) and the New Worlds Observer (Lo et al. Reference Lo, Glassman, Dailey, Sterk, Green, Cash and Soummer2010). However, to be useful these imaging studies would need to be sensitive to planets at kiloparsec distances.
Gravitational microlensing is another method well suited to detecting exoplanets at larger radii (e.g. Tsapras Reference Tsapras2018) and hence is highly complementary. As a foreground lens object transits the LoS between a background source and an observer, the magnification of the source will smoothly increase, peak, and then decrease symmetrically in time. If the lens object hosts an exoplanet, it too can focus the light of the background star, adding a secondary off-centre peak in the magnification curve (Mao & Paczynski Reference Mao and Paczynski1991; Gould & Loeb Reference Gould and Loeb1992). Microlensing requires an element of luck as the lens star and planet have to cross sufficiently close to the LoS to the background source, but known lenses could provide an interesting set of target stars.
Astrometric detection could also break the 
$M_p \sin i$ degeneracy. This method relies on the same reflex motion of the star as the radial velocity technique, but is sensitive to the orthogonal velocity component, which induces 2D proper motion variations on the plane of the sky. Combining with radial velocities to yield a full 3D determination of the velocity variation allows for a determination of the inclination of the orbital plane, breaking the 
$M_p \sin i$ degeneracy (Wright & Howard Reference Wright and Howard2009; Wright & Gaudi Reference Wright and Gaudi2013). For a fixed perpendicular velocity, the measured proper motion decreases with distance, meaning that this technique is most effective with nearby stars. For example, a perpendicular velocity of 
$1~ {\rm m\,s}^{-1}$ would translate to a proper motion of 
$\sim\! 150~ \mu {\rm as\,yr}^{-1}$ at 1.3 pc (i.e. the distance to Alpha Centauri), but 
$\sim\! 200~{\rm nas\,yr}^{-1}$ at 1 kpc. Changes in proper motions between the Hipparcos catalogue and Gaia DR2 measurement have been used to measure 
$\mathcal{O}({\rm mas})$ accelerations of a handful of G- and K-dwarf stars and so detect the presence of brown, white, and M-dwarf companions (Brandt Reference Brandt2018; Brandt, Dupuy & Bowler Reference Brandt, Dupuy and Bowler2018).
These complementary techniques could be linked by constructing models for each exoplanetary system, potentially with additional priors from planetary formation and dynamics models, forward modelling into the data space of each experiment, and simultaneously fitting to all available data. Such a fit would yield a measure of the maximum radial acceleration attributable to unseen exoplanets, and the residual acceleration would be attributable to the Milky Way potential.
5.2.3. Multi-star techniques
The most obvious approach is to combine observations of multiple stars, under the assumption that the alignments of their exoplanetary orbital planes and phases, and hence the resultant acceleration vectors, are random and uncorrelated, whereas the acceleration from the Milky Way potential is coherent over large distances. A rudimentary implementation would be to simply stack radial velocity measurements for adjacent stars, in combination with models for the likely exoplanetary background.
A more sophisticated route would be to add a common acceleration term to the multi-experiment models described above, and conduct a simultaneous fit across multiple stellar systems. A correlated residual acceleration across many stars would be more likely to come from the Milky Way acceleration field than a conspiracy of perfectly aligned exoplanets sitting below the detection threshold of the complementary experiments. We save a full treatment of this topic for future work, but note that the next decade(s) will see an explosion of information about the global properties of exoplanetary systems.
6. Galactic gravitational field
Having considered the challenges one would encounter in performing these measurements, we finally turn to examining specific investigations of the galactic gravitational field that would be facilitated by high-precision radial velocity measurements. Beyond the possibility of gaining a new set of constrains on models of the galactic gravitational field and the underlying mass distribution in the Milky Way, we discuss two opportunities; searches for localised over-densities in the (assumed) dark matter distribution and tests of modified gravitational dynamics.
6.1. Subhalos and over-densities
Many scenarios admit the possibility of large inhomogeneities in the dark matter distribution. These include the remnants of halos associated with mergers that have yet been fully disrupted and assimilated, or exotica such as ultracompact minihalos (Ricotti & Gould Reference Ricotti and Gould2009; Bringmann, Scott & Akrami Reference Bringmann, Scott and Akrami2012) and the local inhomogeneities that might arise in models such ultralight (or fuzzy) dark matter (Hui et al. Reference Hui, Ostriker, Tremaine and Witten2017). A large, localised acceleration that was not matched by a similar anomaly is the distribution of baryonic matter would be clear evidence for dark matter—but the existence of such inhomogeneities depends on both the detailed properties of the dark matter and the dynamical history of the Sun’s neighbourhood.
To investigate this, consider an idealised spherical over-density. At a distance r from (and external to) a spherical object of mass M the acceleration it induces is
\begin{align}a = 1.4 \times 10^{-13} \frac{M}{{\rm M}_\odot} \left( \frac{1 \textrm{pc}}{r} \right)^2 ~{\rm m\,s}^{-2} .\end{align}The obvious observational strategies are simple; stars in the vicinity of an isolated over-density will have an additional, anomalous acceleration. For stars on an LoS through the centre of the over-density, the total range of the modulation is twice the value given by Equation (7); the acceleration has the opposite sign on opposite sides of the over-density.
Standard cold dark matter (CDM) cosmology suggests that the Milky Way halo includes many subhaloes. Typical masses and sizes for these objects are in the range of 
$10^7~ {\rm M}_\odot$ and a radius of 250 pc to 
$10^8~ {\rm M}_\odot$ with a radius of 625 pc (Erkal & Belokurov Reference Erkal and Belokurov2015). These would produce a 10-yr velocity change of 1.4 and 
$2.3\ {\rm cm\,s}^{-1}$, respectively, across the object, which is likely to be experimentally accessible. The variation as a function of angular direction would be much larger than that of the background field if the object is at a distance from the sun at least several times larger than its radius. Much denser objects are also frequently discussed—to give one example, Bonaca et al. (Reference Bonaca, Hogg, Price-Whelan and Conroy2018) analysed the dynamics of stellar streams, positing the existence of objects with radii of just 20 pc and masses between 
$10^{5.5}$ and 
$10^{8}~{\rm M}_\odot$, with corresponding 10-yr velocity changes from 7 to 
$\sim\!2\,000~{\rm cm\,s}^{-1}$. The latter signal is well with the reach of present-day instruments, albeit for stars in a relatively small spatial volume in the vicinity of the object.
This analysis is at the proof-of-concept level, but it is sufficient to demonstrate that plausible compact subhalo objects would be detectable via their contributions to stellar accelerations.
6.2. MOND vs. Newtonian dynamics
Theories of MOND (Milgrom Reference Milgrom1983, Reference Milgrom2010) postulate that galaxies are composed of baryonic matter and that their global dynamics are explained by modifications to Newton’s laws that become manifest at large distances and/or very small accelerations. While 
$\Lambda\text{CDM}$ (e.g. Planck Collaboration et al. 2018) is well-established as the standard model for cosmology, modified gravity theories can remain viable in the absence of a definitive identification of the physical basis for dark matter.
Because MONDian theories are designed to mimic the rotational dynamics of spiral galaxies such as the Milky Way, their predictions will typically overlap with Newtonian gravity in presence of dark matter for motion in the galactic plane. However, the baryonic matter has a very different distribution from the putative dark component—the former is roughly axisymmetric, whereas the dark halo is (approximately) spherical.
Given these differing source geometries, it would be surprising if dynamical predictions of the two theories overlapped globally, and in fact they do not. The effective gravitational potential in a MONDian Milky Way is discussed by Wu et al. (Reference Wu, Famaey, Gentile, Perets and Zhao2008) and Bienaymé et al. (Reference Bienaymé, Famaey, Wu, Zhao and Aubert2009). In these scenarios, the Milky Way acquires an apparent ‘phantom disc’ if the observed dynamics in the z direction are solved assuming that the underlying dynamics are Newtonian, and the direction of the acceleration vector likewise differs between these scenarios. This leads to differing predictions for the vertical accelerations as a function of height, allowing observations to break the degeneracy between the two scenarios. Jeans equation-based techniques have already been used to fit acceleration maps to position and velocity data, which were then used to test MONDian against Newtonian dynamics (Loebman et al. Reference Loebman2014).
One potential observational strategy would be to first select stars that were co-rotating with the Sun via a proper motion cut to reduce perspective effects, and then measure the accelerations vertically above and below the plane (z direction). This would mean 
$a_z \approx a_{\rm LoS}$, with a reduced degeneracy in proper motion accelerations, allowing for a simple test of the vertical acceleration prediction derived from the MONDian ‘phantom disc’ against the spherical CDM halo.
With full radial velocity and proper motion acceleration data, the potential of the Milky Way could be derived using the Poisson equation alone in Newtonian gravity. However, models could be fit to even sparse radial acceleration information while marginalising over the unknown proper motion accelerations, and with the use of highly constraining acceleration data these models could be very general and have few assumptions.
7. Conclusions and discussion
Our investigation of a typical model of the galactic potential shows that gravitational accelerations will lead to changes in the LoS velocities of stars on the order of centimetres per second over a decade. Given that the advent of 
${\rm cm\,s}^{-1}$ spectrographic precision is approaching, it may soon be possible to directly measure the acceleration field of the Milky Way.
That said, many challenges need to be overcome to realise such a measurement. Perhaps the biggest unknown is our ability to exclude jitter—induced by motions in stellar atmospheres—which can generate spurious radial velocity signals of 
$\mathcal{O}(100~{\rm cm\,s}^{-1})$ for dwarf stars, and 
$\mathcal{O}(1\,000~{\rm cm\,s}^{-1})$ for giants. The processes responsible for jitter modify the detailed spectra, allowing it to be distinguished from the pure Doppler shifts associated with the mean stellar velocity (e.g. Davis et al. Reference Davis, Cisewski, Dumusque, Fischer and Ford2017), and very high dispersion and resolution spectrographs coupled to large telescopes may permit further improvements. The need to exclude jitter may guide the selection of target stars. Here we face a trade-off: dwarf stars have lower levels of jitter but are inherently fainter, and observations of them may be limited to regions of the galaxy relatively close to the Sun where accelerations are expected to be lower. Conversely, giant stars have higher jitter but are brighter and thus observable further into the high-acceleration regions close to the galactic centre.
As explained in Section 4, precise astrometric observations will be crucial to measurements of the galactic acceleration field. This is due to perspective acceleration, where purely geometric effects associated with proper motion also change stellar radial velocities. However, while this apparent acceleration is larger than the gravitational acceleration, it can be calculated and subtracted if the parallax and proper motion are known to sufficient precision. Excluding perspective acceleration to the point that it is less than 1% of the characteristic Milky Way acceleration with Gaia-level astrometric precision restricts us to stars with relative transverse velocity of less than 
$55~{\rm km\,s}^{-1}$, but roughly half of all local stars survive this cut. An improvement of an order of magnitude in astrometric precision would make most local stars usable as probes of the galactic acceleration field. If astrometric precisions reach the nano-arcsecond level, perspective accelerations would be trivially removable. Further, as we saw in Section 2 such measurements would yield the transverse component of the gravitational acceleration induced by the Milky Way, allowing for direct mapping of the Milky Way’s gravitational field and mass distribution when combined with radial accelerations.
Exoplanet searches have motivated the improvement of spectrographic precisions, and ironically exoplanets will be a key source of noise when measuring the Milky Way’s acceleration field. Accelerations from short-period planets will be modulated over time, making it possible to extract them. However, large planets with long orbital periods (‘Neptunes’) would yield accelerations of similar magnitude to those of the Milky Way, with a roughly linear change in velocity over a period of a decade. We discussed several possible strategies for dealing with this issue. These include working with stars that have been accelerated to high velocities, as the process that accelerated them is likely to have removed any long-period planets, or prevented their formation in the first place. Similarly, binary systems may be less likely to host long-period planets, although they would present a greater spectroscopic challenge. Next, complementary detection techniques such as direct imaging or microlensing may reveal additional information which can better constrain the possible exoplanetary contribution to any residual change in the radial velocity. Moreover, a better general understanding of the possible ‘architectures’ of extrasolar planetary systems may further facilitate this analysis. Finally, we can combine measurements of multiple systems, on the assumption that exoplanet orbital planes and phases will be randomly oriented, while the Milky Way’s acceleration field will vary smoothly throughout space.
While substantial advances in radial velocity and astrometric precision, jitter modelling, and exoplanetary background reduction are needed to detect accelerations arising from the Milky Way’s gravitational potential, the scientific reward of these measurement would be significant. The potential and mass distribution of the Milky Way could be studied without relying on the assumptions required to extract these quantities from positions and velocities with methods such as distribution function or Jeans modelling.
There are many possible dark matter scenarios that are effectively impossible to verify directly. In fact, it could plausibly be argued that such scenarios are actually the norm, given that vast range of models that can be imagined and the absence of any truly meaningful guidance from fundamental theory as to the nature of dark matter. Consequently, despite the large, ongoing experimental effort for specific models, ‘blind tests’ could be critical to definitively testing the dark matter hypothesis. These would need to be searches for fingerprints of dark matter that are independent of its specific properties, and focus instead on the dynamics of stars and galaxies and gravitational physics at kiloparsec length scales.
Acceleration measurements could definitively establish the presence of the extended, roughly spherical dark matter halo predicted by the 
$\Lambda\text{CMD}$ model. Furthermore, the presence of dark matter subhaloes is a key signature of the hierarchical structure formation that is natural in 
$\Lambda$CMD models. These and other types of dark matter over-densities would generate additional accelerations above that of the total Milky Way potential, and would be visible in the radial acceleration distribution as crescents perpendicular to the LoS.
A direct measurement of acceleration at galactic scales would also be a powerful test of modified gravity theories. Typically, these theories find their strongest support in galactic dynamics (Famaey & McGaugh Reference Famaey and McGaugh2012), rather than phenomena with larger characteristic comoving scales, such as the cosmic microwave background and large-scale structure (e.g. Clowe, Gonzalez & Markevitch Reference Clowe, Gonzalez and Markevitch2004; Slosar, Melchiorri & Silk Reference Slosar, Melchiorri and Silk2005; Clowe et al. Reference Clowe, Bradač, Gonzalez, Markevitch, Randall, Jones and Zaritsky2006; Dodelson Reference Dodelson2011; Xu, Wang & Zhang Reference Xu, Wang and Zhang2015). However, galaxies are strongly driven by non-linear dynamics and baryonic feedback, so determining the detailed predictions of competing theories in this regime is inherently difficult. However, direct measurements of accelerations inside our galaxy have a unique capacity to test modified gravity theories in the settings where their proponents find them the most persuasive.
So what would an acceleration measurement campaign look like? The spectrograph and telescope system would need to have a limiting magnitude of 
$m_V \gtrsim 15$ in order to observe dwarf stars in regions with 
$\Delta V_{\rm LoS} \gtrsim 1 {\rm cm\,s}^{-1}$. The ESPRESSO spectrograph at the VLT has a limiting magnitude of 
$m_V = 17$, and any effort to measure the Milky Way acceleration would require a telescope of a similar class or larger. Given that multiple observations of multiple stars will be needed over a decade, a well-designed campaign would be critical to success. In the shorter term, we could imagine carefully investigating improved techniques for understanding jitter with an instrument like ESPRESSO, as a warm-up for a bigger survey. The issues associated with the extraction of planetary motion are discussed in Ravi et al. (Reference Ravi, Langellier, Phillips, Buschmann, Safdi and Walsworth2018). Carefully designed observing strategies that take full advantage of complementary information could make this sort of survey vastly more tractable than a ‘brute force’ analysis, given that accuracy gains will ultimately run as the square root of the total number of stars. Similarly, there will be trade-offs between making high cadence observations of a small number of stars, versus lower cadence observations of many, and survey design will depend strongly on specific science goals.
Our analysis has focussed on high-precision measurements of stellar velocities from which acceleration data can be extracted for individual stars with tolerable signal to noise. However, an alternative strategy would be to observe many stars with 
$\mathcal{{O}}(10-100)\ \text{cm}\ \text{s}^{-1}$ precision and rely on 
$1/\sqrt{N}$ statistics to obtain mean accelerations. Next generation multifibre spectrographs could potentially yield absolute measurements of this quality for up to one million stars, but determining an optimal observing strategy would require careful analysis. For example, planetary reflex velocities can easily be 
$10^{4}$ times larger than the intrinsic acceleration (see Figure 5) but ‘Jupiters’ can be identified and removed by making repeat observations at moderate cadence; ‘filtering’ that reduces the typical apparent radial velocity by one order of magnitude naively reduces the sample size needed to obtain a given precision by two orders of magnitude. Moreover, observations of individual stars may not be fully independent, as effects associated with the Sun’s motion may correlate velocity estimates for stars in a given spatial volume. However, this strategy could yield a coarse-grained view of the acceleration field above and below the plane, or in the directions of the galactic centre and anti-centre and warrants further careful study.
The proposal discussed here is an ambitious one, and its full realisation would take decades rather than years. However, our analysis shows that obtaining direct measurements of accelerations and net gravitational forces at galactic distances is achievable with a reasonable extrapolation of present-day technology, and would create new and exciting opportunities for galactic astronomy, astrophysics, and fundamental physics.
Acknowledgements
HS acknowledges support from the MINECO (Spanish Ministry of Economy)-FEDER project through grant MDM-2014-0369 of ICCUB (Unidad de Excelencia ‘Maria de Maeztu’). RE acknowledges support from the Marsden Fund of the Royal Society of New Zealand. HS would like to thank Monica Silverwood and Rob Dixon for their hospitality during the early development of this project. We thank the Referee for their helpful comments, and the authors of Ravi et al. (2018) for productive discussions. We thank Tristan Cantat-Gaudin, JJ Eldridge, Claus Fabricius, Anna Ferré Mateu, Francesca Figueras, Debra Fischer, Mark Gieles, Kathryn Johnston, Jordi Miralda-Escude, Roger Mor, Pau Ramos, Nicholas Rattenbury, Aakash Ravi, Justin Read, Mercé Romero-Gómez, Daniel del Ser Badia, and Willem van Straten for useful conversations and comments. This work has made use of data from the European Space Agency (ESA) mission Gaia (https://www.cosmos.esa.int/gaia), processed by the Gaia Data Processing and Analysis Consortium (DPAC, https://www.cosmos.esa.int/web/gaia/dpac/consortium). Funding for the DPAC has been provided by national institutions, in particular the institutions participating in the Gaia Multilateral Agreement.
Appendix A. Spectrographic precisions through history
The following spectrographs and measurements were included in Figure 2: Campbell & Walker (Campbell, Walker & Yang Reference Campbell, Walker and Yang1988); Mt Wilson/Hamilton 1989 (Marcy & Benitz Reference Marcy and Benitz1989); CfA/CORAVEL (Latham et al. Reference Latham, Mazeh, Stefanik, Mayor and Burki1989); Hamilton 1992 (Marcy & Butler Reference Marcy and Butler1992); McDonald Observatory Planetary Search (Hatzes & Cochran Reference Hatzes and Cochran1993); ELODIE (Mayor & Queloz Reference Mayor and Queloz1995; Baranne et al. Reference Baranne1996); Hamilton 1996 (Butler et al. Reference Butler, Marcy, Williams, McCarthy, Dosanjh and Vogt1996); HIRES (Vogt et al. Reference Vogt, Crawford and Craine1994); UCLES-AAPS (Butler et al. Reference Butler, Tinney, Marcy, Jones, Penny and Apps2001); Tull (Tull et al. Reference Tull, MacQueen, Sneden and Lambert1995); HRS (Tull Reference Tull and D’Odorico1998); HERCULES (Hearnshaw et al. Reference Hearnshaw, Barnes, Kershaw, Frost, Graham, Ritchie and Nankivell2002; Endl et al. Reference Endl2015); HARPS (Mayor et al. Reference Mayor2003; Pepe et al. Reference Pepe2002); PFS (Crane, Shectman & Butler Reference Crane, Shectman and Butler2006); SOPHIE+ (Perruchot et al. Reference Perruchot2008; Bouchy et al. Reference Bouchy2013); CHIRON (Tokovinin et al. Reference Tokovinin2013); HARPS-N (Cosentino et al. Reference Cosentino2014; Covino et al. Reference Covino2013); PARAS (Chakraborty et al. Reference Chakraborty2010, Reference Chakraborty2014); Levy CPS (Radovan et al. Reference Radovan, Cabak, Laiterman, Lockwood and Vogt2010); Levy LCPS (Radovan et al. Reference Radovan, Cabak, Laiterman, Lockwood and Vogt2010; Vogt et al. Reference Vogt2014); SONG (Grundahl et al. Reference Grundahl, Christensen-Dalsgaard, Gråe Jørgensen, Frandsen, Kjeldsen and Kjærgaard Rasmussen2011); ESPRESSO (Pepe et al. Reference Pepe2010, Reference Pepe2014; ESPRESSO Collaboration 2018); EXPRES (Fischer et al. Reference Fischer, Jurgenson, McCracken, Sawyer, Blackman and Szymkowiak2017); G-CLEF (Szentgyorgyi et al. Reference Szentgyorgyi2012, Reference Szentgyorgyi2018); CODEX, (Pasquini et al. Reference Pasquini2010b,a).
