Quantifying the Stellar Halo's Response to the LMC's Infall with Spherical Harmonics
Emily C. Cunningham, Nicolas Garavito-Camargo, Alis J. Deason, Kathryn V. Johnston, Denis Erkal, Chervin F. P. Laporte, Gurtina Besla, Rodrigo Luger, Robyn E. Sanderson
DD RAFT VERSION J UNE
17, 2020Typeset using L A TEX twocolumn style in AASTeX63
Quantifying the Stellar Halo’s Response to the LMC’s Infall with Spherical Harmonics E MILY
C. C
UNNINGHAM , N ICOLAS G ARAVITO -C AMARGO , A LIS
J. D
EASON , K ATHRYN
V. J
OHNSTON ,
4, 1 D ENIS E RKAL , C HERVIN
F. P. L
APORTE , G URTINA B ESLA , R ODRIGO L UGER , AND R OBYN
E. S
ANDERSON
7, 11
Center for Computational Astrophysics, Flatiron Institute, 162 5th Ave., New York, NY 10010, USA Steward Observatory, University of Arizona, 933 North Cherry Avenue, Tucson, AZ 85721, USA Institute for Computational Cosmology, Department of Physics, University of Durham, South Road, Durham DH1 3LE, UK Department of Astronomy, Columbia University, 550 West 120th Street, New York, NY, 10027, USA Department of Physics, University of Surrey, Guildford GU2 7XH, UK Kavli IPMU (WPI), UTIAS, The University of Tokyo, Kashiwa, Chiba 277-8583, Japan Department of Physics and Astronomy, University of Pennsylvania, 209 S 33rd St., Philadelphia, PA 19104, USA
ABSTRACTThe vast majority of the mass in the Milky Way (MW) is in dark matter (DM); we therefore cannot directlyobserve the MW mass distribution, and have to use tracer populations in order to infer properties of the MWDM halo. However, MW halo tracers do not only feel the gravitational influence of the MW itself. Tracers canalso be affected by MW satellites; Garavito-Camargo et al. (2019) (hereafter GC19) demonstrate that the LargeMagellanic Cloud (LMC) induces a density wake in the MW DM, resulting in large scale kinematic patterns inthe MW stellar halo. In this work, we use spherical harmonic expansion (SHE) of the velocity fields of simulatedstellar halos in an effort to disentangle perturbations on large scales (e.g., due to the LMC itself as well as theLMC-induced DM wake) and small scales (due to substructure). Using the GC19 simulations, we demonstratehow the different terms in the SHE of the stellar velocity field reflect the different wake components, and showthat these signatures are a strong function of the LMC mass. An exploration of model halos built from accreteddwarfs (Bullock & Johnston 2005) suggests that stellar debris from massive, recent accretion events can producemuch more power in the velocity angular power spectra than the perturbation from the LMC-induced wake. Wetherefore consider two models for the Sagittarius (Sgr) stream — the most recent, massive accretion event inthe MW apart from the LMC — and find that the angular power on large scales is generally dominated by theLMC-induced wake, even when Sgr is included. We conclude that SHE of the MW stellar halo velocity fieldmay therefore be a useful tool in quantifying the response of the MW DM halo to the LMC’s infall.
Keywords:
Milky Way stellar halo, Milky Way dark matter halo, Milky Way dynamics INTRODUCTIONA fundamental challenge facing the study of the MilkyWay (MW) galaxy is that most of its mass is in dark mat-ter (DM). Because we cannot directly observe the MW’s DMhalo, we must use tracer populations, such as halo stars, glob-ular clusters, and MW satellites, to study the MW’s DM haloindirectly. In particular, the velocity distributions of tracerpopulations can be used to derive estimates of the MW’smass, utilizing methods derived from Jeans (1915) modeling(e.g., Dehnen et al. 2006, Watkins et al. 2009, Gnedin et al.2010, Deason et al. 2012, Eadie et al. 2017, Sohn et al. 2018,Watkins et al. 2019, Wegg et al. 2019). Prior to the era of
Gaia , underpinning essentially all studies of the global kine-matic structure of the MW stellar halo, as well as estimatesof the mass using Jeans modeling, are three key assumptions:that the halo is in equilibrium, isotropic, and phase-mixed.In our current era with access to full phase space information for halo tracers and detailed high resolution simulations, wecan confront the ways in which these assumptions are vio-lated, and use this information to understand our Galaxy ona deeper level.One major source of disequilibrium in the MW is its mostmassive satellite, the Large Magellanic Cloud (LMC). Theclassical picture of the LMC is as a relatively low mass( ∼ M (cid:12) ) satellite orbiting the MW on a T ∼ Gyr or-bit (e.g., Avner & King 1967, Hunter & Toomre 1969, Mu-rai & Fujimoto 1980, Lin & Lynden-Bell 1982, Lin et al.1995, Bekki & Chiba 2005, Mastropietro et al. 2005, Con-nors et al. 2006). However, proper motion measurements ofthe LMC using the Hubble Space Telescope (
HST ) revealedthat the total velocity of the LMC is much higher than pre-viously measured ( v ∼ km s − ; Kallivayalil et al. 2006,2013). This high velocity, near the escape speed of the MW,indicates that the LMC is likely on its first infall, based on a r X i v : . [ a s t r o - ph . GA ] J un backward orbital integrations (Besla et al. 2007, Kallivay-alil et al. 2013) and statistical predictions from cosmologicalsimulations (Boylan-Kolchin et al. 2011, Busha et al. 2011,Gonz´alez et al. 2013, Patel et al. 2017). In addition, thereis mounting evidence that the LMC is more massive thanpreviously thought ( ∼ M (cid:12) ), including arguments basedon models of the Magellanic system (e.g., Besla et al. 2010,2012, Pardy et al. 2018); abundance matching (e.g., Behrooziet al. 2010, Guo et al. 2010, Moster et al. 2010, 2013); themeasured rotation curve of the LMC (van der Marel & Kalli-vayalil 2014); the presence of satellites around the LMC, in-cluding the Small Magellanic Cloud (e.g., Kallivayalil et al.2018, Erkal & Belokurov 2020, Pardy et al. 2019, Patel et al.2020); the timing argument (Pe˜narrubia et al. 2016); and per-turbations in the Orphan Stream (Koposov et al. 2019, Erkalet al. 2019).While concerns about the LMC’s influence on dynamics inthe MW were raised as early as Avner & King (1967), the re-vised picture of the LMC as a massive ( ∼ M (cid:12) ) satelliteapproaching the MW for the first time is increasingly wor-risome for estimates of the MW gravitational potential thatneglect the LMC, as the LMC mass is a significant fractionof the MW halo mass. Several studies of simulations of theLMC’s infall demonstrate that a massive LMC invalidates theassumption of an inertial Galactocentric reference frame, asthe center-of-mass (COM) can be substantially displaced (byas much as 30 kpc; G´omez et al. 2015) from the center ofthe Galaxy, resulting in net motion of the halo with respectto the MW disk. This net COM motion is predicted to be ∼
40 km s − (GC19, Erkal et al. 2019, Petersen & Pe˜narrubia2020); G´omez et al. (2015) find that the net motion could beas high as 75 km s − . Erkal et al. (2020) show that ignoringthe influence of the LMC leads to systematic overestimates(as high as ) of the MW mass when using equilibriummodels. Therefore, when considering the motions of tracerpopulations that we use to study the MW DM halo, we mustaccount for the influence of the LMC in our models.In addition to perturbations as a result of the COM motionof the LMC, MW halo tracers are also predicted to be per-turbed by the LMC-induced DM wake (Garavito-Camargoet al. 2019; hereafter GC19). In the Λ CDM cosmologicalparadigm, host halos are predicted to respond to the infall ofsatellites; this response can be thought of as a gravitational ordensity wake. One component of this wake arises due to localinteractions of particles with the satellite. The satellite trans-fers kinetic energy to nearby resonant particles, creating anoverdensity trailing in its orbit and causing an effective dragforce on the satellite (i.e., dynamical friction; e.g., Chan-drasekhar 1943, White 1983, Tremaine & Weinberg 1984).GC19 refer to this component of the wake as the
Transientresponse , as it is expected to weaken over time. In additionto the Transient response, there is also a global response in the DM halo, resulting in large scale over and under densi-ties in the DM halo (e.g., Weinberg 1989) that can potentiallyeven excite structure in the disk (Weinberg 1998, Weinberg& Blitz 2006 Laporte et al. 2018a). GC19 refer to this com-ponent of the wake as the
Collective response . For the benefitof the reader, we have included these definitions in Table 1.Using detailed N -body simulations, GC19 demonstrated thatthe density wake induced by the infall of the LMC gives riseto distinct, correlated kinematic patterns in the MW stellarhalo.GC19 explored these wake signatures in the context oftwo MW halo models, one isotropic halo and one radiallyvarying, radially anisotropic halo. They find that while thereare similarities in the wake morphology for the two models,there are also key differences: the Transient response is muchstronger for the model with the radially anisotropic halo,whereas the Collective response is stronger for the isotropichalo. Therefore, understanding how the velocity anisotropy β = 1 − σ T /σ R behaves in the MW is important for thepredicted morphology of the LMC-induced DM wake. Untilrelatively recently, our knowledge of the motions of halo trac-ers was limited to one component of motion, the line-of-sight(LOS) velocity; given this major observational constraint, itwas necessary to make assumptions about the tangential mo-tions of stars, and isotropy ( β = 0 ) was the most commonassumption. However, simulations predict that β should be-come increasingly radially biased as a function of radius (see,e.g., Rashkov et al. 2013, Loebman et al. 2018), and β in thesolar neighborhood is radially biased ( β ∼ . − . ; Smithet al. 2009, Bond et al. 2010).We now have full phase space information for distant halotracers, from the Gaia mission and
HST proper motion (PM)studies, and can measure β outside the solar neighborhooddirectly. Estimates of β outside of the solar neighborhoodhave generally found radially biased β , using GCs (Watkinset al. 2019, Sohn et al. 2018) and halo stars (Bird et al. 2019,Lancaster et al. 2019b, Cunningham et al. 2019b).Detecting the halo response to the LMC-induced DM wakewould be an exciting advancement in testing our assumptionsabout the properties of dark matter, as well as providing keyconstraints on the potential of the MW and the mass and or-bital history of the LMC. However, the GC19 simulationsgive predictions for the response in the context of smoothMW DM and stellar halos. In reality, the MW stellar halocontains a wealth of substructure that is not yet phase-mixed,in the form of stellar streams (e.g., Odenkirchen et al. 2001,Newberg et al. 2002, Belokurov et al. 2006, Grillmair 2006,Shipp et al. 2018; also see Newberg & Carlin 2016 for a re-cent review) and stellar clouds (e.g., Newberg et al. 2002,Rocha-Pinto et al. 2004, Juri´c et al. 2008, Li et al. 2016). Inaddition, using a sample of MW halo main sequence turnoffstars from the HALO7D survey (Cunningham et al. 2019a),Cunningham et al. (2019b) observed that the estimated pa-rameters of the velocity ellipsoid (i.e., (cid:104) v φ (cid:105) , σ φ , σ R , σ θ ) weredifferent in the different survey fields; these differing esti-mates could be interpreted as evidence that the halo is notphase-mixed over the survey range ( (cid:104) r (cid:105) = 23 kpc). Theyalso showed maps of the halo velocity anisotropy β in two ha-los from the Latte suite of FIRE-2 simulations (introduced inWetzel et al. 2016), finding that the anisotropy can vary overthe range [ − , across the sky. Some of the variation in the β estimates appeared to correlate with stellar overdensitiesin the halos, indicating that galactic substructure is at leastin part responsible for the different velocity distributions.While some substructure in the halo can be clearly identifiedas overdensities in phase-space and removed from analysis,the presence of velocity substructure in the halo could com-plicate attempts to detect signatures of the LMC-induced DMwake. For example, Belokurov et al. (2019) recently arguedthat the Pisces Overdensity (Sesar et al. 2007, Watkins et al.2009, Nie et al. 2015) might be stars in the wake trailing theLMC in its orbit, because of their net negative radial veloci-ties. However, it remains difficult to conclusively argue thisscenario given that these stars could also be in Galactic sub-structure (or, perhaps, stars that are in substructure and havebeen perturbed by the DM wake).In summary, there is substantial observational evidencethat MW stellar halo is in disequilibrium, on average radi-ally anisotropic, and rife with unphase-mixed substructure,in clear violation of the three central assumptions of equi-librium models. The velocity field of the halo contains infor-mation about the potential of the MW, the dwarf galaxies thatwere consumed as the MW assembled its mass, and the prop-erties of its current most massive perturber, the LMC. How-ever, separating out the different origins of the features inthe MW halo velocity field remains a formidable challenge.One way forward is to consider the spatial scale of perturba-tions: we expect substructure to cause velocity variation onrelatively small spatial scales, as opposed to the large scaleperturbations from the LMC-induced DM wake. Therefore,to disentangle these effects, we seek a quantitative descrip-tion of the kinematic structure of the halo that incorporatesvariation on different spatial scales. Spherical harmonic ex-pansion is a natural tool to address this problem.While ideally we would embark on a full basis function ex-pansion (BFE) of the phase space structure of the halo, in thiswork, we focus on the spherical harmonic expansion of thethree components of motion in spherical coordinates, overdifferent distance ranges in the halo (as a complement to thiswork, Garavito-Camargo et al. 2020, in prep, will presentfull BFEs of the spatial distributions for these simulations).GC19 explored the density structure and the properties of thevelocity dispersions in addition to the mean velocities; wechoose to focus on the mean velocities here because of the challenges of estimating densities (i.e., deeply understand-ing completeness and survey selection functions) and the factthat estimates of the mean of a distribution require fewer trac-ers than dispersion estimates.This paper is organized as follows. In Section 2, we presenta brief overview of spherical harmonic expansion and definethe notation used in the rest of the paper. We then show theresults of using spherical harmonic expansion on the velocityfields from the GC19 simulations of the LMC’s infall into theMW in Section 3. We demonstrate how the different compo-nents of the wake are described in terms of spherical harmon-ics. In Section 4, we investigate how Galactic substructuremight complicate our ability to measure perturbations to thevelocity field as a result of the LMC-induced DM wake, bystudying the Bullock & Johnston (2005) purely accreted stel-lar halos. In Section 5, we use two models of the Sagittariusstream to estimate how the MW’s most massive stream mightinfluence the angular power spectrum of the MW halo veloc-ity field and interfere with signatures from the wake. Wesummarize our conclusions in Section 6. SPHERICAL HARMONICSWe seek to describe the variation on different spatial scalesin halo velocity fields by using spherical harmonic expan-sion. In this section, we define the notation we use through-out the paper for spherical harmonics. As a reference, wehave also included many of these definitions in Table 1.Laplace’s spherical harmonics of order (cid:96) and degree m aredefined as: Y m(cid:96) ( θ, φ ) = (cid:115) (cid:96) + 14 π ( (cid:96) − m )!( l + m )! P m(cid:96) (cos θ ) e imφ , (1)where θ is the colatitudinal angle (i.e., the polar angle mea-sured downward from the north pole) and φ is the azimuthalangle (i.e., the angle in the x − y plane measured from the x − axis), and P m(cid:96) are the associated Legendre polynomials.Spherical harmonics comprise an orthogonal basis for anyfunction f ( θ, φ ) defined on the surface of the sphere: f ( θ, φ ) = ∞ (cid:88) (cid:96) =0 m = (cid:96) (cid:88) m = − (cid:96) a (cid:96)m Y m(cid:96) ( θ, φ ) , (2)where a (cid:96)m are the spherical harmonic coefficients are givenby: a (cid:96)m = (cid:90) Ω f ( θ, φ ) Y m ∗ (cid:96) ( θ, φ )dΩ . (3)When the spherical harmonics are complex, the coeffi-cients are also complex; we define the phase ϕ (cid:96)m of a spher-ical harmonic coefficient as: ϕ (cid:96)m = tan − (cid:18) Im[ a (cid:96)m ]Re[ a (cid:96)m ] (cid:19) , (4) Table 1.
Useful Definitions
LMC-Induced Wake Components from GC19
Transient Response The overdensity trailing the LMC in its orbit, that arises due to local scattering. This component can bethought of in terms of classical dynamical friction.Collective Response Refers both to the overdensity in the north (which arises due to particles in resonance with the LMC’s orbit)and the motion of the MW disk with respect to the new MW-LMC barycenter. The global response of theMW halo to the LMC’s infall.
Spherical Harmonics θ Colatitudinal angle, measured downward from the z − axis; cos θ = z/Rφ Azimuthal angle, angle in the x − y plane measured from the x-axis; tan φ = y/x(cid:96) Order of Spherical Harmonic Y (cid:96)m m Degree of Spherical Harmonic Y (cid:96)m a (cid:96)m Spherical harmonic coefficient for mode (cid:96)mϕ (cid:96)m
Phase of spherical harmonic coefficient a (cid:96)m C (cid:96) Average power in mode of order (cid:96) ; total power is (2 (cid:96) + 1) × C (cid:96) Zonal Spherical Harmonics Spherical harmonics with (cid:96) = 0 ; rotational symmetry about z − axisSectorial Spherical Harmonics Spherical harmonics with (cid:96) = | m | N OTE —We use the function arctan2 in NumPy (Oliphant 2006) to compute azimuthal angle φ and phase ϕ , such that both angles take valuesover the range [ − π, π ] . Figure 1.
The real spherical harmonics, plotted in Mollweide projection, evaluated from (cid:96) = 0 to (cid:96) = 4 , for each − (cid:96) ≤ m ≤ l . In thisprojection, the z − axis is oriented upwards, with colatitudinal angle θ = 0 at the north pole and θ = π at the south pole. The azimuthal angleruns from [ π, − π ] from left to right. The zonal spherical harmonics ( m = 0 ), which have rotational symmetry about the z − axis, are plotted inthe central column. The sectorial harmonics ( l = | m | ) are shown in the outermost panels of each row. The only difference between modes with ± m is a phase shift of 90 degrees. where we use the function arctan2 implemented inNumPy (Oliphant 2006) to compute the inverse tangent, tak-ing into account the quadrant in which a (cid:96)m lies in the com-plex plane.The angular power spectrum C (cid:96) can be computed from thespherical harmonic coefficients a (cid:96)m : C (cid:96) = 12 (cid:96) + 1 (cid:88) m | a (cid:96)m | . (5)The total power in a given order (cid:96) is thus (2 (cid:96) + 1) × C (cid:96) , asthere are (2 (cid:96) + 1) total m values for a given (cid:96) value. There-fore, in this paper, power spectra will always have the quan-tity (2 (cid:96) +1) × C (cid:96) (in units of (km s − ) ), as we are expandingthe velocity field) plotted on the y-axis.In this work we use the Python package healpy (Zoncaet al. 2019) , based on the Healpix scheme (G´orski et al.2005), to perform all of our analysis relating to sphericalharmonics. All maps are made using the healpy plottingroutine mollview ; power spectra and spherical harmoniccoefficients are computed using the function anafast ; andsynthetic maps are generated using synfast .While healpy works with the spherical harmonics incomplex form, the real spherical harmonics are defined as: Y (cid:96)m ( θ, φ ) = P (cid:96)m (cos θ ) cos( mφ ) m ≥ P (cid:96) | m | (cos θ ) sin( | m | φ ) m < . (6)For illustrative purposes, we show the real spherical har-monics from (cid:96) = 0 to (cid:96) = 4 in Figure 1. For a given Y (cid:96)m , thedegree m corresponds to the number of waves along a lineof constant latitude. The order (cid:96) , in conjunction with degree m , determines how many times zero is crossed along a lineof constant longitude: there are (cid:96) − | m | zero crossings alonga meridian. In the case of (cid:96) = | m | , there are no zero cross-ings along the meridian (outer column in each row of Figure1), and a total of (cid:96) complete waves along the equator; thesemodes are referred to as the sectorial spherical harmonics.When m = 0 , there are (cid:96) zero crossings along the merid-ian (central column of Figure 1), and no change in amplitudewith longitude; these are known as the zonal spherical har-monics, and have symmetry about the z − axis. Modes with m < are of sine type, while modes with m > are of co-sine type; as demonstrated by Figure 1, for the real sphericalharmonics, changing the sign of m results in a 90 ◦ rotationabout the z − axis.Spherical harmonic expansion and angular power spectrahave many applications in astrophysics, most famously instudies of the Cosmic Microwave Background (e.g., Planck https://healpy.readthedocs.io/en/latest/ Collaboration et al. 2018, Planck Collaboration et al. 2019,and references therein). While spherical harmonics are com-monly used to describe the angular dependence in full ba-sis function expansions of the potential of dark matter ha-los (e.g., Hernquist & Ostriker 1992; Weinberg 1996, 1999;Lowing et al. 2011), they have not generally been used to de-scribe velocity fields. In the following sections, we discussspherical harmonic expansion of the velocity fields of severaldifferent types of simulations. THE LMC-INDUCED DARK MATTER WAKEIn this section, we perform spherical harmonic expansionof the velocity fields on the high resolution, N -body simu-lations of the LMC’s infall into the MW from GC19. Thekinematic patterns (for mean velocities in all components aswell as the densities and velocity dispersions) are discussedin detail in GC19. Here, we discuss how spherical harmonicexpansion of the mean velocities can be used to characterizethe MW halo response to the LMC’s infall.This section is organized as follows. We summarize keyGC19 simulation properties in Section 3.1. In Section 3.2,we discuss the spherical harmonic expansion of the velocitymaps at 45 kpc in detail. In Section 3.3, we discuss the ra-dial evolution of the power spectrum, for both the isotropicand radially anisotropic MW models. The dependence of thepower spectra on the LMC mass is explored in Section 3.4.3.1. GC19 Simulation Details
For the full description of the numerical methods em-ployed in these simulations, we refer the reader to Section3 of GC19. However, we summarize some of the key detailshere.The GC19 N -body simulations were carried out withTree Smoothed Particle Hydrodynamics code Gadget-3 (Springel et al. 2008), with initial conditions specified by thepublicly available code
GalIC (Yurin & Springel 2014). TheMW model has a virial mass of M MW , vir = 1 . × M (cid:12) ,with a DM halo represented by a Hernquist profile with par-ticle masses m p = 1 . × M (cid:12) . The simulations includea disk and bulge component as well, in order to create a re-alistic potential in the inner halo. GC19 presents results forboth an isotropic halo as well as a halo with radially biasedvelocity anisotropy. In this work, we focus primarily on theradially anisotropic halo, given that simulations and observa-tions agree that the MW halo should be radially biased (see,e.g., Loebman et al. 2018, Bird et al. 2019, Cunningham et al.2019b). However, we do discuss the isotropic MW model inSection 3.3.For the LMC, they construct four models, with virialmasses (prior to infall) of M LMC , vir = 0 . , . , . , . × M (cid:12) . They focus on their fiducial model with M LMC , vir = 1 . × M (cid:12) , which is consistent with LMCmass estimates from abundance matching as well as a firstinfall scenario (see Section 1).When considering these simulations, it is important to keepin mind that only the DM is simulated in time, with thestellar halo constructed in post-processing using a weight-ing scheme (as in Laporte et al. 2013a, 2013b; a generalizedversion of the scheme used in Bullock & Johnston 2005).The stellar halo is constructed to be in equilibrium with theDM halo, given a specified stellar density and velocity dis-persion profile. While in GC19 they construct two stellar ha-los (one using the K-Giant density profile measured in Xueet al. 2015, and one using the density profile as measuredfrom RR Lyrae in Hernitschek et al. 2018), for the purposesof this work, we consider only the stellar halo constructedwith the Xue et al. (2015) density profile.GC19 identify two main components of the wake: theTransient and Collective responses. As discussed in the Sec-tion 1, the Transient response refers to the DM overdensitytrailing the LMC in its orbit, corresponding to the classicalChandrasekhar (1943) wake. The Collective response refersto the global response of the halo to the LMC’s infall (Wein-berg 1989), which results in an extended overdensity in thenorth, as well as the motion of the MW about the new MW-LMC barycenter. As we refer to these two components ofthe wake frequently throughout the remainder of the paper,we have included these definitions in Table 1 as a referencefor the reader.3.2. The Velocity Field Near R LMC
We first discuss the velocity field from the GC19 fiducialLMC model ( M LMC = 1 . × M (cid:12) ) with the radiallyanisotropic MW model (referred to as “Model 2” in GC19)at 45 kpc (in the Galactocentric frame), very near the present-day position of the LMC ( R LMC = 50 kpc). We expectthe velocity maps over this radial range to be sensitive tothe Transient response, given that the LMC passed throughvery recently, and the Transient response arises due to lo-cal scattering. The mean velocity maps in three componentsof spherical motion (with respect to the Galactic center) areshown in the top panels of Figure 2.At z = 0 in these simulations, the LMC is located at 50kpc, and its angular position indicated by the star symbol inthe top panels of Figure 2. The color of the star symbol indi-cates the sign of the motion of the LMC in each component ofmotion: ( v R , v θ , v φ ) LMC = (99 , − , − km s − . At 45kpc, the motions of the stars in the Transient Response nearthe LMC’s position trace the COM motion of the LMC. Thenet motion outwards in v R and the converging motions v θ and v φ near the LMC’s position result from stars in the tran-sient response, accelerating towards the LMC. In addition,stars at 45 kpc are being accelerated towards the overdensityin the North (i.e., the Collective response): this is reflected in the large areas in the north of net positive v R and net negative v θ .The middle panels of Figure 2 show the (cid:96) max = 5 ex-pansion of each component of velocity. Because the kine-matic variation induced by the LMC occurs on large scales,the dominant features in the velocity maps are effectivelycaptured by a low-order spherical harmonic expansion. Thelower panels of Figure 2 show the magnitudes of the spheri-cal harmonic expansion coefficients ( | a (cid:96)m | ), color coded by m value. In the radial velocity map (left panels), the domi-nant term is the (cid:96) = 2 , m = ± mode; this mode capturesthe outward radial motion in the upper left quadrant and thelower right quadrant (near the LMC), as well as the inwardradial motion in the lower left quadrant and upper right quad-rant. In the v θ maps, the monopole term ( (cid:96) = m = 0 ) isdominant, reflecting the net upwards motion of the halo withrespect to the disk. In the v φ maps, the (cid:96) = 2 , m = ± modedominates; the sectorial (cid:96) = 2 mode captures the convergingmotions of stars in the Transient Response near the locationof the LMC. The fact that there is no power at (cid:96) = 0 in v φ is indicative of the fact that the GC19 MW models have nonet rotation. If the MW does have any net rotation (whichhas been observed, but not with high statistical significance;Deason et al. 2017), this would result in power at (cid:96) = 0 in v φ ,which would not interfere with the predicted wake signature.We note that we have not included error bars on the spheri-cal harmonic coefficients in Figure 2, nor do we include errorbars on our estimates of the power spectra in subsequent fig-ures. This is because the simulations are very high resolution,so the statistical errors on the coefficients are very small; thedominant sources of uncertainty here are in the models, notin the noise in the simulations.3.3. Evolution with Distance
Figure 3 shows the velocity maps in the three componentsof spherical motion for this simulation at 45 kpc, 70 kpc, and100 kpc, in the Galactocentric frame. Radial velocity mapsare shown in the left panels, polar motion v θ is plotted in themiddle panels, and azimuthal motion is shown in the rightpanels. As a result of the Collective Response, there is a ra-dial velocity dipole that increases in strength as a functionof distance out into the halo. In addition, the net motion v θ becomes increasingly negative at larger Galactocentric radii.The behavior in v θ and v R can also be represented in termsof v z : in these simulations, while the net motion in the planeof the disk is fairly stable ( (cid:104) v x (cid:105) = (cid:104) v y (cid:105) ∼ km s − over all It is worth noting that while the total power in order (cid:96) is invariant underrotation, the amount of power in a given ± m value is only invariant underrotations about the z − axis. Therefore, our choice to orient these simula-tions in Galactocentric coordinates aligned with the disk is important tokeep in mind, and the dominant m values will be very sensitive to the or-bital history of the LMC. | a m | [ k m / s ] v R v m=0m= ± ± ± ± ± v Figure 2.
Top panels: average velocity maps, computed in a 5 kpc shell centered on R = 45 kpc, from the GC19 fiducial simulation( M LMC = 1 . × M (cid:12) ) with the radially anisotropic MW model. The average radial velocity map (cid:104) v R (cid:105) is shown on the left; average polarvelocity (cid:104) v θ (cid:105) is in the middle panel; and average azimuthal velocity (cid:104) v φ (cid:105) is shown on the right. The angular position of the LMC (located at R = 50 kpc) is indicated by the star. The star is color coded to indicate the sign of the LMC’s velocity in each component of motion; themagnitude of the LMC’s velocity in all components is greater than the range shown by colorobar ( ( v R , v θ , v φ ) = (99 , − , − km s − ).All velocities and positions are computed with respect to the Galactic center. Middle panel: the (cid:96) max = 5 spherical harmonic expansion of thesemaps. A low order spherical harmonic expansion expansion effectively captures the salient features in these velocity maps. Bottom panels:magnitudes of the spherical harmonic coefficients, color coded by degree m . For the radial velocity, the dominant mode is (cid:96) = 2 , m = ± ; thismode captures the net outward motions of particles in the Transient response (near the LMC) as well as the outward motions of particles in theCollective response (in the north). In v θ , the monopole term is dominant ( (cid:96) = m = 0 ), reflecting the net motion of the halo with respect to theMW disk, as a result of the new MW-LMC barycenter. In v φ , the (cid:96) = 2 , m = ± mode dominates; this sectorial mode captures the convergingmotions of particles trapped in the Transient response, moving towards the orbit of the LMC. radii), there is net upwards motion in the halo ( (cid:104) v z (cid:105) > ),which increases as a function of distance from the MW’sdisk. Erkal et al. (2020) show that the MW globular clustersand dwarf satellites also show net motion in v z (and no netmotion in v x , v y ); however, they note the caveat that thesetracers may not be phase mixed in the MW potential. Wenote that the velocity shifts in these simulations result fromtwo sources: the DM overdensity in the north (which getsstronger with distance, because the LMC spends more timeat larger radii) as well as the net acceleration of the MW disktowards the LMC (e.g., Petersen & Pe˜narrubia 2020). Disen-tangling the relative contributions to the overall velocity shiftfrom these two sources is beyond the scope of this work. The power spectra for these maps are shown by the boldlines in Figure 4. The kinematic patterns are concisely sum-marized by the angular power spectra. The radial veloc-ity dipole that increases in magnitude with Galactocentricdistance is reflected by the increasing power in the (cid:96) = 1 modes; the increasing mean polar velocity as a function ofdistance is captured by the increasing power in (cid:96) = 0 (i.e.,the monopole). The power in v φ is strongest at 45 kpc, wherethe stars in the Transient response are closest to the present-day position of the LMC and are accelerated by its COMmotion.The faded lines in Figure 4 are the resulting power spectrafor the GC19 simulation using an isotropic MW halo. As Figure 3.
Mean velocity maps in the three components of motion in spherical coordinates ( v R , v θ , v φ ) for the fiducial GC19 simulation withthe radially anisotropic MW model. Mean velocity maps are computed in 5 kpc shells. Top panels show the velocity maps at 45 kpc; middlepanels show the maps at 70 kpc, and the lower panels show the maps at 100 kpc. As in Figure 2, the angular position of the LMC is indicatedby the star in the top panels, color coded by the sign of the LMC’s velocity in each component of motion. discussed in GC19, the Collective response is stronger in theisotropic simulations. This is reflected in the resulting powerspectra. In the radial velocity power spectrum, at large radii,the magnitude of the velocity dipole is larger for the isotropichalo, as is the magnitude of the v θ monopole, reflecting thestronger Collective response.3.4. LMC Mass Dependence
Because the mass of the LMC is still very uncertain,GC19 simulated the LMC’s infall at four different masses: M LMC , vir = 0 . , . , . , . × M (cid:12) . Figure 5 shows theresulting power spectra for the mean velocities in the threespherical components of motion at 45 kpc, 70 kpc, and 100kpc (top, middle and bottom panels, respectively). The dif-ferent linestyles in each figure represent the different LMCmasses.While the shape of these power spectra at a given distancefor a given component of motion are overall very similar toone another (highlighted by the bottom panel, which showsthe power spectra on a logarithmic scale), the total power ata given (cid:96) value clearly trends with LMC mass. In GC19, by quantifying the “strength” of the wake as the magnitudeof the density fluctuations, they find that the strength of thewake is comparable at 45 kpc for all LMC masses (see theirFigure 25; though this is for the isotropic MW model, whichhas a very weak Transient response). Based on the top panelsof Figure 5, we can see that the power at different (cid:96) valuesincreases strongly with LMC mass even at 45 kpc. We em-phasize that the y − axis labels are different in each panel, toemphasize the effect of changing the LMC mass; we notethat the peak of the power spectrum is largest in v R at all dis-tances, and the v φ power spectrum has the least amount ofpower at all distances.The sensitivity of these signals to the LMC mass shown inFigure 5 emphasizes the significance of the paradigm shiftfrom a M (cid:12) mass LMC to a favored ∼ M (cid:12) massLMC. If the LMC were only M (cid:12) , this would be a factorof eight less massive than the least massive LMC simulatedby GC19; based on the power spectra in Figure 5, we can seethat the signatures of the DM wake in this scenario would bevery weak. Because the LMC is favored to be ∼ of theMW mass, as opposed to ∼ (like the other MW satel- ( + ) × C [ ( k m / s ) ] v R Power Spectrum
Radial MWIsotropic MW 0 1 2 3 4 50100020003000400050006000 v Power Spectrum
D=45 kpcD=70 kpcD=100 kpc 0 1 2 3 4 5050100150200250300350 v Power Spectrum
Figure 4.
Corresponding power spectra for the velocity maps shown in Figure 3. Bold lines show the power spectra for the radially anisotropichalo; faded lines show the power spectra for the isotropic halo (maps not shown). The Collective response causes power in (cid:96) = 1 in v R and (cid:96) = 0 in v θ , both of which increase as a function of distance. The power spectra illustrate that the Collective response is stronger in the isotropichalo. The Transient response is captured by (cid:96) = 2 in v R and v φ . We note that the y − axis ranges are different for each component of motion,to highlight the differences within each component as a function of distance; in particular, the power in v φ is much, much lower than the othertwo components of motion (except at 45 kpc). lites), we cannot ignore its gravitational influence on MWhalo tracers. 3.5. Summary
In summary, low order spherical harmonic expansion isable to capture the salient features in the kinematic patternsthat are predicted to arise as a result of the LMC’s infall. Wehave shown that the shape and magnitude of the power spec-tra depend on the kinematic state of the halo, because therelative strengths of the different wake components dependon the kinematic state of the halo. The overall power is astrong function of the mass of the LMC.However, one key simplification of the GC19 simulationsis that the MW halo model is smooth. The MW stellar halois known observationally to contain substructure: remnantsfrom disrupted dwarf galaxies, consumed by the MW dur-ing its hierarchical formation. In the subsequent section, weinvestigate how substructure due to accreted dwarf galaxiesmight obscure the phenomena described in GC19. SPHERICAL HARMONIC EXPANSION OFACCRETED SUBSTRUCTUREThe GC19 simulations model the MW DM halo (and, asresult, their stellar halo) as smooth; however, we know thatthe MW stellar halo is structured. In this section, we investi-gate how Galactic substructure might complicate our abilityto characterize the LMC-induced DM wake using the spher-ical harmonic expansion of the velocity field, using the Bul-lock & Johnston (2005) suite of simulations of purely ac-creted stellar halos (hereafter BJ05). The publicly available BJ05 simulations are N -body sim-ulations of accreted dwarf galaxies onto a MW-like parentgalaxy. The full suite of simulations consists of 1515 individ-ual accretion events, with a variety of masses, orbital param-eters, and accretion times, that together make up the eleventraditional BJ05 halos. Each disrupted satellite is modeledwith DM particles. The parent galaxy is represented bya time-evolving potential with disk, halo and bulge compo-nents. Johnston et al. (2008) explore in depth how the observ-able properties of the substructure in these simulations arerelated to the properties of their satellite progenitors. Here,we explore the links between a galaxy’s accretion history, theresulting velocity maps, and the angular power spectra of thevelocity field. We note that because there are no DM particlesin the parent galaxy halo, there are no density wakes inducedin the BJ05 simulations. In this section we are only con-cerned with spatially varying mean velocities arising fromdebris from accreted dwarfs.Specifically, we consider the six halos with “artificiallyconstructed” accretion histories discussed in Johnston et al.(2008). While the standard eleven BJ05 halos have accretionhistories from merger trees constructed in Λ CDM cosmolog-ical context, the artificially constructed halos contain debrisfrom accretion events selected from the full BJ05 library thathave the desired properties. While all six halos end up with atotal luminosity L ∼ L (cid:12) , they assemble their halos verydifferently. The six artificial halos are: • Circular (Radial) Orbits: halo assembled only from ac-cretion events with J sat /J circ > . ( J sat /J circ < . )0 ( + ) × C [ ( k m / s ) ] v R v M LMC = 8.0 × 10 MM LMC = 1.0 × 10 MM LMC = 1.8 × 10 MM LMC = 2.5 × 10 M 0 1 2 3 4 50100200300400
45 kpc v ( + ) × C [ ( k m / s ) ]
100 kpc ( + ) × C [ ( k m / s ) ] D=45 kpcD=100 kpc
Figure 5.
Power spectra for the mean velocities for the GC19 simulations with the radially anisotropic MW model and different LMC masses.Left panels show the angular power spectra for the radial velocity (cid:104) v R (cid:105) ; middle panels show the polar velocity (cid:104) v θ (cid:105) ; and right panels show (cid:104) v φ (cid:105) . The first and second rows show power spectra (plotted linearly) for velocities at computed at 45 kpc and 100 kpc, respectively. Differentlinestyles show the range of LMC masses simulated in GC19: M LMC = 0 . , . , . , . × M (cid:12) . We emphasize that the y − axis rangesare different in each panel, to highlight the differences in the power spectra for the different LMC masses. The spherical harmonic coefficientsclearly increase with mass of the LMC. Bottom panel: same as first and second rows, but power spectra are plotted on a logarithmic scale. Theshape of the power spectra are broadly the same for all the different LMC masses, but scale with LMC mass. • Recent (Early) Accretion: halo assembled only fromaccretion events with t acc < Gyr ( t acc > Gyr) • High L sat (Low L sat ): halo assembled only from ac-cretion events with L sat > L (cid:12) ( L sat < L (cid:12) )Figure 6 shows velocity maps of the six artificially con-structed BJ05 halos for stars in the distance range of − kpc (with respect to the Galactic center), projected in Moll-weide coordinates. As in previous figures, left panels areradial velocity ( v R ) maps, polar velocity ( v θ ) maps are in themiddle panels, and right hand panels are azimuthal velocity ( v φ ) maps. We emphasize that the colorbars for these mapsrange from [ − , km s − , a much larger range thanshown for the GC19 simulations; the amplitudes of the ve-locity fluctuations in these maps are much greater than thosedue to the LMC-induced DM wake as seen in GC19. As a re-sult, the signatures from substructure are difficult to comparewith the wake signatures using the maps alone; the angularpower spectra corresponding to these maps, along with theGC19 power spectra, are plotted in Figure 7. We also notethat due to the resolution of the BJ05 simulations, there arepixels that contain no star particles; these pixels are assigned1 Figure 6.
Velocity maps for the six artificially constructed BJ05 halos, for stars with kpc < r < kpc. Lefthand panels show radialvelocity v R ; middle panels show polar velocity v θ ; and righthand panels show azimuthal velocity v φ . Pixels in each map are colored by theaverage velocity. Because of the drastically different accretion histories experienced by these halos, their velocity maps look very different: thehalos that experienced only circular, high- L sat and recent accretion events have many more features than the halos that experienced only radial,low- L sat and early accretion events. (cid:104) v (cid:105) = 0 km s − , consistent with the assumptions ofa equilibrium model.As a result of the different (and extreme) accretion his-tories experienced by these halos, their velocity maps lookvery different. Unsurprisingly, the halo that experienced onlyearly accretion has no features in any components of mo-tion in its velocity maps (bottom panels of Figure 6), giventhat all of its accreted material has had sufficient time to be-come phase mixed. In contrast, the halo that accreted mas-sive, high luminosity satellites has large scale features in allcomponents of motion (second row in Figure 6). It is alsoworth noting that this halo accreted ∼ of its mass withinthe last 8 Gyr (see Figure 7 of Johnston et al. 2008); it hastherefore experienced recent accretion in addition to massiveaccretion. The halo that experienced mostly circular accre-tion events (top row of Figure 6) has low levels of variationin its radial velocity map, with many thin features of nearlyconstant (and approximately km s − ) radial velocity, cor-responding to streams. These streams have more energy intangential than radial motion: they appear bands of nearlyconstant but high velocity in the v θ and v φ maps. The halothat experienced only recent accretion (third row of maps inFigure 6) contains both debris from massive accretion events(as seen by large patches of stars at common velocities) aswell as kinematically cold streams (from recent, lower mass,circular events).The halos built from radial accretion events and low-luminosity accretion events also do not have large scale fea-tures in any components of motion; however, it is importantto keep in mind that these halos also have had relatively qui-escent recent accretion histories. Both halos had assembled ∼ of their mass 8 Gyr ago (see again Figure 7 of John-ston et al. 2008); therefore, any features in these maps aredue to relatively recent, low mass events. A radial stream ap-pears as a small bright spot in the v R maps for both halos;circular streams can be seen as thin bands in the v θ and v φ maps for the low- L sat halo.The angular power spectra corresponding to these veloc-ity maps are plotted in Figure 7. The halos with the mostpower at all (cid:96) values are the halo built from recent accretionand the halo built from high luminosity satellites. The halobuilt from circular accretion events also has high power in v θ and v φ . The thin streams in the low- L sat halo also result insubstantial power ( > km s − ) over many (cid:96) in the threecomponents of motion, though not as much power as the re-cently accreted, high- L sat and circular halos.The halos that experienced only radial accretion events andonly early accretion events both have nearly featureless mapsin v θ and v φ ; while with few to no features in radial velocity,we note that the v R maps appear noticeably noisier than theother components of motion. This is a result of the fact thatthese halos have radially biased velocity anisotropy: their ra- dial velocity dispersions are much greater than their tangen-tial velocity dispersions. The resulting radial velocity mapshave greater fluctuations from pixel to pixel than their tan-gential velocity counterparts. In addition, because there aremany fewer particles in these simulations than in GC19 (eventhough we are looking at a larger radial range in BJ05), thepixel to pixel variation is higher for the BJ05 maps than theGC19 maps. This results in some power in the radial velocitypower spectrum, albeit with less power than the other halos,and with hardly any power in the tangential velocity compo-nents.The purple shaded in region in Figure 7 shows the range ofpower at 45 kpc for the GC19 simulations, from the lowest-mass to the highest-mass LMC. The power from the haloformed through recent accretion and the high- L sat halo ismuch greater than the power from the LMC-induced DMwake at nearly all (cid:96) in all components of motion. Even thecircular and low- L sat halos can have comparable signals tothe wake in the tangential components of motion. However,the shape of the power spectrum is very different for Galac-tic substructure than for the LMC-induced DM wake. Thepower spectra are characterized by a sawtooth pattern, withpeaks at odd (cid:96) values for v R and v θ and peaks at even (cid:96) val-ues for v φ . This sawtooth pattern is indicative of the fact thatspherical harmonics are not the ideal basis for the velocitiesof stars in substructure; we explore why the power spectrahave these features in the Appendix.Based on the power spectra plotted in Figure 7, if the MWstellar halo is dominated by debris from recent, massive ac-cretion events, the signal from the LMC-induced DM wakein the velocity field would be overwhelmed by the Galac-tic substructure. This signal is from debris that has not yetphase-mixed; based on the velocity maps, we can see that thissubstructure is clearly visible as overdensities in phase-space,not just in velocity space. Therefore, one could take advan-tage of the fact that many of these features would be clearlyidentifiable observationally as overdensities, and could be re-moved from the analysis relatively easily.In addition, while it is likely that debris from an early, mas-sive accretion event dominates the inner halo (i.e., the GaiaSausage/Gaia-Enceladus, ∼ Gyr ago; e.g., Belokurovet al. 2018, Helmi et al. 2018) the current consensus is thatthe MW has had a fairly quiescent recent accretion history.This consensus has emerged based on numerous studies, in-cluding studies of the structure and kinematics of stars of theGalactic disk plane (e.g., Gilmore et al. 2002, Hammer et al.2007, Ruchti et al. 2015); the steep stellar density profile be-yond ∼ kpc in the halo (e.g., Deason et al. 2013, Pillepichet al. 2014, Deason et al. 2018); and the amount of substruc-ture in the halo relative to predictions from simulations (e.g.,Lancaster et al. 2019a). An alternative scenario is one inwhich the MW has experienced more recent low luminosity3or radial accretion events with debris that is harder to find ob-servationally; for example, Donlon et al. (2019) suggest thata recent ( ∼ Gyr), radial merger (with M ∼ M (cid:12) ), thatmixes efficiently, could explain the Virgo Overdensity (Vivaset al. 2001). Regardless, the MW’s accretion history is notbelieved to be dominated by recent massive accretion.The major exceptions to the picture of the MW as havinga quiescent (massive) recent merger history are the relativelyrecent accretion of Sagittarius (Sgr; ∼ Gyr ago) and theLMC ( ∼ Gyr ago). In the following section, we explorehow the presence of debris from Sgr might impact the spher-ical harmonic expansion of the halo velocity field. SAGITTARIUSIn Section 3, we discussed in detail the spherical harmonicexpansion of the kinematic variation that arises due to the in-fall of the LMC. In Section 4, we saw that recent, massive ac-cretion events can cause kinematic variation on large scales,which in turn can result in substantial power over many (cid:96) val-ues in the power spectrum. In the MW, the debris from themost recent, massive accretion event (aside from the LMC)is found in the Sgr stream (Ibata et al. 2001). The progenitorof Sgr was a relatively massive, luminous satellite ( L ∼ , M > M (cid:12) ; e.g., Pe˜narrubia et al. 2010, Niederste-Ostholtet al. 2012, Deason et al. 2019; Gibbons et al. 2017 suggestan even higher mass M > × ). Debris from the Sgraccretion event is found all over the sky over Galactocentricdistances ranging from ∼ kpc to ∼ kpc (e.g., Ma-jewski et al. 2003, Belokurov et al. 2014, Hernitschek et al.2017, Sesar et al. 2017). In this section, we investigate howthe presence of Sgr stars might obscure the signal from theLMC-induced DM wake in the velocity field.We consider two different models for the Sgr stream. Thefirst Sgr model we consider is a fit of the Sagittarius streamin the presence of the LMC which we dub the Erkal model.This model uses the same stream fitting machinery of Erkalet al. (2019) which accounts for the reflex motion of theMilky Way due to the LMC. This technique rapidly gener-ates streams using the modified Lagrange Cloud strippingtechnique from Gibbons et al. (2014). For this model, we fitthe radial velocity and distances from Belokurov et al. (2014)and on-sky positions from Belokurov et al. (2006); Koposovet al. (2012) for the bright stream.Motivated by the results of Law & Majewski (2010), wemodel Sgr as a . × M (cid:12) Plummer sphere with a scaleradius of . kpc. The progenitor is rewound for 5 Gyrin the combined presence of the Milky Way and LMC andthen disrupted to the present to form the Sgr stream. For theMilky Way potential, we take the triaxial NFW generaliza-tion from Bowden et al. (2013) which allow for different in-ner and outer density flattenings. We fix the concentration to c = 15 . As a further generalization, we allow for an arbitrary rotation of this triaxial halo so that its axes are not necessarilyaligned with the galactic Cartesian coordinates. We also in-clude a similar disk and bulge to the MWPotential2014 from Bovy (2015): a Miyamoto-Nagai disk (Miyamoto &Nagai 1975) with a mass of . × M (cid:12) with a scale ra-dius of kpc and a scale height of . kpc, and a Hernquistbulge (Hernquist 1990) with a mass of × M (cid:12) and ascale radius of 0.5 kpc. We use the dynamical friction pre-scription of Jethwa et al. (2016) both for the dynamical fric-tion from the Milky Way on the LMC and from the MilkyWay on Sgr. The distance, radial velocity, and proper mo-tions of the Sgr are left as free parameters with priors set byobservations (McConnachie 2012). For the LMC, we give ita fixed position and velocity based on its mean observed dis-tance (Pietrzy´nski et al. 2013), radial velocity (van der Marelet al. 2002), and proper motion (Kallivayalil et al. 2013). Wemodel the LMC as a Hernquist profile (Hernquist 1990) witha scale radius of 25 kpc and a free mass with a uniform priorfrom − × M (cid:12) . In order to account for the fact that Sgrwas initially more massive and had a substantial dark mattercomponent which would have experienced more dynamicalfriction, we have an additional free parameter that increasesthe mass of Sgr by λ DF when computing its dynamical fric-tion. This has a uniform prior between 0 and 20. Thus, alltogether we have 15 free parameters: the mass and scale ra-dius of the NFW profile ( M NFW , r s NFW ), an inner and outerminor and intermediate axis flattening ( q , p , q ∞ , p ∞ ), threeangles to describe the rotation of the triaxial halo, the massof the LMC ( M LMC ), the mass multiplier λ DF , and finallythe proper motions, radial velocity, and distance of Sgr pro-genitor at the present day.We use the MCMC package from Foreman-Mackey et al.(2013) to estimate the parameters of Sgr. We use 100 walkersfor 2000 steps with a 1000 step burn in. The best-fit param-eters require an LMC mass of . × M (cid:12) , flattenings of q = 0 . , p = 0 . , q ∞ = 0 . , p ∞ = 0 . . The massmultiplier λ DF = 9 . suggesting that fitting the Sgr streamrequires more dynamical friction than the low mass we haveassumed would provide. The best-fit Milky Way mass is . × M (cid:12) with a scale radius of . kpc. Althoughthis Milky Way mass is relatively modest, we note that thescale radius is also quite small. Despite the flexibility of thismodel, we note that it does not perfectly match the distancebut most importantly for this work, it gives a good match tothe radial velocity across the sky. A comparison of the modelwith radial velocities from Belokurov et al. (2014) is shownin Figure 13. The positions of the stars for this Sgr modelin the x − z plane, color coded by heliocentric line-of-sightvelocity, are shown in the top left panel of Figure 8.4 ( + ) × C [( k m / s )] v R Power Spectrum
Circular OrbitsRadial OrbitsRecent AccretionEarly AccretionHigh L
Sat
Low L
Sat v Power Spectrum v Power Spectrum M LMC = 0.8 2.5 × 10 M at 45 kpc
Figure 7.
Corresponding angular power spectra for the velocity maps from the BJ05 halos with artificially constructed accretion histories.The purple shaded region indicates the range of power spectra at 45 kpc from GC19, for the full range of simulated LMC masses ( M LMC =0 . − .. × M (cid:12) ). The halos that experienced recent and high- L sat accretion events have more power than the GC19 simulations fornearly all (cid:96) . The halo that experienced only circular accretion events also has more power in the tangential components of motion than theGC19 simulations. The second Sgr model we consider is the publicly avail-able model from Dierickx & Loeb (2017) (hereafter DL17). The x − z positions of stars from this model, again color-coded by heliocentric line-of-sight velocity, are shown inthe righthand panel of Figure 8. In DL17, they first utilizea semi-analytic approach to derive initial conditions for theSgr progenitor, by integrating the equations of motion for-ward in time over 7-8 Gyr and comparing the resulting po-sition and velocity vector to the observed properties of theSgr remnant. They assume virial masses M Sgr = 10 M (cid:12) and M MW = 10 M (cid:12) . To then model the disruption ofSgr, they run an N -body simulation using the derived initialconditions from the semi-analytic approach, modelling botha live MW and Sgr. While this model does reproduce manyof the features of the Sgr stream, including the positions ofstars observed in 2MASS (Majewski et al. 2004) and SDSS(Belokurov et al. 2014), and the large apocentric distancesobserved in Sesar et al. (2017), we emphasize that the N -body simulation is not tuned to fit the observations of the Sgrstream. As a result, certain properties of the stream (for ex-ample, the LOS velocities along the leading arm) are not wellmatched by the data.To investigate how stars from Sgr might impact the powerspectrum of the MW halo’s velocity field, we overlay the twomodels for the Sgr stream on to the fiducial anisotropic GC19simulation (with M LMC = 1 . × M (cid:12) ). To combine thetwo independent simulations, we assign the total Sgr stellarmass to be 10 % of the total stellar mass of the GC19 halo.This ratio is consistent with current estimates of the total stel-lar mass of the MW ( ∼ M (cid:12) ; e.g., Deason et al. 2019, https://mdierick.github.io/project2.html Mackereth & Bovy 2020) and Sgr ( M Sgr , ∗ ∼ M (cid:12) ; e.g.,Deason et al. 2019, Niederste-Ostholt et al. 2012).The resulting velocity maps at 45 kpc of the two Sgr mod-els overlaid on the GC19 simulations are shown in the mid-dle panels (for the Erkal model) and lower panels (for theDL17 model) of Figure 8. The corresponding power spec-tra for the velocity maps in Figure 8 are shown by the thicksolid lines in Figure 9. We only show the full power spec-tra at 45 kpc, as Sgr does not substantially contribute to theoverall power spectrum at larger radii (with the exception ofthe DL17 model at 70 kpc; this results from stars acceler-ating towards and away from the stream apocenters, whichare at larger distances for the DL17 model than for the Erkalmodel). From left to right, power spectra for v R , v θ , v φ areplotted; the thick solid lines are the power spectra from thecombination of the GC19 simulation with the Sgr models(top panels show the results when using the Erkal model;lower panels show the results from the DL17 model) at 45kpc. The dashed line shows the power spectrum from theGC19 simulation (excluding Sgr). Dotted dashed lines showthe difference between the halo including Sgr and exclud-ing Sgr (i.e., the contribution of Sgr to the overall powerspectrum), at 45 kpc (purple), 70 kpc (orange), and 100 kpc(blue), computed in 5 kpc shells.The dot-dashed lines in Figure 9, representing the contri-bution to the power spectrum due to Sgr, have a similar mor-phology to the power spectra discussed in Section 4 and inthe Appendix: they are characterized by a saw-tooth pattern,with peaks at odd (cid:96) values in v R and v θ . Figure 9 showsthat the two Sgr models result in different signatures. Forthe Erkal model, including Sgr increases the peak at (cid:96) = 1 in v R , while the v φ power spectrum is mostly unaffected.The DL17 model hardly affects the low (cid:96) power in v R , while5 H e li o c e n t r i c L O S V e l o c i t y ( k m / s ) Figure 8.
Top panels: positions in the x − z plane, color coded by heliocentric LOS velocity, for stars in the Erkal Sgr model (left) and theDL17 Sgr Model (right). Dashed lines indicate 45 kpc, 70 kpc, and 100 kpc. Middle panels: velocity maps for the Erkal Sgr model overlaidonto the GC19 stellar halo, in a 5 kpc shell centered at 45 kpc. The angular position of the LMC is indicated by the star symbol, color codedto show the sign of the LMC’s velocity in each component of motion. We note that we have restricted the velocity colorbar to ± km s − , sothat the wake signatures are still visible by eye; as can be seen in the colorbar in the top panels, the range of velocities of stars in Sgr is muchgreater than the range of mean velocities in the GC19 halo at 45 kpc. In addition, we note that the overdense region in the North in these mapsis not the Sgr progenitor, but rather part of the leading arm. Lower panels: same as middle panels, but for the DL17 Sgr model overlaid on theGC19 halo. the power in v φ is slightly enhanced. Both models result inhigher power in v θ at all (cid:96) ; at 45 kpc, v θ is the only compo-nent of motion for which the signal from Sgr is comparableto the signal of the LMC-induced DM wake.While including Sgr does increase the overall power in or-ders (cid:96) that contain signatures of the LMC-induced DM wake,the features sensitive to the Collective response (the (cid:96) = 1 peak in v R as well as the monopole (cid:96) = 0 peak in v θ ) arestronger at all distances than the power due to both Sgr mod-els alone (the dot-dashed lines in Figure 9; though we notethat the DL17 model does substantially increase the power ofthe monopole in v θ ). In addition, the signatures from the Col-lective response increase as a function of distance; the overall power from Sgr generally decreases as a function of distance(with the exception of the increase in signal at 70 kpc in theDL17 model in v θ ). While the signal from Sgr in these key (cid:96) values does not overwhelm the signal from the LMC, itdoes contribute to the overall power in the orders sensitiveto the wake signatures at 45 kpc. Sgr also contributes powerin specific l, m modes that are sensitive to the LMC-inducedDM wake (e.g., l = 1 , m = 0 in v R ); however, the phasesof the different coefficients are not strongly affected by Sgrin the low (cid:96) values. Modeling the influence of the inclusionof Sgr will be essential in quantifying the strengths of theLMC-induced wake components observationally at smallerGalactocentric distances.6In summary, while Sgr stars will affect the power spectraof the MW halo velocity field at smaller Galactocentric dis-tances, the signatures from the LMC-induced DM wake aspredicted by GC19 should still be distinguishable from theSgr signatures. The primary signature of the Transient re-sponse (power in (cid:96) = 2 in v φ ) is largely unaffected by theinclusion of Sgr for both models. In v R and v φ , the domi-nant features in the power spectra that arise due to the LMC-induced DM wake are stronger than the features arising dueto the inclusion of Sgr, and the power spectra have differ-ent morphologies. Because Sgr stars are likely to contributepower at (cid:96) = 1 in v R and (cid:96) = 0 in v θ , modeling Sgr will beimportant in characterizing the Collective response at closerdistances in the halo ( ∼ < kpc).In addition, Sgr has been extensively studied, and ourknowledge of its velocity structure is better known than everbefore with the release of Gaia DR2 (e.g., Antoja et al. 2020,Ramos et al. 2020, Ibata et al. 2020). Like the substructurestudied in BJ05, the debris from Sgr is also well known tobe overdense on the sky, and many halo studies remove starsbelieved to be associated with the stream. Therefore, giventhat we have shown that the wake signatures should still beidentifiable even if all Sgr stars are included in the analysis,the prospects for observing these signatures only improve ifdebris from Sgr is subtracted, even if this subtraction is im-perfect.However, we note that we have only considered the effectsof including Sgr stars in the analysis of the stellar halo, andnot the effects the infall of the Sgr dwarf may have had onthe MW disk and halo over the last ∼ − Gyr. Dependingon the assumed mass of Sgr, it may have resulted in a sub-stantial shift in the MW barycenter and also caused a wakein the MW DM halo (Laporte et al. 2018b). GC19 comparethe magnitude of the density perturbations that arise due toSgr and the LMC using the Laporte et al. (2018b) simula-tions, and find that the contribution from Sgr is negligible(see GC19’s Figure 26); however, they did not discuss theperturbations to the velocity field. We leave the explorationof these effects to future work. CONCLUSIONSIn this paper, we use spherical harmonic expansion to de-scribe the perturbed velocity fields of the MW as a result ofthe LMC’s infall using the simulations from GC19. We ex-plore the ways in which Galactic substructure might obscurethe signatures from the wake in the power spectrum, using theBJ05 simulations as well as two models for the Sgr stream.We summarize our primary findings as follows:1. We study the perturbation to the velocity field causedby the LMC-induced DM wake using the simulationsfrom GC19. We find that low-order spherical har-monic expansion of the velocity field in these simula- tions usefully captures the salient features of the LMC-induced DM wake. We found that increasing powerwith Galactocentric radius in (cid:96) = 1 in v R and (cid:96) = 0 in v θ are signatures of the Collective response. At 45kpc (near the LMC), power in (cid:96) = 2 in v R and v φ aresignatures of the Transient response. We find that theamplitude of the power spectra scale with LMC mass.2. We investigate how Galactic substructure might affectthe angular power spectrum of the MW’s velocity field,using the BJ05 simulations with artificially constructedaccretion histories. We find that massive, recent ac-cretion causes large scale, high amplitude fluctuationsin the velocity field. Velocity substructure arising dueto debris from recent, massive satellites creates muchmore power in the power spectrum of the velocity fieldthan the perturbation due to the LMC. However, theMW is not believed to have experienced much recent,massive accretion, with the exceptions of Sgr and theLMC itself.3. Given that Sgr is the most recent, massive accretionevent experienced by the MW (with the exception ofthe LMC), we investigate how Sgr stars could impactmeasurements of the overall MW power spectra andour ability to measure the signatures associated withthe LMC-induced DM wake. The power spectrum onlarge scales (i.e., low (cid:96) ) remains generally dominatedby the signatures from the LMC-induced DM wake;this result complements the GC19 findings that the am-plitude of the density wake induced by the LMC’s in-fall is much greater than the density wake induced bySgr. In addition, overall power due to Sgr decreases asa function of distance, in contrast to the Collective re-sponse signatures. However, including Sgr stars doesincrease the power in modes that are sensitive to theCollective response, especially at 45 kpc. Care shouldtherefore be taken to model the impact of Sgr stars onthe power spectrum in studies attempting to use thismethod for detecting and quantifying the Collective re-sponse.Based on our findings, performing spherical harmonic ex-pansion in the MW velocity field could be a method foridentifying and characterizing the LMC-induced DM wake,which would in turn provide constraints on the mass and or-bital history of the LMC. There remain technical challengesassociated with implementing this technique with observa-tional data (e.g., incorporating measurement uncertainties,limited sky coverage of spectroscopic programs, combiningdata from different surveys); we leave a detailed explorationof how to estimate the spherical harmonic expansion coeffi-cients from realistic data to future work. However, the future7 ( + ) × C [ ( k m / s ) ] v R Power Spectrum v Power Spectrum v Power Spectrum
GC19+Erkal SgrGC19 OnlyDifference, 45 kpcDifference, 70 kpcDifference, 100 kpc0 2 4 6 8 100200400600800 ( + ) × C [ ( k m / s ) ] v R Power Spectrum v Power Spectrum v Power Spectrum
GC19+DL17 SgrGC19 OnlyDifference, 45 kpcDifference, 70 kpcDifference, 100 kpc
Figure 9.
Power spectra for the Erkal Sgr model (top panels) and the DL17 Sgr model (lower panels) overlaid onto the GC19 anisotropicsimulation with M LMC = 1 . × M (cid:12) . Solid lines show the resulting power spectra at 45 kpc when the two simulations are combined;dashed lines show the power spectra from GC19 simulation alone at 45 kpc. The difference between the resulting power spectra are plottedas dot-dashed lines, at 45 kpc (purple), 70 kpc (orange), and 100 kpc (blue). The power at low (cid:96) (i.e., large spatial scales), remains generallydominated by the signatures of the LMC-induced DM wake, especially at larger distances (though at 45 kpc, the power due to Sgr in v θ iscomparable to the power due to the wake). However, Sgr does substantially contribute to modes that are sensitive to the LMC-induced wake(e.g., (cid:96) = 1 in v R , (cid:96) = 0 in v θ ); the influence of Sgr stars should therefore be modeled if quantifying the strength of the wake using SHE. is bright given the upcoming observational programs that willmap our Galaxy’s phase-space structure. The first two Gaiadata releases have already transformed our understanding ofthe MW’s kinematic structure; with future releases from Gaiamission in conjunction with the Gaia spectroscopic follow-up programs (e.g., DESI, 4MOST, WEAVE), as well as fu-ture astrometric (e.g., Rubin Observatory Legacy Survey ofSpace and Time, WFIRST) and spectroscopic programs (e.g.,SDSS-V MW Mapper), the halo velocity field will be betterknown than ever before.While the GC19 only include the MW and the LMC, theyreveal the complex behavior that arises in the phase spacestructure of the MW halo due to the infall of the LMC. How-ever, many complications remain to be addressed. While weexplored how Galactic substructure might obscure the signalsfrom the wake, the assumption of smoothness of the MW halo in GC19, as well as their mapping of the stellar haloon to the DM halo, remain important to keep in mind. Whilemuch of the stellar halo is phase-mixed in the inner regions ofthe MW, at larger distances (e.g., r > kpc), we may needto rely on debris that is not phase mixed in order to see wakesignatures. In addition, the stellar halo may not be in equi-librium with the DM halo for several reasons. First of all,simulations show that accretion is not the only mechanismby which stellar halos form: simulated halos show substan-tial fractions of stars that formed in-situ (e.g., Zolotov et al.2009, 2010). Yu et al. (2020) show that in-situ stars (formedin outflows of the host galaxy) can be ejected to large dis-tances in the halo, and can comprise a substantial fraction (5-40%) of the stars in outer halos (50-300 kpc). Bonaca et al.(2017) use the kinematics of stars from the Gaia data releaseto argue that the MW halo has an in-situ component. Second,8because stars are more concentrated within halos than DM,the DM is preferentially stripped initially (e.g., Smith et al.2016). Third, the dynamical mass-to-light ratios of dwarfgalaxies have been observed to vary over orders of mag-nitude (e.g., McConnachie 2012): the relative amounts ofmass accreted in stars and DM is therefore not constant. Fi-nally, a significant fraction of the DM accretion is “smooth”(e.g. Angulo & White 2010; Genel et al. 2010), in contrastto the relatively “lumpy” accretion that builds up the stellarhalo. While methods have been developed using cosmologi-cal simulations to determine DM halo distributions based onstellar halo distributions (Necib et al. 2019), these complexeffects were not included in the creation of stellar halos fromDM halos in GC19. In addition, while GC19 varies the massof the LMC, they assume a single mass for the MW and donot simulate a range of orbital histories.Many of these complications can be addressed by study-ing DM wakes more generally in cosmological simulations.These simulations contain both in-situ and accreted halocomponents and have experienced a variety of accretion his-tories. By applying this technique to encounters betweenMW-like galaxies and massive satellites in cosmological sim-ulations, we can identify wake signatures for a range of massratios and orbital histories. In the present era of wide-field 3Dkinematic datasets coupled with detailed high resolution sim-ulations, we can move beyond equilibrium models and de-velop new methods for characterizing the complex processesthat have shaped our Galaxy’s formation. ACKNOWLEDGMENTSWe thank the anonymous referee for their helpful com-ments and suggestions. ECC and RL are supported by Flat-iron Research Fellowships at the Flatiron Institute. The Flat-iron Institute is supported by the Simons Foundation. ADis supported by a Royal Society University Research Fel-lowship. KVJ’s contributions were supported by NSF grantAST-1715582. NGC and GB acknowledge support fromHST grant AR 15004 and NASA ATP grant 17- ATP17-0006.CFPL’s contributions are supported by world premier inter-national research center initiative (WPI), MEXT, Japan. Thisproject was developed in part at the 2019 Santa Barbara GaiaSprint, hosted by the Kavli Institute for Theoretical Physicsat the University of California, Santa Barbara. This researchwas supported in part at KITP by the Heising-Simons Foun-dation and the National Science Foundation under Grant No.NSF PHY-1748958. ECC would like to thank Andrew Wet-zel, David Hogg, Adrian Price-Whelan and David Spergelfor helpful scientific discussions.
Software:
Astropy (Astropy Collaboration et al. 2013,2018), Healpy (Zonca et al. 2019), IPython (Perez & Granger2007), Matplotlib (Hunter 2007), Numpy (Oliphant 2006,van der Walt et al. 2011), Pandas (McKinney 2010), Scipy(Virtanen et al. 2020), Starry (Luger et al. 2019)REFERENCES
Angulo, R. E., & White, S. D. M. 2010, MNRAS, 401, 1796,doi: 10.1111/j.1365-2966.2009.15742.xAntoja, T., Ramos, P., Mateu, C., et al. 2020, A&A, 635, L3,doi: 10.1051/0004-6361/201937145Astropy Collaboration, Robitaille, T. P., Tollerud, E. J., et al. 2013,A&A, 558, A33, doi: 10.1051/0004-6361/201322068Avner, E. S., & King, I. R. 1967, AJ, 72, 650, doi: 10.1086/110288Behroozi, P. S., Conroy, C., & Wechsler, R. H. 2010, ApJ, 717,379, doi: 10.1088/0004-637X/717/1/379Bekki, K., & Chiba, M. 2005, MNRAS, 356, 680,doi: 10.1111/j.1365-2966.2004.08510.xBelokurov, V., Deason, A. J., Erkal, D., et al. 2019, MNRAS, 488,L47, doi: 10.1093/mnrasl/slz101Belokurov, V., Erkal, D., Evans, N. W., Koposov, S. E., & Deason,A. J. 2018, MNRAS, 478, 611, doi: 10.1093/mnras/sty982Belokurov, V., Zucker, D. B., Evans, N. W., et al. 2006, ApJL, 642,L137, doi: 10.1086/504797Belokurov, V., Koposov, S. E., Evans, N. W., et al. 2014, MNRAS,437, 116, doi: 10.1093/mnras/stt1862Besla, G., Kallivayalil, N., Hernquist, L., et al. 2007, ApJ, 668,949, doi: 10.1086/521385 —. 2010, ApJL, 721, L97, doi: 10.1088/2041-8205/721/2/L97—. 2012, MNRAS, 421, 2109,doi: 10.1111/j.1365-2966.2012.20466.xBird, S. A., Xue, X.-X., Liu, C., et al. 2019, AJ, 157, 104,doi: 10.3847/1538-3881/aafd2eBonaca, A., Conroy, C., Wetzel, A., Hopkins, P. F., & Kereˇs, D.2017, ApJ, 845, 101, doi: 10.3847/1538-4357/aa7d0cBond, N. A., Ivezi´c, ˇZ., Sesar, B., et al. 2010, ApJ, 716, 1,doi: 10.1088/0004-637X/716/1/1Bovy, J. 2015, The Astrophysical Journal Supplement Series, 216,29, doi: 10.1088/0067-0049/216/2/29Bowden, A., Evans, N. W., & Belokurov, V. 2013, MNRAS, 435,928, doi: 10.1093/mnras/stt1253Boylan-Kolchin, M., Besla, G., & Hernquist, L. 2011, MNRAS,414, 1560, doi: 10.1111/j.1365-2966.2011.18495.xBullock, J. S., & Johnston, K. V. 2005, ApJ, 635, 931,doi: 10.1086/497422Busha, M. T., Marshall, P. J., Wechsler, R. H., Klypin, A., &Primack, J. 2011, ApJ, 743, 40,doi: 10.1088/0004-637X/743/1/40 Chandrasekhar, S. 1943, ApJ, 97, 255, doi: 10.1086/144517Connors, T. W., Kawata, D., & Gibson, B. K. 2006, MNRAS, 371,108, doi: 10.1111/j.1365-2966.2006.10659.xCunningham, E. C., Deason, A. J., Rockosi, C. M., et al. 2019a,ApJ, 876, 124, doi: 10.3847/1538-4357/ab16cbCunningham, E. C., Deason, A. J., Sanderson, R. E., et al. 2019b,ApJ, 879, 120, doi: 10.3847/1538-4357/ab24cdDeason, A. J., Belokurov, V., Evans, N. W., & An, J. 2012,MNRAS, 424, L44, doi: 10.1111/j.1745-3933.2012.01283.xDeason, A. J., Belokurov, V., Evans, N. W., & Johnston, K. V.2013, ApJ, 763, 113, doi: 10.1088/0004-637X/763/2/113Deason, A. J., Belokurov, V., & Koposov, S. E. 2018, ApJ, 852,118, doi: 10.3847/1538-4357/aa9d19Deason, A. J., Belokurov, V., Koposov, S. E., et al. 2017, MNRAS,470, 1259, doi: 10.1093/mnras/stx1301Deason, A. J., Belokurov, V., & Sanders, J. L. 2019, MNRAS, 490,3426, doi: 10.1093/mnras/stz2793Dehnen, W., McLaughlin, D. E., & Sachania, J. 2006, MNRAS,369, 1688, doi: 10.1111/j.1365-2966.2006.10404.xDierickx, M. I. P., & Loeb, A. 2017, ApJ, 836, 92,doi: 10.3847/1538-4357/836/1/92Donlon, Thomas, I., Newberg, H. J., Weiss, J., Amy, P., &Thompson, J. 2019, ApJ, 886, 76,doi: 10.3847/1538-4357/ab4f72Eadie, G. M., Springford, A., & Harris, W. E. 2017, ApJ, 835, 167,doi: 10.3847/1538-4357/835/2/167Erkal, D., Belokurov, V., & Parkin, D. L. 2020, arXiv e-prints,arXiv:2001.11030. https://arxiv.org/abs/2001.11030Erkal, D., & Belokurov, V. A. 2020, MNRAS, 495, 2554,doi: 10.1093/mnras/staa1238Erkal, D., Belokurov, V., Laporte, C. F. P., et al. 2019, MNRAS,487, 2685, doi: 10.1093/mnras/stz1371Foreman-Mackey, D., Hogg, D. W., Lang, D., & Goodman, J.2013, PASP, 125, 306, doi: 10.1086/670067Garavito-Camargo, N., Besla, G., Laporte, C. F. P., et al. 2019,ApJ, 884, 51, doi: 10.3847/1538-4357/ab32ebGenel, S., Bouch´e, N., Naab, T., Sternberg, A., & Genzel, R. 2010,ApJ, 719, 229, doi: 10.1088/0004-637X/719/1/229Gibbons, S. L. J., Belokurov, V., & Evans, N. W. 2014, MNRAS,445, 3788, doi: 10.1093/mnras/stu1986—. 2017, MNRAS, 464, 794, doi: 10.1093/mnras/stw2328Gilmore, G., Wyse, R. F. G., & Norris, J. E. 2002, ApJL, 574, L39,doi: 10.1086/342363Gnedin, O. Y., Brown, W. R., Geller, M. J., & Kenyon, S. J. 2010,ApJL, 720, L108, doi: 10.1088/2041-8205/720/1/L108G´omez, F. A., Besla, G., Carpintero, D. D., et al. 2015, ApJ, 802,128, doi: 10.1088/0004-637X/802/2/128Gonz´alez, R. E., Kravtsov, A. V., & Gnedin, N. Y. 2013, ApJ, 770,96, doi: 10.1088/0004-637X/770/2/96 G´orski, K. M., Hivon, E., Banday, A. J., et al. 2005, ApJ, 622, 759,doi: 10.1086/427976Grillmair, C. J. 2006, ApJL, 645, L37, doi: 10.1086/505863Guo, Q., White, S., Li, C., & Boylan-Kolchin, M. 2010, MNRAS,404, 1111, doi: 10.1111/j.1365-2966.2010.16341.xHammer, F., Puech, M., Chemin, L., Flores, H., & Lehnert, M. D.2007, ApJ, 662, 322, doi: 10.1086/516727Helmi, A., Babusiaux, C., Koppelman, H. H., et al. 2018, Nature,563, 85, doi: 10.1038/s41586-018-0625-xHernitschek, N., Sesar, B., Rix, H.-W., et al. 2017, ApJ, 850, 96,doi: 10.3847/1538-4357/aa960cHernitschek, N., Cohen, J. G., Rix, H.-W., et al. 2018, ApJ, 859,31, doi: 10.3847/1538-4357/aabfbbHernquist, L. 1990, ApJ, 356, 359, doi: 10.1086/168845Hernquist, L., & Ostriker, J. P. 1992, ApJ, 386, 375,doi: 10.1086/171025Hunter, C., & Toomre, A. 1969, ApJ, 155, 747,doi: 10.1086/149908Hunter, J. D. 2007, Computing in Science and Engineering, 9, 90,doi: 10.1109/MCSE.2007.55Ibata, R., Bellazzini, M., Thomas, G., et al. 2020, ApJL, 891, L19,doi: 10.3847/2041-8213/ab77c7Ibata, R., Irwin, M., Lewis, G. F., & Stolte, A. 2001, ApJL, 547,L133, doi: 10.1086/318894Jeans, J. H. 1915, MNRAS, 76, 70, doi: 10.1093/mnras/76.2.70Jethwa, P., Erkal, D., & Belokurov, V. 2016, MNRAS, 461, 2212,doi: 10.1093/mnras/stw1343Johnston, K. V., Bullock, J. S., Sharma, S., et al. 2008, ApJ, 689,936, doi: 10.1086/592228Juri´c, M., Ivezi´c, ˇZ., Brooks, A., et al. 2008, ApJ, 673, 864,doi: 10.1086/523619Kallivayalil, N., van der Marel, R. P., Alcock, C., et al. 2006, ApJ,638, 772, doi: 10.1086/498972Kallivayalil, N., van der Marel, R. P., Besla, G., Anderson, J., &Alcock, C. 2013, ApJ, 764, 161,doi: 10.1088/0004-637X/764/2/161Kallivayalil, N., Sales, L. V., Zivick, P., et al. 2018, ApJ, 867, 19,doi: 10.3847/1538-4357/aadfeeKoposov, S. E., Belokurov, V., Evans, N. W., et al. 2012, ApJ, 750,80, doi: 10.1088/0004-637X/750/1/80Koposov, S. E., Belokurov, V., Li, T. S., et al. 2019, MNRAS, 485,4726, doi: 10.1093/mnras/stz457Lancaster, L., Belokurov, V., & Evans, N. W. 2019a, MNRAS,484, 2556, doi: 10.1093/mnras/stz124Lancaster, L., Koposov, S. E., Belokurov, V., Evans, N. W., &Deason, A. J. 2019b, MNRAS, 486, 378,doi: 10.1093/mnras/stz853Laporte, C. F. P., G´omez, F. A., Besla, G., Johnston, K. V., &Garavito-Camargo, N. 2018a, MNRAS, 473, 1218,doi: 10.1093/mnras/stx2146 Laporte, C. F. P., Johnston, K. V., G´omez, F. A.,Garavito-Camargo, N., & Besla, G. 2018b, MNRAS, 481, 286,doi: 10.1093/mnras/sty1574Laporte, C. F. P., Walker, M. G., & Penarrubia, J. 2013a, MNRAS,433, L54, doi: 10.1093/mnrasl/slt057Laporte, C. F. P., White, S. D. M., Naab, T., & Gao, L. 2013b,MNRAS, 435, 901, doi: 10.1093/mnras/stt912Law, D. R., & Majewski, S. R. 2010, ApJ, 714, 229,doi: 10.1088/0004-637X/714/1/229Li, T. S., Balbinot, E., Mondrik, N., et al. 2016, ApJ, 817, 135,doi: 10.3847/0004-637X/817/2/135Lin, D. N. C., Jones, B. F., & Klemola, A. R. 1995, ApJ, 439, 652,doi: 10.1086/175205Lin, D. N. C., & Lynden-Bell, D. 1982, MNRAS, 198, 707,doi: 10.1093/mnras/198.3.707Loebman, S. R., Valluri, M., Hattori, K., et al. 2018, ApJ, 853, 196,doi: 10.3847/1538-4357/aaa0d6Lowing, B., Jenkins, A., Eke, V., & Frenk, C. 2011, MNRAS, 416,2697, doi: 10.1111/j.1365-2966.2011.19222.xLuger, R., Agol, E., Foreman-Mackey, D., et al. 2019, AJ, 157, 64,doi: 10.3847/1538-3881/aae8e5Mackereth, J. T., & Bovy, J. 2020, MNRAS, 139,doi: 10.1093/mnras/staa047Majewski, S. R., Skrutskie, M. F., Weinberg, M. D., & Ostheimer,J. C. 2003, ApJ, 599, 1082, doi: 10.1086/379504Majewski, S. R., Kunkel, W. E., Law, D. R., et al. 2004, AJ, 128,245, doi: 10.1086/421372Mastropietro, C., Moore, B., Mayer, L., Wadsley, J., & Stadel, J.2005, MNRAS, 363, 509,doi: 10.1111/j.1365-2966.2005.09435.xMcConnachie, A. W. 2012, AJ, 144, 4,doi: 10.1088/0004-6256/144/1/4McKinney, W. 2010, in Proceedings of the 9th Python in ScienceConference, ed. S. van der Walt & J. Millman, 51 – 56Miyamoto, M., & Nagai, R. 1975, PASJ, 27, 533Moster, B. P., Naab, T., & White, S. D. M. 2013, MNRAS, 428,3121, doi: 10.1093/mnras/sts261Moster, B. P., Somerville, R. S., Maulbetsch, C., et al. 2010, ApJ,710, 903, doi: 10.1088/0004-637X/710/2/903Murai, T., & Fujimoto, M. 1980, PASJ, 32, 581Necib, L., Lisanti, M., Garrison-Kimmel, S., et al. 2019, ApJ, 883,27, doi: 10.3847/1538-4357/ab3afcNewberg, H. J., & Carlin, J. L. 2016, Tidal Streams in the LocalGroup and Beyond, Vol. 420, doi: 10.1007/978-3-319-19336-6Newberg, H. J., Yanny, B., Rockosi, C., et al. 2002, ApJ, 569, 245,doi: 10.1086/338983Nie, J. D., Smith, M. C., Belokurov, V., et al. 2015, ApJ, 810, 153,doi: 10.1088/0004-637X/810/2/153Niederste-Ostholt, M., Belokurov, V., & Evans, N. W. 2012,MNRAS, 422, 207, doi: 10.1111/j.1365-2966.2012.20602.x Odenkirchen, M., Grebel, E. K., Rockosi, C. M., et al. 2001, ApJL,548, L165, doi: 10.1086/319095Oliphant, T. E. 2006, A guide to NumPy, Vol. 1 (Trelgol PublishingUSA)Pardy, S. A., D’Onghia, E., & Fox, A. J. 2018, ApJ, 857, 101,doi: 10.3847/1538-4357/aab95bPardy, S. A., D’Onghia, E., Navarro, J. F., et al. 2019, MNRAS,2789, doi: 10.1093/mnras/stz3192Patel, E., Besla, G., & Sohn, S. T. 2017, MNRAS, 464, 3825,doi: 10.1093/mnras/stw2616Patel, E., Kallivayalil, N., Garavito-Camargo, N., et al. 2020, ApJ,893, 121, doi: 10.3847/1538-4357/ab7b75Pe˜narrubia, J., Belokurov, V., Evans, N. W., et al. 2010, MNRAS,408, L26, doi: 10.1111/j.1745-3933.2010.00921.xPe˜narrubia, J., G´omez, F. A., Besla, G., Erkal, D., & Ma, Y.-Z.2016, MNRAS, 456, L54, doi: 10.1093/mnrasl/slv160Perez, F., & Granger, B. E. 2007, Computing in Science andEngineering, 9, 21, doi: 10.1109/MCSE.2007.53Petersen, M. S., & Pe˜narrubia, J. 2020, MNRAS, 494, L11,doi: 10.1093/mnrasl/slaa029Pietrzy´nski, G., Graczyk, D., Gieren, W., et al. 2013, Nature, 495,76, doi: 10.1038/nature11878Pillepich, A., Vogelsberger, M., Deason, A., et al. 2014, MNRAS,444, 237, doi: 10.1093/mnras/stu1408Planck Collaboration, Akrami, Y., Arroja, F., et al. 2018, arXive-prints, arXiv:1807.06205. https://arxiv.org/abs/1807.06205Planck Collaboration, Aghanim, N., Akrami, Y., et al. 2019, arXive-prints, arXiv:1907.12875. https://arxiv.org/abs/1907.12875Price-Whelan, A. M., Sip˝ocz, B. M., G¨unther, H. M., et al. 2018,AJ, 156, 123, doi: 10.3847/1538-3881/aabc4fRamos, P., Mateu, C., Antoja, T., et al. 2020, arXiv e-prints,arXiv:2002.11142. https://arxiv.org/abs/2002.11142Rashkov, V., Pillepich, A., Deason, A. J., et al. 2013, ApJ, 773,L32, doi: 10.1088/2041-8205/773/2/L32Rocha-Pinto, H. J., Majewski, S. R., Skrutskie, M. F., Crane, J. D.,& Patterson, R. J. 2004, ApJ, 615, 732, doi: 10.1086/424585Ruchti, G. R., Read, J. I., Feltzing, S., et al. 2015, MNRAS, 450,2874, doi: 10.1093/mnras/stv807Sesar, B., Hernitschek, N., Dierickx, M. I. P., Fardal, M. A., & Rix,H.-W. 2017, ApJL, 844, L4, doi: 10.3847/2041-8213/aa7c61Sesar, B., Ivezi´c, ˇZ., Lupton, R. H., et al. 2007, AJ, 134, 2236,doi: 10.1086/521819Shipp, N., Drlica-Wagner, A., Balbinot, E., et al. 2018, ApJ, 862,114, doi: 10.3847/1538-4357/aacdabSmith, M. C., Evans, N. W., Belokurov, V., et al. 2009, MNRAS,399, 1223, doi: 10.1111/j.1365-2966.2009.15391.xSmith, R., Choi, H., Lee, J., et al. 2016, ApJ, 833, 109,doi: 10.3847/1538-4357/833/1/109Sohn, S. T., Watkins, L. L., Fardal, M. A., et al. 2018, ApJ, 862,52, doi: 10.3847/1538-4357/aacd0b Springel, V., Wang, J., Vogelsberger, M., et al. 2008, MNRAS,391, 1685, doi: 10.1111/j.1365-2966.2008.14066.xTremaine, S., & Weinberg, M. D. 1984, MNRAS, 209, 729,doi: 10.1093/mnras/209.4.729van der Marel, R. P., Alves, D. R., Hardy, E., & Suntzeff, N. B.2002, AJ, 124, 2639, doi: 10.1086/343775van der Marel, R. P., & Kallivayalil, N. 2014, ApJ, 781, 121,doi: 10.1088/0004-637X/781/2/121van der Walt, S., Colbert, S. C., & Varoquaux, G. 2011, Computingin Science and Engineering, 13, 22, doi: 10.1109/MCSE.2011.37Virtanen, P., Gommers, R., Oliphant, T. E., et al. 2020, NatureMethods, doi: https://doi.org/10.1038/s41592-019-0686-2Vivas, A. K., Zinn, R., Andrews, P., et al. 2001, ApJL, 554, L33,doi: 10.1086/320915Watkins, L. L., van der Marel, R. P., Sohn, S. T., & Evans, N. W.2019, ApJ, 873, 118, doi: 10.3847/1538-4357/ab089fWatkins, L. L., Evans, N. W., Belokurov, V., et al. 2009, MNRAS,398, 1757, doi: 10.1111/j.1365-2966.2009.15242.xWegg, C., Gerhard, O., & Bieth, M. 2019, MNRAS, 485, 3296,doi: 10.1093/mnras/stz572Weinberg, M. D. 1989, MNRAS, 239, 549,doi: 10.1093/mnras/239.2.549 —. 1996, ApJ, 470, 715, doi: 10.1086/177902—. 1998, MNRAS, 299, 499,doi: 10.1046/j.1365-8711.1998.01790.x—. 1999, AJ, 117, 629, doi: 10.1086/300669Weinberg, M. D., & Blitz, L. 2006, ApJL, 641, L33,doi: 10.1086/503607Wetzel, A. R., Hopkins, P. F., Kim, J.-h., et al. 2016, ApJL, 827,L23, doi: 10.3847/2041-8205/827/2/L23White, S. D. M. 1983, ApJ, 274, 53, doi: 10.1086/161425Xue, X.-X., Rix, H.-W., Ma, Z., et al. 2015, ApJ, 809, 144,doi: 10.1088/0004-637X/809/2/144Yu, S., Bullock, J. S., Wetzel, A., et al. 2020, MNRAS,doi: 10.1093/mnras/staa522Yurin, D., & Springel, V. 2014, MNRAS, 444, 62,doi: 10.1093/mnras/stu1421Zolotov, A., Willman, B., Brooks, A. M., et al. 2009, ApJ, 702,1058, doi: 10.1088/0004-637X/702/2/1058—. 2010, ApJ, 721, 738, doi: 10.1088/0004-637X/721/1/738Zonca, A., Singer, L., Lenz, D., et al. 2019, The Journal of OpenSource Software, 4, 1298, doi: 10.21105/joss.01298 (cid:96) values for v R and v θ and peaks at even (cid:96) values for v φ . To understand why these featuresin the power spectra are arising, we show the resulting power spectra from several simple intensity maps using the software starry (Luger et al. 2019) in Figure 10. The top panels of Figure 10 show intensity maps, while the bottom panels show thecorresponding power spectra. While a single spot gives rise to a rather smooth power spectrum, with power in even as well asodd (cid:96) values, two spots of opposing sign (located 180 ◦ apart from one another) result in a power spectrum with peaks at odd (cid:96) values. As seen in the right panels, banding can also give rise to the sawtooth pattern: banding at constant intensity yields peaksat even (cid:96) values, whereas a band weighted by a dipole creates peaks at odd (cid:96) values.To understand why the halos have these features in their velocity maps and power spectra, we explore the properties of twospecific accretion events. First we consider BJ05 satellite 1066 ( M Sat = 9 . × M (cid:12) , t Sat = 4 .
58 Gyr , J Sat / J circ = 0 . ).This halo is accreted both by the recently accreted BJ05 halo and the high- L sat halo. The top panels of Figure 11 show thepositions of all stars from this satellite in Cartesian coordinates, color coded by their radial velocities. Dashed lines indicatethe radial range of − kpc; velocity maps in spherical coordinates are shown in the second row of panels. This massive,recently accreted satellite has formed large shell-type features. When we only consider stars from this satellite that are locatedbetween 30-50 kpc from the center, we see that our cross section does not include the apocenters of these shells (where stars have v R ∼ km s − ), but rather stars that are moving quickly either towards an apocenter (yellower spots) or away from apocenterback towards the center of the galaxy (bluer points). In the radial velocity map, the patches of stars with opposite velocity areseparated by approximately 180 ◦ from each other in this projection. This accretion event has a much narrower range of velocitiesin v θ and v φ ; the v θ map has a fairly smooth continuum, going from positive v θ as the stars moves towards and away the firstapocenter to negative v θ as the stars pass through the second apocenter. This results in banding weighted by a dipole, a versionof the map shown to the rightmost panel in Figure 10. In v φ , the velocities are always negative; the velocity map appears as arelatively constant band.The power spectra for these velocity maps are shown as the purple curves are shown in Figure 12. This radial accretion eventhas the most power in the radial velocity power spectrum; the power spectra peak at odd (cid:96) values for v R and v θ , while the v φ power spectrum is dominated by peaks at even (cid:96) values. P o w e r ( + ) × C Figure 10.
Top panels show intensity maps generated by starry for four different patterns, plotted in Mollweide projection. Bottom panelsshow the corresponding power spectra. In the power spectra, the power at odd (cid:96) values is shown in red. While a single spot generates a rathersmooth spectrum as a function of (cid:96) , when there are two spots of opposite signs across from one another in the map, the resulting power spectrumshows spikes at odd (cid:96) . Banding can also cause spikes, but at even (cid:96) (third panel); if the band is weighted by a dipole, then the spikes occur atodd (cid:96) (rightmost panels). M Sat = 5 . × M (cid:12) ), recently accreted ( t Sat = 6 .
55 Gyr) , circular J Sat /J circ = 0 . ) accretion event leaves aw-wrapped kinematically cold stream that lies approximately in the ( x, y ) plane. When we take the cross section of stars from − kpc, we capture many of the stream stars, which form nearly continuous bands in the Mollweide projection. Because thestream is not perfectly circular, at opposite points in the orbit, the radial and polar velocities will be of similar magnitude but ofopposite sign, while the tangential velocity is approximately constant along the stream.The power spectra are shown as the green curves in Figure 12. This circular accretion event has more power in v θ and v φ than v R . The v R and v θ power spectra again shown peaks at odd (cid:96) values, while the v φ power spectrum is peaked at odd (cid:96) values.The fact that power spectra have these sawtooth patterns is indicative that spherical harmonic expansion of the different velocitycomponents is not the best basis for Galactic substructure, given that the signatures in the angular power spectra are complicated.Potentially using vectorized spherical harmonics or a different coordinate transformation could provide ways forward to addressthis; we leave this to future work.4 Figure 11.
Top panels: ( x, y, z ) positions for all stars originating from BJ05 satellite 1066 ( M Sat = 9 . × M (cid:12) , t Sat =4 .
58 Gyr , J Sat / J circ = 0 . ). Points are color coded by their radial velocities, over the range [ − , km s − . Dashed lines indi-cate − kpc. Second row: velocity maps, for stars in the radial range of − kpc, in the three components of motion in sphericalcoordinates. This massive, recent accretion events creates shells in the halo with strong velocity gradients. When we take a cross section of thissubstructure, we see “spots” in the v R maps of opposite signs that are approximately separated by ◦ . Lower panels: same as top panels, butfor BJ05 satellite 1228 ( M Sat = 5 . × M (cid:12) , t Sat = 6 .
55 Gyr , J Sat / J circ = 0 . ). Points are color coded by their radial velocities, overthe range [ − , km s − . This circular, recent accretion event creates a stream, which create banding features in the velocity maps. ( + ) × C [( k m / s )] v R v v Satellite 1066Satellite 1228
Figure 12.
Power spectra for the accretion events mapped in Figure 11. Both accretion events result in power spectra characterized by asawtooth pattern, with peaks at odd (cid:96) in v R and v θ and even (cid:96) in v φ . v g s r ( k m / s ) ModelBelokurov+2014,trailingBelokurov+2014,leading