OPEN ACCESS

<figure><img src="image_1.png" /></figure>

CrossMark

# All-sky Kinematics of the Distant Halo: The Reflex Response to the LMC

Vedant Chandra<fnref n="1" /> ID, Rohan P. Naidu<fnref n="2" />,<fnref n="14" /> ID, Charlie Conroy<fnref n="1" /> ID, Nicolas Garavito-Camargo<fnref n="3" /> ID, Chervin Laporte<fnref n="4" />,<fnref n="5" />,<fnref n="6" /> ID, Ana Bonaca<fnref n="7" /> ID, Phillip A. Cargile<fnref n="1" /> ID, Emily Cunningham<fnref n="8" />,<fnref n="14" /> ID, Jiwon Jesse Han<fnref n="1" /> ID, Benjamin D. Johnson<fnref n="1" /> ID, Hans-Walter Rix<fnref n="9" /> ID, Yuan-Sen Ting (丁源森)<fnref n="10" />,¹¹,<fnref n="12" /> ID, Rebecca Woody<fnref n="1" /> ID, and Dennis Zaritsky<fnref n="13" /> ID

Received 2024 June 3; revised 2025 May 14; accepted 2025 May 18; published 2025 July 21

## Abstract

The infall of the Large Magellanic Cloud (LMC) is predicted to displace the inner Milky Way (MW), imprinting an apparent “reflex motion” on the observed velocities of distant halo stars. We construct the largest all-sky spectroscopic data set of luminous red giant stars from 50–160 kpc, including a new survey of the southern celestial hemisphere. We fit the full 6D kinematics of our data to measure the amplitude and direction of the inner MW’s motion toward the outer halo. The observed velocity grows with distance such that, relative to halo stars at 100 kpc, the inner MW is lurching at $\approx 40 \text{ km s}^{-1}$ toward a recent location along the LMC’s past orbit. Our measurements align with N-body simulations of the halo’s response to a $1.8 \times 10^{11} M_{\odot}$ LMC on first infall, suggesting that the LMC is at least 15% as massive as the MW. Our findings highlight the dramatic disequilibrium of the MW outskirts, and will enable more accurate measurements of the total mass of our Galaxy.

*Unified Astronomy Thesaurus concepts: Milky Way stellar halo (1060); Stellar kinematics (1608); Large Magellanic Cloud (903); Milky Way dynamics (1051)*

*Materials only available in the online version of record: animation*

## 1. Introduction

The LMC and SMC (hereafter the “Clouds”) are by far the Milky Way’s (MW’s) most prominent satellite galaxies, and were likely the first external galaxies discovered (e.g., M. Dennefeld2020). There is growing evidence that the Clouds form a binary pair—with the LMC being $\approx 10 \times$ more massive than the SMC—and are on their first orbital infall toward the Milky Way (N. Kallivayalil et al. 2006; G. Besla et al. 2007, 2010; N. Kallivayalil et al. 2013; but see E. Vasiliev 2024 about the plausibility of a second-infall scenario). Multiple independent measurements estimate the total mass of the LMC to be between 1.0 and $2.5 \times 10^{11} M_{\odot}$ (J. Peñarrubia et al. 2016; D. Erkal et al. 2019b, 2020; N. Shipp et al. 2021; E. Vasiliev et al. 2021; L. Correa Magnus & E. Vasiliev 2022; K. J. Fushimi et al. 2024; S. E. Koposov et al. 2023; E. Vasiliev 2023; L. L. Watkins et al. 2024), 10%–25% the mass of the Milky Way itself (e.g., J. Bland-Hawthorn & O. Gerhard 2016; W. Wang et al. 2020; J. Shen et al. 2022).

Given its large mass, the LMC is expected to significantly perturb our own Galaxy (see E. Vasiliev 2023 for a thorough review). On a global scale, the LMC exerts a gravitational acceleration on various MW components as it approaches, and is predicted to shift the center-of-mass (COM) of the Galaxy significantly. While the disk and inner halo can quickly react to the changing COM by moving toward the LMC, stars in the outer halo are left behind due to their large distances and long orbital timescales ($\gtrsim 1$ Gyr). This results in an apparent “reflex motion” imprinted in the velocities of distant halo stars, which should become more prominent at larger distances (M. D. Weinberg 1995; F. A. Gómez et al. 2015; D. Erkal et al. 2019a; A. Kravtsov & S. Winney 2024; Y. Sheng et al. 2024). The dominant observational effect is that distant halo stars should have apparent “upward” (positive z) velocities in the Galactocentric frame (C. F. Laporte et al. 2018; D. Erkal et al. 2019a; N. Garavito-Camargo et al. 2019, 2021; E. C. Cunningham et al. 2020; M. S. Petersen & J. Peñarrubia 2020; Y.-N. Cao et al. 2023). D. Erkal et al. (2020) demonstrated that this effect differentially distorts the MW halo, and must be accounted for in dynamical analyses; neglecting the influence of the LMC can inflate mass estimates of the MW by up to 50% (see also L. Correa Magnus & E. Vasiliev 2022; K. Chamberlain et al. 2023; A. Kravtsov & S. Winney 2024). On a more local scale, the LMC affects dark matter (DM) in the Milky Way’s halo, imprinting a dynamical

<figure><img src="image_2.png" /></figure>

<footnote n="1">¹Center for Astrophysics | Harvard & Smithsonian, 60 Garden Street, Cambridge, MA 02138, USA; vedant.chandra@cfa.harvard.edu</footnote>

<footnote n="2">²MIT Kavli Institute for Astrophysics and Space Research, 77 Massachusetts Avenue, Cambridge, MA 02139, USA</footnote>

<footnote n="3">³Center for Computational Astrophysics, Flatiron Institute, 162 5th Avenue, New York, NY 10010, USA</footnote>

<footnote n="4">⁴LIRA, Observatoire de Paris, Université PSL, Sorbonne Université, Université Paris Cité, CY Cergy Paris Université, CNRS, 92190 Meudon, France</footnote>

<footnote n="5">⁵Institut de Ciencies del Cosmos (ICCUB), Universitat de Barcelona (IEEC-UB), Martí i Franquès 1, E-08028 Barcelona, Spain</footnote>

<footnote n="6">⁶Kavli IPMU (WPI), UTIAS, The University of Tokyo, Kashiwa, Chiba 277-8583, Japan</footnote>

<footnote n="7">⁷The Observatories of the Carnegie Institution for Science, 813 Santa Barbara Street, Pasadena, CA 91101, USA</footnote>

<footnote n="8">⁸Department of Astronomy, Columbia University, 550 West 120th Street, New York, NY, 10027, USA</footnote>

<footnote n="9">⁹Max-Planck-Institut für Astronomie, Königstuhl 17, D-69117 Heidelberg, Germany</footnote>

<footnote n="10">¹⁰Research School of Astronomy & Astrophysics, Australian National University, Cotter Road, Weston, ACT 2611, Australia
¹¹School of Computing, Australian National University, Acton ACT 2601, Australia</footnote>

<footnote n="12">¹²Department of Astronomy, The Ohio State University, Columbus, OH 43210, USA</footnote>

<footnote n="13">¹³Steward Observatory and Department of Astronomy, University of Arizona, Tucson, AZ 85721, USA</footnote>

<footnote n="14">¹⁴NASA Hubble Fellow.</footnote>

<footnote>Original content from this work may be used under the terms of the Creative Commons Attribution 4.0 licence. Any further distribution of this work must maintain attribution to the author(s) and the title of the work, journal citation and DOI.</footnote>
friction “wake” in the DM density along its past orbit (S. Chandrasekhar 1943; N. Garavito-Camargo et al. 2019; T. Tamfal et al. 2021; H. R. Foote et al. 2023). These effects are predicted to leave a prominent signature in the density and kinematics of MW halo stars, with the perturbations being most prominent in the outer halo beyond the LMC at $\gtrsim$50 kpc.

In recent years, excellent progress has been made toward observationally detecting the influence of the LMC on MW halo stars. V. Belokurov et al. (2019) used RR Lyrae and blue horizontal branch (BHB) stars to study a prominent overdensity in the southern sky, the Pisces Plume. This overdensity is elongated along the past orbit of the Clouds, and is systematically offset toward negative Galactocentric velocities, as predicted by their simulations of the LMC’s wake. C. Conroy et al. (2021) constructed an all-sky sample of luminous red giant stars, and observed a wake-like overdensity in the southern Galactic hemisphere, along with a more diffuse overdensity in the northern hemisphere. Both of these features are predicted by the simulations of N. Garavito-Camargo et al. (2019, hereafter GC19), although the density contrast in the data is twice as strong as the simulations. J. A. S. Amarante et al. (2024) recently characterized the Pisces Plume with a large sample of BHB stars, although they do not detect a corresponding diffuse overdensity in the north.

Attempts to measure the LMC’s effects on the “bulk” MW halo can be significantly confounded by the presence of undetected coherent substructure in the outer halo. The most prominent of these structures in the outer halo is the Sagittarius stream (Sgr), although 6D kinematics provide a powerful way to excise Sgr stars on the basis of their angular momenta (B. D. Johnson et al. 2020; J. Peñarrubia & M. S. Petersen 2021). The LMC and SMC are themselves predicted to contribute stripped stars to the Galactic halo, forming a stellar counterpart to the well-known gaseous Magellanic Stream (MS; e.g., G. Besla et al. 2013; M. S. Petersen et al. 2022). D. Zaritsky et al. (2020) discovered coherent debris in the H3 Survey (C. Conroy et al. 2019) at $\approx$50 kpc, and argued that these might have a Magellanic origin. V. Chandra et al. (2023b) subsequently discovered a population of stars that appear to be stripped from the Clouds (the Magellanic Stellar Stream, or MSS), spanning the entire extent of the trailing MS from 50–100 kpc. The MSS and other debris could partially constitute the Pisces Plume overdensity, but there remains a wake-like structure in the outer halo once the MSS is removed (V. Chandra et al. 2023b; J. A. S. Amarante et al. 2024).

Given the overlapping substructures present in the halo, full 6D kinematics are vital to disentangle coherent structures from the bulk MW halo. Recently, wide-field spectroscopic surveys—combined with astrometry from the Gaia space observatory—have enabled the kinematic detection of the LMC’s global influence on the MW halo. D. Erkal et al. (2021) compiled Galactocentric radial velocities of halo stars beyond 30 kpc, and found that stars in the northern Galactic hemisphere are systematically redshifted relative to those in the south, as expected for the LMC-induced reflex motion. The amplitude of their measured velocity dipole was broadly consistent with a $1.5 \times 10^{11} M_{\odot}$ LMC. Stellar streams are also emerging as powerful tools to map the gravitational field of the Galaxy, including the LMC’s influence (e.g., D. Erkal et al. 2019b; S. E. Koposov et al. 2023; S. Lilleengen et al. 2023; R. A. N. Brooks et al. 2025). E. Vasiliev et al. (2021) comprehensively modeled the formation of the Sgr stream in the presence of the LMC, finding that the motion of the inner

MW toward the LMC is a crucial ingredient to reproduce the observed properties of Sgr.

M. S. Petersen & J. Peñarrubia (2021, hereafter PP21) fit a physical model of the reflex velocity to the 6D phase-space positions of RGB and BHB stars beyond 40 kpc, and measured a global velocity vector $v_{\text{travel}} = 32 \pm 4$ km s$^{-1}$ with an apex (maximum velocity) direction that points toward a location along the past orbit of the LMC. While the amplitude of their measured signal is consistent with simulations of a $\gtrsim 10^{11} M_{\odot}$ LMC, the measured apex direction is somewhat different from most simulations of the MW–LMC interaction (E. Vasiliev 2023). Whereas PP21 measured an apex direction in the opposite southern Galactic quadrant as the LMC today ($\ell_{\text{apex}} > 0^{\circ}$), many simulations predict that the velocity dipole should point toward a more recent point on the LMC’s orbit ($\ell_{\text{apex}} < 0^{\circ}$; e.g., N. Garavito-Camargo et al. 2019; E. Vasiliev et al. 2021). More recent work by R. Yaaqib et al. (2024, hereafter YPP24) found an increasing $v_{\text{travel}}$ velocity toward larger distances using a sample spanning $20 < r_{\text{gal}}/\text{kpc} < 60$, with an apex direction that varies by tens of degrees as a function of distance.

One major limitation of existing outer halo data sets is that they lack coverage at southern declinations, since all spectroscopic surveys probing the distant halo (e.g., the SEGUE survey; B. Yanny et al. 2009, used by PP21 and YPP24) have been conducted from the northern hemisphere. The H3 Spectroscopic Survey has significantly increased the number of distant halo stars with full 6D information, but this too is limited to declinations $\delta \gtrsim -20^{\circ}$ (C. Conroy et al. 2019). Kinematic samples of the outer halo therefore have a large hole in the Galactic south, in precisely the quadrant where the LMC resides and exerts its maximal dynamical influence. It is therefore desirable to construct an all-sky data set of outer halo stars, to obtain a global perspective on the Galaxy’s kinematic response to the LMC.

Over the past 2 yr, we have been conducting a tailored spectroscopic survey of the most distant giants in the southern hemisphere, using the Magellan Echellette Spectrograph (MagE; J. L. Marshall et al. 2008) on the 6.5 m Magellan Baade telescope. Our survey has currently observed $\approx 290$ giants with heliocentric distances $d_{\text{helio}} > 50$ kpc, including $\approx 100$ giants beyond $d_{\text{helio}} > 100$ kpc. In combination with past northern surveys like SEGUE and H3, this is the premier all-sky data set to kinematically map out the farthest outskirts of our Galaxy. Previous results from our survey include identifying distant remnants of our Galaxy’s last major merger (V. Chandra et al. 2023a) and the discovery of the Magellanic Stellar Stream (V. Chandra et al. 2023b).

In this work, we combine our MagE survey data with existing data from the H3 and SEGUE surveys to construct the largest all-sky data set of red giants with full 6D kinematics beyond 40 kpc, with all spectra homogeneously analyzed using the same pipeline (see Section 2). We illustrate the all-sky kinematics of the outer halo in Section 3, demonstrating the dramatic dynamical disequilibrium of the MW outskirts. Subsequently in Section 4 we construct a framework to measure the LMC-induced reflex motion with a distance-dependent continuity model for the 6D kinematics of our sample. We present our fitted model in Section 4.2 and compare our results to N-body simulations and past measurements. We simultaneously measure the reflex-corrected bulk kinematic properties of the outer halo in Section 4.3, and
examine effects of the local dynamical wake of the LMC in Section 4.4. The broad implications of our findings are discussed in Section 5, and our results are summarized in Section 6.

## 2. Data

This section describes the construction of an all-sky data set of RGB stars with full 6D phase-space information beyond 40 kpc. Throughout this work, we adopt a right-handed Galactocentric frame with a solar position $\mathbf{x}_{\odot} = (-8.12, 0.00, 0.02)$ kpc, and solar velocity $\mathbf{v}_{\odot} = (12.9, 245.6, 7.8)$ km s$^{-1}$ (M. J. Reid & A. Brunthaler 2004; R. Drimmel & E. Poggio 2018; GRAVITY Collaboration et al. 2018). Coordinate transformation are carried out using *astropy* (Astropy Collaboration et al. 2013, 2018, 2022) and *gala* (A. M. Price-Whelan 2017; A. Price-Whelan et al. 2020).

### 2.1. MagE Survey

Since 2022, we have been executing a tailored spectroscopic survey of luminous red giants in the outer halo. The selection procedure for our parent sample of red giant stars is detailed in V. Chandra et al. (2023a), and our spectroscopic survey is further described in V. Chandra et al. (2023b). Briefly, we utilize Gaia DR3 astrometry and unWISE photometry to filter out nearby dwarf stars and isolate giants off the Galactic plane ($|b| > 20^{\circ}$) with isochrone-predicted distances $\gtrsim 80$ kpc (A. Mainzer et al. 2014; E. F. Schlafly et al. 2019; A. G. Brown et al. 2021; L. Lindegren et al. 2021; Gaia Collaboration et al. 2023).

We obtained spectroscopic observations of a subset of the above selected giants with MagE (J. L. Marshall et al. 2008) on the 6.5 m Magellan Baade Telescope at Las Campanas Observatory (PIs: Chandra & Naidu). The stars are typically between $17 \lesssim G \lesssim 18.6$, and we utilize magnitude-dependent exposure times ranging from 10–40 minutes per star to achieve signal-to-noise ratio (SNR) $\approx 10$ pixel$^{-1}$ at 5100 Å. Spectra are reduced with a fully automated pipeline$^{15}$ built around the *PypeIt* utility (J. X. Prochaska et al. 2020).

Stellar parameters are estimated with the Bayesian MINE-Sweeper code (P. A. Cargile et al. 2020). We fit a region of the MagE spectra (from 4800–5500 Å) that contains the Mg I triplet, where our line lists are best-calibrated. We also include archival optical-infrared broadband photometry in the likelihood, along with the Gaia parallax (M. Fukugita et al. 1996; J. E. Gunn et al. 1998; M. F. Skrutskie et al. 2006; A. Mainzer et al. 2014; K. C. Chambers et al. 2016; A. G. Brown et al. 2021; Gaia Collaboration et al. 2023). MINE-Sweeper compares these observables to synthetic Kurucz model spectra, and additionally constrains the models to lie on MIST isochrones (R. L. Kurucz 1970; R. L. Kurucz & E. H. Avrett 1981; J. Choi et al. 2016). The posterior distribution of stellar parameters is sampled with *dynesty* (J. S. Speagle 2020), producing measurements of the radial velocity $v_r$, effective temperature $T_{\text{eff}}$, surface gravity log g, metallicity [Fe/H], [α/Fe] abundance, and heliocentric distance. Repeat observations of bright radial velocity standard stars—HIP4148 and HIP22787—over multiple nights indicate that our systematic RV accuracy floor is $\approx 1$ km s$^{-1}$. For more details on the MagE outer halo survey, we refer to V. Chandra et al. (2023a, 2023b).

<figure><img src="image_3.png" /><figcaption>**Figure 1.** Validating the spectroscopic distances of stars in our MagE survey beyond 30 kpc. Black points show giant stars from the H3 Survey, which were fitted with a similar analysis pipeline as our MagE survey. The two colored points show mean distances of observed stars in the Crater 2 and Sextans dwarf galaxies, with their literature distances sourced from G. Torrealba et al. (2016) and R. R. Muñoz et al. (2018), respectively. The lower panel shows the fractional residual, the difference between the MagE and literature distance divided by the literature distance.</figcaption></figure>

Figure 1 shows a cross-validation of the spectroscopic distance scale of our MagE survey. Apart from observing stars in common with the H3 Survey (C. Conroy et al. 2019), we observed giants from our parent giant catalog belonging to the Sextans (86 ± 1 kpc; M. J. Irwin et al. 1990; R. R. Muñoz et al. 2018) and Crater 2 (117 ± 1 kpc; G. Torrealba et al. 2016) dwarf galaxies. Figure 1 illustrates that we obtain reliable distance estimates out to $\gtrsim 100$ kpc. The typical statistical distance uncertainties are smaller in the MagE data set than H3—for stars in common—since we aimed for higher SNR in the MagE survey, leading to more precise metallicities and consequently more precise isochrone distances. Using the scatter in the relation shown in Figure 1 as a guide, we adopt an uncertainty floor of 10% on our distance measurements. This is consistent with the standard deviation of the measured distances in Crater 2 and Sextans, which are both $\approx 10\%$.

There is a 5% systematic offset between the H3 and MagE spectroscopic distance scales (see bottom panel of Figure 1), most likely caused by the MagE metallicity scale being $\approx 0.1$ dex lower than H3. Although we have not yet determined the cause of this offset, it is smaller than the typical statistical uncertainties on these distances. We have verified that artificially shrinking the MagE distances by 5% to the H3 scale has a negligible effect on the results presented here.

<footnote n="15">15 https://github.com/vedantchandra/merlin/releases/tag/v0.2</footnote>
Furthermore, there is a ≈3 km s⁻¹ mean offset between MagE and H3 RVs of the same stars. We have likewise verified that shifting the MagE RVs to the H3 scale has a negligible effect on our results.

As of 2024 May, we have observed 400 stars in this survey, of which 291 are spectroscopically confirmed to lie at heliocentric distances beyond 50 kpc, and 106 are beyond 100 kpc. This is already the largest data set of stars with precise metallicities and 6D phase-space information beyond 100 kpc, and by far the largest data set of outer halo stars beyond ≈50 kpc observed from the southern hemisphere. Our sample selection has a very low contamination rate from dwarfs (≈10%). All data from this survey will eventually be made publicly available, once our primary scientific projects have been completed.

## 2.2. H3 and SEGUE Data

The H3 Spectroscopic Survey (C. Conroy et al. 2019) has been conducting a spectroscopic survey of halo stars with the Hectochelle instrument (A. Szentgyorgyi et al. 2011) on the 6.5 m MMT telescope at the Whipple Observatory in Arizona. We select stars observed by H3 up to 2024 January that have reliable stellar parameter fits from *MINESweeper*, and are not associated with cold substructures (stars with on-sky positions and velocities consistent with membership in dwarf galaxies, star clusters, or known stellar streams). Before applying any further selection cuts, there are ≈300 giants beyond $d_{helio} > 50$ kpc from H3. Typical RV uncertainties for these data are of the order of 0.5–1 km s⁻¹ in absolute terms (C. Conroy et al. 2019; P. A. Cargile et al. 2020).

The Sloan Extension for Galactic Understanding and Exploration (SEGUE; B. Yanny et al. 2009) survey observed a quarter of a million stars with the low-resolution BOSS spectrograph as a part of the Sloan Digital Sky Survey (SDSS; D. G. York et al. 2000; K. N. Abazajian et al. 2009). We independently fit these spectra with the *MINESweeper* routine, from which ≈300 stars with reliable parameters and distances $d_{helio} > 50$ kpc are used in this work. Based on a comparison to literature parameters, we determine that SEGUE stars with SNR < 10 per pixel do not produce reliable parameters, so we excise them from the data set. Due to the low resolution of BOSS, the SEGUE stars have typical RV uncertainties ≈5 km s⁻¹.

The data set used here is not an exhaustive compilation of distant halo stars; rather, we have prioritized building a pure sample of stars with homogeneous stellar parameters derived with the *MINESweeper* pipeline. In cases where the same star has been observed by multiple surveys, we follow an order of precedence to select a single data point: MagE > H3 > SEGUE. This order is chosen to maximize RV accuracy; although the H3 Survey operates at a higher resolution than our MagE survey, the distant giants studied here typically have much higher SNR in our MagE data. Future work could incorporate the LAMOST survey, which has not been analyzed with *MINESweeper* yet (X.-Q. Cui et al. 2012; J. Li et al. 2023; L. Zhang et al. 2023).

## 2.3. Sample Selection

Put together, our MagE+H3+SEGUE data set contains ≈2600 stars with $40 < r_{gal}/kpc < 160$. However, in order to reliably measure the halo's response to the LMC, it is vital to select a high-fidelity subset of these data with known

substructures removed. To start with, we exclude ≈200 stars with reported RV uncertainties >15 km s⁻¹ and fractional distance uncertainties >25%.

The most prominent stellar structure at these distances is the Sagittarius Stream (Sgr; e.g., S. R. Majewski et al. 2003; E. Vasiliev et al. 2021). We kinematically exclude Sgr stars using the angular momentum cuts proposed by B. D. Johnson et al. 2020, removing stars with $|b_{Sgr}| < 25°$ and $L_Y < -2.5 - 0.3 \times L_Z$, where the angular momenta are in units of kpc km s⁻¹. Within 25° of the Sgr orbital plane (as defined by D. R. Law & S. R. Majewski 2010a), an additional excision is performed. Using the models of E. Vasiliev et al. (2021) as a guide at each location of the stream, stars with distances within 10 kpc of the stream and radial velocities within 50 km s⁻¹ of the stream (at the 3σ level given reported uncertainties) are excised.

Stars with unphysical radial or tangential velocities >1000 km s⁻¹ are removed, as are stars within 1° of known dwarf galaxies and globular clusters (A. W. McConnachie 2012; E. Vasiliev & H. Baumgardt 2021). The 13 stars in our MagE data set identified by V. Chandra et al. (2023b) as belonging to the Magellanic Stellar Stream are excised. Finally, although stars in our MagE survey are selected to lie far from the LMC and SMC (see Figure 2), we search for any stars within 50° of the LMC having proper motions and RVs consistent with belonging in an extended LMC halo. No stars confidently associated with the LMC are found.

Our final sample consists of 851 field halo stars from SEGUE, H3, and MagE with $40 < r_{gal}/kpc < 160$. The top panel of Figure 2 shows their on-sky spatial distribution in the outer halo, and the bottom panel shows their cumulative distance distributions. The bottom panel of Figure 2 illustrates the high efficiency of our MagE survey in pushing to the farthest Galactic outskirts, far outnumbering the number of $r_{gal} \gtrsim 60$ kpc stars compared to the H3 and SEGUE surveys. We adopt a systematic uncertainty floor of 1 km s⁻¹ on RVs and 10% on distances. Since most stars have reported catalog uncertainties smaller than this, the median uncertainties in our final catalog are 1 km s⁻¹ for RVs and 10% for distances. The median absolute magnitude of stars in our sample is $M_G = -1.3$, the median metallicity is [Fe/H] = -1.3, and the median alpha abundance is [α/Fe] = +0.2.

# 3. The Velocity Structure of the Outer Halo

## 3.1. Radial Velocities

We begin by considering the overall radial velocity structure of the halo. Figure 3 shows the on-sky distribution of Galactocentric radial velocities $v_{GSR}$ from our data set, for stars with $r_{gal} = 50-100$ kpc. We use this limited distance range for this figure instead of our full 40–160 kpc sample to minimize the effect of the different selection functions of the surveys combined here. As Figure 2 illustrates, most stars beyond 100 kpc are from our MagE survey and lie only in the southern hemisphere. The $r_{gal} = 50-100$ kpc range ensures relatively even depth across the sky. In Appendix D, we further investigate the distance-dependence of these smoothed velocity maps, finding consistent behavior across all distances.

The top panel of Figure 3 shows $v_{GSR}$ measurements for individual stars, and the middle panel shows the smoothed mean velocity map. At each pixel on this map, the mean velocity of surrounding stars is computed, weighted by an
<figure><img src="image_4.png" /><figcaption>**Figure 2.** Top panel: Galactic map of our cleaned all-sky data set of red giant stars beyond 40 kpc, colored by Galactocentric distance. Stars with spectra from our MagE survey are outlined in black. Bottom panel: cumulative number of stars in our sample beyond a given Galactocentric distance, split by spectroscopic survey. Our MagE survey fills out the southern celestial hemisphere, and pushes to larger distances than previous surveys.</figcaption></figure>

angular Gaussian kernel with a radius of 45°. This effectively produces a smoothed velocity map in a continuous rather than binned manner. The maps are qualitatively unchanged but noisier if we use a smaller 15° or 30° kernel. The sky-averaged mean velocity is subtracted out and indicated in the colorbar label. The present-day location and past orbit of the LMC are overlaid, utilizing the orbital trajectory simulated by GC19 (model LMC3 in their simulated suite, with a halo mass $1.8 \times 10^{11} M_{\odot}$).

It is immediately apparent from Figure 3 that the outer halo has a structured velocity distribution with variations across the sky. The Galactocentric radial velocity exhibits a net negative ≈−12 km s⁻¹ signal. Furthermore, there is an unambiguous dipole signature in $v_{GSR}$, with the southern halo systematically blueshifted, and the northern halo redshifted. Although this feature is most apparent in the smoothed map in the middle panel of Figure 3, it is visible even in the raw data in the top panel. This $v_{GSR}$ dipole does not appear to be aligned purely in the north-south direction, but is rather tilted along an axis that approximately points toward the LMC. We emphasize that these are subtle mean trends in $v_{GSR}$ brought out by averaging over large on-sky scales: the intrinsic radial velocity dispersion is 70–110 km s⁻¹ depending on the sky location. In Appendix D, we demonstrate that this observation is robust to the distance

<figure><img src="image_5.png" /><figcaption>**Figure 3.** Comparing the radial velocity structure of our MW data to a numerical simulation of the halo's response to the LMC. Top panel: individual $v_{GSR}$ measurements for stars in our data set from $r_{gal} = 50–100$ kpc. Dashed lines demarcate the quadrants used for the tomographic $v_{GSR}$ measurements shown in Figure 4. Middle panel: average $v_{GSR}$ map of our data, smoothed with a 45° kernel. Bottom panel: smoothed $v_{GSR}$ computed using the MW+LMC3 simulation from GC19, which assumed a $1.8 \times 10^{11} M_{\odot}$ halo mass for the LMC.</figcaption></figure>

range considered, and rule out this effect being an artifact due to the more distant coverage of the MagE data.

From this $v_{GSR}$ map alone, the reflex velocity imparted by the LMC is clearly visible. Furthermore, the orientation of the signal in the data is consistent with numerical simulations of the interaction. The bottom panel of Figure 3 shows the mean $v_{GSR}$ map from the MW+LMC simulations of GC19, specifically their LMC3 model with a $1.8 \times 10^{11} M_{\odot}$ LMC virial mass. The data are quite similar to the results from this simulation—particularly the tilted axis of the dipole—although
<figure><img src="image_6.png" /><figcaption>**Figure 4.** Mean $v_{\text{GSR}}$ as a function of Galactocentric distance for the quadrants demarcated in the top panel of Figure 3: the quadrant containing the present-day location of the LMC, and the one diagonally opposing it. $1\sigma$ error bars are obtained by bootstrapping, and the number of stars in each bin is indicated on the markers. Dashed lines show the same measurements made with the MW+LMC simulations of GC19, with the line width being proportional to the LMC halo mass in each simulation: from thinnest to thickest corresponding to $[0.8, 1.0, 1.8, 2.5] \times 10^{11} M_{\odot}$.</figcaption></figure>

the observed dipole is slightly stronger than the simulated one. The simulation also exhibits a net negative sky-averaged $v_{\text{GSR}} \approx -10 \text{ km s}^{-1}$, hinting that this signal may be LMC-induced as well (as recently argued by YPP24).

Dashed lines in the top panel of Figure 3 show the selection of two diagonally opposing quadrants that contain the regions of maximum mean radial velocity in both the data and simulations. The mean radial velocity in these quadrants is shown as a function of Galactocentric distance $r_{\text{gal}}$ in Figure 4: the quadrant containing the LMC in orange, and the diagonally opposing quadrant in blue. In each uncorrelated $r_{\text{gal}}$ bin, the $3\sigma$-clipped mean $v_{\text{GSR}}$ is shown, with $1\sigma$ error bars computed via bootstrapping. Note that the measurement in the southern quadrant goes much farther than that in the north, due to the unprecedented depth of our MagE survey (see Figure 2). An identical measurement is made on star particles in the MW+LMC simulation of GC19, and from thinnest to thickest the dashed lines correspond to LMC halo masses of $[0.8, 1.0, 1.8, 2.5] \times 10^{11} M_{\odot}$. These stars particles are “tagged on” to the dark matter particles simulated by GC19, using the formalism of C. F. P. Laporte et al. (2013).

Figure 4 emphasizes what is apparent in the maps of Figure 3: these quadrants have drastically different radial velocities over a wide range of distances. Note that the bins used to compute the average velocities in Figure 4 are entirely independent. The 80–100 kpc bin in the LMC quadrant exhibits markedly negative radial velocities and warrants future investigation. This could be due to unresolved substructure at these distances, for example, stars stripped from the LMC or SMC (e.g., V. Belokurov & S. E. Koposov 2016; V. Chandra et al. 2023b). In Appendix D, we split the all-sky $v_{\text{GSR}}$ maps by distance, showing that the overall behavior of these maps is consistent across distance bins.

<figure><img src="image_7.png" /><figcaption>**Figure 5.** All-sky mean tangential velocity along the Galactic latitude coordinate as a function of Galactocentric distance, with $1\sigma$ error bars from bootstrapping. The recent measurement from Y. Sheng et al. (2024) using RR Lyrae stars is shown in gold. Dashed lines show the same measurements made with the MW+LMC simulations of GC19, with the line width being proportional to the LMC halo mass in each simulation: from thinnest to thickest $[0.8, 1.0, 1.8, 2.5] \times 10^{11} M_{\odot}$.</figcaption></figure>

Furthermore, we have verified that removing the 26 stars in the discrepant bin of Figure 4 does not significantly alter the results from our 6D model fit in Section 4.2. The observed behavior matches the MW+LMC simulations shown in Figure 4: the amplitude of the negative radial velocity signal in the southern quadrant reaches $\approx -45 \text{ km s}^{-1}$ at $\approx 130 \text{ kpc}$, most consistent with the simulations containing an LMC halo $\gtrsim 15\%$ the mass of the MW’s.

## 3.2. Tangential Velocities

Typical tangential velocity uncertainties are more than a factor of 10 larger than radial velocity uncertainties for the giants in our data set. Furthermore, any random or systematic uncertainties in the distance estimates propagate to the calculated tangential velocities. Therefore, it is challenging to straightforwardly interpret tangential velocity maps as we have done for radial velocities in the preceding section. However, the mean tangential velocity can be a powerful tracer of the LMC’s influence, particularly in the Galactic latitude ($v_{\text{t,b}}^*$) direction (e.g., D. Erkal et al. 2021; Y. Sheng et al. 2024). This component of the tangential velocity should capture the apparent “upward” (Galactocentric $z > 0$) reflex motion of halo stars as the inner Galaxy lurches “downward.”

Figure 5 shows the all-sky mean tangential velocity along the Galactic latitude $\langle v_{\text{t,b}}^* \rangle$ as a function of Galactocentric radius $r_{\text{gal}}$. These tangential velocities have been corrected for the solar reflex motion. There is a striking trend of increasing $\langle v_{\text{t,b}}^* \rangle$ as a function of $r_{\text{gal}}$, saturating at a maximum of $\langle v_{\text{t,b}}^* \rangle \approx 50 \text{ km s}^{-1}$ beyond 100 kpc. As in Figure 4, the dashed lines show simulation-based predictions from GC19. The increasing $\langle v_{\text{t,b}}^* \rangle$ over distance is a tell-tale signature of the LMC’s reflex motion in proper motions, and this is the first measurement of the growing tangential velocity to date (E. Vasiliev 2023). Y. Sheng et al. (2024) recently measured
the mean tangential velocity of RR Lyrae stars from 50–100 kpc, and that measurement is shown for comparison in Figure 4. The observed velocities shown in Figure 5 are consistent with the most massive LMCs in the GC19 simulations, with halo masses >$10^{11}$ M$_{\odot}$.

## 4. Modeling the Global Velocity of the Outer Halo

The preceding section demonstrated spatial variations in the mean velocities of halo stars—both on-sky and as a function of distance—that can be attributed to the LMC-induced reflex motion. These trends already point to a relatively massive $\gtrsim 10^{11}$ M$_{\odot}$ LMC (see Figures 4 and 5). In this section, we develop a physically motivated model of the inner Galaxy’s travel velocity $v_{\text{travel}}$ toward the outer halo. As Figures 4 and 5 demonstrate, the apparent reflex motion of halo stars should sharply grow as a function of distance. Therefore, we develop a flexible continuity model that allows the travel velocity $v_{\text{travel}}$ and all associated parameters to linearly vary across the distance range covered by our sample.

### 4.1. Global Velocity Model with Reflex Motion

Our goal is to fit the full 6D kinematics of observed stars to measure the amplitude and direction of the apparent velocity $v_{\text{travel}}$ of the inner galaxy relative to its distant outskirts. PP21 developed an elegant formalism that models $v_{\text{travel}}$ as a velocity vector in a rotated Galactocentric reference frame, and fit this model with a Gaussian likelihood for the observed parameters. We adapt the basic ingredients of the PP21 model formalism here, developing a likelihood for the observed on-sky position, distance, radial velocity, and proper motion of each star in terms of the model parameters.

Briefly, at a given $r_{\text{gal}}$, the stellar halo is assumed to have a mean 3D velocity vector in Galactocentric coordinates

$$ \langle \mathbf{v} \rangle = \mathbf{v}_{\text{travel}} + \langle v_r \rangle + \langle v_{\phi} \rangle + \langle v_{\theta} \rangle. \quad (1) $$

The conversion from heliocentric to Galactocentric coordinates is described at the start of Section 2. The reflex motion is encoded in the $v_{\text{travel}}$ parameter, which is expressed by a scalar amplitude $v_{\text{travel}}$ and two angular coordinates encoding the apex (maximum velocity) direction on the sky. Note that the travel velocity $v_{\text{travel}}$ represents the velocity of the inner MW toward the distant halo. The “reflex” velocity is the negative of the $v_{\text{travel}}$ vector, i.e., the apparent motion of distant halo stars toward observers on Earth.

The direction of $v_{\text{travel}}$ can be encoded by any on-sky angular quantity. Although PP21 used the Galactic coordinates $\ell$ and $b$ directly, we find that this leads to inefficient convergence of our fitting routine near the pole at $b \approx -90^{\circ}$, resulting in posterior distributions “pinched” at the southern pole. We therefore use the rotated Magellanic Stream coordinate frame from D. L. Nidever et al. (2008). The $B_{\text{MS}} = 0^{\circ}$ latitude of this frame closely follows the LMC orbit shown in Figure 3, and the present-day LMC is located at $L_{\text{MS}} = 0^{\circ}$. We have verified that this choice of frame does not influence the best-fit model parameters in any way, but enables more efficient sampling (away from any poles) and well-behaved posterior distributions in our probabilistic model-fitting routine.

Apart from the amplitude and direction of the reflex motion, the model in Equation (1) includes mean velocities in each of the spherical Galactocentric components—$\langle v_r \rangle$, $\langle v_{\phi} \rangle$, $\langle v_{\theta} \rangle$. This

allows for any departures in the bulk motion of the halo from the simple $v_{\text{travel}}$ model, for example, the halo’s overall velocity prior to the arrival of the LMC.

Finally, we include intrinsic dispersion parameters for each Galactic velocity component—$\sigma_{v,r}$, $\sigma_{v,l}$, $\sigma_{v,b}$—that are added in quadrature to the observed uncertainties in the model likelihood. Apart from measuring the intrinsic dispersion of the halo stars, these parameters can also absorb any under-estimation in the measurement uncertainties.

We depart from the PP21 model in that we explicitly expand all model parameters in terms of Galactocentric distance, making a tomographic measurement of the reflex motion. We adopt a simple linear expansion of each model parameter, such that the number of parameters is increased by a factor of 2. Each parameter can now have different values at $r_{\text{gal}} = 40$ and at $r_{\text{gal}} = 120$, with a linear trend between those two values. This results in two independent parameter vectors $S_{40}$ and $S_{120}$, each containing nine parameters: $v_{\text{reflex}}$, $\ell_{\text{MS,apex}}$, $b_{\text{MS,apex}}$, $\langle v_r \rangle$, $\langle v_{\phi} \rangle$, $\langle v_{\theta} \rangle$, $\sigma_{v,r}$, $\sigma_{v,l}$, and $\sigma_{v,b}$. We only fit stars out to $r_{\text{gal}} = 120$ with this model, reducing the sample size to 821 stars—373 from H3, 261 from SEGUE, and 187 from MagE. Restricting the fit to $r_{\text{gal}} < 120$ ensures relatively even coverage across the sky at the upper distance range, since almost all stars beyond this distance are from our southern MagE survey. We have tested three other choices of the fitted distance range (40–100 kpc, 30–120 kpc, 30–100 kpc), and find that they do not alter the main conclusions of this work. In the 30–120 kpc fit, the overall amplitude of the fitted $v_{\text{travel}}$ is lower at all distances. We speculate that this is due to the linearity assumption breaking down at $r_{\text{gal}} < 40$ kpc, since the $v_{\text{travel}}$ signature is expected to flatten at these distances (GC19, YPP24). The 40–120 kpc sample has the most homogeneous sky coverage across all distances, and is adopted as our fiducial sample for the rest of this work.

To prevent outliers from biasing our model fit, we clip velocity outliers from our sample. Specifically, we split the sample into three distance bins bounded by the following edges: [40, 60, 90, 120] kpc. In each distance bin, we compute the sigma-clipped mean and standard deviation of each velocity component: one radial and two tangential. We remove stars that are $\ge 3.5\sigma$ away from the mean in either velocity component. This results in 34 stars being removed from the sample, the vast majority of which are within $r_{\text{gal}} < 60$ kpc. These stars could belong to undetected substructures (or known structures like Sagittarius, missed by our angular momentum selection for Sgr members), or could simply have spurious stellar parameters or velocity measurements. Our fitted sample of stars from 40–120 kpc therefore has 772 stars.

For each star, the velocity model is calculated using the model parameters at that distance—linearly evaluated using a given parameter set {$S_{40}$, $S_{120}$}—and the likelihood is computed. This continuity model avoids the pitfalls of binning the data in distance. The substantial distance uncertainties are bound to place stars in incorrect bins, smearing out the signal by placing (intrinsically more numerous) nearby stars at larger distances, damping the inferred reflex velocity. The continuity model utilizes statistical information from the entire sample at once to uncover the distance-dependence of the model parameters, under the assumption that the true parameters vary linearly and continuously over Galactocentric distance.

Conceptually, the continuity model reduces the variance at any given distance, while potentially increasing the bias—for
<table><caption>**Table 1**
Summary of the Best-fit Parameters from Our $v_{\text{travel}}$ Continuity Model</caption><thead><tr><th>Parameter</th><th>r<sub>gal</sub> = 40 kpc</th><th>r<sub>gal</sub> = 120 kpc</th></tr></thead><tbody><tr><td>v<sub>travel</sub> [km s<sup>-1</sup>]</td><td>6<sup>+5</sup><sub>-4</sub></td><td>47<sup>+12</sup><sub>-12</sub></td></tr><tr><td>l<sub>MS,apex</sub> [deg]</td><td>-77<sup>+21</sup><sub>-32</sub></td><td>-10<sup>+17</sup><sub>-15</sub></td></tr><tr><td>b<sub>MS,apex</sub> [deg]</td><td>12<sup>+20</sup><sub>-19</sub></td><td>3<sup>+14</sup><sub>-14</sub></td></tr><tr><td>⟨v<sub>r</sub>⟩ [km s<sup>-1</sup>]</td><td>-8<sup>+5</sup><sub>-5</sub></td><td>-25<sup>+9</sup><sub>-10</sub></td></tr><tr><td>⟨v<sub>φ</sub>⟩ [km s<sup>-1</sup>]</td><td>14<sup>+5</sup><sub>-5</sub></td><td>-14<sup>+8</sup><sub>-8</sub></td></tr><tr><td>⟨v<sub>θ</sub>⟩ [km s<sup>-1</sup>]</td><td>-1<sup>+5</sup><sub>-5</sub></td><td>-3<sup>+12</sup><sub>-12</sub></td></tr><tr><td>σ<sub>v,r</sub> [km s<sup>-1</sup>]</td><td>106<sup>+4</sup><sub>-4</sub></td><td>69<sup>+6</sup><sub>-5</sub></td></tr><tr><td>σ<sub>v,l</sub> [km s<sup>-1</sup>]</td><td>78<sup>+4</sup><sub>-4</sub></td><td>20<sup>+5</sup><sub>-4</sub></td></tr><tr><td>σ<sub>v,b</sub> [km s<sup>-1</sup>]</td><td>56<sup>+3</sup><sub>-3</sub></td><td>33<sup>+6</sup><sub>-5</sub></td></tr></tbody></table>

**Notes.** Each parameter has a free value at $r_{\text{gal}} = 40$ kpc and at $r_{\text{gal}} = 120$ kpc, and is assumed to vary linearly between those values as a function of distance. The marginalized posterior distribution of each parameter is summarized using the 16th, 50th, and 84th quantiles, corresponding to $1\sigma$ uncertainties for a Gaussian distribution.

example, the parameters at nearby distances are, by definition, influenced by data at larger distances. By applying our model to the simulations from GC19, we have verified that the linear continuity model is sufficient to capture the overall trend of the $v_{\text{travel}}$ parameters, as well as the velocity dispersion parameters. In Appendix C we fit a simulation that matches the size and spatial sampling of our data set, affirming that our methodology is robust against uneven spatial sampling.

To compute the likelihood, the mean velocity model (Equation (1)) is projected into the “observed” heliocentric coordinates: the radial velocity and proper motions. The conversion from velocity to proper motion uses the heliocentric distance to each star, and is governed by the proportionality constant $k \equiv 4.74 \text{ km s}^{-1} \text{ kpc}^{-1}$ (mas yr$^{-1}$)$^{-1}$. The mean velocity model defines the Gaussian likelihood for our model fit:

$$ \mathcal{L}(\ell, b, d_{\text{helio}}, v_r, \mu_\ell, \mu_b | S_{40}, S_{120}). \quad (2) $$

The likelihood utilizes the full proper-motion covariance matrix provided by Gaia, and propagates our measured distance uncertainties into the tangential velocity uncertainties (see Equations (5)–(10) in PP21). The uncertainties of the on-sky positions are assumed to be negligible. For more details on the model formalism and likelihood, we refer to PP21 and YPP24.

The $v_{\text{travel}}$ continuity model has 18 free parameters: nine describing the model parameters at $r_{\text{gal}} = 40$ kpc, and nine at $r_{\text{gal}} = 120$ kpc. We adopt uninformative priors on all parameters. For the $B_{\text{MS,apex}}$ latitude, a uniform prior on $\cos[B_{\text{MS,apex}}]$ is used, and for the intrinsic dispersion parameters, a uniform prior on the logarithm of the dispersion is used. An initial maximum-likelihood fit is run using `scipy optimize` (P. Virtanen et al. 2020). The posterior distribution of model parameters is sampled with the nested sampling routine `dynesty` (J. S. Speagle 2020). The full corner plot of posterior distributions is shown in Appendix F, and the best-fit parameters are summarized in Table 1. The best-fit parameters from the linear continuity model are also evaluated at a reference distance of $r_{\text{gal}} = 100$ kpc, and listed in Table 2.

To ensure that our model implementation is consistent with previous measurements from PP21, we fit their data set using

<table><caption>**Table 2**
Summary of the Best-fit Parameters from Our $v_{\text{travel}}$ Continuity Model (See Table 1), Linearly Evaluated at a Reference Distance of 100 kpc</caption><thead><tr><th>Parameter</th><th>r<sub>gal</sub> = 100 kpc</th></tr></thead><tbody><tr><td>v<sub>travel</sub> [km s<sup>-1</sup>]</td><td>37<sup>+9</sup><sub>-9</sub></td></tr><tr><td>l<sub>MS,apex</sub> [deg]</td><td>-27<sup>+11</sup><sub>-10</sub></td></tr><tr><td>b<sub>MS,apex</sub> [deg]</td><td>6<sup>+9</sup><sub>-9</sub></td></tr><tr><td>⟨v<sub>r</sub>⟩ [km s<sup>-1</sup>]</td><td>-21<sup>+6</sup><sub>-7</sub></td></tr><tr><td>⟨v<sub>φ</sub>⟩ [km s<sup>-1</sup>]</td><td>-7<sup>+6</sup><sub>-6</sub></td></tr><tr><td>⟨v<sub>θ</sub>⟩ [km s<sup>-1</sup>]</td><td>-2<sup>+9</sup><sub>-9</sub></td></tr><tr><td>σ<sub>v,r</sub> [km s<sup>-1</sup>]</td><td>77<sup>+5</sup><sub>-4</sub></td></tr><tr><td>σ<sub>v,l</sub> [km s<sup>-1</sup>]</td><td>28<sup>+5</sup><sub>-4</sub></td></tr><tr><td>σ<sub>v,b</sub> [km s<sup>-1</sup>]</td><td>38<sup>+4</sup><sub>-4</sub></td></tr><tr><td>l<sub>apex</sub> [deg]</td><td>-80<sup>+18</sup><sub>-18</sub></td></tr><tr><td>b<sub>apex</sub> [deg]</td><td>-59<sup>+11</sup><sub>-9</sub></td></tr><tr><td>v<sub>travel,x</sub> [km s<sup>-1</sup>]</td><td>3<sup>+5</sup><sub>-5</sub></td></tr><tr><td>v<sub>travel,y</sub> [km s<sup>-1</sup>]</td><td>-17<sup>+6</sup><sub>-5</sub></td></tr><tr><td>v<sub>travel,z</sub> [km s<sup>-1</sup>]</td><td>-31<sup>+9</sup><sub>-9</sub></td></tr></tbody></table>

**Notes.** For convenience and ease of comparison, we also provide the apex direction transformed to Galactic coordinates, and the $v_{\text{travel}}$ velocity in Galactocentric coordinates (see Figure 9). The marginalized posterior distribution of each parameter is summarized using the 16th, 50th, and 84th quantiles, corresponding to $1\sigma$ uncertainties for a Gaussian distribution.

our model in Appendix A, reproducing their results. As an additional check to verify our results, we fit all stars in our data set with a distance-independent nine-parameter $v_{\text{travel}}$ model, directly analogous to that used by PP21 and YPP24. These results are presented in Appendix E, and are broadly consistent with the 18-parameter continuity model results presented in the rest of this section.

## 4.2. Reflex Response Results

Figure 6 summarizes the key results from the continuity model measurement of $v_{\text{travel}}$. The top panel illustrates the posterior distribution of the $v_{\text{travel}}$ direction in Galactic coordinates, evaluated at the reference distance of $r_{\text{gal}} = 100$ kpc. As a reminder, the continuity model evaluated at this reference distance utilizes information from stars spanning $r_{\text{gal}} = 40$ to $120$ kpc. Our measured travel velocity relative to the Galactic outskirts points directly at a past location along the LMC’s orbit, close to the present-day location of the LMC.

The bottom panel of Figure 6 shows the posterior distribution of the measured amplitude of $v_{\text{travel}}$ as a function of Galactocentric distance $r_{\text{gal}}$, with the $1\sigma$ credible interval shaded. The $v_{\text{travel}}$ amplitude significantly grows as a function of distance, from $\approx 6 \text{ km s}^{-1}$ at $r_{\text{gal}} = 40$ kpc to a maximum of $v_{\text{travel}} = 47 \pm 12$ at $r_{\text{gal}} = 120$ kpc. Along with the recent measurements of YPP24, this is the first measurement of the growing reflex velocity as a function of distance. Whereas YPP24 have a single bin containing all of their stars beyond 50 kpc, our larger data set and continuity model allow us to model the $v_{\text{travel}}$ variation out to 120 kpc.

The MW+LMC simulations from GC19 are also shown in the bottom panel of Figure 6. Since we only require the maximum-likelihood parameters (and not their uncertainties) from the simulations, we fit the simulations with `scipy optimize` (P. Virtanen et al. 2020) instead of `dynesty`. We
<figure><img src="image_8.png" /></figure>

**Figure 6.** The travel velocity of the Milky Way relative to the outer halo, as measured by our $v_{\text{travel}}$ continuity model as a function of Galactocentric distance $r_{\text{gal}}$. The MW+LMC simulations from GC19 are shown in blue, with line widths in the bottom panel corresponding to LMC halo masses from $[0.8, 1.0, 1.8, 2.5] \times 10^{11} M_{\odot}$. Top panel: the apex direction of the apparent reflex motion in Galactic coordinates, evaluated from our best-fit continuity model at a reference distance of $r_{\text{gal}} = 100$ kpc. The LMC's past orbit is overlaid in gray, and the background image shows the density of stars in Gaia DR3. Contours are placed at the 40th, 70th, and 90th quantiles, corresponding to Gaussian-equivalent levels of $0.5\sigma$, $1\sigma$, and $1.5\sigma$. The blue marker shows the apex direction inferred from the $1.8 \times 10^{11} M_{\odot}$ simulation of GC19. Bottom panel: amplitude of the travel velocity as a function of Galactocentric distance. The mean of posterior draws from our continuity model is shown, along with the $1\sigma$ credible interval. The reference distance used in the top panel is indicated with a tick mark along the top axis of the bottom panel.

have verified that measuring $v_{\text{travel}}$ directly from the simulation particles—by averaging particle velocities in spherical shells, as done by YPP24—produces a similar slope, and the linear model is sufficient to capture the overall growth of $v_{\text{travel}}$. However, fitting the simulation with the full model is important to disentangle the mean halo velocity terms from the $v_{\text{travel}}$ amplitude. The bottom panel of Figure 8 demonstrates that the amplitude of the travel velocity $v_{\text{travel}}$ is most consistent with simulations of a $\gtrsim 1.8 \times 10^{11} M_{\odot}$ LMC from GC19, or an LMC that is 15% the mass of the Milky Way.

Figure 7 shows the apex direction of $v_{\text{travel}}$ as a function of Galactocentric radius. The $1.8 \times 10^{11} M_{\odot}$ LMC simulation from GC19 is shown for comparison. In both the data and simulations, there is little evolution in the apex direction as a function of distance.

For the remainder of this section, we focus on the $v_{\text{travel}}$ continuity model evaluated at a reference distance of $r_{\text{gal}} = 100$ kpc. This distance provides a balance between being well-sampled by our data set while having a confidently detected reflex motion amplitude $\approx 40 \text{ km s}^{-1}$. Figure 8 contextualizes our measurement by comparing it to past measurements of the reflex velocity from PP21 and YPP24.

<figure><img src="image_9.png" /></figure>

**Figure 7.** The apex direction of $v_{\text{travel}}$ measured in this work (red lines) as a function of Galactocentric radius $r_{\text{gal}}$ with $1\sigma$ uncertainties shaded. The dashed blue line shows the same quantities from the $1.8 \times 10^{11} M_{\odot}$ LMC3 simulation of GC19.

Predictions from the self-consistent MW+LMC simulations of GC19 are also shown for the four LMC masses considered in that work.

The top panel of Figure 8 shows that the GC19 simulations predict a reflex velocity pointing toward a relatively recent position of the LMC's orbit, consistent with our measurement. On the other hand, the measurements of PP21 and YPP24 point toward the opposing Galactic quadrant. The difference between those measurements and the GC19 simulations is acknowledged in those works (see also E. Vasiliev 2023). The measurements in YPP24 get closer to the present-day location of the LMC at larger distances, and their appendices explore how different selection effects can nontrivially affect the inferred apex direction. Furthermore, the predicted $v_{\text{travel}}$ amplitude is quite small ($\lesssim 20 \text{ km s}^{-1}$) for the distances covered by the majority of the PP21 and YPP24 samples, making it more challenging to disentangle from other bulk motions in the halo than at larger distances.

Figure 9 offers an alternative perspective of our measurements by showing the Galactocentric components of the $v_{\text{travel}}$ velocity. As a reminder, $v_{\text{travel}}$ is the inferred velocity of the inner galaxy relative to the outer halo, which sources the apparent reflex velocity of those outer halo stars as measured by us. The corresponding $v_{\text{travel}}$ values for each MW+LMC simulation from N. Garavito-Camargo et al. (2019) are also evaluated at the reference distance of $r_{\text{gal}} = 100$ kpc, along with past measurements of $v_{\text{travel}}$ from E. Vasiliev et al. (2021), PP21, and YPP24. Note that E. Vasiliev et al. (2021) measured $v_{\text{travel}}$ by fitting the orbit of the Sagittarius stream, whereas PP21 and YPP24 directly fit the kinematics of halo stars. The E. Vasiliev et al. (2021) measurement in Figure 9 was provided by E. Vasiliev (2025, private communication).
<figure><img src="image_10.png" /><figcaption>**Figure 8.** Contextualizing the results of our $v_{\text{travel}}$ model with N-body simulations and past measurements. The MW+LMC simulations from GC19 are shown in blue, with increasing marker sizes and line widths corresponding to LMC halo masses from [0.8, 1.0, 1.8, 2.5] × 10¹¹ $M_{\odot}$. Previous measurements from PP21 and YPP24 are shown for comparison. For the YPP24 measurements, the marker size is proportional to distance bin they fitted, ranging from 25–60 kpc.</figcaption></figure>

Figure 9 demonstrates that although all measurements and simulations broadly agree regarding the $v_Z$ velocity (unsurprisingly, since the LMC remains in the southern hemisphere throughout its recent orbit), there is a difference in the $v_Y$ velocity. Specifically, we measure $v_Y \lesssim 0$ in this work (as do E. Vasiliev et al. 2021), whereas PP21 and YPP24 measured $v_Y > 0$. This difference is modest relative to the uncertainties in both the past and present measurements. In Appendix A, we investigate this discrepancy in detail, finding that the measurements in PP21 and YPP24 were likely influenced by the presence of Sagittarius Stream stars in their data sets (see also A. Byström et al. 2024). We conclude that the difference between our measurement and past measurements from PP21 and YPP24 is primarily due to the different selection cuts used to remove Sagittarius stars from the data sets.

Figure 10 provides a 3D illustration of the reflex motion measured in this work. A model of the MW disk is shown in the X-Z (left) and Y-Z (right) Galactocentric frame, with halo stars from our data set shown at their true 3D locations, and the 3D orbit of the LMC reproduced from GC19. A red arrow denotes our measured $v_{\text{travel}}$ vector evaluated at 100 kpc, arbitrarily scaled to reflect 500 Myr of motion. The stars in our data set are plotted at their measured positions, and colored by their corresponding Galactocentric radial velocity in the best-fit $v_{\text{travel}}$ model. This illustration conveys the main results of this work: the inner MW exhibits a significant net motion toward a the recent orbit of the LMC—relative to distant stars

<figure><img src="image_11.png" /><figcaption>**Figure 9.** The mean velocity $v_{\text{travel}}$ of the inner MW relative to the Galactic outskirts, expressed in Galactocentric coordinates. The posterior distribution of our $v_{\text{travel}}$ continuity model is evaluated here at a reference distance of 100 kpc. For comparison, we show the MW+LMC simulations of GC19 evaluated at the same distance, with increasing marker sizes corresponding to LMC halo masses from [0.8, 1.0, 1.8, 2.5] × 10¹¹ $M_{\odot}$. Previous measurements from E. Vasiliev et al. (2021, also evaluated at $r_{\text{gal}} = 100$ kpc), PP21, and YPP24 are also shown.</figcaption></figure>

in the outer halo—and the amplitude of this signal grows stronger as a function of distance.

### 4.3. The Global Velocities of the Outer Halo

In addition to the reflex velocity of outer halo stars, our $v_{\text{travel}}$ model also measures the bulk motion (and intrinsic dispersion) of the outer halo as a function of distance. Figure 11 shows the best-fit bulk halo velocity parameters from our model, the mean velocity components of the stellar halo (top panel), and their corresponding intrinsic dispersions (bottom panel). These measurements represent the sky-averaged motion of the MW halo, once the reflex response to the LMC has been accounted for.

The top panel of Figure 11 shows significant trends in the average radial and azimuthal velocity of the outer halo. The halo exhibits a net “compression” signal ($\langle v_r \rangle < 0$) in radial
<figure><img src="image_12.png" /><figcaption>**Figure 10.** Visualization of our best-fit reflex velocity model in 3D Galactocentric coordinates. A mock MW disk is shown, along with the measured 3D positions of stars in our data set. Stars are colored by their corresponding $v_{GSR}$ velocity evaluated from the best-fit $v_{travel}$ continuity model. The red arrow shows the measured $v_{travel}$ vector of the inner galaxy relative to the halo at 100 kpc, arbitrarily scaled up to 500 Myr of motion. The past orbit of the LMC from GC19 is shown in blue. An animated version of this figure is available, showing a full rotation about the Z-axis. The real-time duration of the animation is 15 s. A copy of this animation is available at https://bit.ly/mw-vtravel.</figcaption></figure>

(An animation of this figure is available in the online article.)

<figure><img src="image_13.png" /><figcaption>**Figure 11.** Spherical velocity components of the MW halo, fit to our all-sky data set simultaneously with the reflex velocity model. The top panel shows mean velocities in spherical Galactocentric coordinates, and the bottom panel shows intrinsic dispersions in Galactic velocities. Additionally, $1\sigma$ credible intervals are shaded.</figcaption></figure>

velocity at 40 kpc (e.g., YPP24), growing to ≈−20 km s⁻¹ at ≈120 kpc. Simulations of the LMC infall predict a similar dynamical compression—at the level of −20 km s⁻¹ at these distances (YPP24)—suggesting that the LMC’s gravitational pull may be responsible for this signature.

Another interesting signature in the top panel of Figure 11 is the nonzero azimuthal rotation measured by the $\langle v_\phi \rangle$ component. Whereas PP21 and YPP24 measured a net negative $\langle v_\phi \rangle$ at all distances, our best-fit model has net positive $\langle v_\phi \rangle$ out to ≈80 kpc, after which it switches to being negative. Our relative trend matches YPP24, with $\langle v_\phi \rangle$ becoming more negative as a function of distance. In our Galactocentric convention, this corresponds to a net retrograde halo that switches to prograde in the outskirts, although the uncertainties are large (see A. J. Deason et al. 2017).

The bottom panel of Figure 11 shows the radial trends of velocity dispersion fitted alongside our $v_{travel}$ model. Note that the continuity model assumed a linear trend of the log-dispersion as a function of $r_{gal}$, producing the apparent curvature of the lines in the bottom panel of Figure 11. We emphasize that these are the intrinsic dispersions from our $v_{travel}$ model, and are therefore inherently de-convolved from the measurement uncertainties in radial velocities and proper motions. All components of the halo velocity dispersion fall as a function of distance, as has been seen in prior studies mapping the outer halo (e.g., S. A. Bird et al. 2021). It is apparent that the radial dispersion is significantly higher than the tangential dispersions across all Galactocentric radii. We have verified that directly calculating the velocity dispersion in radial bins—and subtracting the median measurement error in each bin, in quadrature—produces quantitatively similar profiles. We do note that it is challenging to measure the intrinsic tangential dispersions at large distances, since they are sensitive to the adopted proper-motion and distance measurement uncertainties.

The radial bias (or lack thereof) of the stellar velocity ellipsoid is compactly summarized in the velocity anisotropy parameter (J. Binney & S. Tremaine 2008):

$$ \beta \equiv 1 - \frac{\langle v_{\phi}^2 \rangle + \langle v_{\theta}^2 \rangle}{2 \langle v_r^2 \rangle} $$

This is a crucial parameter for dynamical mass estimates of the MW that use the Jeans formalism (e.g., J. Binney & S. Tremaine 2008; W. Wang et al. 2020). We compute $\beta$ from 40–120 kpc
<figure><img src="image_14.png" /><figcaption>**Figure 12.** Velocity anisotropy parameter $\beta$ as a function of Galactocentric distance, as derived from the velocities shown in Figure 11. Previous measurements are shown for comparison: E. C. Cunningham et al. (2019) from the HALO7D survey, S. A. Bird et al. (2021) using SDSS/LAMOST K giants, and J. J. Han et al. (2024) using the H3 Survey. Our sample constrains $\beta$ twice as far as previous measurements, suggesting a radially anisotropic halo out to $\gtrsim 100$ kpc.</figcaption></figure>

using the bulk halo velocities (and their intrinsic dispersions) shown in Figure 11. For example, $\langle v_r^2 \rangle = \langle v_r \rangle^2 + \sigma_{v,r}^2$. Since our intrinsic dispersions are computed in the observed rather than Galactocentric frame, we make the simplifying assumption that $\sigma_{v,r} \approx \sigma_{v,los}$, $\sigma_{v,\phi} \approx \sigma_{v,l}$, and $\sigma_{v,\theta} \approx \sigma_{v,b}$, which is suitable for the distant stars in our sample. We have verified by fitting the simulations of N. Garavito-Camargo et al. (2019) that we can recover their input $\beta(r_{\text{gal}})$ profile using our fitting methodology. Our measured Galactocentric profile of $\beta$ is shown in Figure 12, with some corresponding measurements from the literature displayed for comparison.

Simulations of hierarchical galaxy formation predict that at large distances, the stellar halo should have radially biased orbits with $\beta > 0$ (e.g., J. Diemand et al. 2005; M. G. Abadi et al. 2006; L. V. Sales et al. 2007; V. P. Debattista et al. 2008; N. C. Amorisco 2017; S. R. Loebman et al. 2018). We measure a significant radial bias in the motions of stars out to $r_{\text{gal}} \approx 120$ kpc, with $\beta$ ranging from $\approx 0.6$ to $0.8$. These are the most distant all-sky measurements of the halo's velocity ellipsoid to date, supporting the idea that bulk of the MW halo was hierarchically assembled even in its most distant outskirts. Importantly, along with the reflex motion described earlier, these $\beta$ values should be incorporated in future studies that measure the dynamical mass of the MW using the Jeans equations (see W. Wang et al. 2020, for a recent review).

## 4.4. Local Effects of the LMC

Up to this point, we have mainly focused on the global response of the MW to the LMC, which is dominated by the reflex motion of the MW disk toward the LMC. As explored in detail by the simulations of GC19 and GC21, the LMC is also predicted to impart significant local perturbations along its past

orbit. This “transient response” is a direct consequence of the LMC’s gravitational influence and dynamical friction wake. Although C. Conroy et al. (2021) presented evidence of a wake-like feature in the density of stars, the tell-tale signature of the transient response should be imprinted in the kinematic structure of distant halo stars (see GC19 for a detailed exploration). These kinematic features should appear as higher-order (in angular scale) structures in velocity maps of the outer halo, above the dominant lower-order structure imprinted by the overall travel velocity (see Figure 3).

To more closely examine the local velocity perturbations due to the LMC, we compare selected observed velocity trends to the simulations of GC19 in Figure 13. The observed data points in Figure 13 are directly measured from our stellar data set, i.e., they are independent of the $v_{\text{travel}}$ model described above. An alternative and valuable approach would be to examine the residuals of the data compared to the $v_{\text{travel}}$ model, but we elect to show directly observed quantities here to highlight the simplest observable signatures that future studies should search for. The top panel shows the mean Galactocentric radial velocity along the $L_{\text{MS}}$ coordinate, which is roughly aligned with the past orbit of the LMC (D. L. Nidever et al. 2008). The extent of this plot only contains stars in the southern Galactic hemisphere. The overall trend of decreasing $v_{\text{GSR}}$ versus $L_{\text{MS}}$ is a signature of the global reflex response (red line), but there are higher-order “wiggles” in the trend from N-body MW+LMC simulations (blue lines) that signal local effects from the LMC. It is clear that even along the relatively precisely measured radial velocity component, a much larger sample of stars will be required to disentangle the kinematic signature of the local wake from the global reflex motion.

The middle and bottom panels of Figure 13 are similar to the top panel, but for the tangential velocity components. The bottom panel uses $B_{\text{MS}}$ as the x-coordinate, since the predicted $v_l^*$ motions produce a “converging” S-shaped trend in this space. Even for luminous giants, the typical tangential velocity uncertainties are an order of magnitude larger than the radial velocity uncertainties. Consequently, our data are entirely unable to detect the local LMC signature in tangential velocities. Regardless, we include these panels to demonstrate the strongest tangential velocity trends that future works should search for.

# 5. Discussion

## 5.1. Context and Limitations

Before discussing the conclusions of our work, it is important to highlight the limitations of our data set and methodology. As Figure 2 illustrates, our MagE survey has expanded spectroscopic coverage of the outer halo to the southern hemisphere, a major improvement over past samples. However, the bottom panel of Figure 2 also shows that the MagE survey is much more efficient at identifying stars beyond $\gtrsim 60$ kpc than corresponding northern surveys (H3 and SEGUE). Therefore, there is a strong need to observe more $\gtrsim 60$ kpc stars in the northern hemisphere, to achieve homogeneous sky coverage at these distances. Regardless, we have demonstrated (see Figure 4 and Appendix D) that the deeper southern coverage should not negatively affect the main measurements in our paper.
<figure><img src="image_15.png" /><figcaption>**Figure 13.** Kinematic signatures of the LMC-induced dynamical friction wake along the Galactocentric radial velocity (top panel), Galactic latitude (middle panel), and Galactic longitude (bottom panel) components. The tangential velocities have been corrected for the solar reflex motion. Binned mean measurements from our data set are shown in black with $1\sigma$ uncertainties, and predictions from the N-body simulations of GC19 are shown in blue. For the latter, thicker lines correspond to more-massive LMC halo masses. Predictions from the pure-reflex $v_{\text{travel}}$ model (with no local LMC perturbations) are shown in red.</figcaption></figure>

Although the dipole nature of the outer halo's Galactocentric radial velocity structure is apparent from the raw measurements (Figure 3), eventually we fit a kinematic model (motivated by PP21) to the 6D phase-space positions of our sample. Since the MW travel velocity $v_{\text{travel}}$ is elucidated in Galactocentric coordinates, any uncertainties in the solar

position and velocity would propagate into our measurements. As noted by PP21, increasing the solar velocity around the MW results in a lower $v_{\text{travel}}$ measurement, and vice versa.

The reflex velocity measured in this work is most sensitive to the mass ratio between the MW and LMC, not the absolute mass of the LMC itself. Therefore, our key result is perhaps best expressed as evidence that the LMC is $\gtrsim$15% the mass of the MW, since the GC19 simulations referred to in this work assume an MW virial mass $M_{\text{vir,MW}} = 1.2 \times 10^{12} M_{\odot}$. Simulations with a proportionally more-/less-massive MW and LMC should broadly match our observed signature too. Furthermore, although our measurements suggest a massive LMC on first infall, E. Vasiliev (2024) argued that the LMC may have had another (much more distant) pericenter passage in the past. Our observations cannot rule out this particular second-infall scenario, since the older pericenter passage has a negligible influence on the present-day reflex motion of the MW halo (see E. Vasiliev 2023; Y. Sheng et al. 2024; E. Vasiliev 2024).

In our kinematic $v_{\text{travel}}$ model, we retain the parameterization of the bulk halo velocities—$\langle v_r \rangle$, $\langle v_\phi \rangle$, $\langle v_\theta \rangle$—from PP21 for a straightforward comparison to these previous measurements. However, we note that future work should ideally utilize a more physically motivated model for the bulk motions. The $\langle v_\phi \rangle$ and $\langle v_\theta \rangle$ motions are challenging to interpret, and it might be preferable to model a single mean rotation component of the halo, with the amplitude and axis of rotation left as free parameters.

The $v_{\text{travel}}$ measured here differs in the sign of $v_Y$ compared to the past measurements of PP21 and YPP24. In Appendix A, we investigate this discrepancy in detail. Briefly, we find that the PP21 and YPP24 measurements appear to be significantly influenced by the presence of Sagittarius Stream (Sgr) stars in their data sets (as first noticed by A. Byström, private communication). These Sgr stars have an on-sky velocity distribution that closely mimics a dipole pointing in the quadrant opposing the LMC. In Appendix A, we demonstrate that after removing this Sgr contamination, the PP21 data set does not have a detectable $v_{\text{travel}}$ signature (Figure 15). Therefore, the most plausible explanation for the difference between our measurements and those from PP21 and YPP24 is the presence of Sgr contamination in the earlier data sets.

Another possible reason our measurement could differ from previous ones is the larger distance and sky coverage of our data set. In Appendix B we test the effect of sky coverage by fitting limited subsamples of our data, finding that our SEGUE+H3 +MagE data set can recover $v_{\text{travel}}$ even if only stars with $\delta > -20^\circ$—observable from the northern hemisphere—are fitted, although the uncertainties are larger. However, if we only use SEGUE data, $v_{\text{travel}}$ becomes undetected. Therefore, we concur with YPP24 that $v_{\text{travel}}$ can in principle be measured using only data observed from the northern hemisphere, but conclude that the increased sample size (and larger distance coverage) of the H3 and MagE data sets are vital contributors to our measurement.

Our agreement with simulations is by no means perfect, with the amplitude of our measured $v_Y$ being larger than the simulations of GC19 and the Sgr measurements of E. Vasiliev et al. (2021). One caveat is that our inclusion of southern data could make our measurements more susceptible to local halo perturbations of the LMC. Although the simulations of GC19 suggest that the all-sky velocity structure is quite smooth and dominated by the reflex motion itself, any deviations from the simulation (or unresolved substructure from the LMC halo
itself) could bias our measurement. The LMC's local influence, including the presence of stripped stars from the Clouds, can be explicitly modeled in the future.

Finally, an implicit assumption in our measurements—and the comparison to simulations—is that there is a relatively smooth component of the MW halo, relative to which a reflex velocity can be measured. The $v_{\text{travel}}$ measurement uncertainties quoted in this work are purely statistical, assuming a particular kinematic model. We have carefully removed known substructures including dwarf galaxies and the Sagittarius Stream from our data set. However, additional unknown substructures may be present in the data. When fitting our $v_{\text{travel}}$ model, we iteratively clip stars that disagree with the model by >3$\sigma$. These cuts should prevent small-scale undetected substructures from biasing our overall measurement. Regardless, given the emerging picture of the various constituents of the outer halo, it would be desirable for future simulations of the MW+LMC interaction to incorporate more realistic and structured models of the MW halo. In addition to including the Sagittarius Stream (e.g., E. Vasiliev et al. 2021), simulations could utilize the latest observations of the Galactic halo's spatially varying density profile (e.g., J. J. Han et al. 2022; J. A. S. Amarante et al. 2024) and velocity ellipsoid (e.g., S. A. Bird et al. 2021; see also Figures 11 and 12 in this work).

## 5.2. The Dynamical Disequilibrium of the Outer Halo

The vast majority of dynamical analyses of our Galaxy—for example, most measurements of the total mass of the MW—assume that the stellar halo is in a state of equilibrium. However, it remains an open question to what extent the MW has any “smooth” stellar halo component at these distances. Recent surveys have revealed that the inner halo is dominated by substructure from accreted dwarf galaxies, some well-mixed in phase space, and some retaining coherent structure to this day (e.g., R. P. Naidu et al. 2020; K. Malhan et al. 2022). Such a paradigm likely extends to the outer halo, where the fraction of surviving debris from smaller dwarf galaxies is higher (e.g., N. C. Amorisco 2017).

There is a growing picture that the density distribution of stars beyond ≈50 kpc is highly anisotropic and structured (see Figure 11 in V. Chandra et al. 2023b; also see C. Conroy et al. 2021; J. A. S. Amarante et al. 2024). Apart from the dominant contribution of the Sagittarius stream, there exist prominent stellar overdensities including the Outer Virgo Overdensity, the Pisces Overdensity, and the elongated Pisces Plume. V. Chandra et al. (2023a) attributed the former two overdensities to the MW's most recent major merger, Gaia-Sausage Enceladus (A. Helmi et al. 2017; V. Belokurov et al. 2018). The Pisces Plume appears to have a Magellanic origin, perhaps comprising both the dynamical friction wake of the infalling LMC, and the stripped Magellanic Stellar Stream (V. Belokurov et al. 2019; C. Conroy et al. 2021; V. Chandra et al. 2023b; J. A. S. Amarante et al. 2024).

Even if there exists a relatively smooth and well-mixed component of the stellar halo beyond ≈50 kpc, it is a priori unlikely that the outer halo retains dynamical equilibrium, since the orbital timescales of these stars are longer (≈1 Gyr) than the timescales over which the gravitational potential of the MW is changing due to the arrival of the Clouds (e.g., D. Erkal et al. 2019b). However, the precise amplitude and nature of this disequilibrium has been challenging to measure due to the paucity of all-sky stellar kinematic samples at these large distances.

Prior to any detailed modeling, Section 3 and Figures 3–5 in this work clearly illustrate the dramatic departure from any notion of “dynamical equilibrium” beyond ≈50 kpc. In Section 4 we use a physical model to demonstrate that the kinematic disequilibrium of the Galactic outskirts can be explained—to first order—by the gravitational influence of the LMC. The kinematic influence of the LMC on our sample of stars is notably strong, suggesting an LMC dark matter halo on the massive side of recent estimates, ≳15% the mass of the Milky Way.

Neglecting the LMC-induced disequilibrium of the halo can inflate dynamical mass estimates of the MW itself by up to 50% (D. Erkal et al. 2019b; L. Correa Magnus & E. Vasiliev 2022; A. Kravtsov & S. Winney 2024). Even mass measurements from the so-called “timing argument” can be inflated by 10% if the travel velocity is ignored (K. Chamberlain et al. 2023). Although we have made strides here toward observationally measuring this disequilibrium, the physical interpretation depends not only on the mass of the MW–LMC system, but also the past orbit of the LMC. Although the SMC is often ignored in simulations due to its low mass (≈10% of the LMC), it can substantially alter the implied past orbit of the LMC (E. Vasiliev 2023). There is therefore a growing need for a comprehensive dynamical analysis that includes all four protagonists in the outer halo—the MW, LMC, SMC, and Sagittarius—a “tango for four” in the same spirit as E. Vasiliev et al. (2021). Our survey of the outer halo has contributed two strong observational constraints for this problem: the reflex motion measured in this work, and the Magellanic Stellar Stream (V. Chandra et al. 2023b).

## 5.3. The Elusive Dark Matter Wake of the LMC

The dynamical friction wake of the LMC is predicted to leave imprints in both the density and kinematics of halo stars along the past orbit of the LMC (GC19; V. Belokurov et al. 2019; C. Conroy et al. 2021; M. Cavieres et al. 2025). Although this wake has been detected in the density of stars—forming the so-called “Pisces Plume” overdensity—a confident detection of the wake’s kinematic signature remains elusive. In Figure 13 we search for the LMC’s local kinematic perturbation in radial and tangential velocities, over and above the bulk reflex motion. However, even in the relatively well-measured radial velocity direction, the detection is weak to nonexistent given current observational errors.

The middle and bottom panels of Figure 13 show predicted and observed tangential velocity signatures along the past orbit of the LMC. These signatures are more straightforward to disentangle from the reflex-only model (compared to radial velocities), as indicated by the separation between the blue and red curves in Figure 13. Despite the size and fidelity of our giant sample, the observational errors are currently far too large to detect these local variations in the tangential velocity—the typical uncertainty in each bin is $30 \text{ km s}^{-1}$, of the same order as the predicted peak-to-peak signal.

Based on our results, detecting the tangential velocity signature of the LMC’s dark matter wake at the $2\sigma$ level will require ≳3× smaller velocity uncertainties, or a ≈10× larger sample size (≈8000 stars beyond 40 kpc). Updated proper motions from the next Gaia data release should also shrink proper-motion uncertainties by almost a factor of 2. A mixture of the data set presented here and ongoing surveys like SDSS-V (J. A. Kollmeier et al. 2017) and DESI (A. P. Cooper et al. 2023) should eventually be sensitive to the velocity signature of the LMC’s dark matter wake.
## 6. Conclusions

We have assembled the largest all-sky data set of 851 luminous red giant stars from 40–160 kpc with full 6D kinematics, including ≈300 new stars in the southern hemisphere from our own tailored spectroscopic survey. In this work, we utilize this data set to measure the motion of the inner galaxy relative to the outer halo, and find that it is induced by a massive and recently infalling LMC. Our main conclusions are summarized below.

1. We construct an all-sky data set of 851 red giant stars beyond 40 kpc using spectroscopy from H3, SEGUE, and our own MagE outer halo survey. All spectra are homogeneously analyzed with the same stellar parameter pipeline, resulting in the largest all-sky catalog of distant halo stars with 6D kinematic information (Figures 1 and 2).

2. The radial velocity distribution of distant halo stars has significant spatial structure (Figure 3). We measure an all-sky dipole radial velocity signature with an apex aligned toward the present-day location of the LMC. The velocity dipole grows as a function of distance, a hallmark signature of the expected “reflex response” of the MW halo to a massive LMC (Figure 4).

3. We measure—for the first time—the growth of the tangential velocity signature of the reflex motion as a function of distance (Figure 5). The tangential velocity along the Galactic latitude $⟨v_{t,b}^*⟩$ rises to ≈50 km s⁻¹ beyond 100 kpc, consistent with the response to a massive $\gtrsim 10^{11} M_{\odot}$ LMC.

4. We fit the 6D positions and velocities of 821 stars from 40 $< r_{gal} < 120$ to infer the amplitude and direction of the inner MW’s motion relative to distant halo stars. The apex direction of this travel velocity points toward the past orbit of the LMC, close to its present-day location (Figure 6, top panel). The amplitude of the travel velocity grows from ~0 at $r_{gal} = 40$ up to 50 ± 20 km s⁻¹ at $r_{gal} = 120$ kpc (Figure 6, bottom panel). A 3D schematic of the measurement is shown in Figure 10 (see animated version here<fnref n="16" />).

5. The $v_{travel}$ direction and amplitude measured in this work align with predictions from N-body simulations of the MW halo’s response to the LMC. The measured amplitude of the reflex velocity is most consistent with an LMC that is ≈15% the mass of the Milky Way, based on the simulations of GC19 (Figure 6).

6. Along with $v_{travel}$ itself, the bulk velocity and velocity dispersion components of the halo are simultaneously measured (Figure 11). We detect a strong radial “compression” signal ($⟨v_r⟩ \lesssim 20$ km s⁻¹), as well as some net azimuthal rotation in the outer halo.

7. The intrinsic velocity dispersion of the MW halo declines as a function of distance along all three spherical coordinates. After accounting for the LMC-induced reflex motion, outer halo stars exhibit significantly radially biased motions, with a rising profile of the velocity anisotropy parameter $\beta \approx 0.6–0.8$ from to 40–120 kpc (Figure 12).

8. We search for the local influence of the LMC on the mean radial velocities of MW halo stars along its past orbit (Figure 13). A much larger data set will be required

to definitively measure these effects, at least twice as large to detect the radial velocity signature, and ≈10× larger for the tangential velocity signature. Future spectroscopic surveys—combined with Gaia DR4— should be sensitive to these signatures.

## Acknowledgments

We thank the referee, Michael Petersen, for thorough feedback that significantly improved this work. V.C. gratefully acknowledges a Peirce Fellowship from Harvard University. R.P.N. acknowledges support for this work provided by NASA through the NASA Hubble Fellowship grant HST-HF2-51515.001-A awarded by the Space Telescope Science Institute, which is operated by the Association of Universities for Research in Astronomy, Incorporated, under NASA contract NAS5-26555. C. L. acknowledges funding from the European Research Council (ERC) under the European Union’s Horizon 2020 research and innovation program (grant agreement No. 852839). It is a pleasure to thank Amanda Byström, Nitya Kallivayalil, David Nidever, Theo O’Neill, Ekta Patel, Michael Petersen, Yanjun Sheng, Eugene Vasiliev, and Rashid Yaaqib for insightful conversations and feedback. This work benefited from discussions at the 2024 “Milky Clouds over Manhattan” Workshop at the Flatiron Institute Center for Computational Astrophysics.

We are indebted to the CfA and MIT telescope time allocation committees for enabling our long-term survey of the Galactic outskirts. We thank the staff at Las Campanas Observatory—including Yuri Beletsky, Carla Fuentes, Jorge Araya, Hugo Rivera, Alberto Pastén, Mauricio Martínez, Roger Leiton, Matías Díaz, Carlos Contreras, and Gabriel Prieto—for their invaluable assistance. The H3 Survey is funded in part by NSF grant NSF AST-2107253. The computations in this paper were run on the FASRC Cannon cluster supported by the FAS Division of Science Research Computing Group at Harvard University.

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 DPAC has been provided by national institutions, in particular the institutions participating in the Gaia Multilateral Agreement. This research has made extensive use of NASA’s Astrophysics Data System Bibliographic Services.

**Facilities:** Magellan:Baade (MagE), Gaia, Sloan, PS1, CTIO:2MASS, FLWO:2MASS, WISE.

**Software:** numpy (C. R. Harris et al. 2020), scipy (P. Virtanen et al. 2020), matplotlib (J. D. Hunter 2007), astropy (Astropy Collaboration et al. 2013, 2018, 2022), gala (A. M. Price-Whelan 2017; A. Price-Whelan et al. 2020), MINESweeper (P. A. Cargile et al. 2020), Reflex-Motion (PP21).

## Appendix A
The PP21 Data Set

A key conclusion from our work is that the apex direction of the MW’s motion points toward a very recent location along the past orbit of the LMC (Figure 6). Although this apex direction agrees well with the results from A. Byström et al. (2024), it is somewhat different from the direction reported in PP21 and YPP24. In particular, the latter two works report

<footnote n="16">¹⁶ https://bit.ly/mw-vtravel</footnote>
<figure><img src="image_16.png" /><figcaption>**Figure 14.** Comparing different angular momentum selections to identify and remove debris from the Sagittarius Stream. The entire sample from PP21 is shown in black. Stars that they identify as Sgr members using a spherical angular momentum selection are shown in orange, and stars that our linear angular momentum cut would identify as Sgr debris are shown in red. For reference, the Sgr simulation from E. Vasiliev et al. (2021) is shown in green.</figcaption></figure>

an apex direction that points toward the opposing Galactic quadrant as the LMC (Figure 8). In this section, we examine the data set used by PP21 to investigate the source of this difference. We use the K giant data set from PP21—although that work also fitted BHB stars and satellite galaxies, the K giants were by far the most constraining due to their low proper-motion uncertainties. We find that the most plausible reason for our differing measurements is the presence of Sagittarius Stream stars in the PP21 and YPP24 data sets.

PP21 and YPP24 utilized the angular momentum-based cuts from J. Peñarrubia & M. S. Petersen (2021) to identify Sgr stars, selecting stars within $3000 \text{ kpc km s}^{-1}$ of the Sgr progenitor's angular momentum: $[L_x = 605, L_y = -4515, L_z = -1267] \text{ kpc km s}^{-1}$. Figure 14 shows the K giant data set from M. S. Petersen & J. Peñarrubia (2021), with stars satisfying their Sgr selection highlighted in orange. We only show stars within 20° of the Sgr orbital plane, as defined by the D. R. Law & S. R. Majewski (2010b) coordinate frame. In red we show stars from their data set that would fail our broader angular momentum selection drawn from B. D. Johnson et al. (2021), which removes all stars with $L_Y < -2.5-0.3 \times L_Z$, where the angular momentum is in units of $10^3 \text{ kpc km s}^{-1}$. In the background, green points show Sgr stars from the simulations of E. Vasiliev et al. (2021), which were carefully fit to observations of the stream.

Figure 14 shows that there is a marked overdensity of Sgr-like stars left in the PP21 data set after their Sgr cuts are applied. Note that the observed stars span a wider swath of angular momentum space than the E. Vasiliev et al. (2021) simulations due to observational uncertainties, chiefly the distance uncertainty. To further illustrate this point, we show the clean PP21 data set in Figure 15 after their Sgr cut has been applied. The top panel shows radial velocities as a function of Galactic longitude, and the middle panel shows an all-sky map

<figure><img src="image_17.png" /><figcaption>**Figure 15.** Illustrating the presence and impact of Sgr debris in the sample from PP21. Top panel: Galactocentric radial velocity as a function of Galactic longitude. Middle panel: on-sky distribution in Galactic coordinates. We show the clean sample from PP21 as black circles, highlighting with red crosses stars that our cuts would identify as Sgr debris. For reference, the Sgr simulation from E. Vasiliev et al. (2021) is shown (green). Bottom panel: stars from PP21 that we would classify as Sgr debris, colored by Galactocentric radial velocity. The E. Vasiliev et al. (2021) simulation of Sgr stars is also shown with the same color scale.</figcaption></figure>

in Galactic coordinates. In both panels, it is clear that a significant fraction of the K giant data set fit by PP21 is coincident with Sgr debris. In particular, of the 548 stars in the clean PP21 data set, our angular momentum cuts would have
<figure><img src="image_18.png" /><figcaption>**Figure 16.** Posterior distribution of $v_{\text{travel}}$ in Galactocentric coordinates inferred by fitting the PP21 data set, before (top panel) and after (bottom panel) applying our angular momentum cuts to remove the ≈30% residual Sgr contamination in that data set. The contours show our measurements, and the diamonds show the reported measurements in PP21. Blue contours/markers correspond to $v_{\text{travel},x}$ on the x-axis and $v_{\text{travel},y}$ on the y-axis, and orange ones correspond the same in the x-z plane.</figcaption></figure>

removed a further 165 stars as Sgr contaminants, 30% of the data set fit by PP21.

The bottom panel of Figure 15 illustrates how the 30% Sgr contamination could bias the reflex motion apex direction inferred by PP21 (and by extension, YPP24). We show the 165 stars in their final sample that fall into the Sgr selection used in this work, colored by Galactocentric radial velocity. In the background are Sgr particles from the simulations of E. Vasieliev et al. (2021). The Sgr debris have a strong signature in radial velocities, which broadly appears like a dipole on the sky: moderately positive in the northern hemisphere, and strongly negative in the southern hemisphere. Furthermore, this dipole is oriented along the Sgr stream, with an apex direction (region of maximum negative Galactocentric radial velocity) in the Galactic quadrant opposing the LMC. This is

the same quadrant as the apex direction inferred by PP21 (indicated in Figure 15; see also Figure 8). Clearly, the presence of these Sgr stars could significantly influence the apex direction inferred by PP21.

To further demonstrate the influence of Sgr contamination, we fit the stars in the PP21 data set with our $v_{\text{travel}}$ model, and also fit their data set after applying our more broad angular momentum selection (Figure 14). Following PP21, we use the distance-independent $v_{\text{travel}}$ model here, fitting all stars in the data set with a single set of $v_{\text{travel}}$ parameters (see Appendix E for the analogous measurement on our own data set). The resulting posterior distributions of $v_{\text{travel}}$ are shown in Figure 16. The top panel shows fits to the PP21 data set after applying their Sgr cuts, and the bottom panel shows the same after our Sgr cuts are applied to their data set. Once our Sgr cuts are applied, the apparent $v_{\text{travel}}$ signature in the PP21 data set goes away, with $v_{\text{travel}}$ being within $1\sigma$ of zero in all components. The apparent $v_{\text{travel}}$ measured by PP21 (and by extension YPP24) therefore appears to be primarily driven by the kinematics of Sgr contaminants in those data sets, although some residual $v_{\text{travel}}$ signal may remain.

The top panel of Figure 16 also shows that any differences with PP21 are not due to a mis-specification of our $v_{\text{travel}}$ model, since we recover their quoted parameters within $1\sigma$—the small disagreement in $v_{\text{travel},x}$ is likely due to subtle differences in the sample selection or sigma-clipping methodology. We conclude that our modeling framework is consistent with that of PP21 and YPP24, and the differences in our results are primarily due to Sgr contamination in the earlier data sets (see also A. Byström et al. 2024, who arrive at the same conclusion).

## Appendix B: The Importance of All-sky Coverage

In Appendix A we argue that the primary reason our results differ from past measurements made by PP21 and YPP24 is our different strategies to remove contamination from the Sagittarius Stream, with apparent Sgr contamination visible in the earlier data sets (see Figure 15). However, our data set differs from past data in two other ways: the inclusion of the H3 Survey adds a significant number of outer halo stars in the northern hemisphere, and the inclusion of our MagE survey fills in a missing southern Galactic quadrant (see Figure 2). In this section, we test the effect of sky coverage on our $v_{\text{travel}}$ measurement.

We repeat the $v_{\text{travel}}$ measurement described in Section 4 (incorporating the linear continuity model from 40 to 120 kpc) on two subsets of our data. The first subset emulates the selection function of surveys observed from the northern hemisphere (e.g., H3 and SEGUE): $\delta > -20^\circ$. The resulting subset is predominantly composed of H3 and SEGUE stars, although some MagE stars at equatorial latitudes are also present. The second subset emulates the SEGUE-only data set used by PP21, albeit with our updated *MINESweeper* stellar parameters and broad cuts to remove Sagittarius contamination.

The resulting measurements of $v_{\text{travel}}$ in Galactocentric coordinates are shown in Figure 17, evaluated at $r_{\text{gal}} = 100$ kpc analogous to Figure 9. Using north-only data—including our SEGUE, H3, and MagE data—produces results somewhat consistent with our full data set, albeit with larger uncertainties. Furthermore, the $v_{\text{travel},x}$ component is
<figure><img src="image_19.png" /><figcaption>**Figure 17.** Posterior distribution of $v_{\text{travel}}$ in Galactocentric coordinates, inferred by fitting limited subsamples of our data set. We show the measurement from using our full data set (solid contour, identical to Figure 9), from using only stars observable from the northern hemisphere (dashed), and from using only data from the SEGUE survey (dotted).</figcaption></figure>

mildly biased relative to the full fit, and the $v_{\text{travel,y}}$ signature is slightly weaker.

On the other hand, using SEGUE-only data—after appropriately removing Sagittarius contamination using our angular momentum cuts (see Figure 14)—results in a nondetection of the $v_{\text{travel}}$ signature. Note that this is analogous to the measurement shown in Figure 16 in that it only uses SEGUE data, but with some differences. We use our own MINE-Sweeper stellar parameters rather than the X. Xue et al. (2015) parameters used by PP21. Furthermore, in Figure 16 we fit a single $v_{\text{travel}}$ at all distances for consistency with PP21, whereas in Figure 17 we use the full linear continuity model. There are fewer stars in our SEGUE data set than in PP21 since we apply more stringent SNR cuts (SNR > 10 per pixel) to ensure reliable stellar parameters. We have verified that there are no major systematic differences between our radial velocities and distances compared to the parameters from

X. X. Xue et al. (2015). Despite these differences, the conclusions from both of these tests are similar: after Sagittarius contamination is removed, there does not appear to be a strong $v_{\text{travel}}$ signature in the SEGUE-only data set.

In summary, the tests shown in Figure 17 demonstrate that although it is feasible to detect $v_{\text{travel}}$—albeit with mildly biased orientations—using only data observed from the northern hemisphere (as argued by YPP24, and recently shown by A. Byström et al. 2024 using DESI data), our SEGUE data set alone appears insufficient to reliably measure $v_{\text{travel}}$. The inclusion of H3 and MagE data substantially improves the north-only measurement, and the best measurement (unsurprisingly) comes from our all-sky data set that incorporates all three surveys. These tests also demonstrate that the reflex response is a sensitive and challenging measurement that can be affected by a variety of factors including sky coverage, distance coverage, and contamination by substructure.

## Appendix C: Fitting Data-matched MW+LMC Simulation

When fitting any model to data, it is important to consider the influence of the data’s selection function. Our data set has the advantage of sampling the entire sky, an improvement over past measurements. Furthermore, our MagE data set systematically spans larger distances than past halo data sets. However, as the top panel in Figure 2 illustrates, this produces an on-sky asymmetry in the typical distance toward stars in our data set. Note that for our $v_{\text{travel}}$ measurement, we restrict the data to 40–120 kpc, which somewhat mitigates this asymmetry. However, it is important to verify that this selection effect does not affect the conclusions of this work.

Throughout this work, we have shown for comparison $v_{\text{travel}}$ measurements made by fitting the simulations of GC19 (e.g., see Figure 6). For this Appendix, we subsample the simulation particles to match the spatial distribution and size of our final data set. Specifically, for each of the 821 stars in our final observed sample, we select the nearest particle (in 3D Galactocentric coordinates) from the GC19 LMC3 simulation, which has a $1.8 \times 10^{11} M_{\odot}$ LMC. This produces a mock data set that exactly matches the spatial selection function of our observed data set. Since distances are expected to be our dominant source of measurement uncertainty, we assume 10% distance uncertainties on this mock data set.

We fit the data-matched simulation using the $v_{\text{travel}}$ model described in Section 4, and show the resulting posterior contours in Figure 18. The $v_{\text{travel}}$ inferred from fitting the simulation without any imposed selection function is also shown. We also verify that the distance-dependence of $v_{\text{travel}}$ is well captured, following the intrinsic behavior of the simulation well within $1\sigma$ uncertainties. Finally, we verify that our result is not affected by our broad angular momentum cut to remove Sgr contamination (based on B. D. Johnson et al. 2020). Although B. D. Johnson et al. (2020) applied this cut to stars across the sky, we only apply it to stars within 25° of the Sgr orbital plane, and remove ≈15% of field stars from the simulated data set. However, removing these stars does not meaningfully alter the best-fit $v_{\text{travel}}$ parameters. These tests support our claim that the size and spatial selection function of our data set—and our strategy to remove Sgr contamination—do not limit or bias its ability to measure $v_{\text{travel}}$.
<figure><img src="image_20.png" /><figcaption>**Figure 18.** Posterior distribution of $v_{\text{travel}}$ in Galactocentric coordinates from our fit to the data-matched GC19 LMC3 simulation. $v_{\text{travel}}$ is evaluated at the reference distance of $r_{\text{gal}} = 100$ kpc, as in Figure 9. The “true” $v_{\text{travel}}$ inferred without any selection function applied is shown with a pentagon. After the imposition our selection function and sample size, $v_{\text{travel}}$ is recovered to within $1.5\sigma$ (the outermost contour).</figcaption></figure>

For ease of comparison to future works, we include in Table 3 our best-fit parameters to the LMC3 simulation from N. Garavito-Camargo et al. (2019). To determine these

<table><caption>**Table 3**
Summary of the Best-fit Parameters from Applying Our $v_{\text{travel}}$ Continuity Model to a Random Subsample of 10,000 Stars from the LMC3 Simulation of GC19</caption><thead><tr><th>Parameter</th><th>r<sub>gal</sub> = 40 kpc</th><th>r<sub>gal</sub> = 120 kpc</th></tr></thead><tbody><tr><td><i>v</i><sub>travel</sub> [km s<sup>-1</sup>]</td><td>11</td><td>51</td></tr><tr><td><i>ℓ</i><sub>MS,apex</sub> [deg]</td><td>-25</td><td>-48</td></tr><tr><td><i>b</i><sub>MS,apex</sub> [deg]</td><td>8</td><td>4</td></tr><tr><td>⟨<i>v</i><sub>r</sub>⟩ [km s<sup>-1</sup>]</td><td>-5</td><td>-20</td></tr><tr><td>⟨<i>v</i><sub>φ</sub>⟩ [km s<sup>-1</sup>]</td><td>-1</td><td>-3</td></tr><tr><td>⟨<i>v</i><sub>θ</sub>⟩ [km s<sup>-1</sup>]</td><td>8</td><td>9</td></tr><tr><td>σ<sub>v,r</sub> [km s<sup>-1</sup>]</td><td>115</td><td>87</td></tr><tr><td>σ<sub>v,f</sub> [km s<sup>-1</sup>]</td><td>90</td><td>62</td></tr><tr><td>σ<sub>v,b</sub> [km s<sup>-1</sup>]</td><td>85</td><td>66</td></tr></tbody></table>

parameters, we randomly subsampled 10,000 particles from the simulation, with no selection function or uncertainties applied. These parameters (and their analogs for the other three simulations from GC19) are used as the reference comparisons throughout this work.

## Appendix D: All-sky Velocity Maps versus Distance

As Figure 2 illustrates, our MagE survey in the southern hemisphere systematically pushes to farther distances than the northern H3 and SEGUE samples used here. In the context of the all-sky velocity maps of Figure 3, one possible concern is that the deeper southern coverage could bias these maps and produce artificial structures. To investigate this further, in Figure 19 we reproduce the mean velocity map of Figure 3 in four disjoint $r_{\text{gal}}$ bins. The stars in each bin are shown as black circles, demonstrating the underlying spatial distribution used to compute each velocity map.

Although the maps shown in Figure 19 are naturally noisier than Figure 3—which uses our entire sample—the features are broadly consistent. Figure 19 therefore provides a reassuring check that the overall velocity structure seen in Figure 3 is not dominated by the inhomogeneous distance ranges spanned by the surveys combined here, but rather represents the underlying velocity structure of the halo. It is naturally desirable to follow up the results presented here with more spatially homogeneous data sets in the future. As Figures 2 and 19 illustrate, there is now a strong need for more distant halo tracers in the northern hemisphere to match the advances made by our MagE survey in the south.
<figure><img src="image_21.png" /><figcaption>**Figure 19.** All-sky maps of the mean Galactocentric radial velocity $v_{\text{GSR}}$ using our data set, smoothed on 45° scales. These figures are computed identically to Figure 3, but are split into four disjoint $r_{\text{gal}}$ bins. On-sky positions of stars in each distance bin are shown as black circles.</figcaption></figure>
## Appendix E

### Distance-independent Reflex Model

In Section 4, we built on the work of PP21 to develop a distance-dependent continuity model to measure the increasing reflex velocity as a function of $r_{gal}$. Although this model is

<figure><img src="image_22.png" /><figcaption>**Figure 20.** Results from the distance-independent $v_{travel}$ fit to all data from 40–120 kpc. The other elements are described in the caption for Figure 6. In the bottom panel, the horizontal error bar denotes the 16th and 84th quantiles of our data set’s $r_{gal}$ distribution. This distance-independent fit provides consistent results with the full distance-dependent model described in Section 4, most importantly regarding the apex direction of the reflex motion.</figcaption></figure>

<table><caption>**Table 4**

Summary of the Best-fit Parameters from Our $v_{travel}$ Model That Uses All Stars from 40–120 kpc</caption><thead><tr><td>Parameter</td><td>Value</td></tr></thead><tbody><tr><td><i>v</i><sub>travel</sub> [km s<sup>-1</sup>]</td><td>17<sup>+4</sup><sub>-4</sub></td></tr><tr><td><i>ℓ</i><sub>apex</sub> [deg]</td><td>-86<sup>+35</sup><sub>-49</sub></td></tr><tr><td><i>b</i><sub>apex</sub> [deg]</td><td>-73<sup>+7</sup><sub>-10</sub></td></tr><tr><td>⟨<i>v</i><sub>r</sub>⟩ [km s<sup>-1</sup>]</td><td>-16<sup>+3</sup><sub>-3</sub></td></tr><tr><td>⟨<i>v</i><sub>φ</sub>⟩ [km s<sup>-1</sup>]</td><td>6<sup>+2</sup><sub>-2</sub></td></tr><tr><td>⟨<i>v</i><sub>θ</sub>⟩ [km s<sup>-1</sup>]</td><td>1<sup>+4</sup><sub>-4</sub></td></tr><tr><td>σ<sub>v,r</sub> [km s<sup>-1</sup>]</td><td>95<sup>+2</sup><sub>-2</sub></td></tr><tr><td>σ<sub>v,l</sub> [km s<sup>-1</sup>]</td><td>58<sup>+2</sup><sub>-2</sub></td></tr><tr><td>σ<sub>v,b</sub> [km s<sup>-1</sup>]</td><td>50<sup>+1</sup><sub>-1</sub></td></tr></tbody></table>

**Note.** The marginalized posterior distribution of each parameter is summarized using the 16th, 50th, and 84th quantiles, corresponding to $1\sigma$ uncertainties for a Gaussian distribution.

highly flexible, it also doubles the number of free parameters, and consequently the uncertainties on the measured parameters are larger.

To ensure that we have robustly measured the basic direction and amplitude of the MW’s LMC-induced travel velocity $v_{travel}$, we fit all stars in our sample from 40–120 kpc with a distance-independent $v_{travel}$ model. This model is directly analogous to that used by PP21 and YPP24, in that it uses a single value for each model parameter, instead of a linearly varying distance-dependent value (see Section 4).

The results of the distance-independent fit to all of our data is shown in Figure 20, and the parameters are summarized in Table 4. As expected, the results are broadly consistent in direction and amplitude with the full continuity model fit shown in Figure 6. It is also clear that ignoring the distance-dependence of the reflex motion makes it very challenging to constrain the MW–LMC mass ratio with these measurements. Most stars in the sample are at the lower end of the distance range, but the underlying reflex motion is strongest for the most distant stars. This further motivates our use of a distance-dependent continuity model in Section 4.
Appendix F

Posterior Distribution of Distance-dependent Reflex Model
Parameters

Figure 21 shows the posterior distributions of all 18 parameters from our distance-dependent $v_{\text{travel}}$ continuity model, described in Section 4. These posterior distributions are quantitatively summarized in Table 1.

Figure 22 shows the information that is missing from Figure 21, namely the covariance of each parameter at $r_{\text{gal}} = 40$ kpc versus its counterpart at $r_{\text{gal}} = 120$ kpc. There is modest correlation in all panels, which is to be expected for the same physical parameter at different distances. However, all parameters are well constrained relative to their prior distributions at both distances, confirming that the linear continuity model is well measured over the distance range considered.

<figure><img src="image_23.png" /><figcaption>**Figure 21.** Posterior distribution of 18 model parameters from our distance-dependent $v_{\text{travel}}$ continuity model (see Section 4). Each parameter has a free value at $r_{\text{gal}} = 40$ kpc (blue contours) and at $r_{\text{gal}} = 120$ kpc (red contours), with the parameter assumed to vary linearly between those two values as a function of $r_{\text{gal}}$. Contours correspond to 0.5$\sigma$, 1$\sigma$, and 1.5$\sigma$ Gaussian equivalents.</figcaption></figure>
<figure><img src="image_24.png" /><figcaption>**Figure 22.** Posterior distributions for all parameters in the $v_{\text{travel}}$ continuity model, comparing the covariance between each parameter at $r_{\text{gal}} = 40$ kpc and $r_{\text{gal}} = 120$ kpc. The axes are scaled to span the prior range of each parameter, and the axis limits are indicated in each panel.</figcaption></figure>

## ORCID iDs

Vedant Chandra https://orcid.org/0000-0002-0572-8012

Rohan P. Naidu https://orcid.org/0000-0003-3997-5705

Charlie Conroy https://orcid.org/0000-0002-1590-8551

Nicolas Garavito-Camargo https://orcid.org/0000-0001-7107-1744

Chervin Laporte https://orcid.org/0000-0003-3922-7336

Ana Bonaca https://orcid.org/0000-0002-7846-9787

Phillip A. Cargile https://orcid.org/0000-0002-1617-8917

Emily Cunningham https://orcid.org/0000-0002-6993-0826

Jiwon Jesse Han https://orcid.org/0000-0002-6800-5778

Benjamin D. Johnson https://orcid.org/0000-0002-9280-7594

Hans-Walter Rix https://orcid.org/0000-0003-4996-9069

Yuan-Sen Ting (丁源森) https://orcid.org/0000-0001-5082-9536

Rebecca Woody https://orcid.org/0000-0002-0721-6715

Dennis Zaritsky https://orcid.org/0000-0002-5177-727X

## References

Abadi, M. G., Navarro, J. F., & Steinmetz, M. 2006, MNRAS, 365, 747

Abazajian, K. N., Adelman-McCarthy, J. K., Agüeros, M. A., et al. 2009, ApJS, 182, 543

Amarante, J. A. S., Koposov, S. E., & Laporte, C. F. P. 2024, A&A, 690, A166

Amorisco, N. C. 2017, MNRAS, 464, 2882

Astropy Collaboration, Price-Whelan, A. M., Lim, P. L., et al. 2022, ApJ, 935, 167

Astropy Collaboration, Price-Whelan, A. M., Sipőcz, B. M., et al. 2018, AJ, 156, 123

Astropy Collaboration, Robitaille, T. P., Tollerud, E. J., et al. 2013, A&A, 558, A33

Belokurov, V., Deason, A. J., Erkal, D., et al. 2019, MNRAS: Letters, 488, L47

Belokurov, V., Erkal, D., Evans, N. W., Koposov, S. E., & Deason, A. J. 2018, MNRAS, 478, 611

Belokurov, V., & Koposov, S. E. 2016, MNRAS, 456, 602

Besla, G., Hernquist, L., & Loeb, A. 2013, MNRAS, 428, 2342

Besla, G., Kallivayalil, N., Hernquist, L., et al. 2007, ApJ, 668, 949

Besla, G., Kallivayalil, N., Hernquist, L., et al. 2010, ApJL, 721, 97

Binney, J., & Tremaine, S. 2008, Galactic Dynamics: Second Edition (Princeton, NJ: Princeton Univ. Press)

Bird, S. A., Xue, X.-X., Liu, C., et al. 2021, ApJ, 919, 66

Bland-Hawthorn, J., & Gerhard, O. 2016, ARA&A, 54, 529

Brooks, R. A. N., Garavito-Camargo, N., Johnston, K. V., et al. 2025, ApJ, 978, 79

Brown, A. G., Vallenari, A., Prusti, T., et al. 2021, A&A, 649, 1

Byström, A., Koposov, S. E., Lilleengen, S., et al. 2024, arXiv:2410.09149

Cao, Y.-N., Tian, H., Shao, S., et al. 2023, RAA, 23, 085016

Cargile, P. A., Conroy, C., Johnson, B. D., et al. 2020, ApJ, 900, 28

Cavieres, M., Chanamé, J., Navarrete, C., et al. 2025, ApJ, 983, 83

Chamberlain, K., Price-Whelan, A. M., Besla, G., et al. 2023, ApJ, 942, 18

Chambers, K. C., Magnier, E. A., Metcalfe, N., et al. 2016, arXiv:1612.05560

Chandra, V., Naidu, R. P., Conroy, C., et al. 2023a, ApJ, 951, 26

Chandra, V., Naidu, R. P., Conroy, C., et al. 2023b, ApJ, 956, 110

Chandrasekhar, S. 1943, ApJ, 97, 255

Choi, J., Dotter, A., Conroy, C., et al. 2016, ApJ, 823, 102

Conroy, C., Bonaca, A., Cargile, P., et al. 2019, ApJ, 883, 107

Conroy, C., Naidu, R. P., Garavito-Camargo, N., et al. 2021, Natur, 592, 534

Cooper, A. P., Koposov, S. E., Allende Prieto, C., et al. 2023, ApJ, 947, 37

Correa Magnus, L., & Vasiliev, E. 2022, MNRAS, 511, 2610

Cui, X.-Q., Zhao, Y.-H., Chu, Y.-Q., & Parkin, D. L. 2020, RAA, 12, 1197

Cunningham, E. C., Deason, A. J., Sanderson, R. E., et al. 2019, ApJ, 879, 120

Cunningham, E. C., Garavito-Camargo, N., Deason, A. J., et al. 2020, ApJ, 898, 4

Deason, A. J., Belokurov, V., Koposov, S. E., et al. 2017, MNRAS, 470, 1259

Debattista, V. P., Moore, B., Quinn, T., et al. 2008, ApJ, 681, 1076

Dennefeld, M. 2020, Msngr, 181, 37

Diemand, J., Madau, P., & Moore, B. 2005, MNRAS, 364, 367

Drimmel, R., & Poggio, E. 2018, RNAAS, 2, 210

Erkal, D., Belokurov, V. A., & Parkin, D. L. 2020, MNRAS, 498, 5574

Erkal, D., Boubert, D., Gualandris, A., Evans, N. W., & Antonini, F. 2019a, MNRAS, 483, 2007

Erkal, D., Belokurov, V., Laporte, C. F., et al. 2019b, MNRAS, 487, 2685

Erkal, D., Deason, A. J., Belokurov, V., et al. 2021, MNRAS, 506, 2677

Foote, H. R., Besla, G., Mocz, P., et al. 2023, ApJ, 954, 163

Fukugita, M., Ichikawa, T., Gunn, J. E., et al. 1996, AJ, 111, 1748

Fushimi, K. J., Mosquera, M. E., & Dominguez, M. 2024, A&A, 688, A147

Gaia Collaboration, Vallenari, A., Brown, A. G. A., et al. 2023, A&A, 674, A1

Garavito-Camargo, N., Besla, G., Laporte, C. F. P., et al. 2019, ApJ, 884, 51

Garavito-Camargo, N., Besla, G., Laporte, C. F. P., et al. 2021, ApJ, 919, 109

Gómez, F. A., Besla, G., Carpintero, D. D., et al. 2015, ApJ, 802, 128

GRAVITY Collaboration, Abuter, R., Amorim, A., et al. 2018, A&A, 615, L15

Gunn, J. E., Carr, M., Rockosi, C., et al. 1998, AJ, 116, 3040

Han, J. J., Conroy, C., Johnson, B. D., et al. 2022, AJ, 164, 249

Han, J. J., Conroy, C., Zaritsky, D., et al. 2024, arXiv:2406.12969

Harris, C. R., Millman, K. J., van der Walt, S. J., et al. 2020, Natur, 585, 357

Helmi, A., Veljanoski, J., Breddels, M. A., Tian, H., & Sales, L. V. 2017, A&A, 598, 1

Hunter, J. D. 2007, CSE, 9, 90

Irwin, M. J., Demers, S., & Kunkel, W. E. 1990, AJ, 99, 191

Johnson, B. D., Conroy, C., Naidu, R. P., et al. 2020, ApJ, 900, 103

Johnson, B. D., Leja, J., Conroy, C., & Speagle, J. S. 2021, ApJS, 254, 22

Kallivayalil, N., van der Marel, R. P., Alcock, C., et al. 2006, ApJ, 638, 772

Kallivayalil, N., van der Marel, R. P., Besla, G., Anderson, J., & Alcock, C. 2013, ApJ, 764, 161

Kollmeier, J. A., Zasowski, G., Rix, H.-W., et al. 2017, arXiv:1711.03234

Koposov, S. E., Erkal, D., Li, T. S., et al. 2023, MNRAS, 521, 4936

Kravtsov, A., & Winney, S. 2024, OJAp, 7, 50

Kurucz, R. L. 1970, Atlas: a Computer Program for Calculating Model Stellar Atmospheres 309, SAOSR

Kurucz, R. L., & Avrett, E. H. 1981, Solar Spectrum Synthesis. I. A Sample Atlas from 224 to 300 nm 391, SAOSR

Laporte, C. F., Gómez, F. A., Besla, G., Johnston, K. V., & Garavito-Camargo, N. 2018, MNRAS, 473, 1218

Laporte, C. F. P., Walker, M. G., & Penarrubia, J. 2013, MNRAS, 433, L54

Law, D. R., & Majewski, S. R. 2010a, ApJ, 714, 229

Law, D. R., & Majewski, S. R. 2010b, ApJ, 718, 1128

Li, J., Long, L., Zhong, J., et al. 2023, ApJS, 266, 4

Lilleengen, S., Petersen, M. S., Erkal, D., et al. 2023, MNRAS, 518, 774

Lindegren, L., Klioner, S. A., Hernández, J., et al. 2021, A&A, 649, A2

Loebman, S. R., Valluri, M., Hattori, K., et al. 2018, ApJ, 853, 196

Mainzer, A., Bauer, J., Cutri, R. M., et al. 2014, ApJ, 792, 30
Majewski, S. R., Skrutskie, M. F., Weinberg, M. D., & Ostheimer, J. C. 2003, ApJ, 599, 1082

Malhan, K., Ibata, R. A., Sharma, S., et al. 2022, ApJ, 926, 107

Marshall, J. L., Burles, S., Thompson, I. B., et al. 2008, Proc. SPIE, 7014, 701454

McConnachie, A. W. 2012, AJ, 144, 4

Muñoz, R. R., Côté, P., Santana, F. A., et al. 2018, ApJ, 860, 66

Naidu, R. P., Conroy, C., Bonaca, A., et al. 2020, ApJ, 901, 48

Nidever, D. L., Majewski, S. R., & Butler Burton, W. 2008, ApJ, 679, 432

Peñarrubia, J., Gómez, F. A., Besla, G., Erkal, D., & Ma, Y. Z. 2016, MNRAS: Letters, 456, L54

Peñarrubia, J., & Petersen, M. S. 2021, MNRAS: Letters, 508, L26

Petersen, M. S., Peñarrubia, J., & Jones, E. 2022, MNRAS, 514, 1266

Petersen, M. S., & Peñarrubia, J. 2020, MNRAS: Letters, 494, L11

Petersen, M. S., & Peñarrubia, J. 2021, NatAs, 5, 251

Price-Whelan, A., Sipőcz, B., Lenz, D., et al. 2020, adrn/gala: v1.3, Zenodo, doi:10.5281/zenodo.4159870

Price-Whelan, A. M. 2017, JOSS, 2, 388

Prochaska, J. X., Hennawi, J. F., Westfall, K. B., et al. 2020, JOSS, 5, 2308

Reid, M. J., & Brunthaler, A. 2004, ApJ, 616, 872

Sales, L. V., Navarro, J. F., Lambas, D. G., White, S. D. M., & Croton, D. J. 2007, MNRAS, 382, 1901

Schlafly, E. F., Meisner, A. M., & Green, G. M. 2019, ApJS, 240, 30

Shen, J., Eadie, G. M., Murray, N., et al. 2022, ApJ, 925, 1

Sheng, Y., Ting, Y.-S., Xue, X.-X., Chang, J., & Tian, H. 2024, MNRAS, 534, 2697

Shipp, N., Erkal, D., Drlica-Wagner, A., et al. 2021, ApJ, 923, 149

Skrutskie, M. F., Cutri, R. M., Stiening, R., et al. 2006, AJ, 131, 1163

Speagle, J. S. 2020, MNRAS, 493, 3132

Szentgyorgyi, A., Furesz, G., Cheimets, P., et al. 2011, PASP, 123, 1188

Tamfal, T., Mayer, L., Quinn, T. R., et al. 2021, ApJ, 916, 55

Torrealba, G., Koposov, S. E., Belokurov, V., & Irwin, M. 2016, MNRAS, 459, 2370

Vasiliev, E. 2023, Galax, 11, 59

Vasiliev, E. 2024, MNRAS, 527, 437

Vasiliev, E., & Baumgardt, H. 2021, MNRAS, 505, 5978

Vasiliev, E., Belokurov, V., & Erkal, D. 2021, MNRAS, 501, 2279

Virtanen, P., Gommers, R., Oliphant, T. E., et al. 2020, NatMe, 17, 261

Wang, W., Han, J., Cautun, M., Li, Z., & Ishigaki, M. N. 2020, SCPMA, 63, 109801

Watkins, L. L., van der Marel, R. P., & Bennet, P. 2024, ApJ, 963, 84

Weinberg, M. D. 1995, ApJL, 455, L31

Xue, X. X., Rix, H. W., Ma, Z., et al. 2015, ApJ, 809, 144

Yaaqib, R., Petersen, M. S., & Peñarrubia, J. 2024, MNRAS, 531, 3524

Yanny, B., Rockosi, C., Newberg, H. J., et al. 2009, AJ, 137, 4377

York, D. G., Adelman, J., Anderson, J. E., Jr., et al. 2000, AJ, 120, 1579

Zaritsky, D., Conroy, C., Naidu, R. P., et al. 2020, ApJ, 905, L3

Zhang, L., Xue, X.-X., Yang, C., et al. 2023, AJ, 165, 224